转录组分析工具大比拼 （完整翻译版）
本文是完整翻译版，生信宝典公众号同一期还有导读版
文献阅读与翻译 - 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。下面研究人员使用短读长的Illumina hESC、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参考转录组注释用作金标准。使用短读长和长读长转录本预测的组合方法(用”+”标记)稍微提升了短读长转录本预测方案的性能。R语言学习 - 散点图绘制
对于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样品之间表达差异百分比。R语言学习 - 箱线图一步法 R语言学习 - 箱线图（小提琴图、抖动图、区域散点图）R语言学习 - 热图简化 R语言学习 - 线图绘制
差异表达
不同样品和条件下差异表达基因的识别是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分析。
为了评价我们提出的分析流程在功能富集分析(看到了富集分析 R语言学习 - 富集分析泡泡图 （文末有彩蛋） & 富集分析DotPlot，可以服)的性能，相对于正常样品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，这些基因在胚胎发育过程中发挥重要作用。