知識的學習沒有一蹴而就缠捌，沒有捷近，扎實的學習是唯一的捷近苍姜。
一篇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）吝沫。讀長較長的序列及雙端序列更有利于注釋信息較差的轉(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ì)控以及文庫大小矯正。 若測序量較大踏烙，降低技術(shù)誤差：文庫準備時不同批次及l(fā)ane的樣本完全隨機师骗，或每個樣本單獨進行barcoding，然后在多個illumina lane中讨惩，加入所有的樣本進行測序辟癌。
RNA-seq數(shù)據(jù)分析
數(shù)據(jù)分析的主要步驟：指控，比對（有參考基因組及無參考基因組）荐捻，獲得基因及轉(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ù)重復系數(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)體含有共同的外顯子而更顯著摔蓝，結(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讀長短坚芜，難跨越剪接區(qū)不能直接的到轉(zhuǎn)錄本全長览芳；轉(zhuǎn)錄本的起始和終止位點難確定。PE reads和該覆蓋率有利于低表達轉(zhuǎn)錄本的發(fā)現(xiàn)鸿竖，重復有利于解決假陽性率（false-positive call）沧竟。Cufflinks, iReckon , SLIDE和StringTie與注釋相結(jié)合將其加到可能的異構(gòu)體中铸敏，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)錄本片段，長測序技術(shù)如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. 可變剪接分析
同一基因轉(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ù)時比被，大量不相關(guān)聯(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嗡载，然后檢測基因型和表達水平的關(guān)系，用來解釋復雜性狀基因偏好性仍稀。大量的eQTL研究表明基因變異影響多數(shù)基因的表達
RNA-seq在檢測eQTL方面有兩個優(yōu)勢：發(fā)現(xiàn)影響轉(zhuǎn)錄過程的變異洼滚；雜合性SNP可以分布比對到父本和母本上，對個體內(nèi)等位基因特異性表達進行定量分析技潘。
2. DNA甲基化
DEGs和甲基化模型的相關(guān)分析遥巴，然而通過線性相關(guān)性，貝葉斯相關(guān)性享幽，邏輯相關(guān)性模型得出兩者的相關(guān)性較低铲掐。網(wǎng)絡互作分析RNA-seq與DNA甲基化之間的關(guān)系，發(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é)合有爭議由于兩者的相關(guān)性低（～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只需極少量的供試材料户誓，可通過單個細胞的擴增得到∧幌溃可用于發(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技術(shù)準備文庫時復合和限制DNA分子長度颗祝，將這些特定長度的鏈分開標記然后重新融合測序浊闪。 PacBio最常見。
缺點：測序錯誤率高螺戳，不能用于從頭合成需要參考基因組搁宾；SMRT細胞數(shù)量較低阻礙了轉(zhuǎn)錄本定量分析。
