RNA-seq分析流程-概述4836

RNA-seq分析流程——概述前?接下來我們要介紹的是RNA-seq數(shù)據(jù)的處理分析流程，根據(jù)RNA-seq測(cè)序技術(shù)的不同，可以分為三種：Starketal.NatRevGenet(2019)1.short-read2.long-read3.directRNA-seq?我們?般的RNA-seq測(cè)序數(shù)據(jù)分析流程算法，基本上都是基于short-read（短讀長）技術(shù)所產(chǎn)?的數(shù)據(jù)?件?前，我們可以從ShortReadArchive(SRA)數(shù)據(jù)庫獲取的RNA-seq數(shù)據(jù)中，有超過95%的數(shù)據(jù)是由Illumina公司的shortread測(cè)序技術(shù)所產(chǎn)?的其分析過程可以?下?的路線圖表?Conesaetal.GenomeBiology(2016)該路線圖?致分為三個(gè)部分：1.數(shù)據(jù)獲?。喊▽?shí)驗(yàn)設(shè)計(jì)、測(cè)序設(shè)計(jì)以及數(shù)據(jù)下機(jī)后的rawreads數(shù)據(jù)的質(zhì)控2.數(shù)據(jù)分析在獲取到?凈的數(shù)據(jù)之后，可以進(jìn)?reads的?對(duì)，然后進(jìn)?基因表達(dá)的量化、差異表達(dá)分析、功能富集分析等3.?級(jí)分析包括數(shù)據(jù)的可視化，其他?分?RNA分析、融合分析以及與其他類型的數(shù)據(jù)進(jìn)?整合分析等?我們分析的起始點(diǎn)，是從原始數(shù)據(jù)開始的，也就是獲取rawreads數(shù)據(jù)。通常這種?通量測(cè)序數(shù)據(jù)會(huì)保存為FASTQ格式的?件。FASTQ格式是?種以ASCII碼字符的形式保存?物序列及其對(duì)應(yīng)的每個(gè)堿基的質(zhì)量的?本?件。FASTQ?件中每條序列（通常是?條read）是由4?組成，其中：第??以@字符開頭，之后的字符為序列的標(biāo)識(shí)符和描述信息第??為具體的序列第三?以+符號(hào)開頭，之后可以可選地加上與第???樣的序列標(biāo)識(shí)或描述信息第四?為堿基質(zhì)量分?jǐn)?shù)（Phred），其字符數(shù)量與第??相等，每個(gè)字符表?對(duì)應(yīng)堿基的質(zhì)量得分，例如@SEQ_IDGATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT+!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65其中，堿基質(zhì)量值的編碼?式為1.先將堿基錯(cuò)誤率P進(jìn)?負(fù)對(duì)數(shù)轉(zhuǎn)換，得到Q值2.然后將Q值加上33或64得到的值所對(duì)應(yīng)的ASCII碼即為堿基質(zhì)量分?jǐn)?shù)例如，錯(cuò)誤率P=0.01，則Q=20，如果是Phred33則對(duì)應(yīng)的質(zhì)量為字符5（53），如果是Phred64則對(duì)應(yīng)的字符為T(84)分析流程1.數(shù)據(jù)獲取?般情況下，如果??有送樣檢測(cè)數(shù)據(jù)的話，測(cè)序公司會(huì)提供原始的FASTQ格式的數(shù)據(jù)。如果我們要使?別??章中發(fā)表的公開數(shù)據(jù)，還需要從數(shù)據(jù)庫中下載對(duì)應(yīng)的數(shù)據(jù)例如，我們從SRA數(shù)據(jù)中下載的原始測(cè)序?件是sra格式，我們需要先使??具將其轉(zhuǎn)換為FASTQ格式2.質(zhì)量控制主要在三個(gè)地需要對(duì)數(shù)據(jù)的質(zhì)量進(jìn)?監(jiān)控獲取原始數(shù)據(jù)之后?對(duì)完之后表達(dá)定量之后2.1rawread對(duì)rawreads數(shù)據(jù)進(jìn)?質(zhì)量控制，需要分析序列的質(zhì)量、GC含量、是否存在接頭、短重復(fù)序列的分布、測(cè)序錯(cuò)誤以及PCR重復(fù)和污染質(zhì)控軟件：FastQC：?于分析Illumina測(cè)序平臺(tái)的數(shù)據(jù)NGSQC：可應(yīng)?于所有平臺(tái)?般來說，reads的質(zhì)量會(huì)朝著3'端遞減，如果堿基的質(zhì)量太低，我們需要?jiǎng)h除它以提??對(duì)率FASTX-Toolkit和Trimmomatic兩個(gè)軟件可以?于切除低質(zhì)量的堿基和接頭序列2.2?對(duì)后reads通常需要?對(duì)到?個(gè)參考基因組或轉(zhuǎn)錄組，??對(duì)的質(zhì)量是評(píng)估測(cè)序準(zhǔn)確率和是否存在DNA污染的?個(gè)重要指標(biāo)?對(duì)質(zhì)量通常為?對(duì)到的reads數(shù)占總reads數(shù)的?例。例如，?對(duì)到?類參考基因組的?對(duì)質(zhì)量通常需要在70-90%，且有?量的reads映射到?個(gè)相同的區(qū)間內(nèi)。如果是?對(duì)到轉(zhuǎn)錄本上，由于可變剪切的影響，可以適當(dāng)放寬?對(duì)質(zhì)量在外顯?和?對(duì)?向上的read覆蓋率的均?性，也是評(píng)估質(zhì)量的重要指標(biāo)。如果reads主要聚集在轉(zhuǎn)錄本的3'端，可能表明原始樣本的RNA質(zhì)量較低?對(duì)上的reads的GC含量，可能揭?了PCR的錯(cuò)誤率主要軟件有：Picard、RSeQC和Qualimap2.3定量后在計(jì)算完表達(dá)的量化值之后，可以計(jì)算GC含量和基因長度的誤差，在必要時(shí)可以使?標(biāo)準(zhǔn)化?法來進(jìn)?校正如果參考轉(zhuǎn)錄組注釋得很好，則可以分析樣本的?物構(gòu)成，來評(píng)估RNA純化步驟的質(zhì)量。例如，rRNA和smallRNA不能出現(xiàn)在polyAlongRNA的制備中NOISeq和EDASeq等R包可以使?圖形來展?count數(shù)據(jù)的質(zhì)量控制2.4可重復(fù)性上?的質(zhì)量控制都只是針對(duì)單個(gè)樣本的，此外，不同樣本之間的可重復(fù)性評(píng)估，對(duì)于評(píng)價(jià)整個(gè)數(shù)據(jù)集的質(zhì)量也是?關(guān)重要的技術(shù)重復(fù)樣本的可重現(xiàn)性?般很?（spearman），但是?物學(xué)重復(fù)樣本之間并沒有明確的標(biāo)準(zhǔn)，取決于實(shí)驗(yàn)系統(tǒng)的異質(zhì)性。如果不同實(shí)驗(yàn)系統(tǒng)之間存在差異基因，則同?條件下的?物學(xué)重復(fù)在主成分分析（PCA）中會(huì)被聚類在?起。3.序列?對(duì)在對(duì)樣本的rawreads進(jìn)?質(zhì)控之后，就可以進(jìn)?序列?對(duì)了，序列?對(duì)主要有三種策略，如下圖Conesaetal.GenomeBiology(2016)如果有參考序列，根據(jù)參考序列的不同，可以分為?對(duì)到基因組：使?間隔?對(duì)算法，如TopHat、STAR等，然后根據(jù)是否提供了注釋?件（GFF格式?件，包含轉(zhuǎn)錄本位置信息），?可以分為轉(zhuǎn)錄本識(shí)別和轉(zhuǎn)錄本發(fā)現(xiàn)并進(jìn)?定量分析?對(duì)到轉(zhuǎn)錄組：使??間隔?對(duì)算法，如Bowtie等，然后使?RSEM或Kallisto?法識(shí)別轉(zhuǎn)錄本并計(jì)算定量信息如果沒有參考序列，則需要先把序列組裝成轉(zhuǎn)錄本，再將reads?對(duì)到組裝后的參考轉(zhuǎn)錄本上，然后使?HTseq-count等算法對(duì)轉(zhuǎn)錄本進(jìn)?定量3.1轉(zhuǎn)錄本發(fā)現(xiàn)使?Illumina技術(shù)檢測(cè)的shortreads來發(fā)現(xiàn)新的轉(zhuǎn)錄本是RNA-seq分析中的?個(gè)挑戰(zhàn)。通常來說，短reads很少會(huì)跨越多個(gè)剪切位點(diǎn)，這就很難直接推斷出?個(gè)轉(zhuǎn)錄本的整體長度。此外，轉(zhuǎn)錄的起始和終?位置也?較難識(shí)別，?些像GRIT的?具，通過合并5'端的信息可以提?異構(gòu)體識(shí)別的準(zhǔn)確性。其他如Cufflinks、iReckon、SLIDE和StringTie等?法，通過結(jié)合現(xiàn)有的注釋信息，作為?個(gè)可能的異構(gòu)體列表?些尋找基因的?具，如Augustus，結(jié)合RNA-seq數(shù)據(jù)，可以更好的注釋蛋?編碼轉(zhuǎn)錄本，但是對(duì)?編碼轉(zhuǎn)錄本的性能更差。3.2Denovo轉(zhuǎn)錄本重構(gòu)在沒有轉(zhuǎn)錄本或轉(zhuǎn)錄本不全的情況下，可以對(duì)reads進(jìn)?組裝來重構(gòu)?份轉(zhuǎn)錄本?？蛇x的?法很多，如SOAPdenovoTrans、Oases、Trans-ABySS或Trinity通常來說，使?雙端鏈特異性測(cè)序和longreads測(cè)序包含更多的信息，會(huì)有更好的效果雖然，對(duì)于低表達(dá)的轉(zhuǎn)錄本進(jìn)?組裝的可靠性較低，但是reads太多也會(huì)導(dǎo)致潛在的組裝錯(cuò)誤和較長的時(shí)間消耗等問題。因此，在深度測(cè)序的樣本中，可以適當(dāng)減少reads的數(shù)量對(duì)于多樣本的?較，可以將所有樣本作為?個(gè)輸?來構(gòu)建參考轉(zhuǎn)錄本，然后分別對(duì)每個(gè)樣本的reads進(jìn)??對(duì)?論是使?參考序列還是從頭開始組裝，使?短reads的Illumina技術(shù)來完全重構(gòu)轉(zhuǎn)錄組仍然是?個(gè)具有挑戰(zhàn)性的問題4.轉(zhuǎn)錄組定量RNA-seq最?泛的應(yīng)?就是?來評(píng)估基因和轉(zhuǎn)錄本的表達(dá)，這?應(yīng)?主要是基于?對(duì)到轉(zhuǎn)錄組區(qū)間內(nèi)的reads的數(shù)量最簡單的?法是，使?HTSeq-count或featureCounts計(jì)算區(qū)間內(nèi)的reads數(shù)來量化基因的表達(dá)。這種基因?平的（不是轉(zhuǎn)錄本?平）的量化?法使?的是GTF?件，這種?件包含外顯?和基因在基因組上的坐標(biāo)。但?般不能直接使?readcount來?較基因的表達(dá)?平，因?yàn)樵撝禃?huì)受到轉(zhuǎn)錄本長度、reads總數(shù)以及測(cè)序偏差等因素的影響。所以需要先進(jìn)?標(biāo)準(zhǔn)化，標(biāo)準(zhǔn)化?法有1.RPKM/FPKM:每百萬reads每?千堿基對(duì)中包含的reads數(shù)該?法先計(jì)算測(cè)序深度系數(shù)，即總reads數(shù)除以?百萬，然后計(jì)算基因或轉(zhuǎn)錄本的長度（單位為kb），標(biāo)準(zhǔn)化順序?yàn)橄认郎y(cè)序深度的影響，再消除長度的影響：其中x表??個(gè)基因或轉(zhuǎn)錄本，或基因組上?段特定的區(qū)域表??對(duì)到x外顯?區(qū)域的reads數(shù)；R表?當(dāng)前樣本中包含的全部reads數(shù)表?x外顯?區(qū)域包含的堿基數(shù)（長度，bp）FPKM與RPKM的計(jì)算公式?樣，只是RPKM?于單端測(cè)序，F(xiàn)PKM?于雙端測(cè)序2.TPM:其與RPKM最?的區(qū)別是，標(biāo)準(zhǔn)化順序?yàn)橄认蜷L度的影響，再消除測(cè)序深度的影響?先，將readscount除以基因或轉(zhuǎn)錄本的長度（kb）得到RPK(readsperkilobase)，然后將樣本中所有的RPK加起來除以化系數(shù)，最后使?RPK除以標(biāo)注化系數(shù)，得到標(biāo)準(zhǔn)其中x表??個(gè)基因或轉(zhuǎn)錄本，或基因組上?段特定的區(qū)域表??對(duì)到x外顯?區(qū)域的reads數(shù)表?x外顯?區(qū)域包含的堿基數(shù)（kp）N表?基因或轉(zhuǎn)錄本總數(shù)這樣，每個(gè)樣本的TPM總和是?樣的，便于?較樣本間的差異?前，也有許多復(fù)雜的算法通過解決相關(guān)轉(zhuǎn)錄本共享reads的問題來評(píng)估轉(zhuǎn)錄本?平的表達(dá)，例如，Cufflinks使?TopHat的?對(duì)結(jié)果，應(yīng)?期望最?化算法來評(píng)估轉(zhuǎn)錄本的豐度。這??法考慮到長度不同的基因的reads分布并不均勻等因素的影響。還有其他算法也可以量化轉(zhuǎn)錄組的表達(dá)，例如RSEM、eXpress、Sailfish和kallisto等。這些?法允許轉(zhuǎn)錄本之間存在多?對(duì)的reads，并輸出經(jīng)測(cè)序偏差校正的樣本內(nèi)歸?化值。5.差異表達(dá)分析差異表達(dá)分析是對(duì)樣本間基因的表達(dá)值進(jìn)??較，雖然RPKM、FPKM和TPM標(biāo)準(zhǔn)化?法消除了測(cè)序深度和基因或轉(zhuǎn)錄本的長度因素的影響，但這些?法依賴于總的或有效的reads數(shù)，當(dāng)樣本的具有異質(zhì)性轉(zhuǎn)錄本分布或當(dāng)?表達(dá)或差異表達(dá)的特征扭曲了count分布時(shí)，表現(xiàn)?佳?像TMM、DESeq、PoissonSeq和UpperQuartile等?法會(huì)忽略?變異或?表達(dá)的特征。?擾樣本內(nèi)?較的其他因素包括不同樣本的轉(zhuǎn)錄本長度變化、轉(zhuǎn)錄本覆蓋位置的偏差、平均?段??以及基因的GC含量等NOISeq這個(gè)R包提供了多種繪圖，來識(shí)別RNA-seq數(shù)據(jù)中的誤差來源，并應(yīng)?相應(yīng)的?法來標(biāo)準(zhǔn)化這些誤差除了這些樣本內(nèi)特異的標(biāo)準(zhǔn)化?法，還需要解決數(shù)據(jù)集之間的批次效應(yīng)（不同實(shí)驗(yàn)條件下產(chǎn)?的數(shù)據(jù)之間存在的差異），批次矯正?法有COMBAT和ARSyN等，雖然這些?法是針對(duì)芯?數(shù)據(jù)設(shè)計(jì)的，但是在RNA-seq數(shù)據(jù)中也有很好的效果計(jì)算差異表達(dá)的?法有很多，有些?法，如edgeR將原始的readcounts作為輸?，并在統(tǒng)計(jì)模型中加?了標(biāo)準(zhǔn)化，另?些?法，需要先對(duì)數(shù)據(jù)進(jìn)?標(biāo)準(zhǔn)化，如DESeq2使?的是負(fù)?項(xiàng)分布作為參考分布，并提供了??的標(biāo)準(zhǔn)化?法。baySeq和EBSeq是貝葉斯?法，還有?些基于線性模型的?法。最后，?些?參數(shù)?法，如NOISeq和SAMseq對(duì)于?樣本量的研究，負(fù)?項(xiàng)分布會(huì)存在噪?污染，這種情況下，?些簡單點(diǎn)?法，如基于Poisson分布的DEGseq，或者基于經(jīng)驗(yàn)分布的NOISeq可能會(huì)更好些。但是需要強(qiáng)調(diào)的是，在沒有?夠?物學(xué)重復(fù)的情況下，?法進(jìn)?總體的推斷，因此任何p值計(jì)算都是?效的。許多獨(dú)?的研究都已經(jīng)證實(shí)，選擇不同的?法會(huì)對(duì)結(jié)果有?定的影響，?且沒有哪?種?法能夠適?于所有的數(shù)據(jù)，所以，推薦在分析的時(shí)候使?多個(gè)軟件進(jìn)?相互驗(yàn)證。6.可變剪切分析可變剪接(AlternativeSplicing)是指轉(zhuǎn)錄形成的前體RNA通過去除內(nèi)含?、連接外顯??形成成熟RNA的過程，從?實(shí)現(xiàn)?個(gè)基因同時(shí)編碼多種蛋?質(zhì)，實(shí)現(xiàn)?物功能多樣性在不同組織或者發(fā)育的不同階段，可變剪接不是?成不變的，在特定的組織或條件下，通過連接不同的外顯?，會(huì)產(chǎn)?特定的剪接異構(gòu)體（isoform）。有?量的研究發(fā)現(xiàn)，可變剪接的變化與癌癥等多種疾病相關(guān)，所以研究可變剪接在不同組織中的作?是?常有意義的。轉(zhuǎn)錄本?平的差異表達(dá)分析可以潛在地檢測(cè)同?基因的轉(zhuǎn)錄異構(gòu)體表達(dá)的變化，已經(jīng)有?些算法應(yīng)?于RNA-seq數(shù)據(jù)的中進(jìn)?可變剪切分析這些?法主要分為兩?類：1.異構(gòu)體表達(dá)估計(jì)與差異表達(dá)相結(jié)合，來揭?總基因表達(dá)中每種異構(gòu)體的?例變化例如，BASIS?法使?分層貝葉斯模型來直接推斷轉(zhuǎn)錄異構(gòu)體的差異表達(dá)；CuffDiff2?法先評(píng)估異構(gòu)體的表達(dá)，然后?較它們之間的差異；rSeqDiff?法使?分層似然率檢驗(yàn)同時(shí)檢測(cè)?剪接變化的差異基因表達(dá)和差異異構(gòu)體表達(dá)。所有這些?法通常都受限于短讀長測(cè)序的內(nèi)在局限性，?法在異構(gòu)體?平上進(jìn)?準(zhǔn)確識(shí)別2.?種所謂的exon-based的?法，它跳過了對(duì)異構(gòu)體表達(dá)的估計(jì)，通過?較樣本之間基因外顯?和連接點(diǎn)上的reads分布來檢測(cè)可變剪接的信號(hào)其基本假設(shè)為：可以在外顯?及其連接點(diǎn)的信號(hào)中追蹤異構(gòu)體表達(dá)的差異。DEXseq和DSGSeq采?類似的思路，通過檢測(cè)基因的外顯?（和連接點(diǎn)）上readcounts的差異顯著性來識(shí)別不同的異構(gòu)體。rMATS是通過?較?連接點(diǎn)的reads定義的外顯?inclusionlevels表達(dá)?平的差異7.融合分析基因融合是指兩個(gè)基因的全部或?部分的序列相互融合為?個(gè)新的基因的過程。其有可能是染?體易位、中間缺失或染?體倒置所導(dǎo)致的，可在DNA或RNA層?上表達(dá)。融合基因通過基因失調(diào)、融合產(chǎn)?嵌合體蛋?這兩種機(jī)制引發(fā)癌癥的發(fā)?。?前，RNA-seq融合算法100多種，有?對(duì)常?的15中融合檢測(cè)算法進(jìn)?了?較LiuetalNucleicAcidsResearch,2016沒有哪?個(gè)算法具有明顯的優(yōu)勢(shì)，整體來看，SOAPfuse可能會(huì)好?些，F(xiàn)usionCatcher和JAFFA其次。8.功能注釋標(biāo)準(zhǔn)的轉(zhuǎn)錄組分析的最后?步，是使?差異表達(dá)基因來進(jìn)?功能或通路的注釋。最常?的兩類?法是：基于超?何分布的過表達(dá)富集分析GSEA富集分析?些?具，如GOseq考慮了基因長度等因素對(duì)差異表達(dá)結(jié)果的影響，并使?超?何分布進(jìn)?富集分析，GSVA和SeqGSEA使?類似GSEA的?法進(jìn)?功能富集功能富集需要預(yù)先定義的基因集合或通路，包括GO、KEGG、Reactome等數(shù)據(jù)庫。通過在蛋?質(zhì)數(shù)據(jù)庫（例如SwissProt）和包含保守蛋?質(zhì)結(jié)構(gòu)域（例如Pfam和InterPro）的數(shù)據(jù)庫中搜索相似序列，使?直系同源分析對(duì)蛋?質(zhì)