Evaluation of tools for long read RNA-seq splice-aware alignment论文详解

Krešimir Križanović, Amina Echchiki, Julien Roux, Mile Šikić, Evaluation of tools for long read RNA-seq splice-aware alignment, Bioinformatics, Volume 34, Issue 5, 01 March 2018, Pages 748–754, https://doi.org/10.1093/bioinformatics/btx668

前言

三代测序数据的长度变长，错误率变高。本研究探索了目前可用的RNA-seq剪接对齐工具如何应对读段长度和错误率的增加。所有测试工具最初都是为短NGS读取而开发的，但有些工具声称支持长太平洋生物科学(PacBio)甚至牛津纳米孔技术(ONT)的长读段。

1 引言

目的：确定目前可用的RNA-seq剪接定位器是否能够处理第三代测序数据，即更长读段长度和明显更高错误率数据。

这种RNA-seq校准工具和管道的基准测试

• 以前在真实的和合成的Illumina reads (Engstro¨m et al.， 2013)上都进行过，证明非常有帮助。
• 另一种RNA-seq比对工具的基准是对不同错误率和复杂度的合成数据进行比对(Baruzzo et al.， 2017)。
• 然而，据我们所知，没有对第三代测序数据进行测试。

能够识别拼接的RNA-seq比对工具可以分为两组。

• 首先，引导剪接感知配对器使用基因组序列和已知的基因注释来计算基因或转录本丰度，但不能用于识别新的剪接连接。
• 其次，能够识别新剪接的对准器可以将RNA-seq读序列与参考基因组序列进行比对，而无需事先获得基因注释的信息。

工具：

• BBMap：唯一明确声明支持PacBio和ONT读段的工具。使用短k-mers将读段码直接与基因组对齐，跨越内含子寻找新的亚型。使用一个自定义的仿射变换矩阵来生成对齐分数；
• STAR：在未压缩后缀数组中使用序列最大可映射种子搜索，然后进行种子聚类和拼接过程。检测新的经典的、非经典的剪接连接和嵌合序列；
• GMAP：GMAP/GSNAP包的一部分，通过diagonalizetion对角化找到外显子区域，通过链接短k-mers的低聚物来细化他们，并在核苷酸水平上进行动态规划来解决不匹配、indels和内含子边界。
• TopHat2：Illumina reads最流行的对准器，实现两步方法，首先分析初始读段序列发现exon-exon链接，然后在第二步中使用它们来确定最终的序列。
• HISAT2：TopHat2的继承者。使用全局FM索引，以及一大堆小的FM索引（称为局部索引），共同覆盖了整个基因组，这种策略使得跨多个外显子的RNAseq读段能够有效对齐。

高错误率-错误纠正。

2 材料与方法

为什么使用模拟和真实数据集？

• 由于真实数据集中读段的实际来源是未知的，只能通过比对过程来估计，因此真实数据集并不最适合评估对齐工具的性能；
• 对准器的准确度和精度可以再合成数据上进行评估，但是，模拟器无法模拟真实数据集的每一个方面，有可能导致基准测试结果出现偏差。

数据集：

• 所有真实数据集由RNA转化为cDNA并在测序前扩增而成。
  为了进行仿真，使用PacBio读段模拟器PBSIM；
  用不同的参数模拟多个数据集，并使用不同生物体的注释转录组（酿酒酵母、黑腹果蝇和人类。由上图所示）；
• 为了更精确探索测试一些对准器的不佳性能，模拟一个包含很少错误的长读段数据集，可以知道对准器性能差是由于读段长还是错误率高。
• 测试重点是PacBio技术，因此拥有大量真实数据和专用模拟器（PBSIM）。
• 包含一个真实的ONT数据集。为了比对，我们使用PBSIM在黑腹果蝇上模拟了一个ONT MinION数据集，根据统计的ONT MinION R9真实数据设置参数。虽然PacBio模拟器并不完全适合于ONT MinION数据，但我们认为模拟他们的读段长度和错误配置（插入、删除和不匹配率）应该可以提供一些有用的剪接。——在进行模拟实验的时候，没有意识到一个专用的MinION读段模拟器，开始了解NanoSim，但是由于时间限制，不将其纳入基准测试。
• 利用人类染色体19模拟了额外的合成ONT MinION数据集。结果类似于第一个模拟数据集上取得的结果。
• 为了探索读段错误教程对比对的影响，使用最近认可的工具Racon对高质量的真实PacBio数据集进行了误差校正。探讨了使用外部Illumina读段的矫正和自校正。

从table1中可以看出数据集的大小和复杂性各不相同。例如，数据集2和4具有相似的大小，因为他们是使用相同的基因覆盖直方图近似值生成的，而由于MinION ONT读段平均比PacBio读段长，数据集2比数据集4包含更多的读段。所有用于创建测试数据集的数据（以及数据本身）可以通过FigShare获得https://figshare.com/projects/RNAseq_benchmark/24391。

2.1 数据集

为了生成模拟数据集，使用从https://code.google.com/archive/p/pbsim/下载的1.0.3版本的PBSIM。
合成数据集由一下生物体构建：

• 酿酒酵母Sacharomyces cerevisiae S288(baker’s yeast)
• 黑腹果蝇Drosophia melanogaster r6(fruit fly)
• 人类Homo Sapiens GRCh38.p7(human)

所有生物体的参考基因组从 http://www.ncbi.nlm.nih.gov下载。
PBSIM将被用作基因组读段模拟器，以参考序列和一组模拟参数（如覆盖率、读段长度、误差分布）作为输入。为了生成RNA-seq读段，使用 https://genome.ucsc.edu/cgi-bin/hgTables 下载基因注释，将PBSIM应用于由特定基因组生成的一组转录本。为了使数据尽可能真实，分析和使用真实数据集确定仿真参数。使用真实表达数据集选择一组转录本进行模拟（下载地址： http://bowtie-bio.sourceforge.net/recount/，三种物种数据集都被使用）。
附件S1中详细描述了模拟合成数据的整个过程。

本基准测试中使用的真实RNA-seq数据集来自D. melanogaster。
同一样品的技术重复测序采用三种不同的技术:Illumina HiSeq、PacBio RSII和ONT MinION。
Illumina数据用于所有测试工具的基准比较和PacBio读段的误差校正。
PacBio和MinION数据被用来评估对齐性能，并确定错误轮廓，然后用于模拟合成数据。

我们总共使用了:

• 1GB的Illumina读取，从一个更大的数据集随机抽样。读数大小为101 bp。Illumina数据只是作为一个基线被包括进来，以表明所有的工具在Illumina上工作得相当好，并使用它们进行错误校正。正因为如此，我们使用了Illumina read没有匹配的末端信息。
• 超过5GB的PacBio子片段，测序自三种不同大小的转录本(1-2 kb, 2 - 3 kb和3-7 kb，每个大小片段有2个smrt细胞测序)。这对应于从子读段中提取大约2GB的Insert读段。
• 350MB ONT MinION读段使用R9化学。由于一维读取的质量很低，所以在这个基准测试中只使用了二维读段。

2.2 纠错

为了测试使用误差校正是否可以改善对齐结果，我们对质量最高的PacBio数据集(包含ROIs)进行了校正。
使用Racon进行误差校正(Vaser et al.， 2017)。
使用Illumina reads进行校正，并进行自校正测试。
由于自校正数据集被证明有更好的错误配置，只有这个数据集被保留作为基准(数据集统计信息在补充表S1中给出)。
补充表S1显示了所有真实数据集的错误率和读取长度统计数据，包括使用纠错获得的所有数据集。

2.3 评估RNA-seq工具

我们测试了最近更新的5个RNA-seq校准工具，反映出它们仍在维护中。

2.3.1 STAR

从 https://github.com/alexdobin/STAR 下载。使用的版本是2.5.2b。在Illumina公司集星运行使用regular script STAR ,在长期使用STARlong脚本读取数据集星运行参数显示在Bioinfx研究:优化明星对准器Iso-Seq PacBio GitHub的数据页面( https://github.com/PacificBiosciences/cDNA_primer/wiki /Bioinfx-study:-Optimizing-STAR-aligner-for-Iso-Seq-data )。此网址应失效，PacBio公司加了blasr，长读段对准器，还有minimap2。

2.3.2 Tophat2

从 https://ccb.jhu.edu/software/tophat/index.shtml 与Bowtie2一起使用下载二进制文件。使用的版本是2.1.1，带有用于对齐的默认参数。使用SAMTools 1.2版本将Tophat输出从BAM格式转换为SAM格式。

2.3.3 Hisat2

从 https://ccb.jhu.edu/software/hisat2/index.shtml 下载二进制文件。版本2.0.4使用了默认的对齐参数。

2.3.4 BBMap

从 https://sourceforge.net/projects/bbmap/ 下载。使用了脚本mapPacBio.sh。使用的是BBMap 35.92版本。读取首先使用samscripts工具( https://github.com/isovic/samscripts )转换为FASTA格式(最初是FASTQ格式)。然后运行程序，并将fastareadlen选项设置为适合每个数据集的值。

2.3.5 GMAP/GSNAP

源代码可以从 http://research-pub.gene.com/gmap/ 下载。使用版本2016-11-07。GMAP与默认参数一起使用，正如在使用PacBio数据的教程中推荐的那样( https://github.com/pacificbiosciences/cdna_primer/wiki /Aligner-tutorial%3A-GMAP%2C-STAR%2C-BLAT%2Cand-BLASR 网址已失效)。
我们也在Illumina数据集上运行GSNAP(因为它是专为短读而设计的)，但是在默认参数和无匹配端信息的情况下，它比GMAP映射的读数要少一些，所以我们决定不使用它。
用于运行每个工具的确切命令可以在补充说明S2中找到。

2.4 RNAseqEval工具

• 5个RNA序列校准器中的3个校准器对资源使用和校准质量进行了评估。使用Cgmemtime工具的分支计算CPU和内存消耗( https://github.com/isovic/cgmemtime.git )。
• 为了评估每个aligner的质量，我们开发了RNAseqEval ( https://github.com/kkrizanovic/RNAseqEval )这是一个用于评估RNA-seq比对的通用工具。

它是用Python编写的，包含两个主要脚本

• 一个用于评估使用PBSIM模拟的数据；
• 另一个用于评估真实数据或来源未知的数据。

这两个脚本都需要SAM格式的对准器输出，它们将其与基因注释进行比较，如果是模拟数据，则需要以MAF格式的对齐文件来描述每次模拟读取的来源。

2.4.1 评估合成数据

用于评估合成或模拟数据的脚本目前只工作在用PBSIM模拟的数据上，但是将来可以扩展以支持其他模拟器。

除了SAM格式的对准器输出和GTF或BED格式的基因注释外，该脚本还接受一个包含PBSIM生成的文件的文件夹——包含PBSIM数据的文件夹需要有一个特定的结构并遵循一个特定的命名约定，如程序文档中所描述的那样。

对于每一次从对准器读取的输出，脚本将使用PBSIM生成的MAF文件和基因注释在参考基因组中查找其来源，并将其与对准器计算的比对进行比较。
比对序列的开始和结束位置以及读取起点的位置，可以容忍五个核苷酸的错误。
脚本输出与染色体、链和起始位置精确对齐的读取数的摘要信息。

2.4.2 评估真实数据

用于评估真实数据的脚本仅将SAM格式的对准器输出和GTF或BED格式的基因注释作为输入。

由于读段的来源未知，该脚本将检查读段重叠的基因注释，然后评估read对齐与该基因的外显子和内含子匹配的程度。

当对注释中的每个外显子进行比对时，五个核苷酸的错误是可以容忍的。类似地，对齐和外显子注释之间的重叠需要至少有5个碱基对才能被认为是有效的。我们测试了允许误差(和最小重叠)的不同值，将其增加到5个碱基对以上并没有显著改善结果。

3 结果

3.1 基线比较

我们首先检查了对齐工具在Illumina“基线”数据集A上的表现(表2)。我们发现所有对齐工具都成功地对齐了Illumina的大部分读段。
然而，在包含更长的和更多错误读取的数据集(数据集1到数据集8)上，工具之间存在很大的差异。特别地，使用默认参数的Tophat2和Hisat2对所有长时间读取数据集的读取对齐<7%。为了公平起见，必须声明他们没有声称使用长读取，并且为了完整性的考虑而被包含在测试中。因此，我们没有在进一步的分析中考虑这两种工具，而是将重点放在剩下的三个对准器上:BBMap、GMAP和STAR。

• 如果我们看一下结果数据集B(低错误的长读取)，我们会发现Tophat2和Hisat2无法对齐使用默认参数的几乎所有读取(表中的数字被四角化)。我们可以得出结论，Tophat2和Hisat2是专为短NGS读取而设计的，不能处理较长的读段。

• 基于对齐的读段百分比，GMAP获得了最好的结果，它对所有测试数据集中>85%的读数进行了对齐。

• BBMap表现略好 Illumina公司(数据集A)和合成酿酒酵母PacBio数据集(数据集1,它包含很少multi-exon成绩单),但是读取对齐的分数落后GMAP更复杂的合成数据集和真实数据集(例如,只有43%的合成智人PacBio读取数据集4是一致的)。

• STAR成功对齐了Illumina大部分读取值(96.8%)，但其在第三代测序数据集上的表现并不均匀，对齐率从0.1%到67.2%，而且常常对齐不到一半的读取值。STAR似乎受到了数据集复杂性增加的影响，以及错误率增加的影响(Illumina和经过错误校正的PacBio数据集获得了最佳性能)。由于STAR成功地对齐了数据集B的重要部分(长读取低错误)，我们可以得出结论，它可以处理长读段，但在较高的错误率方面存在问题，特别是在更复杂的数据集上。

• 通过对数据集5和数据集7的比较，可以看出，在所有工具中，纠错改进了对齐率。

总之，对于一些对准者来说，第三代测序技术reads的比对百分比与Illumina reads的比对百分比相似。然而，只看每个工具对基因组的读取量并不能可靠地衡量总体比对质量。例如，一种工具可以对大部分读取序列进行排列，但只能对它们读段的一部分进行排列，或者可以将它们排列到基因组不正确的位置。

3.2 合成数据

fig1 合成数据的评估

为了获得更多的质量比对,我们评估了四种对准器合成数据集生成不同复杂性的转录组使用PBSIM工具(材料和方法),并应该反映特征PacBio(数据集1-3)和ONT MinION技术(数据集4)。

在这些数据集,每个读段的确切起源是已知的,通过检查比对位置与基因组中原始位置的匹配程度，可以评估比对质量。
使用RNAseqEval工具对这些数据集的对齐结果进行评估，如图1所示。

所有结果显示为数据集中所有读取的百分比。
读段显示对齐的百分比(没有评估的准确性)：对齐读段的匹配率的读段百分比，开始、末尾和内部的便捷在五个碱基对内被正确的比对（Correct），重叠读段的所有外显子起源(Hit all)和重叠读段至少一个外显子的读段来源(Hit one)。匹配率是按与参考基因组上对应的碱基数相等的对齐碱基的百分比计算的。Hit one和Hit all需要至少五个base。
各合成reads的评价结果见表3。还显示了对剪接读段子集(即对参考基因组上多个非连续位置进行对齐的读取)的评估。剪接读段，如果比对正确，应该在原始转录本上至少重叠一个外显子-外显子连接，从而覆盖两个或更多的外显子。表3所示读段的百分比与输入读段的数量有关;相对于对齐读取数的百分比显示在补充表S2中。

总的来说，GMAP给出的排列最准确，BBMap紧随其后，而STAR的排列比其他两个要差。例外是数据集1，其中BBMap证明略优于GMAP。在数据集2、3和4上，GMAP在校正一般基因组位置(命中所有和命中一个)和正确确定它们的确切起始位置(正确)方面都超过了其他两种工具。

STAR读段大多对齐到正确的一般基因组位置(全部命中和一个命中)，并且显示出非常好的匹配率，然而，总体低读段比对比例(表3和4)不允许该工具与GMAP和BBmap进行良好的比较。此外，STAR在正确对齐读段的开头和结尾方面并没有表现得特别好。

数据集2、3和4包含大量的剪接读段。聚焦于这些数据集的剪接读段统计数据，BBMap的表现明显比GMAP差，有时甚至比STAR差:在数据集3上，它设法从一个读源重叠所有外显子(分裂命中所有外显子)的精确度低于STAR(10.2%对19.4%)。对于STAR, 剪接读段的结果与总体结果一致，但总体对齐读段的数量很少，因此不推荐STAR用于第三代测序RNA-seq reads的对齐。

总的来说，BBMap在更低的复杂性的数据集1的校准精度上比GMAP更好(更少的多外显子基因)，但是在更复杂的数据集上总体校准效率落后，有时差距很大。这表明在使用BBMap时，应该谨慎地对齐拆分的RNA-seq读取。在这种情况下，GMAP的性能最好，应该成为首选，尽管数据集1的结果表明，在处理第三代测序数据的高错误率方面，GMAP仍有一定的改进空间。

3.3 真实数据

对于真实的数据，我们不知道每个读段的来源，因此我们通过将读段的比对位置与给定的一组基因注释进行比较来评估比对结果。还提取了比对匹配率、表达基因数等相关数据(表4)。表4所示reads的百分比与输入reads的数量有关。补充表S3还显示了读取相对于对齐的读取数的百分比。

该表显示了与基因注释(连续的外显子对齐)相对应的对齐读段的百分比(不评估准确性)，重叠至少一个外显子的读段的百分比(外显子命中)和重叠一个或多个外显子的读段的百分比。所有值都显示为数据集中所有读取操作的百分比。该表还显示了表达基因的数量和对齐reads的平均匹配率。匹配率是按与参考基因组上对应的碱基数相等的对齐碱基的百分比计算的。外显子命中统计的重叠至少需要5个碱基。

所有真实数据集都是在不同平台上对相同的D. melanogaster样本进行rna测序的技术复制。有趣的是，这些数据集具有不同的错误特征(补充表S1)。

正如之前的测试所预期的那样，GMAP显示了最好的结果，紧随其后的是BBMap。GMAP在将reads对齐到基因组中注释的外显子位置方面稍好一些。对齐读段的匹配率大致等于每个数据集的确定的错误配置(如补充表S1所示)，因此表明读段是对齐到正确位置的。GMAP甚至能够以合理的准确性对齐所有的ONT MinION数据。有趣的是，根据某些标准，GMAP在低质量数据集7(由子读取组成)上显示的结果比高质量数据集5(由ROI组成)和数据集6(经错误修正的ROI)上显示的结果更好。
BBMap和GMAP都报告了很大比例的ONT MinION reads对齐，但是匹配率和外显子命中率都低于PacBio数据集，这表明这些对齐在错误位置上的比例更大。
STAR显示了最差的对齐结果。成功对齐的读取显示了高匹配率，这可能反映了这样一个事实:STAR无法以最高的错误率对齐读段，或者对齐设置非常保守。

补充表S1显示，纠错在一定程度上改善了错误配置，平均匹配率增加了2-3%。然而，即使是这一微小的改进也导致了数据集6上所有对准器的明显更好的比对结果:更多的读段报告为对齐，更多的外显子击中，更多的基因表达和更高的匹配率。如表4所示，STAR从错误校正中获益最大，BBMap略少，而GMAP获益最少。由此可以得出的结论是，GMAP对错误的容忍度最高，其次是BBMap, STAR的容忍度最低。表2所示的“长读低错误”数据集B的结果支持这一点。
Fig. 2. Aligned read percentage violin plots for GMAP and STAR

最后，我们检查了读段长度的哪一部分被对齐(图2)。结果与其他比对质量的方法一致，与GMAP相比，STAR平均能够对读段长度的更大一部分进行对齐。不显示BBMap结果，因为在测试的设置中，所有的对齐都在读取的整个长度上(全局对齐)。这使得图2中为BBMap绘制的小提琴图毫无用处，因为每次读取都是沿着其长度的100%对齐的。这种行为在报告的结果中有一定的隐含意义，因为读取两端的对齐有时是不正确的，从而导致较低的匹配率。剪接从BBMap产生的比对可能是好的，例如使用’ local '标志，它通过剪辑全局比对转化为本地比对，如果结果是更高的分数。

3.4 资源使用情况

为了评估每个RNA对准器的效率，我们测量了CPU时间和最大内存使用量(常驻集大小- RSS)。所有工具都在多线程环境中运行，尽可能在12个线程上运行，并测量总CPU时间。结果如图S1所示。本次分析中省略了Illumina数据(dataset A)和长读低误差数据(dataset B)，因为本文的重点是第三代测序数据。
运行时间似乎取决于数据集的大小。在所有设置中，GMAP使用的内存最少，运行速度最快。STAR是最慢的，一直使用60-80 GB的RAM。BBMap的内存占用也一直在10-15 GB左右。

4 总结

近年来，第三代测序仪器在基因组学研究中稳步站稳脚跟。这些技术有望解决NGS读取长度短造成的问题。对于RNA-seq分析，更长的读取应该显著改善转录本识别。然而，第三代测序技术也带来了新的生物信息学挑战，主要是由于其高错误率。

在这项研究中，我们试图评估目前可用的RNA-seq比对工具对第三代测序数据的工作能力。利用真实数据集和合成数据集对五种对准工具进行了测试。

• Hisat2和Tophat2几乎无法对齐任何读段。
• STAR只在错误少的数据集上显示可通过的结果，但在高错误易发率的MinION数据上几乎完全失败。
• BBMap的表现相当好，尤其是在PacBio ROI reads(错误率较低)和具有较少多外显子基因的简单生物体上。这似乎表明，虽然BBMap是一种可识别拼接的对准器，但它在连续对准上取得了最佳性能(例如，来自dna序列分析)，可能并不最适合RNA序列分析数据。
• 最后，GMAP显示了最佳的对齐结果。它运行速度最快，占用内存最少，通常产生最高的对齐率，特别是在复杂数据集上。
  BBMap仅在低复杂度的模拟数据集上优于GMAP，其中包含很少的分割读取，这表明尽管GMAP在很大程度上优于其他对准器，但它仍有改进的空间。

GMAP在数据集4中特别突出，该数据集包含基于wine fly基因组的模拟ONT MinION读取。GMAP能将97%以上的数据映射到与读源至少一个外显子重叠的近似正确的位置，而第二好的对准器(BBMap)则能将数据映射到小于50%的位置。在真实的ONT MinION数据集(数据集8)和在补充说明S4中模拟的人类染色体19上的ONT MinION数据集上的映射质量差异要小得多。

总的来说，目前在一些可用的工具(即GMAP和BBMap)上对齐第三代测序RNA reads是可行的，但我们对其定位的低精度感到惊讶。除了数据集1(主要包含单外显子读取)，最佳对准器(GMAP)将读段的20%到31%赋值于它们正确的起始位置(65个碱基)。目前还不清楚这一结果是由于技术的高错误率所固有的，还是由于校准算法(最初不是为这些类型的数据开发的)，还是由于基准测试中使用的特定参数。例如，测试剪切BBMap对齐对其整体性能的影响将是需要注意的。

通过开发新的更复杂、更敏感的算法，或者在读取比对之前在生物信息学管道中加入纠错步骤，可能还有很大的改进空间，因为在我们的测试中，这明显改善了比对结果。