轉錄組分析工具大比拼
文獻閱讀與翻譯 - Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis (Nature Communications DOI: 10.1038/s41467-017-00050-4)
摘要
RNA-sequencing (RNA-seq)是轉錄組研究的重要技術。自從RNA-seq技術問世以來，已經開發了大量的分析工具。雖然有有一些研究嘗試評估最新的分析工具，但是還沒有綜合評價分析流程評估單個或組合分析工具的優劣。在這篇文章中，研究人員深入探索了大量的RNA-seq分析流程。這些流程不僅包括表達分析技術，也包含了RNA variant-calling、RNA編輯和RNA融合檢測技術。具體來說，研究人員檢查了短讀長和長讀長的RNA-seq技術，包括了39個分析工具的120種組合，15個樣品的490種分析，樣品包括生殖細胞、癌細胞和干細胞數據集，最后報道了分析工具的性能，提出了一個綜合RNA-seq分析手冊：RNACocktail。這個手冊包含獲得高准確性的分析流程。通過與實時定量PCR數據的比較驗證，表明作者提出的流程可以幫助研究着獲得更多准確的生物相關的研究信息。
簡介
高通量二代測序(NGS)RNA-seq將轉錄組分析引入一個新的時代。RNA-seq需要各種分析流程以滿足測序技術、樣品類型、基因組獲取以及計算機資源的需求。不同的分析流程具有顯著不同的分析准確度、速度和代價。因此，在受到代價和性能限制的條件下研究RNA-seq分析每一步使用哪個或哪些分析工具獲得最高的准確性顯得尤為重要。通常，整體最佳的流程用於特定樣品分析可能是次優的，找出最優的分析流程更具挑戰性，這需要分析大量不同的數據集。
有很多研究比較了不同RNA-seq分析工具的性能。但是，這些研究主要關注RNA-seq分析的某一步，或者局限於比對和定量。因此，綜合系統分析有助於最大化了解RNA-seq數據。為了解決以前研究的限制，研究人員深入調查了RNA-seq分析的所有主要步驟，評價了不同步驟下分析工具組合的准確性、效率和一致性，提出了一個綜合的RNA-seq分析流程手冊。他們認為RNACocktail分析流程可以獲得高准確性。RNACocktail分析流程在檢測不同樣品生物相關的差異表達基因和臨床重要的轉錄本得到了進一步的驗證。RNACocktail分析流程是開源的並且可以在網站http://bioinform.github.io/rnacocktail/免費下載使用。
結果
數據集
為了進行綜合評估，研究人員分析了不同種類的RNA-seq數據, 包括15個Illumina和Pacific Biosciences (PacBio)數據集。這些數據具體來自正常人類樣品NA12878、人類MCF-7乳腺癌細胞、H1人類胚胎干細胞(hESC)以及測序質量控制聯盟(SEQC)。
RNA-seq分析流程手冊
下圖展示了研究人員提出的綜合分析RNA-seq數據的流程手冊。在流程中的每一步都列出了常用分析方法 (使用的工具名和版本號請看補充表2)。下面會詳細闡述每一個步驟。
RNACocktail分析流程手冊。RNACocktail是一個綜合的RNA-seq數據分析手冊。本圖概括了RNA-seq分析關鍵步驟中廣泛使用的分析工具。研究人員可以使用這個圖展示的工作流程分析RNA-seq數據。
通過短讀長進行轉錄本檢測
識別表達的轉錄本往往是RNA-seq分析的第一步。通常先將reads比對回參考基因組(或者參考轉錄組)，然后根據reads比對結果拼裝轉錄本。比對到參考基因組可以檢測新的轉錄本，但需要消耗大量計算資源。而比對回參考轉錄組要容易得多，但檢測不到新的轉錄本。另外如果沒有可靠的參考基因組或者轉錄組存在，也可以從頭組裝轉錄本。
基於參考的轉錄本識別
比對和結合點預測
RNA-seq reads拼接比對到參考基因組要根據外顯子-內含子邊界拼接reads。下面研究人員使用短讀長的IlluminahESC、NA12878、SEQC、100 bp和300 bp長的MCF7數據集評價了TopHat、STAR和HISAT2的性能，如圖所示。
不同比對方案的性能比較。a.不同方案檢測到的拼接位點之間的overlap以及與dbEST數據庫可靠的拼接點相比的驗證率。可靠的EST拼接點至少包括兩個EST支持的拼接點。圓的大小反映了每個方案檢測到的拼接點數量。圖中展示了每個工具檢測到的拼接點數量和驗證率(位於小括號中)。b.read比對分析。左圖展示測序片段比對狀態的分布(展示了NA12878、MCF7和SEQC樣品的雙端reads比對狀態，對於hESC樣品，展示了單端reads的比對狀態，藍色表示獨一無二比對到基因組，橙色表示基因組上有多個比對位置，紅色表示沒有比對到基因組上)。中圖展示了比對回基因組的片段上鹼基被軟件去除 (soft-clip)的數量分布。右圖展示了比對回基因組的片段錯配鹼基的數量分布。
HISAT2在所有樣品中擁有最高的拼接點驗證率，但是其預測的拼接點數量小於TopHat和STAR，如圖a所示。STAR具有最高比例的獨一無二比對到基因組上的reads，尤其是300讀長的MCF7樣品(b)。與TopHat和HISAT2不同，STAR會把雙端reads比對到基因組，否則移除雙端reads，以避免單端reads的比對。另一方面，STAR獲得了較低質量的比對，具有更多的soft-clipped比對和錯配鹼基(b)。TopHat禁止截斷reads(b)。對長讀長樣品MCF7-300和單端read樣品hESC的分析結果表明，與TopHat和HISAT2相比，STAR具有更高的容忍性，接受鹼基錯配和soft-clipping以將更多的reads比對回參考基因組。在比對速度方面，HISAT2比STAR快2.5倍，比TopHat快大約100倍。
基於比對的轉錄組組裝
拼接比對之后，表達的轉錄本集合通過轉錄組組裝來識別。下面研究人員關注兩個廣泛使用的基於比對的轉錄組發掘工具Cufflinks和StringTie。作為這兩個組裝工具的輸入，我們使用了上述討論的三個比對工具。我們以Ensembl參考轉錄組注釋為指導。
除了短讀長轉錄本預測工具外，也研究了轉錄本檢測和預測工具(IDP)。IDP通過短讀長比對以幫助長讀長的轉錄本檢測，是一個混合的轉錄本預測工具。IDP通過GMAP和STARlong的長讀段比對和TopHat、STAR和HISAT2的短讀段比對進行評價。此外，也分析了長讀段轉錄本預測工具Iso-Seq。
研究人員將預測的轉錄本與GENCODE v19中的參考轉錄組注釋比較，不存在於參考轉錄組中的轉錄本認為是假陽性(FPs)。Cufflinks和StringTie預測了許多單個外顯子轉錄本，大部分是假陽性(圖a)。StringTie預測的轉錄本數量比Cufflinks多50–200%。IDP在不同的樣品中預測最少的外顯子，它預測不出單外顯子的基因。對於多個外顯子轉錄本來說，IDP預測的轉錄本數量和Cufflinks相似(圖a)，而且，其預測的外顯子數量分布與GENCODE更相似。Iso-Seq算法預測的94%的單個外顯子轉錄本和77%的多個外顯子轉錄本在GENCODE上是沒有的，這反映了其組裝准確度不高，但是具有更高的敏感度檢測新轉錄本。
不同轉錄組重構方案的性能比較。a.不同的組裝算法預測的轉錄本外顯子數量分布。標簽標明了組裝工具、長讀段比對工具(IDP)、短讀段比對工具，之間用”-“分開。b.基因和轉錄本水平上不同轉錄組重構方法的敏感度和准確度。GENCODE參考轉錄組注釋用作金標准。使用短讀長和長讀長轉錄本預測的組合方法(用”+”標記)稍微提升了短讀長轉錄本預測方案的性能。
對於MCF7-300樣品，有STAR工具的組合預測出更多的轉錄本(圖a)。IDP與長讀段比對工具GMAP和短讀段比對工具HISAT2聯合使用會預測出更多轉錄本(圖a)。
與短讀段組裝工具不同，IDP會檢測到一個基因的多個轉錄本(圖)。與Cufflinks相比，StringTie預測更多的基因，而且每個基因超過5個轉錄本(圖)。StringTie的預測結果最好地匹配了GENCODE上每個基因轉錄本數量的分布。
不同轉錄組重構算法預測的每個基因轉錄本數量歸一化后的分布比較。標簽標明了組裝工具、長讀段比對工具(IDP)、短讀段比對工具，之間用”-“分開。使用短讀長和長讀長轉錄本預測的組合方法(用”+”標記)稍微提升了短讀長轉錄本預測方案的性能。
IDP在基因水平評價中獲得了最好的准確度和靈敏度(圖b)。Cufflinks比StringTie更敏感和准確。對於MCF7-300樣品，組合工具之間性能差異較大。與StringTie組合使用的TopHat和HISAT2表現比STAR更佳。Iso-Seq敏感度最低，准確度在IDP和Cufflinks、StringTie之間。
在轉錄本水平上，IDP在准確度方面表現最佳(圖b)。IDP的預測局限於多個外顯子的轉錄本，它的靈敏度稍高於Cufflinks，低於StringTie。在短讀長組裝工具中，在轉錄本水平，StringTie的靈敏度和准確度都高於Cufflinks，分別高出25%和11%(圖b)。
不同的工具預測Ensembl沒有注釋但GENCODE v19注釋的3681個新的轉錄本，結果表明StringTie預測到了最多的轉錄本，比Cufflinks多2.5倍，比IDP多6.5倍(圖)。
不同轉錄組重構算法預測新的轉錄本性能比較。新的轉錄本是GENCODE (v19)注釋的但在Ensembl中沒有注釋的多外顯子轉錄本。標簽標明了組裝工具、長讀段比對工具(IDP)、短讀段比對工具，之間用”-“分開。使用短讀長和長讀長轉錄本預測的組合方法(用”+”標記)稍微提升了短讀長轉錄本預測方案的性能。
StringTie是運行最快的工具，完成組裝比Cufflinks快約60倍，比IDP快約50倍(默認輸入是錯誤校正過的數據和比對后的數據)，如下表所示。
轉錄本定量
基於比對的轉錄本定量
經典的比對分析是將reads比對回參考基因組或者參考轉錄組，之后估計轉錄本豐度。如果研究目的是測量已知的和新的轉錄本豐度，比對回參考基因組后使用Cufflinks和StringTie進行組裝和豐度估計。如果使用參考轉錄組是發現不了新的轉錄本的，reads可以直接比對到轉錄組之后使用RSEM和eXpress進行豐度估計。
不經過比對的轉錄本定量
不經過比對的定量方法直接將reads分配給轉錄本，這與拼接比對方法相比需要更少的計算資源。Sailfish、Salmon、quasi-mapping和kallisto四種工具可以解決哪個轉錄本可以生成哪個read，或者尋找部分比對回轉錄組的reads。
下面我們比較了基於基因組比對的工具StringTie和Cufflinks(使用不同的比對工具)、基於轉錄組的比對工具eXpress和Salmon-Aln、不經過比對的工具kallisto、Sailfish (with quasi-mapping)、Salmon-SMEM和Salmon-Quasi、基於長讀段的技術IDP(使用不同的短讀長和長讀長比對工具)的性能。
基於log尺度表達水平之間的Spearman等級相關分析的不同定量方法聚類分析表明采用相似方法的定量工具聚在一起(圖a)。不經過比對的工具與StringTie聚的更近。Salmon-SMEM與基於轉錄組比對的工具eXpress和Salmon-Aln聚在一起，Salmon-SMEM運行速度更快。
基於短讀段的技術中，兩個不經過比對的工具kallisto和Salmon-SMEM對MCF7-100和MCF7-300樣品具有最一致的預測結果(圖b,c)。對於MCF7-100和MCF7-300樣品，IDP表現出較高的一致性(圖b),尤其是刪除低表達的基因后一致性更高(圖c)。
通常，StringTie與高效的比對工具如HISAT2組合是最高效的基於比對方法，但是不進行比對的工具也非常有效(運行速度比比對工具快一個數量級)。之前的研究表明與比對方法相比，定量方法對豐度估計影響更大，我們的研究結果也證實了這一點。
轉錄本豐度估計工具的性能比較。a.基於NA12878樣品log表達的spearman等級相關分析的不同方案的聚類分析比較。b. MCF7-100和MCF7-300樣品之間log2-fold表達變化的分布。每種方法中，虛線表示分布的平均值，點線表示四分位數。c.用不同閾值去掉底表達轉錄本后在 MCF7-100和MCF7-300樣品之間表達差異百分比。
差異表達
不同樣品和條件下差異表達基因的識別是RNA-seq分析的重要目標。有多種方法檢測差異表達基因，包括基於計數技術的DESeq、limma和edgeR、基於組裝技術的Cuffdif和Ballgown、不經過比對定量進行差異分析的sleuth。
首先，通過從SEQC樣品(SEQC-A vs. SEQC-B)1001個基因中檢測差異表達的基因與反轉錄定量PCR(qRT-PCR)測量的表達變化進行比較來評價工具的性能(圖)。與其他工具相比，DESeq2表現最佳。sleuth、edgeR和limma性能較差。Cuffdiff和Ballgown的准確度沒有基於計數的工具准確度高。對於AUC-30的測量，edgeR表現最佳。
我們比較了不同的工具在預測SEQC數據集中的92個外部RNA控制聯盟 spike-in 基因的准確度。對於Spearman等級相關系數和RMSD測量，DESeq2獲得最佳的性能。對於AUC-30測量，Cufflinks和Ballgown表現最佳。基於計數的工具比基於組裝的工具更高效。不經過比對技術如Salmon和kallisto能夠獲得高質量的預測結果。
差異基因表達分析工具不同性能的比較。a.qPCR測量的基因的Spearman等級相關系數、根均方差(RMSD)和AUC-30得分。b.qRT-PCR測量的基因(左)和ERCC基因(右)。
運行時間比較
不經過比的方法運行速度最快，StringTie-HISAT2組合是基於比對方法中運行最快的，但是比不經過比的方法運行速度慢一個數量級。Cufflinks-TopHat組合與基於長讀段的方法比StringTie-HISAT2組合慢兩個數量級(圖)。
不同RNA-seq分析方法運行的總CPU時間。
高准確度的分析流程
沒有一款分析工具在各種條件下表現最佳，基於RNA-seq分析工具的整體性能，我們提出了RNACocktail分析流程，這個流程中每一步都是由高准確度的分析工具組成，可以用於一般目的RNA-seq分析(圖)。當前廣泛使用的Cufflinks-TopHat流程在准確度和計算代價上的表現不如基於比對的工具如StringTie-HISAT2組合和不經過比對工具如Salmon-SMEM。
當前RNACocktail計算流程。這個分析流程每一步由高准確度的工具組成，可以用於一般目的的RNA-seq分析。
為了評價我們提出的分析流程在功能富集分析(看到了富集分析 )的性能，相對於正常樣品NA12878，MCF7和hESC樣品高表達基因使用我們的分析流程和Tuxedo方法進行識別。使用ToppGene工具對前10個高表達的基因進行功能富集分析。對於MCF7-100和MCF7-300樣品，Cufflinks-TopHat預測的基因集沒有與MCF7或者乳腺癌相關的基因富集在一起，而StringTie-HISAT2和Salmon-SMEM預測的前10個表達基因與許多MCF7和乳腺癌細胞系相關的基因集高度富集在一起。對於hESC樣品，也得到了相同的結果。
我們的分析流程是經過廣泛而深入的評價后提出的。這個流程也更具綜合性，包含了其他分析流程如Galaxy和Grape缺少的從頭組裝、變異位點檢測、RNA編輯位點檢測、長讀段RNA-seq分析等。
討論
通過綜合分析RNA-seq分析流程中不同步驟的工具性能發現不同的分析工具和方法對分析結果的准確度和分析時間影響巨大。HISAT2表現出最快的速度和最准確的拼接比對，但是沒有STAR的敏感度高。StringTie在速度和准確度上都優於Cufflinks。長讀段方法如IDP和Iso-Seq會識別許多短讀段技術沒有識別到的多外顯子轉錄本，但是會丟失一些單外顯子轉錄本。通常，在從頭組裝工具中，Oases表現最佳。不經過比對的工具如Salmon-SMEM和kallisto獲得了最好的一致性和最高准確度，因此，如果目標不是發現新的轉錄本，如Salmon-SMEM和kallisto可以作為准確而快速的解決方案。DESeq2和edgeR與不經過比對的工具聯用可以獲得高准確度的差異表達分析結果。GATK是一個准確的變異位點檢測工具，可以與不同的比對工具聯用。當與HISAT2或者STAR比對工具聯用時，GIREMI可以不依賴基因組准確預測RNA編輯位點。長讀段方法如IDP-fusion可以准確預測RNA融合，而短讀段方法如FusionCatcher或者SOAPfuse具有更高的靈敏度。通常情況下，整體最好的分析流程對於特定的數據集特定的研究目的來說可能是次優的。比如，對於比對和轉錄組構建，HISAT2``-StringTie組合具有更高的准確度和更快的速度。但是對於MCF7-300樣品來講，STAR- StringTie組合具有更高的靈敏度(圖a)。
對hESC和MCF7樣品中高表達基因的詳細分析表明新開發的工具表現比標准的Tuxedo手冊好。比如，六個人類胚胎干細胞中常見的上調基因集中的89個基因列表中，StringTie-HISAT2和Salmon-SMEM分別預測了位於89個基因列表中的10個基因的6個和4個，而Cufflinks-TopHat預測的基因都不在上述基因列表中。 StringTie-HISAT2發現的6個高表達基因是TDGF1、CRABP1、SFRP2、GJA1、GAL、LIN28A，這些基因在胚胎發育過程中發揮重要作用。