用于生产MegaBLAST搜索的数据库索引

摘要

1简介

BLAST（高等学校等。,1997)是一个程序包，用于将查询生物序列与数据库进行匹配，并识别与查询的一部分具有统计意义的局部对齐的数据库序列。使BLAST快速的主要启发式方法之一是对查询进行预处理，以构建一个查找表，该表随后用于查找短的高分未映射对齐，这里称为“种子”。种子可以扩展到更长的未映射路线，然后在BLAST项目的后期阶段扩展到全间隙路线。种子的默认长度是蛋白质的3个氨基酸，用BLASTN模块进行高灵敏度核苷酸搜索的11个核苷酸，用MegaBLAST模块进行低灵敏度核苷酸搜索时的28个核苷酸（Zhang等。,2000). BLAST 2.0版中引入的一项重要创新是要求附近存在两个种子，以减少通过种子阶段的候选匹配数（Altschul等。,1997). 对于蛋白质，最近的研究表明，需要一个长度为6的连续种子和更复杂的预处理可以获得更好的性能（Shiryev等。,2007)。

而不是通过搜索派生自查询，可以通过搜索派生自数据库.大量研究，最近在江审查等。(2007)建议对数据库进行预处理并搜索数据库派生的数据结构可以产生更快的搜索时间。

对数据库进行预处理以进行生物序列比较的两个广泛使用的软件包是SSAHA（Ning）等。,2001)和BLAT（肯特，2002). 然而，这两个包都被定义为只找到几乎相同的匹配[例如，BLAT和BLAST的TBLASTN模块的比较，请参见Gertz等。(2006)]. 我们唯一能找到的软件包是：

1. 产生与BLAST的某些模块类似或相同的结果，
2. 预处理数据库而不是查询，为种子搜索阶段构建数据结构，
3. 为数据结构和数据库使用相当数量的内存（排除后缀树解决方案），以及
4. 具有当前在网络上免费提供的可执行文件或源代码

是miBLAST（Kim等。,2005). miBLAST是为100个基数以下的短查询设计的，我们确认，对于较长的查询，其相对性能确实会下降。下面，我们将介绍一些方法，例如miBLAST，这些方法通过集合术语“数据库索引”从数据库构建搜索结构。

由于数据库索引似乎很有前途，但NCBI BLAST的任何模块都没有使用这种范式，因此我们着手设计一种替代（至少）NCBI PLAST的一个模块，并评估是否可以在生产使用中实现声称的性能优势。NCBI的BLAST Web服务中最常用的搜索数据库是数量，这是一个全面的“非冗余”序列集合。索引数量与其他数据库相比，数据库带来了特定的挑战，因为数量很大，每天都会更新。用于网络BLAST搜索的其他流行数据库是人类和小鼠基因组，在较小程度上是其他脊椎动物基因组。

部分由于2007年NCBI的BLAST网页的重新设计，人类和小鼠的单基因组数据库越来越多地被用户选择，他们每周分别提交约10000和3000个人类和小鼠数据库查询。当前的默认设置是使用MegaBLAST搜索基因组特定核苷酸数据库（Zhang等。,2000)模块，从程序内调用序列比对分析从NCBI C++工具包编译而成。

因此，我们着手开发一种新版本的MegaBLAST，它可以使用从基因组特异性数据库中获得的数据结构来搜索初始种子。我们还需要构建一个程序，在这里称为mk索引名，构建索引结构。“索引的MegaBLAST”的基本目标是，它应该找到与我们开始使用的版本相同的匹配序列和对齐，并且对实践中出现的大多数查询执行得更快。与miBLAST不同，索引MegaBLAST的性能应该随着查询时间的延长或匹配数量的增加而适度下降。

在项目早期，我们意识到，对良好性能的一个基本障碍是查询与基因组中多次出现微小变异的DNA序列很好地对齐。为了解决这个困难，我们开发了软件包WindowMasker（Morgulis等。,2006年a)快速屏蔽频繁重复的序列而不使用重复库。我们意识到，大多数用户很少希望看到BLAST的重复匹配输出，因此我们修改了MegaBLAST，使其能够搜索“soft-masked”数据库。Soft-masked意味着种子不能与遮罩间隔相交，但对齐可以延伸到遮罩间隔并穿过遮罩间隔。MegaBLAST中早就提供了查询的软屏蔽功能。

第2节和三和补充材料描述了我们如何设计和测试索引MegaBLAST，以记录潜在的性能改进。关键测试是NCBI的web BLAST服务的用户将如何响应。索引MegaBLAST于2007年10月部署，用于查询人类和小鼠基因组的BLAST核苷酸搜索变体（单独），并在定期发行说明中宣布。默认情况下，使用WindowMasked版本的基因组，但用户可以关闭屏蔽。用户投诉为零，只有一项一般性调查表明索引MegaBLAST在生产中运行良好。

2方法

在本节和补充材料中，我们描述了用于组织数据库索引的数据结构，以及用于在查询和数据库中查找初始相同子字符串的“种子搜索算法”。我们还描述了测试策略以及索引MegaBLAST如何投入生产使用。如下文所述，种子搜索算法已并入NCBI的MegaBLAST的修改版本中，对后续处理种子的算法代码部分几乎没有更改。

2.1指标结构

MegaBLAST数据库索引包含压缩序列数据以及k个-默斯。除了长度k个，另外两个重要参数是最小种子长度w个(≥k个)和步幅秒用于在数据库中移动。确保每个长度的精确匹配w个在查询和主题之间，我们使用关系秒=w个−k个+1.关于k个-mer结束于秒-第个位置存储在索引中，它满足下面指定的条件，确保我们可以找到所有种子。当前实现默认为w个=16,k个=12和秒=5.

数据库索引由几个称为索引卷的文件组成，每个文件对应于底层数据库中的连续序列范围。索引卷文件包含三个部分：标头、序列数据和偏移数据，如中所示图1数据库中核苷酸的数量表示为n个。在补充资料中，我们得出了将索引存储为以下形式所需字节数的估计值：

(1)

字节。例如，对于大小为1GB的未屏蔽数据库和我们的默认值k个=12,秒=5，则指数的估计大小为~1.175 GB。

在单独的窗口中打开

图1。

用于数据库索引的数据结构示意图。

2.1.1标题

索引卷的标头部分主要存储由该卷索引的基础数据库中的序列范围。标题部分还存储了一些用于内务管理的辅助值，例如索引格式版本。在代码附带的外部文档中可以找到索引格式当前版本的确切定义。

2.1.2序列数据

序列数据以压缩格式存储，每个序列基使用NCBI2NA编码2位。在这种编码中，基A、C、G和T被相应地编码为整数0、1、2和3。包含不明确字符的位置编码为0。下一小节中描述的种子搜索算法从不涉及序列中的模糊部分。这些部分仅用作占位符。

种子搜索算法与逻辑序列大小大致一致。这使得映射k个-mer定位到相应序列更有效。如下文所述，逻辑序列可以通过拆分长序列或组合几个短序列来形成。逻辑序列和实际数据库序列之间的映射存储在索引卷的序列数据部分中。

为参与给定索引卷的所有逻辑序列分配一个连续整数序列id。在索引的MegaBLAST和NCBI BLAST软件的其他变体中，长数据库序列被分割为重叠序列大块优化处理。从现在开始，由输入序列我们指的是数据库序列或由这种拆分产生的块，这使我们能够避免将算法描述与处理长序列的细节混淆。

输入序列到逻辑序列的映射。让Ş是BLAST代码中定义的块大小（当前默认值为5 Mbase），并让B=⌈log₂(Ş）。由于分块，我们可以假设所有输入序列最多都是长度ℓ. 只要连续输入序列的组合长度不超过ℓ. 逻辑序列被分配为连续整数逻辑序列ID从0开始。

索引卷的序列数据部分包含一个表，该表将每个逻辑序列id映射到组成它的输入序列。表中的每个元素包含四个整数：

1. 逻辑序列中第一输入序列的序列id；
2. 逻辑序列中最后一个输入序列的序列id；
3. 此逻辑序列的开始与序列存储开始的偏移量（字节），以及
4. 此逻辑序列结束时与序列存储开始时的偏移量（字节）。

序列存储。序列存储是包含压缩序列数据的索引卷的序列数据部分。逻辑序列按增加逻辑序列id的顺序存储在那里。

逻辑序列数据包含级联编码的组成输入序列。如果需要，使用每基2位的输入序列编码在末尾用额外的0位填充，以确保它们以字节边界结束。索引卷存储所有输入序列的真实长度，因此种子搜索算法不会将填充的0位视为核苷酸编码的一部分。

2.1.3偏移数据

索引卷的偏移数据部分由两部分组成：查找表和偏移列表。查找表包含4个^k个条目。这个我-th条目包含指向数据库中名为偏移量列表，其中k个-带有值的mer我发生。

偏移列表位于数据结构中查找表的后面。单个列表由哨兵4字节字编码0终止。我们的偏移量编码方法确保0不是列表项的有效值。的位置k个-mer可以按照其在基础数据库中出现的顺序添加到偏移列表中，但实际上，为了更快地进行搜索，会按照补充材料中的描述排列顺序。

要确定哪些k个-要存储的mer位置种子合格区间为每个输入序列计算。一个种子合格区间是输入序列的子序列，不包含歧义或小写（屏蔽）字母。一个k个-如果以下三个条件成立，则mer存储在索引中：

1. 这个k个-mer完全包含在种子合格区间内；
2. 种子合格间隔至少为w个核苷酸长度和
3. 的最后一个字母的偏移量k个-mer相对于相应逻辑序列的开始是可被步长整除的，秒.

让秒₂=2^{⌈日志₂秒⌉}; 即2的最小幂大于或等于秒。对于k个-mer实际值o（o）存储在偏移列表中的按以下方式计算：

哪里小时是逻辑序列id，⌈2^B/秒⌉是用于编码k个-mer在其逻辑序列中从容不迫秒和第页是的最后一个字母的偏移量k个-mer相对于逻辑序列的开始。o（o）存储为32位宽无符号整数。

特殊偏移量和前缀。对于一些人k个-mer实例中，额外信息与其编码位置一起存储。让d日_我(d日_第页)是距第一个（最后一个）底面的距离k个-mer到相应有效间隔的左（右）端。我们称之为k个-mer实例特殊的如果有d日_我<秒，或d日_第页<秒或两者兼而有之。我们定义u个_我(u个_第页)等于d日_我(d日_第页)如果d日_我<秒(d日_第页<秒)否则为0。

种子搜索算法在偏移列表中使用特殊前缀，以确保报告的种子永远不会跨越种子合格间隔的边界。对于特殊情况k个-mer实例，一个额外的32位前缀值o（o）_秒=秒₂u个_我+u个_第页存储在其编码位置之前的偏移量列表中o（o）通常情况下o（o）_秒<秒₂²和o（o）≥秒₂²，因此前缀很容易与实际编码位置区分开来。

偏移条目示例。考虑来自掩蔽的人类基因组的第一个重叠群NT_019273.18的以下子序列，起始于偏移8 500 000（第一个碱基位于偏移0）：gagaggACAACACTTAAAGGTTCAACTAGCAATA，（默认）值为k个=12,秒=5，块大小为5 Mb，块重叠为100 bp，输入中偏移量为8500000的基数位于第二个逻辑序列中，其逻辑块中偏移量达3500 100。偏移量3500 100（gagaggACACA）和3500 105（gACACTCTAA）处的12 mers部分被掩盖。因此，它们不属于任何符合种子条件的间隔，也不会添加到索引结构中的查找表偏移量列表中。考虑添加到查找表中的下一个12 mer位于偏移量3 500 110处（CACTTAAAGGTT），并且处于种子合格间隔中。这也是一个特殊的偏移量，因为到其符合种子条件的区间边界的左侧距离只有四个。偏移量3500 115处的12 mer（AAAGGTTCAACT）也处于种子合格区间，但它不是特殊偏移量。

上述子序列中CACTTAAAAGGTT的计算结果导致屏蔽人类基因组索引结构中的条目为：（i）使用四个碱基的位编码的查找表条目索引4702 383，（ii）使用公式2的特殊偏移值32^⌈日志₂5⌉·4+0，以及（iii）使用公式（2^{⌈日志₂5 ⌉)2}+ 3 500 110/5 + 1 · ((2^⌈日志₂5 ⌉)²+ ⌈ 2^⌈日志₂(5 000 000)⌉/5 ⌉ + 1). 在查找表条目索引4 702 383的偏移列表中，特殊偏移值32后面是偏移值2 377 873。

3结果

为了使BLAST数据库索引在生产中可行，第一个要求是索引数据结构使用的内存量应该合理。表1显示了索引卷相对于实际数据库长度的大小。这些数据用于中所述的定时测试数据库第2.2.3小节和2.2.4以及更短的果蝇属基因组。索引是使用创建的秒=5和k个=12.表的最后一列显示了公式1提供的上限。对于未屏蔽的数据库，这个界限相当紧。在遮罩的情况下，它不太严密，因为它没有考虑到未遮罩基的实际数量。

第二个要求是，创建索引数据结构所需的预处理时间应该合理。我们测量了创建人类未屏蔽指数、人类屏蔽指数、小鼠未屏蔽指数和小鼠屏蔽指数所需的时间，分别为1319.11s、1000.67s、1176.56s和911.27s。对于每个数据库，前一句中的时间是数据库所有卷的索引卷创建时间之和。

已经表明，用于数据库索引的数据结构可以在合理的内存和时间内创建，我们继续评估搜索时间，特别是查找种子的时间。索引种子搜索的性能优劣取决于种子候选列表的大小（参见补充材料），我们指的是使用索引数据库的BLAST的种子查找阶段。当数据库（或来自相应基因组的查询）被重复和低复杂性区域屏蔽时，种子候选列表可能会很短。然而，对于无掩码数据库和查询，如果查询包含高度重复的区域，与无索引情况相比，索引搜索实际上可能具有更差的性能。

3.3生产环境中索引和非索引MegaBLAST的比较

图2显示生产设置中100个索引搜索和100个非索引搜索的运行时间，但用于此测试的计算集群没有用于其他任务。报告的时间是针对每个查询和每个方法进行的两次搜索中的第二次搜索。索引搜索的平均墙锁时间为0.5秒，非索引搜索的为0.8秒。第75个百分位时间提供了更好的用户体验度量；索引搜索为0.38秒，非索引搜索为1.37秒。根据查询大小的对数绘制时间，以表明索引搜索对于较短的查询通常更快，而非索引搜索对于较长的查询更快。

在单独的窗口中打开

图2。

生产设置中100个索引和非索引搜索的墙锁时间是查询长度对数的函数。考虑到如果时间在0.01秒以内，那么索引搜索会快75倍，无索引搜索会更快19倍，并且它们会绑定6个查询。索引搜索在较短的查询上更快，在最长的查询上较慢。

由于无索引搜索使用10-20个节点，而索引搜索使用一个节点，因此即使无索引墙锁时间较短，无索引搜索也可能比索引搜索使用更多的总CPU时间。只要用于索引搜索的专用机器能够跟上传入的对人类和小鼠基因组的MegaBLAST查询，用户就会感觉到基因组特定查询的响应时间要快得多，这是因为更快的运行和对其他数据库的查询没有竞争的综合影响。量化所有类型查询对吞吐量的好处是很困难的，因为专用和共享机器的工作负载随时间变化。

4讨论

我们介绍了MegaBLAST的种子搜索阶段的新实现（Zhang等。,2000)通过搜索索引结构找到种子k个-mers来源于对数据库的预处理。我们发现，这种“索引MegaBLAST”比“基线MegaBLAST”更快，后者在大多数情况下对查询进行预处理，尤其是对屏蔽数据库。当索引的MegaBLAST由于种子太多而速度较慢时，性能下降的程度足够有限，因此代码可以用于生产。索引MegaBLAST满足输出与基线MegaBLAS相同或几乎相同的要求，只要所有选项设置一致。此外，预处理器程序的使用mk索引名正在NCBI进行其他项目的测试，例如将“下一代测序技术”生成的短序列与参考基因组进行比对。

开展这一项目受到了几篇研究论文的积极影响，这些论文表明，对数据库进行预处理以进行初始查找比预处理查询可以获得更快的搜索时间。我们决定将精力投入到相对简单的数据结构设计、MegaBLAST模块的种子搜索阶段以及对屏蔽数据库的支持中，以便在生产BLAST中更快地使用数据库索引。现在，这一总体目标已经实现，有许多可能性可以扩展和改进我们的实施。

首先，我们的实现可以通过改变步幅和k个-摩尔长度。如后文所述，为了获得良好的性能，可能需要长列表的一些替代数据结构。由于字母表的大小增加，将我们的工作扩展到蛋白质搜索将更具挑战性。Williams和Zobel（2002）实现了一个核苷酸和蛋白质搜索的研究原型，称为CAFE，他们在其中将数据库预处理为一个反向索引。他们在CAFE中表现出了令人印象深刻的性能，但为了获得这一性能，他们大大改变了在同源搜索阶段使用的搜索方法和评分，而不仅仅是最初的单词查找。因此，制作一个替换一些BLAST模块的生产版本将比我们进行的种子搜索阶段的替换更困难。

我们的实现有一些可能的算法改进，可以在不进行重大重新设计的情况下进行尝试。我们总结了其中两种可能的变化。而不是使用统一步幅（默认情况下秒=5）选择时k个-在索引中出现的次数较多的情况下，可以通过不均匀的增量来跨越，这样每个潜在的种子（默认长度28）都可以用一个或多个来充分表示k个-mer子串。可以选择每个位置的步幅，以便k个-出现在索引中的mer尽可能少。这将使k个-mer列出的列表更长，但加快了种子搜索阶段。对人类基因组的初步评估表明，使用不均匀步幅可以同时增加k个-mer列表减少9%，加权平均长度缩短25%加权平均数”是指每个k个-mer根据其在基因组中的出现次数进行加权。使用非均匀步长会使预处理花费更长的时间，但我们不认为数据库预处理时间是性能的主要衡量标准，因为预处理是在MegaBLAST外部完成的。另一个可能的改进是k个-mers可以被更复杂的结构所取代，例如散列表或某种树（吉拉迪等。,2002)。

数据库索引研究实施建议的可能改进包括以下方面，这些改进需要进行重大重新设计。使用“问-gram filtering”，它考虑k个-动态规划对齐矩阵区域内的mer匹配可以为候选匹配提供比单个种子（拉斯穆森等。,2006). 数据库序列可以转换为Lee提出的数据特征向量等。(2007). 本着与k个-我们在已经实现的优化（补充材料）下描述的mer重新排序，“pier”的一些子集k个-mer可以存储在Cao中提出的更小、更快速访问的数据结构中等。(2004). 最后，在MICA中实现的序列表示（Stokes和Glick，2006)可以采用节省空间的方式存储较小值的数据库索引k个[斯托克斯和格利克(2006)状态k个=7作为有效上限k个]. 此类索引可用于BLASTN目前提供的更敏感的搜索。

总之，自2007年10月以来，NCBI的网络BLAST服务的一部分一直在使用数据库索引，而不是种子搜索阶段的查询索引。索引MegaBLAST的生产使用和本文中的测试验证了数据库索引所声称的性能改进。我们的实施可能需要进一步改进，可能需要进行重大的重新设计。

补充材料

致谢

参考文献