前言

NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

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，想象下那個曲線，原文中也有一個示例)。

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

轉錄組研究

• NGS基礎 - FASTQ格式解釋和質量評估

• NGS基礎 - 高通量測序原理

• NGS基礎 - 參考基因組和基因注釋文件

• NGS基礎 - GTF/GFF文件格式解讀和轉換

• NGS基礎 - 測序原始數據下載

高顏值免費在線繪圖

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的39个工具，120种组合深度评估 (转录组分析工具哪家强)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。