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

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

2、分析（105.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ù)進行進行質(zhì)控，保證數(shù)據(jù)質(zhì)量適合下游分析。這里我們使用 fastqc和RNA-SeQC來對數(shù)據(jù)進行質(zhì)量評定。1.1 測序數(shù)據(jù)過濾測序得到的原始下機數(shù)據(jù)往往有許多問題，不能直接使用，通常會經(jīng)過以下過濾 盡量保證測序數(shù)據(jù)的質(zhì)量。a去除帶測序接頭的測序序列（reads;b.去除低質(zhì)重的reads1.2 質(zhì)量值分布 按照現(xiàn)有的測序技術(shù)（ill

3、umina平臺單堿基的錯誤率應(yīng)控制在1%以下，即質(zhì)量值在20以上。q it 0 8 鼻 I bn- I fwigw I kuni- 1 S it1 tiding I4 二.1*5及l(fā) m 3 *、n 了引slh 6 心 x / 熊，州型力函k m h 及 n re w jo v a 阿 用 酊便由6%隴nd Ibp橫坐標為reads的堿基位置，縱坐標為單堿基質(zhì)量值質(zhì)量值與錯誤率的關(guān)系：Q=-10log10（e;其中Q phred為測序堿基質(zhì)量值，e為測phred序錯誤率。1.3 GC含量分布對于RNA測序,鑒于序列通過超聲隨機打斷，所以理論上每個測序循環(huán)上的C、 G及A、T含量應(yīng)分布相等，并且

4、CG-content對于每個物種應(yīng)大致相同。橫坐標為reads的堿基位置，縱坐標為各種堿基的不同比例2 .參考序列比對對于通過質(zhì)量控制的數(shù)據(jù)，可以進行后續(xù)分析。首先需要將 clean reads比對到 參考基因組上。由于測序時reads是隨機的，只有這些reads的堿基信息和質(zhì)量信息 沒有其在基因組上的位置信息，比對這一步就是給所有reads一個在基因組上位置的 信息。在RNA測序中，其實測的是cDNA的序列，由于內(nèi)含子的存在，所以會較常出現(xiàn) 一條read跨內(nèi)含子的情況，tophat2可以較好的處理這種情況，所以我彳門選用tophat2 來做比對。比對率間接反應(yīng)了測序的質(zhì)量和建庫的質(zhì)量，若比對

5、率低，很可能建庫時混入了 其他物種的序列，導(dǎo)致無法比對到研究的物種參考基因組上。Left reads:Input: 11607353Mapped; 11607353of input)q網(wǎng)電(g.S%】 have nultiple aligiments82 howe 26) Rightinput i11607353Hipped :11607352(1加.被of input)of ttese:60103(have nuttplealtgimentshcnreads比對到基因上的位置統(tǒng)計Sample IntragenicRateExonicRateIntronicRateIntergenicRate

6、SplitReadsExpressionProfilingEfficiencyTranscriptsDetectedGenesDetected1BJ 0.885 0.738 0.147 0.114 9,910,010 0.738 32,796 15,434(Sample才羊本名(2IntragenicRate比對到基因內(nèi)的reads比例(3ExonicRate:比對至U外顯子的reads比例(4IntronicRate:比對到內(nèi)含子的reads比例(5IntergenicRate比對到基因間區(qū)的reads比例(6SplitReads:比對至U兩外顯子交接處的 reads數(shù)(7Expressio

7、nProfilingEfficiency:比對到外顯子上的reads占總體的比例(8TranscriptsDetected比對上reads數(shù)大于5的轉(zhuǎn)錄本數(shù)(9GenesDetectect匕對上reads數(shù)大于5的基因數(shù)Left reads;InpMtt； 11 即7353Mapped: 11667353 aee+0K of input)o6010 ( e,SMJ have nultiple olignnertts Cfi2 Iwe 20 Right reads: input ： 11687353Mapped : 11 瑰7352 Q掘I 限 of iwO of tfww; 60103 ( A

8、L5S0 have nulttple a1 ignments Cq haw* 2) IMOK owrall isud nqppin$ rate.AligMd pdifE： L1W3UoF tK&e:601B3 C 0泉D have raultile oligmerrt上8 C or* d1 scorttont atignr帕nt型 lOO.fl% concordant pair allQ-nent rate.3 .基因表達水平3.1 基因表達水平定量在RNA-seq分析中，我們可以通過定位到基因組區(qū)域或基因外顯子區(qū)的reads的計數(shù)來估計基因的表達水平。Reads計數(shù)除了與基因的真實表達水平成

9、正比外，還 與基因的長度和測序深度成正相關(guān)。為了使不同基因、不同實驗間估計的基因表達 水平具有可比性，人們引入了 RPKM的概念,RPKM(Reads Per Kilo bases per Million reads是每百萬reads中來自某一基因每千堿基長度的 reads數(shù)目。RPKM同時考慮 了測序深度和基因長度對reads計數(shù)的影響，是目前最為常用的基因表達水平估算方 法(Mortazavi et al., 2008。Gene_ID Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 ENSG00000000003 49.32 46.94 48

10、.91 22.51 20.60 22.95 ENSG00000000419 35.92 34.58 33.69 32.80 35.65 32.73 ENSG00000000457 1.34 0.94 1.19 2.06 2.13 2.26ENSG00000000460 1.19 1.20 1.22 3.00 3.33 3.06(1 Gene_ID:Ensembl基因 ID(2 Other columns各樣本中該基因的表達水平(RPKM3.2 基因表達水平分步每個樣本所有基因的RPKM盒形圖可以展示出不同實驗條件下基因表達水平 的分布情況。圖3.2.1不同條件下的基因表達水平分布圖3.3 生物

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

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

13、基因。Gene baseMean log2FoldChange pvalue padjFBgn0000370 31324.379200 -1.3665378519 5.6393206e-176 2.9843284e-OI6J 印1psfi曾ie Gi tc4Ci lerOJ tetO5ftiean e(ortssion172FBgn0033913 17544.483454 -1.1571536021 6.3177309e-90 1.3373372e-87(1Gene:基因 ID(2baseMean所有樣本矯正后的平均reads數(shù)(3log2FoldChange:log2后的表達量差異(4pva

14、lue統(tǒng)計學(xué)差異顯著性檢驗指標(5padj:校正后的pvalue。padj越小，表示基因表達差異越顯著4.3 差異基因可視化用火山圖可以推斷差異基因的整體分布情況。圖4.3.1顯著性差異表達基因用紅色點表示；橫坐標表示基因在不同樣本中的表達倍數(shù)變化；縱坐標表示統(tǒng)計學(xué)上基因表達量變化差異的顯著性對于特定基因在不同實驗中的表達情況，和此基因的不同轉(zhuǎn)錄本在不同實驗中 的表達情況中白，U七1圖4.3.2左圖為regucalcin基因在兩個樣本中的表達差異情況；右圖為此基因在不同轉(zhuǎn)錄本中的表達差異情況4.4 差異基因聚類聚類分析用于判斷差異基因在不同實驗條件下的表達模式。通過將表達模式相 同或相近的基因

15、聚集成類，從而識別未知基因的功能或已知基因的未知功能。5 .差異表達基因功能分析5.1 GO富集分析Gene Ontology箭稱 GO, /是基因功能標準分類體 系。研究差異基因在Gene Ontology中的分布情況將闡明差異基因富集的生物學(xué)功 能。M口, Z Scwfl5.2 信號通路富集分析在生物體內(nèi)，不同基因相互協(xié)調(diào)實現(xiàn)其生物學(xué)功能，通過Pathway顯著性富集能 確定差異表達基因參與的最主要信號通路。KEGG(Kyoto Encyclopedia of Genes andGenomes,http:/www.kegg.jp層有關(guān) P

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

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