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

1、北京大學(xué)生科院/CLS生物信息平臺(tái)RNA-Seq測序數(shù)據(jù)分析服務(wù)流程（試運(yùn)行2015.3平臺(tái)聯(lián)系人：李程（lch3000文檔撰寫：張超Table of Contents1 .測序質(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富集

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

3、umina平臺(tái)單堿基的錯(cuò)誤率應(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橫坐標(biāo)為reads的堿基位置，縱坐標(biāo)為單堿基質(zhì)量值質(zhì)量值與錯(cuò)誤率的關(guān)系：Q=-10log10（e;其中Q phred為測序堿基質(zhì)量值，e為測phred序錯(cuò)誤率。1.3 GC含量分布對于RNA測序,鑒于序列通過超聲隨機(jī)打斷，所以理論上每個(gè)測序循環(huán)上的C、 G及A、T含量應(yīng)分布相等，并且

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

5、率低，很可能建庫時(shí)混入了 其他物種的序列，導(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)計(jì)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 .基因表達(dá)水平3.1 基因表達(dá)水平定量在RNA-seq分析中，我們可以通過定位到基因組區(qū)域或基因外顯子區(qū)的reads的計(jì)數(shù)來估計(jì)基因的表達(dá)水平。Reads計(jì)數(shù)除了與基因的真實(shí)表達(dá)水平成

9、正比外，還 與基因的長度和測序深度成正相關(guān)。為了使不同基因、不同實(shí)驗(yàn)間估計(jì)的基因表達(dá) 水平具有可比性，人們引入了 RPKM的概念,RPKM(Reads Per Kilo bases per Million reads是每百萬reads中來自某一基因每千堿基長度的 reads數(shù)目。RPKM同時(shí)考慮 了測序深度和基因長度對reads計(jì)數(shù)的影響，是目前最為常用的基因表達(dá)水平估算方 法(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各樣本中該基因的表達(dá)水平(RPKM3.2 基因表達(dá)水平分步每個(gè)樣本所有基因的RPKM盒形圖可以展示出不同實(shí)驗(yàn)條件下基因表達(dá)水平 的分布情況。圖3.2.1不同條件下的基因表達(dá)水平分布圖3.3 生物

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

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

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后的表達(dá)量差異(4pva

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

15、聚集成類，從而識別未知基因的功能或已知基因的未知功能。5 .差異表達(dá)基因功能分析5.1 GO富集分析Gene Ontology箭稱 GO, /是基因功能標(biāo)準(zhǔn)分類體 系。研究差異基因在Gene Ontology中的分布情況將闡明差異基因富集的生物學(xué)功 能。M口, Z Scwfl5.2 信號通路富集分析在生物體內(nèi)，不同基因相互協(xié)調(diào)實(shí)現(xiàn)其生物學(xué)功能，通過Pathway顯著性富集能 確定差異表達(dá)基因參與的最主要信號通路。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)計(jì)檢驗(yàn)找出差異表達(dá)基因 中顯著性富集的Pathway。Color Ker啊dl HSg.raip5.3 癌基因功能注釋原癌基因(Proto-oncogene是參與細(xì)胞生長、細(xì)胞分裂和細(xì) 胞分化的正 ?；?，當(dāng)其發(fā)生突變后(如基因序列被改變就會(huì)變成致癌基因 (Oncogene通常在月中瘤或惡性細(xì)胞系中某些特異性癌基因會(huì)上調(diào)表達(dá)，通過了解 癌基因在實(shí)驗(yàn)不同組的表達(dá)情況有助于深入認(rèn)識疾病的發(fā)病機(jī)理。Cosmic( https:cancer.sanger.ac.uk/cosmi

17、是英國 Sanger實(shí)驗(yàn)室開發(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á)，表明了有著不同的剪切方式。這時(shí)可以使用 Bioconductor的DEXSeq軟件包。該分析可以給出每個(gè)基因在不同的實(shí)驗(yàn)條件下， 外顯子的使用情況。比 如上圖的10號外顯子在control和knockdown兩組中的表 達(dá)差別較大，此外顯子的表達(dá)量情況，也反映到了在兩組中此基因的剪切形式有 差異。7. SNP 分析 7.1 SNP 檢測 SNP 全稱 Single Nucleotide Polymorphisms.,是指在 基因組上由單個(gè)核甘酸變異形成的遺傳標(biāo)記，其數(shù)量很多，多態(tài)性豐富