基于hash的下一代DNA测序长读序列映射算法剖析

摘要

动机：最近，有人提出了一系列将短读映射到参考基因组的程序。其中许多都针对短读映射进行了大量优化，因此对于较短的查询非常有效，但这使得它们效率低下，或者不适用于长度超过200 bp的读取。然而，许多测序器已经在生成更长的读取，预计还会有更多的读取。对于长读取序列映射，选项有限；BLAT、SSAHA2、FANGS和BWA-SW是最受欢迎的产品。然而，重测序和个性化医学需要更快的软件将这些长测序读数映射到参考基因组，以识别SNP或罕见转录物。

结果：我们提出了AGILE（AliGnIng Long rEads），这是一种基于哈希表的高通量序列映射算法，用于较长的454次读取，该算法使用对角多种子匹配标准、自定义q-gram过滤和动态增量搜索方法等启发式方法来优化映射过程的每一步。在我们的实验中，我们观察到AGILE比BLAT更准确，并且与BWA-SW和SSAHA2相当。对于实际错误率（<5%）和读取长度（200–1000 bp），AGILE明显快于BLAT、SSAHA2和BWA-SW。即使在其他情况下，AGILE也可以与BWA-SW媲美，比BLAT和SSAHA2快几倍。

可利用性： http://www.ece.northwestern.edu/~smi539/agile.html.

联系人： smi539@eecs.northwestern.edu

补充信息： 补充数据可在生物信息学在线。

1简介

新一代测序（NGS）技术的最新进展已导致价格合理的台式测序器，具有低运行成本和高吞吐量。这些测序器产生测序过程中测序的基因组的小片段。通过将这些小片段（读取）映射到参考基因组，我们可以对新个体的DNA进行测序。NGS使得大规模进行这些研究成为可能。这些进展将开创个人基因组学的时代，每个人都可以对自己的DNA进行测序和研究，以开发出更个性化的预测、诊断和治疗疾病的方法(Patrick，2007年).

这种性质的研究已经开始。科学家们发现了导致Charcot–Marie–Tooth等疾病的遗传原因(卢普西语等。, 2010)和米勒综合征(罗奇等。, 2010)通过对患者基因组进行测序。通过降低成本和提高高通量测序的速度，这些研究成为可能。下一代测序器（NGS）通过生成称为reads的DNA小子串来对DNA进行测序。随着测序技术的快速改进，读取的长度不断增加。例如，454个读数的长度从2007年的约250个基点增加到2009年的约500个基点。照明读数的长度从2007年的约30-50个基点增加到2009年的约100个基点。太平洋生物科学公司（Pacific Biosciences）在2009年也宣布了1000个基点的长读数。这些NGS的吞吐量和读取长度的增长速度甚至让摩尔定律都相形见绌。因此，人们越来越需要能够进行更长时间读取并且仍然能够与NGS的速度相匹配的工具。

已经开发了许多工具用于更短的illumina查询。其中包括MAQ(锂等。2008年a)、ELAND、SOAP(锂等。2008年b)、领结(朗米德等。, 2009)、摩赛克、帕斯(坎帕尼亚等。, 2009)和SHRiMP(隆隆声等。, 2009). 然而，大多数这些工具仅适用于读取长度<200的情况。此外，它们允许很少数量的不匹配（通常<2），并且其中许多不允许任何间隙。然而，随着读取长度的迅速增加，我们需要能够处理更长读取长度的工具。此外，这些用于更长时间读取的新工具应该能够处理更多的间隙和不匹配。据我们所知，专门设计用于长时间读取的其他工具只有BWA-SW(Li和Durbin，2010年)和FANGS(米斯拉等。, 2009). BWA-SW是一个基于Burrows–Wheeler Transform（BWT）的软件包。它支持查询的间隙全局对齐，是最快的长读取对齐算法之一，同时还可以查找次优匹配。

哈希表被广泛用于短读映射和许多其他相关问题(阿尔特舒尔等。, 1990；肯特，2002；锂等。2008年b；宁等。2001年；皮尔逊和利普曼，1988年；史密斯等。, 2008). 因此，我们可能已经用尽了用于序列映射的哈希表的所有可能用途。然而，我们会发现，通过适当的启发式，哈希表也可以为较长读取的序列映射提供极好的加速。任何基于哈希表的序列映射算法的一般结构如下：（i）创建基因组的哈希表索引；（ii）使用该索引查找基因组中可能具有同源性的区域；（iii）更详细地检查每个区域并输出确实同源的区域。现有基于哈希表的算法的执行时间主要由阶段（iii）决定。此阶段的处理时间与找到的区域数成正比。FANGS也是一个基于散列表的工具。它使用q-gram过滤和鸽子洞原理过滤出许多区域，以接近100%的灵敏度快速映射454个读取。然而，FANGS在错误率大于1%时效率低下。除了FANGS使用的技术外，使用这些技术快速过滤出不同源的区域将大大减少所花费的时间。我们沿着这条路走。

在本文中，我们将讨论五种不同的技术来过滤出非同源区域及其对高通量长读序列映射问题的适应性。随后，我们提出了另一种基于哈希表的序列映射工具AGILE-yet，我们使用这些过滤技术来优化序列映射过程的每一步。

2高吞吐量序列映射

在NGS的背景下，序列映射问题涉及在参考基因组中搜索允许少量差异的小DNA序列（读取）。参考基因组是从与读数相同物种的生物体中获得的，这意味着读数和参考基因组之间具有高度相似性。允许存在少量差异，以解释单个生物体之间的差异和测序错误。

这就给我们带来了序列映射问题的形式化定义。我们可以将每个基因组序列表示为字母表∑=上的字符串{A类,C,G公司,T型}. 给定基因组数据库G公司主题序列的{S公司₁,S公司₂,···,S公司_我}，查询序列（读取）问长度的|问|，基因组测序问题是找到G公司最小值为E类(α,问)对于所有α。这个问题可以简化为查找查询的最佳匹配α，以便E类(α,问)小于某个界限。我们可以继续增加界限，直到找到匹配的为止。此外，对于给定的排序错误率，与较短的读取相比，较大的读取在对齐方面往往有更多的差异。因此，对对齐中的差异数量设定绝对界限是不合适的，因为读取的长度可能会有所不同。边界应该是读取长度的一小部分。因此，我们定义了ε-匹配序列映射问题：

给定基因组数据库G公司主题序列的{S公司₁,S公司₂,…,S公司_我}，一个查询序列（读取）问长度的|问|错误率ε求的子串αG公司，因此E类(α,问)最小值和E类(α,问) ≤ ε|问|.

3算法

大多数序列映射算法将问题分为两个阶段：搜索阶段和对齐阶段。搜索阶段在基因组中找到可能与读数同源的区域。比对阶段验证这些区域，以检查它们是否真的同源。对齐阶段通常比搜索阶段的计算量更大，所花费的时间与搜索阶段中找到的区域数成正比。因此，快速序列映射算法的最佳策略是在对齐阶段之前过滤出尽可能多的候选区域。各种程序都试图通过针对特定问题范围使用精心设计的过滤器来实现这一点。在这一节中，我们描述了AGILE算法，该算法可以通过一些精心设计的过滤器，在广泛的读取长度和错误率范围内实现这种过滤。其中一些过滤器是通用的，适用于整个读取长度和错误率范围，而其他过滤器则适用于特定的子集。然而，组合这些过滤器可以在很大范围内实现更快的速度。

我们首先将基因组划分为非重叠基因组q粒度并存储q粒度在一个q-gram索引（哈希表）。使用q-gram索引，我们找到了q-位在read和基因组之间。每个q-位可以在任意一侧扩展以创建候选区域。相互重叠较大的区域表示相同的路线，因此可以合并在一起。这为我们提供了大量区域。在下面的小节中，我们描述了我们应用于这些区域的AGILE中使用的过滤技术。

3.1连续完美匹配

3.2多重完美匹配

肯特（2002）)已经讨论了使用多个完全匹配作为过滤器的可能性。每个q-hit都是完美的匹配。匹配项不必是连续的。让q个= 16. 例如，考虑两个q-hits-one起始于read中的位置20和基因组中的位置10 020，第二个q-hits-one开始于read的位置52和基因组中位置10 052。由于q-位在阅读和基因组中，它们很容易成为一个比对的一部分。注意，两次读取的“基因组位置”–“读取位置”=10 000。该值称为对角线。我们可以通过保持对角线相等或稍有不同的q-hits数的最小截点来过滤出大量区域。

詹姆斯·肯特（James Kent）讨论了使用两个完美匹配作为筛选标准的想法，而在这里，我们使用更多的完美匹配作为过滤标准。表2显示的值N个它可以用于序列相似性和读取长度的不同值。

3.3忽略高频率的q-grams

在上述两种优化中，我们将理论约束应用于q-位在读取和基因组之间，以筛选出不需要的区域。然而，一些q粒度比其他许多人出现得更频繁。因此，我们得到了大量q-位用于一些阅读。此外，如果特定q粒度在基因组中出现的频率很高，那么它将在不希望出现的地方产生大量匹配，浪费处理时间。A不太频繁q粒度在精确定位匹配时更有用。因此，我们的启发式忽略了所有q粒度频率大于某一截止频率F类以节省时间。为了确保连续的完美匹配条件仍然有效，我们为所有q粒度频率大于F类正如在FANGS中所做的那样。原则上，我们需要减少N个所以我们仍然能够找到所有的同源区域。除非我们考虑到每个模型的准确概率，否则理论模型可能不准确q粒度，这使得模型非常复杂。因此，我们决定使用已经知道正确答案的合成查询进行实证分析。表3显示了读取长度为1000的此类实验的结果。显然，对于固定值N个，较高的值F类应该会导致在对齐阶段处理更多区域，并生成更正确的映射区域。例如，对于N个=1，偶数F类=4将35 670 70个区域传递到对齐阶段。使用时N个=2和F类=17，AGILE正确映射更多读取并过滤出更多区域。对于N个=1，如果F类减少，这将进一步减少正确映射的区域。另一方面，如果F类增加，这将增加所用的时间。比较N个=2箱N个=3箱，即使有F类= 64,N个=3案例正确映射较少的读取次数，但花费更多时间。因此，对于1000的读取长度和10%的错误率，N个=2最有效。使用N个=2和F类=17，AGILE正确映射了99.8%的查询。与相比F类= 17,F类=29需要额外130秒才能将正确映射的读取次数增加2。因此，在花费的时间和正确映射的查询数量之间存在权衡。在这种特殊情况下，N个=2和F类17似乎是一个相对较好的选择。使用类似分析，对于10000的读取长度和10%的错误率，N个=32和F类=4是一个很好的选择。此过滤条件适用于较大的读取长度。

3.4定制q-gram过滤

我们在AGILE中应用的另一个过滤器是在读取和区域之间保持q-hit数的最小截止。这称为q-gram过滤。为了评估q-gra姆过滤的效果，我们使用合成读取进行了实验。对于每次读取，我们计算了q-位以及该区域是否真的同源。如所示图1和2，同源区域往往有更多q-位而非同源区的q-hits数量较少。因此，经过仔细选择q-位，我们可以过滤掉非同源区域。但是，某些查询可能会更频繁地出现q粒度而不是其他人。在这种情况下，出现频率更高的查询q粒度每个区域可以有多个点击。另一方面，查询频率较低q粒度即使在同源区域也只能有几次命中。因此，一般性截止是不合适的。在AGILE中，我们为每次读取动态选择截止点。对于每次读取，我们发现q-位比如说，在所有地区C.我们将截止值保留为分数（f）属于C因此，如果一个区域≥金融公司q-hits，只有这样才能进一步处理。的最佳价值（f）可以过滤掉最大数量的非同源区域，同时保持正确映射的读数的数量相同。例如，表4演示使用的各种值的效果（f）读取长度1000。在这种情况下（f）为0.5。

3.5选择误差率ε

所有上述优化都假设我们已经知道ε，以过滤出区域。然而，无法预先知道特定读取的误差百分比。不同来源的读取可能具有不同的错误百分比。即使从同一来源读取，错误也可能有所不同。为了说明这一点，我们应用了以下启发式方法。

对于每个查询，我们从一个小的ε值开始，比如3%。我们不断增加数值，直到找到匹配项。可以通过多种方式增加价值。我们可以每次增加1%（例如3、4、5等），或者增加一个更大的常数间隔（例如3，8，13，18等），我们可以每次将增量值增加一个（例如，3，4，6，9，13等），每次可以将增量值加倍（例如，三，四，六，十，十八等）。通过我们的实验，我们发现每次将增量值加倍（指数增量）比上述其他策略更快。

选择合适的ε起始值对于上述步骤至关重要，这取决于读数的误差分布。如果读取中的误差分布在频谱较低一侧的一个很小的频带内，那么选择较高的初始值ε将导致在不必要的处理中浪费大量额外的时间。另一方面，如果读数中的误差分布在频谱较高一侧的一个很小的频带上，则选择较小的ε初始值将意味着我们将无法找到与ε的初始少数值相匹配的值。因此，我们将在不必要且无结果的处理中浪费时间。为了“猜测”一个好的初始值n个首先，我们通过学习前面的查询动态调整ε。我们将所有先前查询的单位长度编辑距离的平均值作为下一个查询的ε的初始值。

4敏捷实施

我们的实现将基因组FASTA文件和查询FASTA文档作为输入。查询文件通常包含一批多个查询。在输出中，我们报告了基因组中读取的完整长度匹配的位置，以及匹配的映射质量和分数。

AGILE的高级工作流如所示图3.AGILE利用了454个定序器的错误率非常小的事实(Rothberg和Leamon，2008年). 因此，我们划分了查找读取的最佳匹配α的问题|问|使其最小化E类(α,问)分为多个问题，允许不同的编辑距离值。现在，每个这样的问题的形式都是一个子串αG公司这样的话E类(α,问) < ε |问|. 我们从允许ε的一个小值开始。如果找到匹配项，则输出该匹配项并转至下一次读取。如果找不到匹配项，则增加ε的值并重试。

AGILE使用q-gram索引要映射查询的基因组。我们逐个处理输入读取。由于读取可以来自任何DNA链，AGILE处理读取和反向补码以搜索匹配。对于任何序列，我们首先确定参考基因组中可能与读序列同源的区域。在下一步中，我们使用上一节中讨论的启发式方法过滤掉许多区域。通过创建编辑距离矩阵，使用动态编程对过滤列表中的每个同源区域进行进一步处理，以检查该区域是否实际具有≤ε的编辑距离|问|. 由于编辑距离以ε为界|问|，我们只需要计算宽度为2ε的对角线带|问|矩阵的+1。如果在任何阶段，编辑距离矩阵中的所有路径的编辑距离都大于ε|问|，我们停止进一步处理并丢弃该区域。此修剪策略大大减少了算法所花费的时间，因为在处理前几行后，许多区域被丢弃。表5演示了通过此优化获得的速度效益。

为了找到不同读取长度和ε值的最有效参数值，我们进行了大量实验。基于这些实验，我们为每个查询自动设置这些参数的最佳值。此外，如中所述第3.5节，我们动态选择最合适的ε起始值和增量ε，直到找到匹配。自动选择参数值使用户更容易使用程序，因为他们不需要担心决定合适的参数设置。此外，为不同的场景定制参数设置值，可以运行具有不同长度和错误率的实际查询。

5结果

5.1测绘质量计算

制图质量的概念是由锂等。（2008年a)估计读取序列是否映射到正确位置的概率。Li和Durbin（2010年)将映射质量近似为250(S公司₁负极S公司₂)/S公司₁，其中S公司₁是最佳对齐的分数S公司₂是次佳排列的得分。我们在实验中采用了相同的方法来计算映射质量。

5.2合成数据结果

为了创建合成查询，我们使用了SAMTOOLS包中提供的WGSIM脚本。修改脚本以适应454个数据。我们创建了总计1000万bp不同读取长度的读取。对于每个读取长度，我们引入了2%、5%和10%的错误。总的来说，20%的错误是indels。我们使用AGILE、BWA-SW、BLAT（选项-fastMap）、SSAHA2（选项-454）、Mosaik和PASS将这些模拟读取与人类基因组hg19对齐。由于这些是合成读取，我们知道它们在基因组中的坐标。因此，我们将对齐的坐标与已知坐标进行比较，以计算对齐误差。SSAHA2、BWA-SW和AGILE报告制图质量。但是，如果工具无法拾取最佳和次最佳路线，则相应的映射质量将不正确。因此，比较不同工具报告的映射质量值可能无效，因为某些工具可能会找到比其他工具更多的次优对齐。为了解决这个问题，我们使用所有工具报告的路线计算每个工具的映射质量。对于每个工具，让S公司₁是该工具找到的最佳对齐α的分数。我们接受S公司₂作为任何工具找到的最佳对齐（α除外）的分数。我们使用以下公式计算映射质量S公司₁和S公司₂以与中所述类似的方式第5.1节为了进行评估，我们在Intel Xeon四核E5430 2.66 GHz处理器上运行了所有实验，该处理器具有2×6 MB缓存和32 GB RAM，运行基于Linux的操作系统。每个工具都以单线程模式运行。

表6显示了不同读取长度和测序错误率值的AGILE（版本0.3.0）、BWA-SW（版本0.5.7）、BLAT（版本34）和SSAHA2（版本2.5.1）的CPU时间、自信（映射质量>20）对齐读取的百分比和错误对齐读取的百分比。我们在补充表SI中报告了AGILE与Mosaik和PASS的比较。除非另有必要，否则我们为每个程序使用默认的命令行选项。仔细选择的命令行选项可能会产生更好的结果。

5.2.1敏捷与BWA-SW

从以下位置观察表6AGILE比BWA-SW更准确，尤其是对于长度更短、错误率更高的读取。然而，对于长度更短、错误率更高的相同案例读取，AGILE也比BWA-SW慢。对于较小的错误率，AGILE比BWA-SW快得多（最多快5倍）。AGILE最有效的读取长度约为1000。随着454个读取长度的增加，1000 bp的读取很快将成为标准。

5.2.2敏捷与爆发

发件人表6很明显，AGILE比BLAT快得多（高达30倍），但在大多数情况下仍然更加准确；最重要的是，错误率较大的情况。即使在AGILE不如BLAT准确的情况下，它也不会差多少。BLAT缺乏准确性也是由于“-fastMap”选项造成的。然而，使用默认参数，BLAT比使用“-fastMap”选项慢一个数量级以上。

5.2.3敏捷与SSAHA2

在所有输入上，AGILE比SSAHA2快几倍。此外，除了读取长度为100且错误率为2%的情况外，对于所有输入，AGILE要么与SSAHA2相当，要么更准确。SSAHA2不适合10 kb的读取，因为它不是为这种读取长度设计的，因此无法在其上进行测试。

5.2.4内存比较

就内存使用而言，为了绘制针对人类基因组的读取，BLAT和BWA-SW使用约4GB内存。SSAHA2的峰值内存需求约为5.6 GB。对于AGILEq-gram索引需要大约3.5 GB的内存，基因组本身需要大约3 GB的内存。每个查询所需的内存对于它们来说微不足道。

5.3实际数据结果

由于缺乏基本事实，很难根据实际数据进行评估。然而，如果两个算法为一个读取输出相同的对齐，则很可能是正确的。如果两个对齐器X和Y为一次读取输出不同的对齐，如果X和Y都报告低映射质量，则对齐不明确，无论哪一个错误都无所谓。如果X报告读取的映射质量很高，而X对齐分数低于或略高于Y，则X映射质量是错误的，X没有意识到同样可能的对齐。此分析与中报告的分析类似Li和Durbin（2010年).

我们对从NCBI短读存档SRR005010获得的实际查询运行了BWA-SW和AGILE。我们过滤了读取集以删除长度小于100 bp的查询，并统一选择平均长度为338 bp的100k读取。AGILE用了283 CPU秒，BWA-SW用了939 CPU秒。这两种工具都发现96 293条常见路线。总共有357个读取未被任何工具映射。表7显示了仅由一个对齐器对齐或映射到不同位置的读取的分解，并且为BWA-SW或AGILE分配的映射质量≥20。总的来说，BWA-SW错过了AGILE很好绘制的203+247=450条路线。AGILE错过了95条（46+49）由BWA-SW很好对齐的对齐。请注意，即使整个读取集的平均读取长度为338，但不一致主要出现在较小的读取长度（平均长度286或更小）的情况下。这与合成读数的结果一致，因为两种校准器在校准较短长度的读数时往往会犯更多错误。BWA-SW在更短的读取长度范围（140-170）中往往会丢失更多对齐，而AGILE在稍长的读取长度（250-286）中会丢失更多的对齐。

6分析

AGILE与BLAT、SSAHA2、FANGS、Mosaik、PASS和BWA-SW在共享种子和扩展策略以获得候选区域方面相似。然而，主要区别在于使用启发式方法来减少要处理的区域数量。BLAT和SSAHA2将短（10–15 bp）精确匹配视为种子。BLAT还提供了使用短的不精确匹配（允许一个不匹配）或两个对角线值略有不同的精确匹配的功能。对于较长的读取，这两种技术都会产生大量候选区域。AGILE通过允许更长的种子（最多28 bp）来筛选候选对象。它还使用最小连续精确种子匹配数的截断，以及对角线值相等或略有不同的最小精确匹配数的较大截断。据我们所知，没有其他工具使用多个完美匹配过滤标准，其截止值超过两个匹配。这类似于使用长间隙的种子。由于类似的原因，BWA-SW也使用长间隙种子。BWA-SW和AGILE的主要区别在于实现长间隙种子的方式。虽然BWA-SW使用前缀Trie和前缀DAWG，但AGILE依赖于更简单的数据结构q-gram索引（散列表）和对角线坐标。BWA-SW和AGILE都进一步使用不同的启发式，以减少搜索空间。

AGILE使用的启发式算法与字符串匹配或序列映射算法中使用的类似。对于每一个启发式，AGILE都会将其应用于长读映射的问题。定制q-gram过滤就是一个例子。Q-gram滤波是一种众所周知的过滤不需要区域的技术。然而，据我们所知，没有任何算法使用q-gram过滤，并为每次读取定制了截止值。AGILE的另一个主要贡献是在找到匹配之前，允许的错误率逐渐增加。大多数工具使用允许的错误率来决定修剪搜索空间的阈值。然而，许多读取都有少量错误。由于这些工具使用的是允许错误率的固定值，因此它们最终会对所有读取使用更大的允许错误率值，以同时映射具有更大错误的读取，从而导致效率损失。因此，逐渐增加映射中允许的错误率会对这些工具的效率产生巨大影响。

7结论

测序技术的进步要求开发高性能、可扩展的算法，以从这些数据集中提取生物相关信息。开发序列映射算法的研究主要集中在映射短读，而对于较长的454读，几乎没有做什么工作。AGILE是一种基于散列的序列映射算法，它使用高效的启发式算法来优化映射过程的不同步骤，快速映射长读取。AGILE可以处理非常大的基因组大小和读取长度。它非常灵活，因为它允许在映射中进行大量不匹配和插入/删除，并提供命令行参数来控制映射过程的每个步骤。当读取时间较长且错误率较小时，AGILE的灵敏度和特异性最佳。考虑到随着测序技术的改进，读取长度将进一步增加，错误率将降低，AGILE在未来将更加有用。

资金：本研究部分得到了国家科学基金会奖号CCF-0621443、SDCI OCI-0724599、CNS-0551639、IIS-0536994和HECURA-0938000的支持。这项工作也得到了能源部拨款DE-FC02-07ER25808和DE-FG02-08ER25848的部分支持。

利益冲突：未声明。

参考文献

作者注释