文獻(xiàn)：Sahraeian S M E, Mohiyuddin M, Sebra R, et al. Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis[J]. Nature Communications, 2017, 8(1):59.

最近在NC上發(fā)這篇RNAseq工具對比的一篇文獻(xiàn)，為咱們小伙伴來分析RNAseq這種應(yīng)用最為廣泛的測序技術(shù)提供了思路。下面小編具體解說一下。

文獻(xiàn)摘要：

RNA-sequencing(RNA-seq)是一個重要的轉(zhuǎn)錄組學(xué)研究技術(shù),數(shù)百款分析工具目前已經(jīng)開發(fā)出來。盡管最近相關(guān)研究評估了最新的可用的RNAseq工具,但他們沒有全面綜合的評估RNAseq分析的工作流。這里我們進(jìn)行廣泛的RNA-seq工作流的研究分析，不僅包括表達(dá)分析,我們的工作還包括了評估的RNA variant-calling,RNA編輯和RNA融合檢測技術(shù)。更為獨特的是我們對二代RNAseq和三代Isoseq技術(shù)都進(jìn)行了研究,39個分析工具，~ 120種組合,涉及15個樣品與各種生殖系、癌癥和干細(xì)胞的數(shù)據(jù)集的~490種分析。我們報告了各流程性能并提出一個全面的，分析準(zhǔn)確性高的RNA-seq分析流程,名字叫做RNACocktail。在不同的樣品中驗證表明,我們提出的流程可以幫助研究人員通過轉(zhuǎn)錄組的分析獲取更多的生物有關(guān)的預(yù)測結(jié)果。

流程下載地址：http://bioinform./rnacocktail/

附錄：39個工具版本號、重要參數(shù)及下載地址：

比對工具

1.TopHat2： –no-coverage-search

http://ccb./software/tophat/index.shtml

2.STAR: -twopassMode Basic –outFilterType BySJout

https://github.com/alexdobin/STAR/releases

3.HISAT2 2.0.1-beta –dta (or –dta-cufflinks)

http://www.ccb./software/hisat/index.shtml

4.RASER 0.52 -b 0.03

https://www.ibp./research/xiao/RASER.html

有參考轉(zhuǎn)錄本組裝工具

1.Cufflinks 2.2.1 –frag-bias-correct

http://cole-trapnell-lab./cufflinks/

2.StringTie 1.2.1 -v -B

http://www.ccb./software/stringtie/

無參考轉(zhuǎn)錄本組裝工具

1.SOAPdenovoTrans 1.04 -K 25

https://github.com/aquaskyline/SOAPdenovo-Trans/

 2.Oases 0.2.09 (Velvetv1.2.10) (velveth haslength: 25) (velvetg options: -read trkg yes)

http://www./~zerbino/oases/

3. Trinity 2.1.1 –normalize reads 

http://trinityrnaseq./

三代長read分析工具

1.LoRDEC 0.6 -k 23 -s 3 

http://atgc./lordec/

2.GMAP 12/31/15 -f 1

http://research-pub./gmap/

3. STARlong 2.5.1b

https://github.com/alexdobin/STAR/releases

 Followed the recommended options : 

–outSAMattributes NH HI NM MD

 –readNameSeparator space 

–outFilterMultimapScoreRange 1 

–outFilterMismatchNmax 2000 

–scoreGapNoncan -20

 –scoreGapGCAG -4

 –scoreGapATAC -8 

–scoreDelOpen -1

 –scoreDelBase -1

 –scoreInsOpen -1

 –scoreInsBase -1

 –alignEndsType Local 

–seedSearchStartLmax 50

 –seedPerReadNmax 100000

 –seedPerWindowNmax 1000

 –alignTranscriptsPerReadNmax 100000

 –alignTranscriptsPerWindowNmax 10000

 –outSAMstrandField intronMotif 

–outSAMunmapped Within

 4. IDP 0.1.9 

https://www.healthcare./labs/au/IDP/

定量工具

1. eXpress 1.5.1 (bowtie2 v2.2.7) (bowtie2 options: -a -X 600 –rdg 6,5 –rfg 6,5 –score-min L,-.6,-.4 –no-discordant –no-mixed)

https://pachterlab./eXpress/index.html

2. kallisto 0.42.4 

http://pachterlab./kallisto/about.html

3. Sailfish 0.9.0 

http://www.cs./~ckingsf/software/sailfish/

4. Salmon-Aln 0.6.1 

https://github.com/COMBINE-lab/salmon

5. Salmon-SMEM 0.6.1

https://github.com/COMBINE-lab/salmon

 index: –type fmd 

quant: -k,19 

6. Salmon-Quasi 0.6.1 

https://github.com/COMBINE-lab/salmon

index: –type quasi -k 31 

7. featureCounts 1.5.0-p1 -p -B -C

http://subread./

差異表達(dá)分析工具

1. DESeq2 1.14.1 

http:///packages/release/bioc/html/DESeq2.html

2. edgeR 3.16.5

http://www./packages/release/bioc/html/edgeR.html

3. limma 3.30.7 

http:///packages/release/bioc/html/limma.html

4. Cuffdiff 2.2.1

 –frag-bias-correct –emit-count-tables

http://cole-trapnell-lab./cufflinks/

5. Ballgown 2.6.0 

https://github.com/alyssafrazee/ballgown

6. sleuth 0.28.1

https://github.com/pachterlab/sleuth

變異分析工具

1. SAMtools 1.2 (bcftools v1.2) 

samtools mpileup -C50 -d 100000

https://github.com/samtools/samtools

2. bcftools filter -s LowQual -e ‘%QUAL<20 —— DP>10000’

https://github.com/samtools/bcftools

 3.GATK v3.5-0-g36282e4 (picard 1.129)

https://software./gatk/download/

 Picard AddOrReplaceReadGroups: SO=coordinate

 Picard MarkDuplicates: CREATE INDEX=true VALIDATION STRINGENCY=SILENTGATK

 SplitNCigarReads: -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60

 -U ALLOW N CIGAR READSGATK 

HaplotypeCaller: -stand call conf 20.0

 -stand emit conf 20.0 -A StrandBiasBySample

 -A StrandAlleleCountsBySampleGATK 

VariantFiltration: -window 35 -cluster 3 -filterName FS -filter

 “FS >30.0” -filterName QD -filter “QD <2.0”

RNA編輯

1. GIREMI 0.2.1

https://github.com/zhqingit/giremi

2.  Varsim 0.5.1

https://github.com/bioinform/varsim

基因融合

1.FusionCatcher 0.99.5a beta 

https://github.com/ndaniel/fusioncatcher

2.JAFFA 1.0.6 

https://github.com/Oshlack/JAFFA

3.SOAPfuse 1.27

http://soap./soapfuse.html

 4.STAR-Fusion 0.7.0 

https://github.com/STAR-Fusion/STAR-Fusion

5.TopHat-Fusion 2.0.14

http://ccb./software/tophat/fusion_index.shtml

其中涉及到幾款常用工具小編有過講解：

轉(zhuǎn)錄組數(shù)據(jù)比對工具：Tophat

新的轉(zhuǎn)錄組組裝方法HISAT，Stringtieand Ballgown(一)

HISAT, StringTie and Ballgown（二）

samtools命令大全

Trinity轉(zhuǎn)錄組組裝軟件介紹

一、數(shù)據(jù)集

    來源于人的 15個Illumina和 Pacific Biosciences (PacBio) 數(shù)據(jù)集 

二、分析結(jié)果

（一）比對工具評價

  不同方案檢測到的剪接點利用與dbEST數(shù)據(jù)庫中鑒定到的可靠的剪接點的一致性衡量各方案的準(zhǔn)確性。 一個可靠的EST剪接點由至少兩個EST支持， 圓圈的大小反映出來每個方案鑒定出的剪接位點數(shù)目。 對于每個工具，顯示出鑒定剪接位點數(shù)和驗證率（括號中）。每個數(shù)據(jù)集的驗證率也在Venn圖上顯示。 b read比對效率分析：測序片段的read映射狀態(tài)的分布（左）（對于NA12878，MCF7和SEQC樣品，顯示配對末端read的映射狀態(tài)，而對于hESC，反映的是唯一映射（藍(lán)色），多映射（橙色）和未映射（紅色）單端read的映射情況），映射片段中soft-clipped的數(shù)目分布（中），映射片段中錯配的數(shù)目的分布（右）

（二）基于比對的轉(zhuǎn)錄組組裝

spliced aligned之后就是轉(zhuǎn)錄本組裝了，有參考二代轉(zhuǎn)錄組數(shù)據(jù)組裝常用的兩個工具：Cufflinks和StringTie。除此之外你還評價了二 三混合組裝工具IDP（分別使用GAMP和STARlong作為比對工具）和Pacbio官方轉(zhuǎn)錄本組裝工具Iso-Seq，準(zhǔn)確性評價采用GENCODE v19中的參考轉(zhuǎn)錄組。

Cufflinks和StringTie報告了更多單一外顯子轉(zhuǎn)錄本（圖3a;補充圖4和5），其主要是假陽性的（補充圖6）。 StringTie比cufflinks多預(yù)測50-200％的轉(zhuǎn)錄本。 IDP在各個樣本中均預(yù)測出外顯子數(shù)目最少，因為它不報告單外顯子基因設(shè)計，在多個外顯子轉(zhuǎn)錄本上，預(yù)測出的數(shù)目與Cufflinks數(shù)量相似

（圖3a;補充圖5）。而且，IDP的預(yù)測出的外顯子數(shù)目分布更好地類似于GENCODE，特別是對于多外顯子轉(zhuǎn)錄本（圖3a）。平均來說，Iso-Seq算法預(yù)測差不多94％的單個外顯子轉(zhuǎn)錄本和77％的多外顯子轉(zhuǎn)錄本在GENCODE缺少。這個可能反映了Iso-Seq方法的組裝準(zhǔn)確性較差，但檢測新的轉(zhuǎn)錄本靈敏度高。對于MCF7-300樣本，STAR預(yù)測的數(shù)量多于其他比對軟件（圖3a;補充圖5），可能是由于它處理更長的read能力。使用長read比對工具GMAP和短read比對工具HISAT2的IDP可以預(yù)測更多的可變剪接。

與短read組裝工具不同，IDP傾向于檢測一個基因的多個轉(zhuǎn)錄本（補充圖7）。和cufflinks相比，StringTie平均預(yù)測基因數(shù)目多50倍以上且每個基因具有超過五種可變剪接。 StringTie的每個基因的可變剪接數(shù)量的分布與GENCODE中觀察到更加一致（補充圖7）。

對于基因水平評估，IDP在所有樣品中達(dá)到最佳精度和靈敏度（圖3b;補充圖8和9）。此外，cufflinks比StringTie更敏感和精確。在MCF7-300樣本上，不同比對工具之間有更多的差異，其中TopHat和HISAT2好與STAR。 Iso-Seq算法敏感度最低，而其精度在IDP和Cufflinks、StringTie之間。

StringTie是最快的工具，組裝速度分別高?60×和?50×比cufflinks和IDP（輸入的是錯誤糾正和對齊數(shù)據(jù)）（補充表4）。我們觀察到，與以前的研究不同，在更多具有挑戰(zhàn)性的例子中，如MCF7-300，STAR報道的更多的轉(zhuǎn)錄本數(shù)量（主要是單個外顯子）但是也有更高的假陽性率（圖3a;補充圖4和5）。

在這里我們分析了三種廣泛應(yīng)用從頭組裝工具Trinity，Oases和SOAPdenovo-Trans。

Trinity傾向于預(yù)測更長的可變剪接，更多的基因和轉(zhuǎn)錄本，但是許多是斷裂的轉(zhuǎn)錄本（圖4a;補充圖16和17）。Oases在所有樣本中產(chǎn)生了最高的N10至N50值（圖4b;補充圖18），表明其檢測長的可變剪接的優(yōu)越性；各軟件檢測到的不同表達(dá)量基因情況如（圖4c;補充圖19）。 SOAPdenovo-Trans最高峰在小的百分位數(shù)上（表達(dá)量從大到小排，類似于基因組N50），表明其傾向于檢測高表達(dá)轉(zhuǎn)錄本。另一方面，Oases擅長檢測低表達(dá)的基因（峰靠近右邊）。

將重建的轉(zhuǎn)錄本與參考轉(zhuǎn)錄本進(jìn)行比較表明，SOAPdenovo-Trans和Trinity在內(nèi)含子水平分別具有最高精度和靈敏度（補充圖21a）。對于內(nèi)含子鏈級水平上，Oase和Trinity優(yōu)于SOAPdenovo-Trans（補充圖21b）。在較低的內(nèi)存和計算要求下，SOAPdenovo-Trans表現(xiàn)最佳（補充表5）。

（四）三代長Read直接獲取轉(zhuǎn)錄本

人類轉(zhuǎn)錄本長度（GENCODE v19注釋）中位數(shù)為783 bp，比目前NGS技術(shù)可以提供的讀長長得多。然而，長讀長測序平臺不用組裝便可以輕松獲得完全跨越大多數(shù)轉(zhuǎn)錄本的Read。 在hESC上樣本，例如，原始PacBio的SubRead的中位數(shù)長度是1164bp，這足以覆蓋大部分轉(zhuǎn)錄本（64％）。 因此，長讀技術(shù)可以方便精確的發(fā)現(xiàn)轉(zhuǎn)錄本，無需外顯子 - 外顯子連接點預(yù)測或者組裝。

我們使用GMAP和STARlong進(jìn)行比對，結(jié)果作為IDP的輸入。平均而言，GMAP的比對率比STARlong高28%（補充表7）。IDP另外一種可選輸入是PacBio的Iso-Seq流程比對MCF7樣品的結(jié)果。

在不同的樣本上，基于長讀技術(shù)的IDP和Iso-Seq預(yù)測了許多新的轉(zhuǎn)錄本或者已知的任何短讀長測序技術(shù)都未檢測到的參考轉(zhuǎn)錄本（補充圖22）。對通過長讀長或短讀長預(yù)測的轉(zhuǎn)錄本統(tǒng)計分析表明只有IDP預(yù)測的轉(zhuǎn)錄本具有廣泛的長度（達(dá)到10,000 bp），而由Iso-Seq預(yù)測的大部分轉(zhuǎn)錄本長度在1000到4000bp之間。

在速度方面，STARlong比GMAP快68倍（補充表8），而IDP每個樣品大約耗時170個CPU小時

（五）轉(zhuǎn)錄本定量

基于比對的轉(zhuǎn)錄本定量。比較傳統(tǒng)方法是將read比對（spliced -aligned）到參考基因組，然后利用Cufflinks和StringTie進(jìn)行轉(zhuǎn)錄本組裝，最后進(jìn)行定量。如果具有參考轉(zhuǎn)錄本序列，reads可以直接跟轉(zhuǎn)錄本序列比對（aligned），然后使用RSEM和eXpress進(jìn)行定量。

不經(jīng)過比對（alignment-free）的轉(zhuǎn)錄本定量。主要提供了四個工具：Sailfish、Salmon、quasi-mapping和kallisto。不經(jīng)過比對就可以確定哪個轉(zhuǎn)錄本生成哪些read或者尋找部分比對回轉(zhuǎn)錄本的reads。

在這里我們比較了基于基因組比對的cufflinks和StringTie（使用不同的比對工具），基于轉(zhuǎn)錄本比對的工具，eXpress和Salmon-Aln，不需要比對的kallisto，Sailfish，Salmon-SMEM和Salmon-Quasi，以及基于長讀長技術(shù)的IDP（使用不同的短讀長和長讀長比對工具）四種方式的性能。

基于不同的定量方法所得表達(dá)值的Spearman相關(guān)性分析表明，具有相似方法的定量方案聚類在一起（圖5a;補充圖23和24）。不經(jīng)過比對的方法各個工具也集中在一起，并且相比Cufflinks更接近于StringTie的位置。 Salmon-SMEM 與基于轉(zhuǎn)錄組比對的各工具聚在一起。鑒于Salmon-SMEM更快的速度，這使得其優(yōu)于eXpress和Salmon-Aln。涉及IDP的組合也聚集在一起，與其他組合的相似性較小，特別是其中的涉及cufflinks的組合（圖5a）。

兩個免比對工具kallisto和Salmon-SMEM對MCF7-100和MCF7-300豐度估計具有最一致的結(jié)果（圖5b，c）。反映出免比對工具在其豐度估計中無樣本特異性和讀長偏好性。 IDP對MCF7-100和MCF7-300豐度估計也表現(xiàn)出高度的一致性，特別是排除低表達(dá)基因（圖5c）。在短讀長比對工具中， HISAT2在不同樣本中豐度估計的一致性最好（圖5c）。

一般來說，免比對工具非常有效（補充表9），而帶有高效比對工具如HISAT2的StringTie在基于對齊的方法中是最為高效的（比免比對工具慢一個數(shù)量級）。以前的研究表明在豐度估計準(zhǔn)確性上估計的方法相對于比對工具而言具有更突出的作用，.我們的結(jié)果（圖5c）清楚地描繪了HISAT2和TopHat相對STAR的優(yōu)越性。

Fig. 5 轉(zhuǎn)錄本豐度估計各方法性能. a Clustering of different schemes based on the Spearman rank correlation of their log expressions on NA12878. b Distribution of log2-fold change of expressions between MCF7-100 and MCF7-300 samples. For each method, dashed line represents the mean of the distribution and the dotted lines represents the quartiles. c Percentage of expression disagreement between MCF7-100 and MCF7-300

samples when low-expressed transcripts are discarded with different thresholds

（六）差異表達(dá)

不同的時空以及不同的條件下差異基因分析是RNAseq分析的重要目標(biāo)。差異表達(dá)分析方法包括：基于Read數(shù)目的DESeq、limma和edgeR；基于組裝技術(shù)的Cuffdif和Ballgown；基于免比對的定量方法sleuth。

通過QPCR對各工具經(jīng)行評價。與其他工具相比，DESeq2表現(xiàn)最佳。sleuth、edgeR和limma性能較差。Cuffdiff和Ballgown的準(zhǔn)確度沒有基于計數(shù)的工具準(zhǔn)確度高。對于AUC-30的測量，edgeR表現(xiàn)最佳。平均而言，DESeq2在不同定量方法中均優(yōu)于其他技術(shù)，而 sleuth,edgeR 和 limma的性能略有下降，這在之前文獻(xiàn)中已經(jīng)證實。Cuffdiff和 Ballgown準(zhǔn)確度均低于基于原始read差異分析的技術(shù)。Salmon-SMEM, Salmon-Aln, kallisto和eXpress與基于原始read差異分析技術(shù)是最佳組合方案。在ROC曲線下低于30％（AUC-30）條件下，edgeR優(yōu)于其他技術(shù)。

作為另一種準(zhǔn)確度量，比較了不同的方案在預(yù)測92個External RNA

Control Consortium (ERCC) spike-in genes in the SEQC數(shù)據(jù)集上的性能（圖6b;補充圖29,35-38）。用Spearman相關(guān)性衡量，edgeR 和 limma明顯著超過其他工具。用Spearman和RMSD同時評估，DESeq2仍然表現(xiàn)最好，而 sleuth優(yōu)于edgeR和limma。然而，在AUC-30測量中，采用 Cufflinks的Ballgown的表現(xiàn)優(yōu)于其他組合?；趓ead計數(shù)的工具比基于組裝的工具更有效率，尤其是采用基于轉(zhuǎn)錄本的比對方法或免方法（補充表10）。Cuffdiff比Ballgown慢四到五倍，是最慢的工具。

總體而言，免對齊工具Salmon 和 kallisto能夠提供高質(zhì)量的差異基因預(yù)測。

Fig. 6 Performance of differential gene expressions analysis tools on SEQC-A vs. SEQC-B samples. a Spearman rank correlation, root-mean-score-deviation (RMSD), and AUC-30 scores for qPCR measured genes. Spearman rank correlation and RMSD scores are measured between the log2-fold change of the qRT-PCR and RNA-seq tools. AUC-30 score represents the area under the ROC curve up to the false positive rate of 30%. b ROC analysis of qRT-PCR measured genes (left) and ERCC (right) genes. For each differential analysis tool the plot reflects average performance when different alignment-based and alignment-free tools are used for abundance estimation and error bar shows the maximum and minimum variations. Results for each tool combination are

shown in Supplementary Figs. 30 and 35

基于短讀長測序技術(shù)，F(xiàn)usionCatcher敏感性和準(zhǔn)確性最高，SOAPfuse也顯示高的敏感性。長讀長技術(shù) IDP fusion融合提供了最高的準(zhǔn)確性（圖7f）。STAR-Fusion是最快的方法（比其他方法快超過10×），而FusionCatcher和TopHat-Fusions具有更高的計算需求（補充表13）。

三、高準(zhǔn)確度的分析流程

作者提出一個新的高準(zhǔn)確度分析流程，RNACocktail，使用的具體軟件如下圖所示。