史上最全 | 39个RNAseq分析工具与对比
文献：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上发表的使用RNAseq工具对比的一篇文献，解读这篇文献对我们使用RNAseq发文提供了思路。下面小编具体解说一下。
文献摘要：
RNA-sequencing(RNA-seq)是一个重要的转录组学研究技术,数百款分析工具目前已经开发出来。尽管最近相关研究评估了最新的可用的RNAseq工具,但他们没有全面综合的评估RNAseq分析的工作流。这里我们进行广泛的RNA-seq工作流的研究分析，不仅包括表达分析,我们的工作还包括了评估的RNA variant-calling,RNA编辑和RNA融合检测技术。更为独特的是我们对二代RNAseq和三代Isoseq技术都进行了研究,39个分析工具，~ 120种组合,涉及15个样品与各种生殖系、癌症和干细胞的数据集的~490种分析。我们报告了各流程性能并提出一个全面的，分析准确性高的RNA-seq分析流程,名字叫做RNACocktail。在不同的样品中验证表明,我们提出的流程可以帮助研究人员通过转录组的分析获取更多的生物有关的预测结果。
流程下载地址：http://bioinform.github.io/rnacocktail/
附录：39个工具版本号、重要参数及下载地址：
比对工具
1.TopHat2： –no-coverage-search
http://ccb.jhu.edu/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.jhu.edu/software/hisat/index.shtml
4.RASER 0.52 -b 0.03
https://www.ibp.ucla.edu/research/xiao/RASER.html
有参考转录本组装工具
1.Cufflinks 2.2.1 –frag-bias-correct
http://cole-trapnell-lab.github.io/cufflinks/
2.StringTie 1.2.1 -v -B
http://www.ccb.jhu.edu/software/stringtie/
无参考转录本组装工具
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.ebi.ac.uk/~zerbino/oases/
3. Trinity 2.1.1 –normalize reads
http://trinityrnaseq.sourceforge.net/
三代长read分析工具
1.LoRDEC 0.6 -k 23 -s 3
http://atgc.lirmm.fr/lordec/
2.GMAP 12/31/15 -f 1
http://research-pub.gene.com/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.uiowa.edu/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.github.io/eXpress/index.html
2. kallisto 0.42.4
http://pachterlab.github.io/kallisto/about.html
3. Sailfish 0.9.0
http://www.cs.cmu.edu/~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.sourceforge.net/
差异表达分析工具
1. DESeq2 1.14.1
http://bioconductor.org/packages/release/bioc/html/DESeq2.html
2. edgeR 3.16.5
http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
3. limma 3.30.7
http://bioconductor.org/packages/release/bioc/html/limma.html
4. Cuffdiff 2.2.1
–frag-bias-correct –emit-count-tables
http://cole-trapnell-lab.github.io/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.broadinstitute.org/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.genomics.org.cn/soapfuse.html
4.STAR-Fusion 0.7.0
https://github.com/STAR-Fusion/STAR-Fusion
5.TopHat-Fusion 2.0.14
http://ccb.jhu.edu/software/tophat/fusion_index.shtml
一、数据集
来源于人的 15个Illumina和 Pacific Biosciences (PacBio) 数据集
二、分析结果
（一）比对工具评价
不同方案检测到的剪接点利用与dbEST数据库中鉴定到的可靠的剪接点的一致性衡量各方案的准确性。 一个可靠的EST剪接点由至少两个EST支持， 圆圈的大小反映出来每个方案鉴定出的剪接位点数目。 对于每个工具，显示出鉴定剪接位点数和验证率（括号中）。每个数据集的验证率也在Venn图上显示。 b read比对效率分析：测序片段的read映射状态的分布（左）（对于NA12878，MCF7和SEQC样品，显示配对末端read的映射状态，而对于hESC，反映的是唯一映射（蓝色），多映射（橙色）和未映射（红色）单端read的映射情况），映射片段中soft-clipped的数目分布（中），映射片段中错配的数目的分布（右）
HISAT2所有样本中剪接位点验证率最高，尽管总数少于TopHat或STAR（图2a;补充图1-3）。STAR唯一映的read比例最高，特别是在MCF7-300上，大概是由于read增加长度（图2b）。 STAR映射只接受双端比对，不像TopHat和HISAT2一样可以接受单端read的比对。另一方面，STAR也产生了较差质量的比对结果，具有更多的soft-clipped和错配数（图2b）。 TopHat不允许截断read（图2b）。虽然这些结果证实了以前的发现，从较长的read的样本（MCF7-300）和单端测序样本（hESC）比对结果中可以看出STAR相对于TopHat和HISAT2具有更高的容忍性，用于接受不匹配和soft-clipped，以调整获得更高的mapping率。（图2b）。平均来说，HISAT2的速度分别比STAR和TopHat快2.5和〜100×（补充表3）
（二）基于比对的转录组组装
spliced aligned之后就是转录本组装了，有参考二代转录组数据组装常用的两个工具：Cufflinks和StringTie。除此之外你还评价了二+三混合组装工具IDP（分别使用GAMP和STARlong作为比对工具）和Pacbio官方转录本组装工具Iso-Seq，准确性评价采用GENCODE v19中的参考转录组。
Cufflinks和StringTie报告了更多单一外显子转录本（图3a;补充图4和5），其主要是假阳性的（补充图6）。 StringTie比cufflinks多预测50-200％的转录本。 IDP在各个样本中均预测出外显子数目最少，因为它不报告单外显子基因设计，在多个外显子转录本上，预测出的数目与Cufflinks数量相似
（图3a;补充图5）。而且，IDP的预测出的外显子数目分布更好地类似于GENCODE，特别是对于多外显子转录本（图3a）。平均来说，Iso-Seq算法预测差不多94％的单个外显子转录本和77％的多外显子转录本在GENCODE缺少。这个可能反映了Iso-Seq方法的组装准确性较差，但检测新的转录本灵敏度高。对于MCF7-300样本，STAR预测的数量多于其他比对软件（图3a;补充图5），可能是由于它处理更长的read能力。使用长read比对工具GMAP和短read比对工具HISAT2的IDP可以预测更多的可变剪接。
与短read组装工具不同，IDP倾向于检测一个基因的多个转录本（补充图7）。和cufflinks相比，StringTie平均预测基因数目多50倍以上且每个基因具有超过五种可变剪接。 StringTie的每个基因的可变剪接数量的分布与GENCODE中观察到更加一致（补充图7）。
对于基因水平评估，IDP在所有样品中达到最佳精度和灵敏度（图3b;补充图8和9）。此外，cufflinks比StringTie更敏感和精确。在MCF7-300样本上，不同比对工具之间有更多的差异，其中TopHat和HISAT2好与STAR。 Iso-Seq算法敏感度最低，而其精度在IDP和Cufflinks、StringTie之间。
转录本水平IDP比其他软件在精度上超过20％以上（图3b）。但是，它的预测结果准确性仅限制在多外显子上，它的敏感度比StringTie低，但是比cufflinks高。在短序列组装工具中，StringTie比cufflinks在转录本水平上高出平均有11％的精确度和25％的敏感性（图3b;补充图8和9）。 Iso-Seq接近零精确度，主要由于其构建转录本较差。对于StringTie和IDP，被预测出更多内含子的基因更有可能代表新的可变剪接，与以前使用长read的研究结果一致（补充图12）
StringTie是最快的工具，组装速度分别高〜60×和〜50×比cufflinks和IDP（输入的是错误纠正和对齐数据）（补充表4）。我们观察到，与以前的研究不同，在更多具有挑战性的例子中，如MCF7-300，STAR报道的更多的转录本数量（主要是单个外显子）但是也有更高的假阳性率（图3a;补充图4和5）。
(三) 从头组装
在这里我们分析了三种广泛应用从头组装工具Trinity，Oases和SOAPdenovo-Trans。
Trinity倾向于预测更长的可变剪接，更多的基因和转录本，但是许多是断裂的转录本（图4a;补充图16和17）。Oases在所有样本中产生了最高的N10至N50值（图4b;补充图18），表明其检测长的可变剪接的优越性；各软件检测到的不同表达量基因情况如（图4c;补充图19）。 SOAPdenovo-Trans最高峰在小的百分位数上（表达量从大到小排，类似于基因组N50），表明其倾向于检测高表达转录本。另一方面，Oases擅长检测低表达的基因（峰靠近右边）。
将重建的转录本与参考转录本进行比较表明，SOAPdenovo-Trans和Trinity在内含子水平分别具有最高精度和灵敏度（补充图21a）。对于内含子链级水平上，Oase和Trinity优于SOAPdenovo-Trans（补充图21b）。在较低的内存和计算要求下，SOAPdenovo-Trans表现最佳（补充表5）。
（四）三代长Read直接获取转录本
人类转录本长度（GENCODE v19注释）中位数为783 bp，比目前NGS技术可以提供的读长长得多。然而，长读长测序平台不用组装便可以轻松获得完全跨越大多数转录本的Read。 在hESC上样本，例如，原始PacBio的SubRead的中位数长度是1164bp，这足以覆盖大部分转录本（64％）。 因此，长读技术可以方便精确的发现转录本，无需外显子 - 外显子连接点预测或者组装。
我们使用GMAP和STARlong进行比对，结果作为IDP的输入。平均而言，GMAP的比对率比STARlong高28%（补充表7）。IDP另外一种可选输入是PacBio的Iso-Seq流程比对MCF7样品的结果。
在不同的样本上，基于长读技术的IDP和Iso-Seq预测了许多新的转录本或者已知的任何短读长测序技术都未检测到的参考转录本（补充图22）。对通过长读长或短读长预测的转录本统计分析表明只有IDP预测的转录本具有广泛的长度（达到10,000 bp），而由Iso-Seq预测的大部分转录本长度在1000到4000bp之间。
在速度方面，STARlong比GMAP快68倍（补充表8），而IDP每个样品大约耗时170个CPU小时
（五）转录本定量
基于比对的转录本定量。比较传统方法是将read比对（spliced -aligned）到参考基因组，然后利用Cufflinks和StringTie进行转录本组装，最后进行定量。如果具有参考转录本序列，reads可以直接跟转录本序列比对（aligned），然后使用RSEM和eXpress进行定量。
不经过比对（alignment-free）的转录本定量。主要提供了四个工具：Sailfish、Salmon、quasi-mapping和kallisto。不经过比对就可以确定哪个转录本生成哪些read或者寻找部分比对回转录本的reads。
在这里我们比较了基于基因组比对的cufflinks和StringTie（使用不同的比对工具），基于转录本比对的工具，eXpress和Salmon-Aln，不需要比对的kallisto，Sailfish，Salmon-SMEM和Salmon-Quasi，以及基于长读长技术的IDP（使用不同的短读长和长读长比对工具）四种方式的性能。
基于不同的定量方法所得表达值的Spearman相关性分析表明，具有相似方法的定量方案聚类在一起（图5a;补充图23和24）。不经过比对的方法各个工具也集中在一起，并且相比Cufflinks更接近于StringTie的位置。 Salmon-SMEM 与基于转录组比对的各工具聚在一起。鉴于Salmon-SMEM更快的速度，这使得其优于eXpress和Salmon-Aln。涉及IDP的组合也聚集在一起，与其他组合的相似性较小，特别是其中的涉及cufflinks的组合（图5a）。
两个免比对工具kallisto和Salmon-SMEM对MCF7-100和MCF7-300丰度估计具有最一致的结果（图5b，c）。反映出免比对工具在其丰度估计中无样本特异性和读长偏好性。 IDP对MCF7-100和MCF7-300丰度估计也表现出高度的一致性，特别是排除低表达基因（图5c）。在短读长比对工具中， HISAT2在不同样本中丰度估计的一致性最好（图5c）。
一般来说，免比对工具非常有效（补充表9），而带有高效比对工具如HISAT2的StringTie在基于对齐的方法中是最为高效的（比免比对工具慢一个数量级）。以前的研究表明在丰度估计准确性上估计的方法相对于比对工具而言具有更突出的作用，.我们的结果（图5c）清楚地描绘了HISAT2和TopHat相对STAR的优越性。
Fig. 5 转录本丰度估计各方法性能. 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
（六）差异表达
不同的时空以及不同的条件下差异基因分析是RNAseq分析的重要目标。差异表达分析方法包括：基于Read数目的DESeq、limma和edgeR；基于组装技术的Cuffdif和Ballgown；基于免比对的定量方法sleuth。
通过QPCR对各工具经行评价。与其他工具相比，DESeq2表现最佳。sleuth、edgeR和limma性能较差。Cuffdiff和Ballgown的准确度没有基于计数的工具准确度高。对于AUC-30的测量，edgeR表现最佳。平均而言，DESeq2在不同定量方法中均优于其他技术，而 sleuth,edgeR 和 limma的性能略有下降，这在之前文献中已经证实。Cuffdiff和 Ballgown准确度均低于基于原始read差异分析的技术。Salmon-SMEM, Salmon-Aln, kallisto和eXpress与基于原始read差异分析技术是最佳组合方案。在ROC曲线下低于30％（AUC-30）条件下，edgeR优于其他技术。
作为另一种准确度量，比较了不同的方案在预测92个External RNA
Control Consortium (ERCC) spike-in genes in the SEQC数据集上的性能（图6b;补充图29,35-38）。用Spearman相关性衡量，edgeR 和 limma明显着超过其他工具。用Spearman和RMSD同时评估，DESeq2仍然表现最好，而 sleuth优于edgeR和limma。然而，在AUC-30测量中，采用 Cufflinks的Ballgown的表现优于其他组合。基于read计数的工具比基于组装的工具更有效率，尤其是采用基于转录本的比对方法或免方法（补充表10）。Cuffdiff比Ballgown慢四到五倍，是最慢的工具。
总体而言，免对齐工具Salmon 和 kallisto能够提供高质量的差异基因预测。
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
（七）基因融合
基于短读长测序技术，FusionCatcher敏感性和准确性最高，SOAPfuse也显示高的敏感性。长读长技术 IDP fusion融合提供了最高的准确性（图7f）。STAR-Fusion是最快的方法（比其他方法快超过10×），而FusionCatcher和TopHat-Fusions具有更高的计算需求（补充表13）。
三、高准确度的分析流程
作者提出一个新的高准确度分析流程，RNACocktail，使用的具体软件如下图所示。
