39個工具，120種組合深度評估 (轉錄組分析工具哪家強)

2019 年 10 月 7 日
筆記

RNA-seq分析工具知多少

RNA-seq是研究轉錄組應用最廣泛，也最重要的技術之一。RNAseq其分析內容包括序列比對、轉錄本拼裝、表達定量、差異分析、融合基因檢測、可變剪接、RNA編輯和突變檢測等，具體流程和常用工具如下圖所示。通常的分析不一定需要走完全部流程，按需進行，某些步驟可以跳過、簡化等。

RNA-seq分析工具最優組合

Nature Communication上一篇文章 Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis對15個樣品 (正常樣品、癌細胞和幹細胞，短讀長和長讀長)的轉錄組數據利用39個分析工具，120種常見組合方式進行的490次深入分析, 並以測序質量控制聯盟(SEQC)的qPCR檢測結果做為正對照，總結出一套普適性流程，如下。

通過綜合分析RNA-seq分析流程中不同步驟的工具性能發現不同的分析工具和方法對分析結果的準確度和分析時間影響巨大。

HISAT2表現出最快的速度和最準確的拼接比對，但是沒有STAR的敏感度高。StringTie在速度和準確度上都優於Cufflinks。

長讀段方法如IDP和Iso-Seq會識別許多短讀段技術沒有識別到的多外顯子轉錄本，但是會丟失一些單外顯子轉錄本。

不經過比對的工具如Salmon-SMEM和kallisto獲得了最好的一致性和最高準確度，因此，如果目標不是發現新的轉錄本，如Salmon-SMEM和kallisto可以作為準確而快速的解決方案。

DESeq2和edgeR與不經過比對的工具聯用可以獲得高準確度的差異表達分析結果。

通常情況下，整體最好的分析流程對於特定的數據集特定的研究目的來說可能是次優的。比如，對於比對和轉錄組構建，HISAT2-StringTie組合具有更高的準確度和更快的速度。但是對於MCF7-300樣品來講，STAR– StringTie組合具有更高的靈敏度。

下面將詳細闡述每部分的評估。

序列比對質量大比拼

STAR具有最高比例的在基因組上有唯一比對位置的reads，尤其是對讀長為300 nt的MCF7樣品也有最高的比對率。

與TopHat和HISAT2不同，STAR只保留雙端reads都比對到基因組的序列，但對低質量的比對 (允許更多的錯配鹼基和soft-clip事件) 容忍度高。這一點在長reads (MCF7-300)樣品中的體現更為明顯。TopHat則不允許soft-clip事件。

soft-clip事件: 即reads末端存在低質量鹼基或接頭導致比對不上的, STAR會自動嘗試截去未比對部分，只保留比對上的部分。

在比對速度方面，HISAT2比STAR快2.5倍，比TopHat快大約100倍。(後續會推出柱狀圖的一步畫法)

Exon-exon junction位點評估

轉錄組reads比對不同於基因組reads比對(如ChIP-seq、WES等)的地方在於比對的reads可能來源於2個被內含子隔開的外顯子區域，導致reads一端比對在第一個外顯子的後面部分，另一端比對在第二個外顯子的前面部分，從而形成exon-exon junction (剪接點)。這些reads又稱為junction reads，對轉錄本的拼接、鑒定和差異分析具有重要的意義。

下面的維恩圖展示了不同比對軟件檢測到的共有和特有的剪接位點的比較 (整數代表每個軟件檢測到的剪接位點的數目，百分數代表每個集合的splice junction被驗證的比例)。可信的剪接點定義為dbEST數據庫中有至少2個表達序列標籤(EST)支持的位點, 做為正對照。

HISAT2在所有樣品中擁有最高的剪接點驗證率 (80%-91%)，TopHat其次 (54%-74%)，STAR最低 (42%-54%)。但是HISAT2預測的剪接點的數量最少，約為TopHat的60%和STAR的50%。

韋恩圖繪製看 R語言學習 – 韋恩圖輕鬆繪製各種Venn圖

基於參考基因組的轉錄組組裝

對於二代測序數據，Cufflinks和StringTie是應用最廣泛的兩個基於比對結果的轉錄本拼裝工具。(比對軟件STAR,HISAT2和TopHat)

對於三代測序數據，PacBio的流程中默認使用軟件Iso-Seq。

二代和三代測序數據雜交拼裝，使用的是IDP (Isoform Detection and Prediction)。(比對軟件GMAP、STAR long)

轉錄本拼裝質量評估的依據是GENCODE v19的參考轉錄組注釋，不存在於這個集合的轉錄本視為假陽性。

每個轉錄本中包含的外顯子的數目是轉錄本拼裝質量的一個評價標準, 通常單外顯子轉錄本可信度最差。Cufflinks的單外顯子轉錄本的數目佔到30%左右，StringTie在15%左右。這些單外顯子轉錄本大約90%為假陽性 (數字為目測附圖的估計)。StringTie拼裝獲得的轉錄本的數目約為Cufflinks的兩倍，其外顯子數目的分佈與GENCODE v19較為相似。

IDP組裝出的都是多外顯子轉錄本，整體數目與Cufflinks排除單外顯子轉錄本後相近，但外顯子數目的分佈與GENCODE v19更一致。與之相比，Iso-Seq的假陽性率較高，但敏感性更強。

堆積柱狀圖的畫法可以參考：是Excel的圖，不！是R的圖

對於基因水平的組裝，IDP的的準確性和靈敏性都是最好的。Cufflinks比StringTie更為準確和靈敏。對於MCF3-300樣品來講，含有STAR的組合拼裝出更多的轉錄本，但拼裝準確性和靈敏性都略低於基於TopHat和HISAT2的結果。IDP和StringTie拼裝出更多的多轉錄本基因。(下圖左)

對於轉錄本水平的組裝，IDP的準確性比其它技術高20%，但其敏感性低於StringTie，高於Cufflinks。相比喻Cufflinks，StringTie轉錄本水平的組裝精確性和敏感性高11%和25%。在預測新的轉錄本上 (ENSEMBL沒有注釋但GENCODE v19有的3681個轉錄本)，StringTie得到的最多，約是Cufflinks和IDP的2.5和6.5倍。(下圖右)

另外StringTie的速度是Cufflinks的50倍，IDP的60倍。

散點圖繪製 R語言學習 – 散點圖繪製

表達定量

傳統的表達分析是將reads比對回參考基因組或者參考轉錄組，然後估計轉錄本丰度。如果研究目的是關注已知的和新的轉錄本的丰度，比對回參考基因組後使用Cufflinks和StringTie進行組裝，然後評估表達丰度。如果只想定量已經注釋的基因，直接比對到參考轉錄組，再使用RSEM和eXpress進行丰度估計。

現在基於轉錄本的定量還有一種方式是不經過比對直接判斷read來源於哪個轉錄本，這比拼接比對定量需要更少的計算資源。Sailfish、Salmon、quasi-mapping和kallisto四種工具是這一計算方式的代表。

對樣品NA12878採用不同方法定量得到的基因表達譜進行log轉換後的Spearman秩和相關性分析表明採用相似方法的定量工具獲得的表達圖譜更相近。Cufflinks的定量結果與其他工具相關性最差，不足0.4. 不需要比對直接定量的工具與StringTie計算的結果更相近 (相關係數0.6-0.8)。Salmon-SMEM與基於轉錄組比對的工具eXpress和Salmon-Aln聚在一起，但Salmon-SMEM運行速度更快。

R語言學習 – 熱圖簡化 R語言學習 – 熱圖美化 R語言學習 – 熱圖繪製 (heatmap)

對於同一個樣品不同測序讀長的數據 (MCF7-100和MCF7-300)的比較分析可以反應比對工具定量的穩定性。兩個不依賴於比對的定量工具kallisto和Salmon-SMEM具有最一致的定量結果。Cufflinks-TopHat組合的結果在基於比對的定量工具組合中表現最優。整體看，基於STAR的比對結果，定量穩定性低於基於HISAT2的比對。

綜上，不基於比對的定量結果效率和穩定性最高。StringTie與HISAT2的組合是基於比對的定量工具中性能最好的, 但也要比不基於比對的工具慢一個數量級。

此圖為小提琴圖 (R語言學習 – 箱線圖（小提琴圖、抖動圖、區域散點圖）R語言學習 – 箱線圖一步法)，展示了數據分佈的密度，越胖的地方數據越集中。縱向表示兩個樣品基因表達變化的幅度，橫向表示變化幅度的集中度，數據越集中於y=0，定量一致性越好。

此圖為線圖（R語言學習 – 線圖一步法 R語言學習 – 線圖繪製），展示的是逐步移除最低表達的部分轉錄本後定量的一致性。線越接近X軸表明一致性越好。

差異表達基因鑒定

不同樣品和條件下差異表達基因的識別是RNA-seq分析的重要目標。有多種方法鑒定差異表達基因，包括基於計數 (reads count)的DESeq、limma和edgeR、基於組裝技術的Cuffdiff和Ballgown、不經過比對定量進行差異分析的sleuth。

SEQC樣品 (SEQC-A vs SEQC-B, SEQC-C vs SEQC-D)中1001個有qRT-PCR定量過的基因作為對照評價工具的性能。

DESeq2在所有組合中表現最佳（DESeq2差異基因分析和批次效應移除），sleuth、edgeR和limma(典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集)略微次之，但差別不大。

Cuffdiff和Ballgown的準確度沒有基於計數的工具準確度高。

對於AUC-30的估計，edgeR表現最佳, DESeq2與之差別不大。

基於來講基於計數的工具比基於組裝的工具更高效, 不經過比對直接定量的工具如Salmon和kallisto能夠獲得高質量的差異分析結果。

以上三個圖都是散點圖，第一個Spearman rank correlation相關性越高越好，第二個RMSD類似於均方差(與對照相比得分偏差的平方和先求均值再開方), 第三個AUC-30表示在假陽性率為30%時ROC曲線下的面積，面積越大表示結果越準確 (縱軸是True positive rate，想像下那個曲線，原文中也有一個示例)。

文獻解讀完了，工具也選擇好了，圖也都可以重複了，就只剩下有人動動手，去實際操作了。