一篇RNA-seq分析流程的綜述，全面而詳細！深度好文，可用來反復閱讀。初學者用于把握RNA-seq真?zhèn)€流程及各個流程選擇上的差異。已經(jīng)開始學習者可用來查缺補漏和發(fā)現(xiàn)新的分析角度。
A survey of best practices for RNA-seq data analysis
摘要：
沒有任何一個RNA-seq分析流程可適用于所有的轉(zhuǎn)錄組分析。討論RNA-seq分析流程主要步驟：實驗設計，質(zhì)控，比對，基因水平和轉(zhuǎn)錄組水平定量，可視化，基因差異表達，可變剪接，功能分析，融合基因檢測，eQTL (expression quantification trait loci,表達數(shù)量性狀位點)。展望轉(zhuǎn)錄組研究存在的問題。
背景：
研究材料基因組信息已知，通過將RNA-seq獲得的序列比對到基因組上獲得轉(zhuǎn)錄信息；研究材料無基因組信息則從頭拼接reads為contigs后將reads比對到轉(zhuǎn)錄組。
基因組注釋已知，基于注釋基因組進行轉(zhuǎn)錄組分析或發(fā)挖掘新的轉(zhuǎn)錄組及其調(diào)控通路。其次研究者可以對感興趣的mRNA亞型表達或microRNA水平或等位變異分析。在此分析過程中可以只進行RNA-seq分析也可以聯(lián)合其他組學一起分析。
不同的RNA-seq分析有不同的轉(zhuǎn)錄組定量，均一化以及差異表達分析，并且質(zhì)控可確保結(jié)果的可重復性和可靠性。圖一為Illumina sequencing實驗設計、分析流程圖。簡單羅列一些數(shù)據(jù)及圖例來說明這些分析中潛在的不足。最后討論single cell RNA-seq（單細胞轉(zhuǎn)錄組）及測序長度比較（3代測序和2代測序）。
實驗設計：
文庫類型、測序深度、重復，準確的實驗操作以確保數(shù)據(jù)未被污染。
首先：RNA提取中去除大量存在的rRNA， 通常占總RNA的90%，mRNA為1-2%。
提取mRNA可選擇用ployA選擇性富集mRNA或刪除rRNA。ployA通過RNA intergrity number （RIN，RNA完整度）來表示mRNA的比例，對于不能產(chǎn)生高質(zhì)量和足夠數(shù)量的材料則用刪除rRNA法來獲得mRNA（例如細菌mRNA無多聚A）。
另一個問題是：是否產(chǎn)生strand-preserving libraries， strand-specific protocols 如dUTP法，通過在第二條cDNA合成時加入UTP，先于接頭連接隨后含有dUTP的鏈被降解。測序長度小于500bp，分單端測序（single end，SE）和雙端測序（paired-end，PE）。讀長較長（long reads）的序列及雙端序列更有利于注釋信息較差的轉(zhuǎn)錄組分析。
其次：測序深度及文庫大小。測序較深的到的轉(zhuǎn)錄組信息及轉(zhuǎn)錄本數(shù)量更加詳細，但不是越深越好。?
5百萬條比對序列對中到高表達基因的量化分析足夠，100萬條序列足以分析低表達基因分析，單細胞轉(zhuǎn)錄組通常為1百萬，高表達基因測序5萬，脾組織只需2萬。
文庫大小取決于目標轉(zhuǎn)錄組的復雜程度，測序深度有利于轉(zhuǎn)錄本的數(shù)量和鑒定，但同時增加了雜質(zhì)信息和脫靶轉(zhuǎn)錄本。飽和曲線可以用來評估給定測序深度下轉(zhuǎn)錄組的覆蓋度。
最后：樣本重復，包括測序時不同批次的差異及樣本的差異。至少3個重復
box2
RNA-seq文庫準備和測序過程中包擴：RNA打段，cDNA合成，接頭，PCR擴增，bar-coding，lane loading，這些過程可能會增加測序結(jié)果的偏好性。
外源參考轉(zhuǎn)錄組（exogenous reference transcripts，‘spike-ins’）可用來作為質(zhì)控以及文庫大小矯正。 若測序量較大，降低技術誤差：文庫準備時不同批次及l(fā)ane的樣本完全隨機，或每個樣本單獨進行barcoding，然后在多個illumina lane中，加入所有的樣本進行測序。
RNA-seq數(shù)據(jù)分析
數(shù)據(jù)分析的主要步驟：質(zhì)控，比對（分：有參考基因組、無參考基因組），獲得基因及轉(zhuǎn)錄本表達矩陣，基因差異分析。也討論可變剪接，轉(zhuǎn)錄本融合，小RNA表達，可視化工具。
1. 質(zhì)控檢測
1.1 原始序列
包括：序列質(zhì)量，GC含量，接頭，過高k-mers，重復reads。同一研究中重復度，k-mer或是GC含量應該已知，不一致性大于30%則剔除。常用FastQC。
準則：3‘末端序列質(zhì)量下降時需要刪除以增加比對率。FASTX-Toolkit 和Trimmomatic用來去除低質(zhì)量序列，去接頭，去掉低質(zhì)量堿基。
1.2 比對
最重要的是比對到 ：基因組或是轉(zhuǎn)錄組上的比對率。
人類基因組的比對率期望值是70-90%，會出現(xiàn)多個序列比對在有限的序列區(qū)稱之為“多重比對序列”（multi-mapping reads）；
轉(zhuǎn)錄組上的比對率較低，由于未注釋的轉(zhuǎn)錄本會被過濾且“多重比對序列”增加，由于同一個基因不同亞型共有外顯子區(qū)。
另一個參數(shù)：序列覆蓋度在外顯子和比對鏈上的均一性。3‘末端轉(zhuǎn)錄本聚集表明序列質(zhì)量差，GC含量可以顯示PCR偏好性，指控工具包括：Picard，RSeQC，Qualimap。
1.3 量化
樣本內(nèi)轉(zhuǎn)錄本定量后需檢測GC含量以及基因長度偏好性來居定是否進行矯正。確認無rRNA，smallRNA（R 包NOISeq或EDASeq 對計數(shù)進行質(zhì)控）。
1.4 重復
整個RNA-seq數(shù)據(jù)的可重復性檢測來排除批次效應（技術重復系數(shù)Spearman R2 > 0.9）。若相同條件下基因表達量有差異則主成分分析（principle component analysis，PCA）應聚在一支。
2. 轉(zhuǎn)錄本
有參分析時將序列比對到參考基因組或是轉(zhuǎn)錄組上獲得表達轉(zhuǎn)錄本。比對到轉(zhuǎn)錄組上會屏蔽新的未注釋的轉(zhuǎn)錄本，只對已知轉(zhuǎn)錄本進行定量分析。
無參時先組裝為長contigs后已contig作為表達轉(zhuǎn)錄組將reads比對上去進行定量分析，或者覆蓋度可用于對轉(zhuǎn)錄本進行定量。區(qū)別在于轉(zhuǎn)錄和定量同時完成還是順序完成。
2.1? 比對
有參比對分兩種：基因組比對和轉(zhuǎn)錄組比對（圖2a，b），一條或多條序列（multireads）都可以比對在特定的位點。
多比對由于重復序列或是有共同結(jié)構(gòu)域的旁系同源基因而導致，在比對在基因組上會產(chǎn)生顯著性的比對結(jié)果，在轉(zhuǎn)錄組為參考基因組時由于基因異構(gòu)體（insoform）含有共同的外顯子而更顯著，結(jié)果保留。在基因表達變化時轉(zhuǎn)錄本的發(fā)現(xiàn)和定量更加困難。
box3 比對到參考序列
比對到參考基因組可發(fā)現(xiàn)新的轉(zhuǎn)錄本和基因，需要gap或剪接map由于序列可能跨越剪接區(qū)。要發(fā)現(xiàn)正確的剪接區(qū)尤其是參考基因組中存在錯誤或差異或者無保守區(qū)和融合轉(zhuǎn)錄本。
?Tophat分兩步進行無剪接序列先比對到外顯子，沒比對的序列被分開比對來尋找外顯子區(qū)。比對時參數(shù)設置取決于文庫，錯配數(shù)，reads的長度和類型及測序長度。
2.1 轉(zhuǎn)錄本發(fā)現(xiàn)
新轉(zhuǎn)錄本的發(fā)現(xiàn)困難在于：Illumina讀長短(short reads)，難跨越剪接區(qū)不能直接的到轉(zhuǎn)錄本全長；轉(zhuǎn)錄本的起始和終止位點難確定。
PE reads(雙端測序）和該覆蓋率有利于發(fā)現(xiàn)低表達轉(zhuǎn)錄本，重復 有利于解決假陽性率（false-positive call）。
Cufflinks, iReckon , SLIDE和StringTie與注釋相結(jié)合將其加到可能的異構(gòu)體（insoform）中，Montebello將異構(gòu)體的發(fā)現(xiàn)與定量用似然法比對，Augustus可講轉(zhuǎn)錄組數(shù)據(jù)與編碼蛋白轉(zhuǎn)錄本注釋很好的結(jié)合，但非編碼轉(zhuǎn)錄本較差。
2.2 從頭合成轉(zhuǎn)錄本重建
無參序列組裝為轉(zhuǎn)錄本，SOAPdenovoTrans， Oases，Trans-ABySS或Trinity。無參轉(zhuǎn)錄組需PE reads和讀長較長的序列。無參分析在計算機分析時測序較深時要降低序列的數(shù)量。樣本間比較分析時，建議將多個樣本的所有序列都合并為一個輸入文件來的到一個穩(wěn)健的contigs（transcripts），然后比對回短序列進行表達量評估。
從頭組裝導致產(chǎn)生十或上百的contigs作為轉(zhuǎn)錄本片段，長測序技術如Bioscience 的SMRT提供讀長可以為多數(shù)基因提供完整的轉(zhuǎn)錄本。
3. 轉(zhuǎn)錄本定量
RNA-seq分析核心為基因和轉(zhuǎn)錄本的定量分析，基于比對到轉(zhuǎn)錄本上的數(shù)量。
最簡單的定量方法是用HTSeq-count或featureCounts累積原始數(shù)量。
基因水平定量使用GTF（genome transfer format ）文件，包含外顯子和基因，通常丟棄很多序列。原始序列數(shù)量不能用于比較樣本與樣本間的表達水平，由于受到轉(zhuǎn)錄本長度，總測序數(shù)以及測序偏好性的影響。
RPKM是樣本內(nèi)均一化方法，用于去除長度和樣本大小的影響（RPKM：reads per kilobases of exon model per millions reads）,FPKM(fragments per kilobase of exon model per million mapped read)與RPKs和TPM（transcripts per million）類似，都用于樣本內(nèi)歸一化，F(xiàn)PKM可以與TPM相互轉(zhuǎn)化。
樣本內(nèi)和樣本間的區(qū)分導致在文章中較為混亂。相同基因在樣本與樣本之間的表達量比較時，其長度不需要矯正。但同一個樣本內(nèi)對基因表達排序時必須的由于較長的序列回累積更多的reads。樣本之間Cufflinks得到基因長度顯著不同不同忽略。？？？？（備注：到底應該怎么辦？）
轉(zhuǎn)錄水平表達計算基于相同的轉(zhuǎn)錄本 共有 多數(shù)序列來進行計算。TopHat用最大期望值來對轉(zhuǎn)錄本的豐富度進行計算。Cufflinks使用GTF信息來發(fā)現(xiàn)轉(zhuǎn)錄本或只從比對序列提供從頭合成的轉(zhuǎn)錄本。從轉(zhuǎn)錄本比對量化表達包括SEM (RNA-Seqby Expectation Maximization)，eXpress，Sailfish，kallisto。
轉(zhuǎn)錄本中容許多比對reads以及將序列偏好性矯正后樣本內(nèi)均一化值輸出。RSEM使用最大期望值并返回TPM值。NURD為SE reads提供轉(zhuǎn)錄組表達評估，占內(nèi)存低。
4. 差異基因表達分析
差異表達分析需要將樣本與樣本之間的基因表達值進行比較。
RPKM，F(xiàn)PKM和TPM在樣本間進行比較時將測序深度進行歸一化，但當樣本有雜合性轉(zhuǎn)錄本分布即高且差異表達特性偏離count分布時結(jié)果較差。NOISeq R包包含大量的分析plots對每種情況進行合適的歸一化步驟。除樣本內(nèi)，樣本間差異，批次效應可能會產(chǎn)生影響，COMBAT或ARSyN可以剔除批次效應。
RNA-seq定量分析基于reads counts絕對或可能匹配到轉(zhuǎn)錄本上（波松或負二項分布）。絕對-離散概率分布-小片段樣本變異不同的表達包括在內(nèi)時不適合。
edgeR將原始輸入reads計數(shù)及可能的偏好性帶入數(shù)據(jù)模型，將歸一化和差異分析同時進行，類似的為DESeq2（負二項分布）。baySeq和EBSeq為貝葉斯法（負二項分布），不同實驗組內(nèi)的差異以及每組內(nèi)每個基因的后驗概率。
無參法NOISeq或SAMseq做最小假設，從真實數(shù)據(jù)中為理論分析做空值分布估算。最小生物學重復為3。不同算法顯著性的影響分析的結(jié)果，因此要表明參數(shù)設置，版本，以及考慮生物學重復。
5. 可變剪接分析：差異異構(gòu)體表達。
同一基因轉(zhuǎn)錄本異構(gòu)體的表達為可變剪接。分析方法分兩類：將異構(gòu)體表達評估與差異表達檢測結(jié)合來對總基因表達中每個異構(gòu)體占比的變化進行計算，兩步結(jié)合后第一步的不確定性考慮在內(nèi)：數(shù)據(jù)分析來尋找差異異構(gòu)體表達。
基于外顯子分析法（exon-based）省略異構(gòu)表達和可變剪接的信號檢測通過比較兩個比對樣本之間基因外顯子和連接區(qū)序列分布DEXseq和 DSGSeq （基因外顯子count），rMATS（連接區(qū)reads），rDiff（可變區(qū)域基因readscounts），DiffSplice用比對圖來發(fā)現(xiàn)可變剪接模型。優(yōu)點：exon或junction法可精準的發(fā)現(xiàn)單個可變剪接；exon-based適合特殊的外顯子和功能結(jié)構(gòu)域，不適合整個異構(gòu)體分析。
6. 可視化
可視化可以在reads水平（ReadXplorer)或在處理深度(read pileup), 未均一化 (總count) 或均一化后（基因組瀏覽器 UCSC browser，Integrative Genomics Viewer (IGV) , Genome Maps 或Savant，RNAseqViewer查看多個RNA-seq樣本，展示風豐富的外顯子，轉(zhuǎn)錄本，連接區(qū)，但比IGV慢。
7. 發(fā)現(xiàn)融合基因
染色體重排產(chǎn)生融合基因與新異構(gòu)體基因鑒定方法類似，但跨度更大。
假的融合基因由于多態(tài)性，同源異記序列錯誤而導致的比對錯誤而產(chǎn)生。過濾多態(tài)性豐富和同源配對基因，也過濾掉不可能參與基因融合的高表達基因如rRNA。另外野生型中在近融合區(qū)存在低頻的二體可能以為著高表達基因的錯配。
若得到正確的chimeric，下一步是得到有生物學功能的融合基因。當融合出現(xiàn)在對照數(shù)據(jù)中時可能會被過濾，當無對照數(shù)據(jù)時，大量不相關聯(lián)的數(shù)據(jù)庫同時出現(xiàn)，且過濾后出現(xiàn)真正的融合時則表明artifacts。
8. Small RNAs
sRNA通常包含18-34堿基，有miRNA, siRNA（小干擾RNA），PIWI-交互RNAs（PIWI-interacting RNA，piRNAs）以及其他類型的調(diào)控分子。由于其復雜度小測序通常為2-10 百萬reads，于RNA-seq分析方法有不同。去接頭：動物中，長度為22和23bp，植物種21和24bp。
sRNA需用Bowtie2，STAR，Burrows-Wheeler Aligner (BWA)比對到參考基因組上。未比對上的潛在的重復序列需要剔除。每個基因組上通常容許5-20個不同的mapping。保證無mRNA降解污染。
下一步的分析步驟包括與已知sRNA比較以及從頭發(fā)現(xiàn)sRNAs。miRDeep用于動物分析，miRDeep-P用于植物，or the trans-acting siRNA預測工具 UEA sRNA Workbench。miRTools 2.0，ShortStack和 iMir能為sRNA文庫綜合注釋，并鑒定多種 sRNAs分類
9. RNA-seq功能注釋
標準轉(zhuǎn)錄組分析最后一步：差異表達基因（differentially expressed genes，DEGs）的功能和通路分析。
兩個主要的方法：比較差異表達基因與剩余基因組，基因富集分析（gene set enrichment analysis, GSEA）基于差異表達轉(zhuǎn)錄本排序。
功能分析需要對研究的材料有可用及豐富的功能注釋。
Gene Ontology，Bioconductor，DAVID或Babelomics包含多數(shù)模式物種的注釋數(shù)據(jù)。
從頭組裝所得到的新轉(zhuǎn)錄本缺乏注釋信息，編碼蛋白注釋可以基于序列相似性用旁系同源功能注釋（SwissProt），以及保守蛋白結(jié)構(gòu)域用Pfam和InterPro。一般有50-80%的轉(zhuǎn)錄本可以被注釋。缺少編碼蛋白的轉(zhuǎn)錄本為長非編碼RNA（long non-coding RNA）,相似性注釋可用于短非編碼RNA，而對于長非編碼RNA還沒有相應的注釋。
與其他數(shù)據(jù)類型相結(jié)合
1. 與DNA測序結(jié)合
RNA與DNA測序相結(jié)合可用來發(fā)現(xiàn)單堿基多態(tài)性（single nucleotide polymorphism, SNP）RNA-編輯，表達數(shù)量性狀位點（expression quantitative trait loci，eQTL）。
經(jīng)典的eQTL研究中，同一類型的組織基因型和轉(zhuǎn)錄組測序數(shù)量大于50，然后檢測基因型和表達水平的關系，用來解釋復雜性狀基因偏好性。大量的eQTL研究表明基因變異影響多數(shù)基因的表達
RNA-seq在檢測eQTL方面有兩個優(yōu)勢：發(fā)現(xiàn)影響轉(zhuǎn)錄過程的變異；雜合性SNP可以分布比對到父本和母本上，對個體內(nèi)等位基因特異性表達進行定量分析。
2. DNA甲基化
DEGs和甲基化模型的相關分析，然而通過線性相關性，貝葉斯相關性，邏輯相關性模型得出兩者的相關性較低。
網(wǎng)絡互作分析RNA-seq與DNA甲基化之間的關系，發(fā)現(xiàn)一個或多個基因有差異表達和差異甲基化的協(xié)同性。
3. 染色質(zhì)特征
RNA-seq與轉(zhuǎn)錄元件（transcription factor，TF）染色質(zhì)免疫沉降測序（ChIP-seq）數(shù)據(jù)用來剔除ChIP-seq中的假陽性和表明目的基因上TF的激活或抑制。
ChIP-seq數(shù)據(jù)組蛋白修飾用來表示表觀修飾對基因表達量的改變。DNase-seq可用于DNA結(jié)合因子的基因組印記，與基因的表達相結(jié)合可用于研究轉(zhuǎn)錄網(wǎng)絡活性。
4. MicroRNAs
兩種數(shù)據(jù)相結(jié)合可能用來解釋轉(zhuǎn)錄穩(wěn)定水平miRNA的調(diào)控作用。
5. 蛋白組及代謝組
與蛋白組數(shù)據(jù)結(jié)合有爭議由于兩者的相關性低（～0.4）。然而仍可以用來發(fā)現(xiàn)新異構(gòu)體。用RNA-seq預測未報道的肽鍵或事轉(zhuǎn)錄后編輯。與代謝組結(jié)合可用來發(fā)現(xiàn)基因表達和代謝水平的調(diào)控通路。
6.多數(shù)據(jù)類型聯(lián)合及可視化
蛋白–蛋白, DNA–蛋白, miRNA–mRNA 互作網(wǎng)絡來發(fā)現(xiàn)miRNA–基因調(diào)控模型。
展望
目前轉(zhuǎn)錄組分析主要方面：少量的供試材料 和長序列中 更好的發(fā)現(xiàn)轉(zhuǎn)錄本
1. 單細胞轉(zhuǎn)錄組（single-cell RNA-seq）
前沿和火熱的研究區(qū)域。Smart-seq和Smart-seq2只需極少量的供試材料，可通過單個細胞的擴增得到?？捎糜诎l(fā)現(xiàn)組織中新的未分類的細胞類型。一類單細胞文庫與細胞群相比，發(fā)現(xiàn)多細胞亞群與表達基因相結(jié)合。
少量的供試材料以及PCR擴增限制了測序深度，因而一般測序少于1百萬reads。scRNA-seq測序深度增加可能有利于同源特異性表達基因的挖掘，但表達量的增加鮮有提高。
scRNA含有3000-8000個表達基因，加入?yún)⒖嫁D(zhuǎn)錄本以及特異性分子標記（uniqe molecule identifiers，UMI）有利于克服偏好性擴增并提高基因定量。
scRNA-seq比對在轉(zhuǎn)錄組參考基因組上不能發(fā)現(xiàn)新的基因，若研究目的未基因表達量則用轉(zhuǎn)錄組未參考基因組來減少工作量。
2.長測序
短序列限制性在于不能精準的沖否完整的轉(zhuǎn)錄本。Pacific-Bioscience（PacBio）SMRT和Oxford Nanopore獲得長序列。PacBio在cDNA分子上加接頭形成一個環(huán)形結(jié)構(gòu)，此單鏈用來多次測序。Nanopore GridION可直接用RNA合成酶和RNA特異性堿基進行測序。Moleculo技術準備文庫時復合和限制DNA分子長度，將這些特定長度的鏈分開標記然后重新融合測序。 PacBio最常見。
缺點：測序錯誤率高，不能用于從頭合成需要參考基因組；SMRT細胞數(shù)量較低阻礙了轉(zhuǎn)錄本定量分析。
