RNA-Seq測序數(shù)據(jù)分析服務(wù)流程(試運行)(精)

1、北京大學(xué)生科院/CLS生物信息平臺RNA-Seq測序數(shù)據(jù)分析服務(wù)流程(試運行2015.3平臺聯(lián)系人:李程(文檔撰寫:張超TableofContents1. 測序質(zhì)量評估(31.1 測序數(shù)據(jù)過濾(31.2 質(zhì)量值分布(31.3 GC含量分布(42. 參考序列比對(43. 基因表達(dá)水平(63.1 基因表達(dá)水平定量(63.2 基因表達(dá)水平分步(63.3生物學(xué)重復(fù)相關(guān)性分析(63.4樣本間層次聚類及PCA分析(74. 差異基因分析(74.1 基因表達(dá)標(biāo)準(zhǔn)化(74.2 差異基因列表(84.3 差異基因可視化(84.4 差異基因聚類(95. 差異表達(dá)基因功能分析(105.1 GO富集分析(105.2信號通

2、路富集分析(105.3 癌基因功能注釋(116. 基因結(jié)構(gòu)差異分析(116.1 可變剪切分析(117. SNP分析(127.1 SNP檢測(127.2 SNP篩選(127.3 GO/KEGG富集(121. 測序質(zhì)量評估通過測序的數(shù)據(jù)進(jìn)行進(jìn)行質(zhì)控,保證數(shù)據(jù)質(zhì)量適合下游分析。這里我們使用fastqc和RNA-SeQC來對數(shù)據(jù)進(jìn)行質(zhì)量評定。1.1測序數(shù)據(jù)過濾測序得到的原始下機(jī)數(shù)據(jù)往往有許多問題,不能直接使用,通常會經(jīng)過以下過濾,盡量保證測序數(shù)據(jù)的質(zhì)量。a.去除帶測序接頭的測序序列（reads;b.去除低質(zhì)量的reads1.2質(zhì)量值分布按照現(xiàn)有的測序技術(shù)（illumina平臺單堿基的錯誤率應(yīng)控制在1%

3、以下，即質(zhì)量值在20以上。橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基質(zhì)量值質(zhì)量值與錯誤率的關(guān)系:Q=-10log10（e淇中Qphred為測序堿基質(zhì)量值,e為測phred序錯誤率1.3 GC含量分布對于RNA測序,鑒于序列通過超聲隨機(jī)打斷，所以理論上每個測序循環(huán)上的C、G及A、T含量應(yīng)分布相等，并且CG-content對于每個物種應(yīng)大致相同。橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為各種堿基的不同比例2. 參考序列比對對于通過質(zhì)量控制的數(shù)據(jù)，可以進(jìn)行后續(xù)分析。首先需要將cleanreads比對到參考基因組上。由于測序時reads是隨機(jī)的，只有這些reads的堿基信息和質(zhì)量信息,沒有其在基因組上的

4、位置信息，比對這一步就是給所有reads一個在基因組上位置的在RNA測序中，其實測的是cDNA的序列，由于內(nèi)含子的存在，所以會較常出現(xiàn)一條read跨內(nèi)含子的情況,tophat2可以較好的處理這種情況,所以我們選用tophat2來做比對。比對率間接反應(yīng)了測序的質(zhì)量和建庫的質(zhì)量，若比對率低，很可能建庫時混入了其他物種的序列，導(dǎo)致無法比對到研究的物種參考基因組上。Leftreods；Input:11607353M叩詛:11607353Ciee.effi-ofinputofthese：國1躬C仏鑒litivfimultipleallgmnEntsC02Itove20)Eightreads;input:

5、n翩吟科Mapped:11闘7352(1轉(zhuǎn).劇ofW幣吐ofthese:C站.531)multipleallprwentshove30)reads比對到基因上的位置統(tǒng)計:SampleIntragenicRateExonicRateIntronicRateIntergenicRateSplitReadsExpressionProfilingEfficiencyTranscriptsDetectedGenesDetected1BJ0.8850.7380.1470.1149,910,0100.73832,79615,434(ISample:樣本名(2IntragenicRate:比對到基因內(nèi)的rea

6、ds比例(3ExonicRate:比對到外顯子的reads比例(4IntronicRate:比對到內(nèi)含子的reads比例(5IntergenicRate:比對到基因間區(qū)的reads比例(6SplitReads:比對到兩外顯子交接處的reads數(shù)(7ExpressionProfilingEfficiency:比對到外顯子上的reads占總體的比例(8TranscriptsDetected:比對上reads數(shù)大于5的轉(zhuǎn)錄本數(shù)(9GenesDetected:比對上reads數(shù)大于5的基因數(shù)Let*wds；Irut:U6073S3Mapped:U607353CL幅”瞬ofnputofthaaei闊1酣

7、C戯菠havflimilttplaalLgnin?itC82haveRightreads;input:11607353Happ記:口翩咱5Zfl軸蔣ofnpirtaftese:曲1虻C枳圖叼叫機(jī)8口噲5gImve曲109JQS9overallreadmopingrale-AKgndpairs：11Q73SJFLIh&b;曲調(diào)C0-590liavenutLipledlLgiwenLs酉CE-0arediscordantaLi.gmerrtscanMi-Jarrtpair口Hgntnarrtata3. 基因表達(dá)水平3.1基因表達(dá)水平定量在RNA-seq分析中，我們可以通過定位到基因組區(qū)域或基因外顯

8、子區(qū)的reads的計數(shù)來估計基因的表達(dá)水平。Reads計數(shù)除了與基因的真實表達(dá)水平成正比外，還與基因的長度和測序深度成正相關(guān)。為了使不同基因、不同實驗間估計的基因表達(dá)水平具有可比性，人們引入了RPKM的概念,RPKM(ReadsPerKilobasesperMillionreads是每百萬reads中來自某一基因每千堿基長度的reads數(shù)目。RPKM同時考慮了測序深度和基因長度對reads計數(shù)的影響，是目前最為常用的基因表達(dá)水平估算方法(Mortazavietal.,2008。Gene_IDSample1Sample2Sample3Sample4Sample5Sample6ENSG000000

9、0000349.3246.9448.9122.5120.6022.95ENSG0000000041935.9234.5833.6932.8035.6532.73ENSG000000004571.340.941.192.062.132.26ENSG0000000046023.003.333.06(1Gene_ID:Ensembl基因ID(2Othercolumns:各樣本中該基因的表達(dá)水平(RPKM3.2基因表達(dá)水平分步每個樣本所有基因的RPKM盒形圖可以展示出不同實驗條件下基因表達(dá)水平的分布情況。日百I亡1片和圖321不同條件下的基因表達(dá)水平分布圖3.3生物學(xué)重復(fù)相關(guān)性分

10、析生物學(xué)重復(fù)主要有兩個用途:一個是證明所涉及的生物學(xué)實驗可重復(fù)性強(qiáng)、差異小，另一個用于估計生物學(xué)變異進(jìn)行差異基因檢測。樣品間基因表達(dá)水平相關(guān)性是檢驗實驗可靠性和樣本選擇是否合理的重要指標(biāo)。相關(guān)系數(shù)越接近1,表明樣品之間表達(dá)模式的相似度越高口an爲(wèi)aILU圖331生物學(xué)重復(fù)散點圖3.4樣本間層次聚類及PCA分析當(dāng)樣本數(shù)目較多時，可以利用基因的表達(dá)量進(jìn)行樣本間聚類分析及PCA分析，對樣本間關(guān)系進(jìn)行探究或者對實驗設(shè)計進(jìn)行驗證。樣本聚類距離或者PCA距離越近,說明樣本越相似。gvi643gvifi44. 差異基因分析4.1基因表達(dá)標(biāo)準(zhǔn)化對于有生物學(xué)重復(fù)的樣品，我們采用DESeq2提出的scalingf

11、actor的方法對原始的readcount進(jìn)行標(biāo)準(zhǔn)化（normalization。以消除非生物學(xué)引起的readcount的差異（最主要消除各個文庫測序數(shù)據(jù)量不同帶來的差異。對于標(biāo)準(zhǔn)化的結(jié)果，我們采用MA-plot或box-plot來評價。圖4.1.1MA-plot橫坐標(biāo)為表達(dá)量，縱坐標(biāo)為log后的表達(dá)差異倍數(shù)基于大部分基因都是非差異表達(dá)的，所以大多點應(yīng)在logfoldchange=0左右，并且不隨表達(dá)量的變化而變化。4.2差異基因列表對于有生物學(xué)重復(fù)的的樣品，我們采用DESeq2來分析差異表達(dá)基因。該方法基于負(fù)二項分布模型(KijNB(“ij，來檢測差異表達(dá)基因。GenebaseMeanlog

12、2FoldChangepvaluepadjFBgn000037031324.379200-1.36653785195.6393206e-1762.9843284e-DESeq2CM!r甲-le-0!telZJW心fElQjnSXpreSSID!172FBgn003391317544.483454-13177309e-901.3373372e-87(1Gene:基因ID(2baseMean:所有樣本矯正后的平均reads數(shù)(31og2FoldChange:log2后的表達(dá)量差異(4pvalue:統(tǒng)計學(xué)差異顯著性檢驗指標(biāo)(5padj:校正后的pvalue。padj越小，表

13、示基因表達(dá)差異越顯著4.3差異基因可視化用火山圖可以推斷差異基因的整體分布情況。sene;3-弓rr-asEI圖4.3.1顯著性差異表達(dá)基因用紅色點表示；橫坐標(biāo)表示基因在不同樣本中的表達(dá)倍數(shù)變化；縱坐標(biāo)表示統(tǒng)計學(xué)上基因表達(dá)量變化差異的顯著性對于特定基因在不同實驗中的表達(dá)情況，和此基因的不同轉(zhuǎn)錄本在不同實驗中的表達(dá)情況lwl-lA-A理FI圖4.3.2左圖為regucalcin基因在兩個樣本中的表達(dá)差異情況;右圖為此基因在不同轉(zhuǎn)錄本中的表達(dá)差異情況4.4差異基因聚類聚類分析用于判斷差異基因在不同實驗條件下的表達(dá)模式。通過將表達(dá)模式相同或相近的基因聚集成類,從而識別未知基因的功能或已知基因的未知功

14、能Clqh-Kn-i-ELL一禹-:H-EL+?L!HJ1二目-iiulwwn-e_f=_ItwE5. 差異表達(dá)基因功能分析5.1GO富集分析GeneOntology（簡稱GO,/是基因功能標(biāo)準(zhǔn)分類體系。研究差異基因在GeneOntology中的分布情況將闡明差異基因富集的生物學(xué)功5.2信號通路富集分析在生物體內(nèi)，不同基因相互協(xié)調(diào)實現(xiàn)其生物學(xué)功能，通過Pathway顯著性富集能確定差異表達(dá)基因參與的最主要信號通路。KEGG(KyotoEncyclopediaofGenesandGenomes,http:/www.kegg.jp/是有關(guān)Path

15、way的主要公共數(shù)據(jù)庫(Kanehisa,2008。Pathway顯著性富集分析以KEGGPathway為單位，應(yīng)用統(tǒng)計檢驗找出差異表達(dá)基因中顯著性富集的PathwayoCoIorKeyAidHiSlflBJWl5.3癌基因功能注釋原癌基因(Proto-oncogene是參與細(xì)胞生長、細(xì)胞分裂和細(xì)胞分化的正常基因，當(dāng)其發(fā)生突變后(如基因序列被改變就會變成致癌基因(Oncogene。通常在腫瘤或惡性細(xì)胞系中某些特異性癌基因會上調(diào)表達(dá)，通過了解癌基因在實驗不同組的表達(dá)情況有助于深入認(rèn)識疾病的發(fā)病機(jī)理。Cosmic(https:/cancer.sanger.ac.uk/cosmic是英國Sanger

16、實驗室開發(fā)并維護(hù)的癌基因及相關(guān)注釋數(shù)據(jù)庫，有較高的權(quán)威性及可信度，通過與數(shù)據(jù)庫進(jìn)行比對，可對差異表達(dá)基因中的癌基因部分進(jìn)行鑒別和注釋。6基因結(jié)構(gòu)差異分析6.1可變剪切分析對于RNA-seq，除了gene水平的差異分析外，還可以進(jìn)行exon水平的差異分析。不用的exon表達(dá)，表明了有著不同的剪切方式。這時可以使用Bioconductor的DEXSeq軟件包。該分析可以給出每個基因在不同的實驗條件下，外顯子的使用情況。比如上圖的10號外顯子在control和knockdown兩組中的表達(dá)差別較大，此外顯子的表達(dá)量情況，也反映到了在兩組中此基因的剪切形式有差異。7.SNP分析7.1SNP檢測SNP全稱SingleNucleotidePolymorphisms，是指在基因組上由單個核苷酸變異形成的遺傳標(biāo)記，其數(shù)量很