CN105138534B - 基于fmd索引和快表的跨越式种子查找算法 - Google Patents

Abstract

本发明公开了一种基于FMD索引和快表的跨越式种子查找算法，包括下述步骤：S0、构建数据库的FMD索引以及快表；S1、从快表中取出查询序列中长度为k的子序列的双区间；S2、这个步骤通过向后搜索算法，逐步找出k种子左边的匹配区域；S3、对步骤S2中缩小前的区间执行向前搜索算法，以找出k种子右边的匹配区域；S4、检查当前检测位置是否位于查询序列的尾部，如果是，则算法终止，否则，执行步骤S5；S5、将当前检测位置向前跳跃w‑k+1个位置，重复执行步骤S2‑S5，其中w是要查找的种子的长度。本发明提出的快表具有占用空间少、访问效率高的特点；在快表和FMD索引的基础上，本发明提出的种子查找算法能够快速地找出所有w种子的双区间。

Description

基于FMD索引和快表的跨越式种子查找算法

技术领域

本发明涉及生物序列数据库搜索的技术领域，特别涉及一种基于FMD索引和快表的跨越式种子查找算法。

背景技术

要搜索在生物序列数据库和查询序列之间长度至少为w的完全匹配的种子，现有的解决方案可以归为三类。

第一类算法对查询序列构建一个快表。该快表是一个Hash表，每个表项是一个线性链表，记录了一个长度为k的序列在查询序列中出现的所有位置。这类算法接着对数据库序列进行跨越式的扫描。跨越式的扫描指的是每隔w-k+1个位置就检测一次长度为k的子序列。检测过程包括通过快表找出这个子序列在查询序列中出现的位置，每个位置对应一个长度为k的种子，然后比较每个k种子的左右两边的区域，检查这个k种子是否包含在一个w种子种。这种种子查找算法被应用于MegaBLAST中。

第二类算法对数据库构建快表，这类快表也是一个Hash表，每个表项对应一个长度为k的短序列，典型的解决方案有两种。第一种解决方案以Indexed MegaBLAST为代表，对数据库每隔w-k+1个位置取出一个长度为k的子序列，将其位置加入相应的Hash表项中。算法然后检测查询序列中所有长度为k的子序列，找到k种子，然后通过方案一的做法，检查这些k种子是否包含在w种子中。第二种解决方案以BLAT为代表，取出数据库序列中所有互不重叠的长度为w的子序列，将其位置记录在相应的Hash表项中，然后检查查询序列的所有长度为w的子序列，在Hash表中找出其位置链表，链表中的每个位置对应一个w种子。

第三类解决方案对数据库建立FM索引或FMD索引，利用该索引查找长度至少为w的最大匹配区域。最大匹配区域指的是不能继续向左右两边延伸的一段完全匹配区域。序列比对软件Cushaw2用FM索引，从查询序列的第一个字符开始，每一步从右边添加一个字符，直到搜索结果为空集为止。算法接着从上一步停止的地方开始，继续上述过程。序列比对软件BWA-MEM利用FMD查找超级最大匹配，超级最大匹配也是最大匹配，但是它们在查询序列上的区段不能被其他最大匹配在查询序列上的区段覆盖。

上述第一种方案属于最初被提出的方案，也是这三类解决方案中效率最低的种子查找方案。第二种方案中，Indexed MegaBLAST在小型数据库和短查询序列上可以运行得很快，然而在大型数据库和长查询序列上，快表将会变得非常庞大，而且性能会急剧下降，以至于比MegaBLAST还要低效。BLAT不是精确的算法，它无法保证找出所有w种子。第三类解决方案虽然性能好，然而它们同样不能保证能够找出所有长度至少为w的种子，从而导致最终搜索精度的下降。

发明内容

本发明的主要目的在于克服现有技术的缺点与不足，提供一种基于FMD索引和快表的跨越式种子查找算法。

为了达到上述目的，本发明采用以下技术方案：

基于FMD索引和快表的跨越式种子查找算法，包括下述步骤：

S0、构建数据库FMD索引的快表，这个快表是一个Hash表，每一个表项对应一个长度为k的短序列，保存的是在FMD索引中搜索这个短序列所得到的双区间；

S1、计算查询序列中k子序列的Hash值，并从快表中取出其相应的长度为k的种子(简称K种子)的双区间；

S2、这个步骤通过向后搜索算法，逐步扩大k种子左边的匹配区域；

S3、对步骤S2中缩小前的区间执行向前搜索算法，以找出k种子右边的匹配区域；

S4、检查当前检测位置是否位于查询序列的尾部，如果是，则算法终止，否则，执行步骤S5；

S5、将当前检测位置向前跳跃w-k+1个位置，重复执行步骤S2-S5，其中w是要查找的种子的长度。

优选的，FMD索引具体为：

核酸序列的一个长度为k的子序列称为k子序列，查询序列与数据库序列之间的一段长度为w的完全匹配区段称为w种子，在数据库的FMD索引中搜索序列P，其搜索结果以双区间的形式表示，而双区间则由三个整数表示，给定核酸序列P的双区间和字符a，a为A、C、G、T中的其中一个，由向后搜索算法，可以得到aP的双区间；由向前搜索算法，可以得到Pa的双区间，双区间中元素的个数称为该区间的大小，它表示P在数据库中出现的次数，如果P的双区间为空区间，则表示P没有在数据库中出现。

优选的，所述快表的大小是固定的，每个字符只能是A、C、G、T中的一个，故该快表一个有4^k个表项。

优选的，所述步骤S2中，在向后搜索过程中，如果出现双区间缩小或者双区间为空的情况，就表明某些k种子向左扩展时遇到了不匹配的字符对，对缩小前的区间，算法还需要通过向前搜索算法找出其对应的种子右边的匹配部分，试图找出长度为w的种子。

优选的，所述步骤S3中，在向前搜索过程中，如果出现搜索区间为空的情况，就说明查询序列中的这段区域在数据库中不存在；否则的话，将得到查询序列中w子序列的双区间，算法将其作为结果输出。

优选的，所述步骤S5中，并不需要检测查询序列中所有的k子序列，相反，它只需每隔w-k+1个位置检测一次即可。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明提出的快表与现有的基于线性链表的快表的区别之处在于，基于线性链表的快表的大小会随着数据库的增长而增加，在大型数据库中，这种快表将占用非常大的存储空间，长的线性链表也会使得种子查找过程非常耗时。与基于FM索引构建的快表相比，本发明提出的快表保存的是短序列的双区间，而基于FM索引构建的快表保存的则是后缀数组区间。

2、本发明提出的种子查找算法与现有的算法相比，具有精度高、效率快的优势。传统的种子查找算法有的无法找出所有的w种子，无法保证搜索精度。有的虽然也采用了跨越式的扫描方式，但在大型数据库和长的查询序列下性能低下，它们只能逐个检查k种子是否包含在w种子中。

3、本发明提出的种子查找算法同样采用跨越式的扫描方式，但是它能够一次检查一批k种子是否包含在w种子中，这使得它具很好的性能优势。

附图说明

图1是本发明基于FMD索引和快表的跨越式种子查找算法的流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

本实施例主要包括两个部分：基于FMD索引的快表和种子查找算法。构建快表针对的是上述背景技术中第二类算法所采用的快表的缺点。在这类算法中，快表的每个表项都是一个线性链表，这使得快表在大型数据库中将会占用非常巨大的存储空间。一些算法例如BLAT等，为了减少快表的大小从而采用索引部分短序列的方式，这会牺牲搜索精度。本发明的快表每个表项均保存一个双区间，每个双区间只需由三个数字表示，因此快表的大小是固定的，不会随着数据库的增长而发生变化。

种子查找算法针对的是上述第二、第三类算法中的缺点。这类算法中，IndexedMegaBLAST能够找出所有的w种子，但其在大型数据库下的效率不高，因为这时候线性链表会非常长，从而需要频繁地检查k种子是否包含在w种子种。其他算法为了提高效率从而采用了牺牲精度的方式，它们只能够找出部分的w种子。本发明的种子查找算法采用了第一类算法中的跨越式扫描方式，能够检查一批k种子是否包含在w种子中，这使得该算法精确的同时具有非常高的执行效率。

FMD索引是bidirectional FM-index的缩写，FM是提出FM索引的两个作者Ferragina Paolo和Manzini Giovanni的名字缩写。

在详细介绍这两个部分内容之前，本实施例先说明与FMD索引相关的信息。核酸序列的一个长度为k的子序列称为k子序列。查询序列与数据库序列之间的一段长度为w的完全匹配区段称为w种子。在数据库的FMD索引中搜索序列P，其搜索结果以双区间的形式表示，而双区间则由三个整数表示。给定核酸序列P的双区间和字符a(a为A、C、G、T中的其中一个)，由向后搜索算法，可以得到aP的双区间；由向前搜索算法，可以得到Pa的双区间。双区间中元素的个数称为该区间的大小，它表示P在数据库中出现的次数。如果P的双区间为空区间，则表示P没有在数据库中出现。

本发明的基本流程如下所述：

S0、构建数据库的FMD索引以及快表；

S1、计算查询序列中长度为k的子序列的Hash值，并从快表中取出其相应的长度为k的种子(简称为k种子)的双区间；

S2、这个步骤通过向后搜索算法，逐步找出k种子左边的匹配区域；

S3、对步骤2中缩小前的区间执行向前搜索算法，以找出k种子右边的匹配区域；

S4、检查当前检测位置是否位于查询序列的尾部，如果是，则算法终止，否则，执行步骤S5；

S5、将当前检测位置向前跳跃w-k+1个位置，重复执行步骤S2-S5，其中w是要查找的种子的长度。

本发明中的快表是在核酸序列数据库的FMD索引的基础上构建的。它也是一个Hash表，每一个表项对应一个长度为k的短序列，保存的是在FMD索引中搜索这个短序列所得到的双区间。这种快表的大小是固定的。由于每个字符只能是A、C、G、T中的一个，因此快表一个有4^k个表项。利用这个快表，可以立刻得到查询序列中k子序列的双区间。

算法的第二部份是种子查找算法，其流程如图1所示，该算法以查询序列最开始的k子序列开始，逐步执行5个主要步骤。

算法的第一个步骤是图1中的S1部分，它计算k子序列的Hash值，并从快表中取出其相应的双区间。利用快表，不仅可以实现跨越式的扫描，还可以一次性得到k子序列的双区间，而不需要逐步利用向前或向后搜索算法，这节省了大量的时间。

算法的第二个主要步骤在图1中是S2部分。这个步骤通过向后搜索算法，逐步找出k种子左边的匹配区域。在向后搜索过程中，如果出现双区间缩小或者双区间为空的情况，就表明某些k种子向左扩展时遇到了不匹配的字符对。对缩小前的区间，算法还需要找出其对应的种子右边的匹配部分，试图找出长度为w的种子。

算法的第三个主要步骤对应图1中的S3部分，它对步骤2中缩小前的区间执行向前搜索算法，以找出k种子右边的匹配区域。在向前搜索过程中，如果出现搜索区间为空的情况，就说明查询序列中的这段区域在数据库中不存在；否则的话，将得到查询序列中w子序列的双区间，算法将其作为结果输出。

算法的第四个步骤为图1中的S4部分，这个步骤检查当前检测位置是否位于查询序列的尾部。如果是，则算法终止，否则的话，执行步骤5。

算法的第五个步骤在图1中是S5部分，它将当前检测位置向前跳跃w-k+1个位置，重复执行步骤二至步骤五，这就是所谓的跨越式扫描。跨越式扫描并不需要检测查询序列中所有的k子序列，相反，它只需每隔w-k+1个位置检测一次即可。这种跨越式的扫描方式在保证能够找到所有w种子的同时，也具有非常高的效率。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.基于FMD索引和快表的跨越式种子查找方法，包括下述步骤：

S0、构建数据库FMD索引的快表，这个快表是一个Hash表，每一个表项对应一个长度为k的子序列，保存的是在FMD索引中搜索这个子序列所得到的双区间；FMD索引具体为：

核酸序列的一个长度为k的子序列称为k子序列，查询序列与数据库序列之间的一段长度为w的完全匹配区段称为w种子，在数据库的FMD索引中搜索序列P，其搜索结果以双区间的形式表示，而双区间则由三个整数表示，给定核酸序列P的双区间和字符a，a为A、C、G、T中的其中一个，由向后搜索算法，可以得到aP的双区间；由向前搜索算法，可以得到Pa的双区间，双区间中元素的个数称为该区间的大小，它表示P在数据库中出现的次数，如果P的双区间为空区间，则表示P没有在数据库中出现；

S1、计算查询序列中k子序列的Hash值，并从快表中取出其相应的长度为k的种子的双区间；

S2、这个步骤通过向后搜索算法，逐步扩大k种子左边的匹配区域；

S3、对步骤S2中缩小前的区间执行向前搜索算法，以找出k种子右边的匹配区域；

S4、检查当前检测位置是否位于查询序列的尾部，如果是，则算法终止，否则，执行步骤S5；

S5、将当前检测位置向前跳跃w-k+1个位置，重复执行步骤S2-S5，其中w是要查找的种子的长度。

2.根据权利要求1所述的基于FMD索引和快表的跨越式种子查找方法，其特征在于，所述快表的大小是固定的，每个字符只能是A、C、G、T中的一个，故该快表一个有4^k个表项。

3.根据权利要求1所述的基于FMD索引和快表的跨越式种子查找方法，其特征在于，所述步骤S2中，在向后搜索过程中，如果出现双区间缩小或者双区间为空的情况，就表明某些k种子向左扩展时遇到了不匹配的字符对，对缩小前的区间，算法还需要通过向前搜索算法找出其对应的种子右边的匹配部分，试图找出长度为w的种子。

4.根据权利要求1所述的基于FMD索引和快表的跨越式种子查找方法，其特征在于，所述步骤S3中，在向前搜索过程中，如果出现搜索区间为空的情况，就说明查询序列中的这段区域在数据库中不存在；否则的话，将得到查询序列中w子序列的双区间，算法将其作为结果输出。

5.根据权利要求1所述的基于FMD索引和快表的跨越式种子查找方法，其特征在于，所述步骤S5中，并不需要检测查询序列中所有的k子序列，相反，它只需每隔w-k+1个位置检测一次即可。