CN116168765A

CN116168765A - 基于改进strobemer的基因序列生成方法及系统

Info

Publication number: CN116168765A
Application number: CN202310449411.8A
Authority: CN
Inventors: 刘卫国; 李金津; 殷泽坤
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2023-04-25
Filing date: 2023-04-25
Publication date: 2023-05-26
Anticipated expiration: 2043-04-25
Also published as: CN116168765B

Abstract

本发明公开了一种基于改进strobemer的基因序列生成方法及系统，涉及生物基因序列数据处理技术领域，该方法包括：针对每条碱基序列，从碱基序列的第一个碱基字符的位置开始，生成strobe序列，并进行合法性判断，若判断为合法，则基于hash函数计算合法strobe序列的hash值并保存；依次对碱基序列中的前n‑k+1个碱基字符循环上述生成序列的过程；以每个合法strobe序列为起始序列，调取所保存的hash值，利用改进的chop_X函数选取strobe序列，构建strobemer。本发明采用多线程并行处理，通过先进行合法性判断及hash值计算，避免冗余的判断和计算，提高处理效率。

Description

基于改进strobemer的基因序列生成方法及系统

技术领域

本发明涉及生物基因序列数据处理技术领域，尤其涉及一种基于改进strobemer的基因序列生成方法及系统。

背景技术

近年来，随着生物测序技术的进步与发展，生物基因数据库的规模越来越大，在计算生物基因组学和基因序列分析中，常采用传统的基于k-mer的方法，即将对基因数据的处理转化为对基因序列中一系列公共子序列片段（k-mer，长度为K的子字符串）的处理，以此实现高效的基因数据处理。然而，基于k-mer的传统方法对于突变十分敏感。现有技术中，虽然已经提出多种改进k-mer算法能够克服这种敏感性，但是仍不能解决基因序列增加和删除的变异情况。目前，strobemer算法作为一种新提出的种子算法，能够很好的克服k-mer算法的不足，具有广阔的研究前景，但是strobemer算法程序的实现却存在计算速度慢、处理效率低等问题。

在生物信息学中，k-mer是生物基因序列中长度为k的子串，如图1所示的一条序列3-mer的示意图。而strobemer可以看做是从生物序列s中按照一定选取规则选出的若干k-mers连接在一起形成基因序列，将strobemer中的每个k-mer称为strobe，即每个strobemer均是由n个strobe组成。采用strobemer算法对基因数据进行处理，生成多个strobemer，其相比于原基因数据要小的多，以此达到数据降维的效果。

现有的strobemer算法通常采用单线程的形式实现，运行速度较慢，处理1GB的生物数据需要上百秒的时间开销，严重影响了测序数据的后续分析，而且构建strobemer的过程中采用的多种构建函数耗时较长，strobemer算法的计算效率低。即，现有的strobemer算法分析处理测序数据的效率低，难以为大规模测序数据提供算力支持。

发明内容

为解决上述现有技术的不足，本发明提供了一种基于改进strobemer的基因序列生成方法及系统，通过改进strobemer构建过程中的数据处理方式，避免冗余的判断和计算，提高单线程处理的效率，并采用多线程并行处理多个碱基序列，提升处理效率，降低耗时。

第一方面，本公开提供了一种基于改进strobemer的基因序列生成方法。

一种基于改进strobemer的基因序列生成方法，包括：

读取多条碱基序列；

分别对每条碱基序列进行strobemer序列的构建，构建过程包括：

从碱基序列中的第一个碱基字符的位置开始，生成一个strobe序列，并对该strobe序列进行合法性判断，若判断为合法，则基于hash函数，计算合法strobe序列的hash值并保存；

循环生成strobe序列的过程，依次对碱基序列中的前n-k+1个碱基字符生成对应的strobe序列，筛选合法strobe序列并计算合法strobe序列的hash值；其中，n为碱基序列的长度，k为strobe序列的长度；

以每个合法strobe序列为起始序列，分别调取所保存的strobe序列的hash值，利用改进的chop_X函数选取strobe序列，构建strobemer序列。

进一步的技术方案，在对每条碱基序列构建strobemer序列的过程中，

保存每个合法strobe序列在碱基序列中的起始位置；

在构建strobemer序列后，调取strobemer序列中每个strobe序列的起始位置；

根据每个strobe序列的起始位置，结合碱基序列，输出所构建的strobemer序列。

进一步的技术方案，构建缓冲区，在缓冲区中保存每个合法strobe序列的hash值，并保存每个合法strobe序列在碱基序列中的起始位置。

进一步的技术方案，采用多线程读取多条碱基序列，每个线程负责处理一条碱基序列，不同线程同时处理不同的碱基序列。

进一步的技术方案，利用生产者线程，读取多条碱基序列，将读取的数据放入缓冲队列中；

利用多个消费者线程，从缓冲队列中分别顺序读取一条碱基序列，对读取的碱基序列构建strobemer序列。

进一步的技术方案，所述改进的chop_X函数是指向量化展开strobemer构建方法中strobe选取过程的内层循环。

进一步的技术方案，所述内层循环向量化展开包括：

将内层循环展开为每次循环同时读入x个strobe的hash值，并初始化为一个向量，确定该向量的最小值；其中，x表示向量的大小；

将获取的向量的最小值与上一次循环得到的最小值进行比较，确定当前循环的窗口最小值。

第二方面，本公开提供了一种基于改进strobemer的基因序列生成系统。

一种基于改进strobemer的基因序列生成系统，包括：

碱基序列读取模块，用于读取多条碱基序列；

strobemer序列构建模块，用于分别对每条碱基序列进行strobemer序列的构建，构建过程包括：

保存每个合法strobe序列在碱基序列中的起始位置；

以上一个或多个技术方案存在以下有益效果：

1、本发明提供了一种基于改进strobemer的基因序列生成方法及系统，通过以碱基序列中每个碱基为起始位置进行划分，得到多个strobe序列，对每个strobe序列先进行合法性判断，并计算合法strobe序列的hash值，在后续针对不同起始序列利用不同函数选取strobe序列时，不必每次都需要进行合法性判断及重新计算hash值，直接调取合法strobe相应的hash值，避免冗余的判断和计算，提高单线程处理的效率；

2、本发明记录每个合法strobe序列在碱基序列中的起始位置，在返回strobemer时仅返回相应strobe序列的起始位置即可，避免多余的拷贝及动态分配需要申请存储空间所造成的耗时，提高处理效率；

3、本发明通过生产者线程负责从fasta文件中读取数据放入缓冲队列，通过消费者线程负责从缓冲队列中读取数据、解析为fasta数据并处理得到strobemers，而且在构建strobemers时，采用多线程并行处理多个碱基序列，提升处理效率，降低耗时。

附图说明

构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1为序列3-mer的示意图；

图2为strobemer窗口的示意图；

图3为基于minstrobes构建strobemer的原理示意图；

图4为基于randstrobes构建strobemer的原理示意图；

图5为基于hybridstrobes构建strobemer的原理示意图；

图6为strobemer构建过程中热点函数的耗时示意图；

图7为现有strobemer构建方法的流程图；

图8为基于minstrobes选取strobe的流程图；

图9为本发明实施例一所述基于改进strobemer的基因序列生成方法的流程图；

图10为本发明实施例一所述方法中单线程处理的优化耗时结果示意图；

图11为本发明实施例一所述方法中多线程处理的加速效果示意图；

图12为本发明实施例一所述方法中向量化minstrobemers构建方法的示意图。

具体实施方式

应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

实施例一

首先，明确strobemer技术。strobemer是指从生物基因序列s中按照一定选取规则选出的若干k-mers连接在一起而形成的基因序列。为了降低基因数据的维度，实现高效的基因数据处理，利用strobemer算法构成strobemer基因序列。每个strobemer由n个strobe组成，对于一个strobemer而言，其包含的第一个strobe是从序列s内的固定位置生成的，其余每个strobe都是从固定位置偏移得到的新窗口中选取。如图2所示，若第1个strobe从位置i开始，那么第2个strobe将从下游窗口中

选取（图示中的wmin即为/>

，wmax即为/>

）；类似的，第3个strobe将从窗口/>

中选取。以此类推，第n个strobe将从窗口/>

选取。因此，将一个strobemer涉及到的参数归纳为/>

。其中，n表示strobes的个数，l表示strobe的长度，即strobe是一个长度为l的l-mer，/>

和/>

分别表示窗口长度的上限和下限。

构建strobemer的重点在于如何将strobe从窗口中选出。根据选取思想的不同，strobemer有三种构建方法，分别是：minstrobes、randstrobes和hybridstrobes。

其中，minstrobes构建方法为：选取窗口内所有的strobe，经过hash函数h计算得到的hash值，hash值最小的strobe即为该窗口中选出的strobe。如图3所示，第一个strobe从位置0开始，第二个为strobe从第一个窗口中选出，分别计算窗口内为TTC、TCT和CTG的strobe序列的hash值（即哈希值），分别为5、2、19，选取hash值最小的strobe序列TTC作为该窗口的strobe。同样的，选出的第三个strobe与此类似。

randstrobes构建方法与minstrobes每个strobe的选取都是独立的不同，其考虑了strobe间的条件依赖性。如图4所示，第一个strobe从位置0开始，第二个为strobe从第一个窗口中选出，分别计算窗口内为TTC、TCT、CTG的strobe序列的条件hash值，即h(TTG|AGT)、h(TCT|AGT)和h(CTG|AGT)分别为7、13、5，选取条件hash值最小的strobe序列CTG作为该窗口的strobe。同样的，选出的第三个strobe与此类似。

上述条件hash值是指：假设窗口内一个候选strobe是子串

，在此之前已经选取了/>

作为strobes，则这个候选strobe对应的条件hash值的计算公式为：

其中，⨁表示子串的拼接。

hybridstrobes构建方法可以看做是前两种方法思想的混合，其选取的strobe依然是窗口内hash值最小的子串，同时还考虑了两个strobe之间的条件依赖性。具体的，将采样窗口划分为x个不相交的子窗口，子窗口的长度

根据/>

计算得到，即第二个strobe的采样窗口被划分为/>

，

，其他strobes的采样窗口也是类似的。用上一个选取到的strobe的hash值取余，来选定当前strobe的采样子窗口，并在子窗口内选hash值最小的作为当前窗口选取到的strobe。如图5所示，第一个strobe从位置0开始，第二个窗口为第二个strobe的采样窗口，其被分为三个不相交的子窗口，即strobe序列TTC所在窗口、strobe序列TCT所在窗口和strobe序列CTG所在窗口，上一个取到的strobe序列即第一个strobe序列AGT对应的hash值为18，18对子窗口个数3取余为0，所以选定strobe序列TTC所在子窗口为采样子窗口，选取该子窗口内hash值最小的strobe为选取结果，该窗口只有一个strobe所以选取TTC为第二个strobe。同样的，选出的第三个strobe与此类似。

strobemer构建过程中，采用上述任一种构建方法将strobe从窗口中选出，再将选出的若干strobe连接在一起，形成的基因序列。针对一条碱基序列完成完整的strobemer构建，现有的strobemer构建方法如图7所示，包括：读入序列，计算strobes的hash值，构建strobemer和strobemer结果保存。具体的，首先，打开一个fasta格式（fasta格式是一种基于文本用于表示核酸序列或多肽序列的格式）的生物数据文件，顺序读出每一条碱基序列；然后，根据碱基序列，利用hash函数，计算每个strobe的hash值，其中，hash函数可以是任何以字符串作为输入而输出一个整数的函数的抽象，其满足均一性和确定性；之后，根据三种构建方法选取strobes组成strobemer；最后，保存得到的strobemer序列。对于一条长度为n、构建的strobemer中strobe长度为k的碱基序列，重复上述过程n-k+1次，获得n-k+1个strobemer基因序列，保证每一个碱基为起始位置的strobe都能作为第一个strobe计算一次strobemer结果。

现有的strobemer算法采用单线程实现，运行速度慢，而且如图6所示，strobemer构建过程中，三种构建方法的耗时都主要出现在chop_X函数（X为minstrobes、randstrobes或hybridstrobes，表示不同的构建方法）和strobemer构建函数的语句中。上述chop_X函数是chop_minstrobes、chop_randstrobes和chop_hybridstrobes三个函数的代称，可以表述为strobemer的切分函数，在程序运行中直接调用该chop_X函数。此外，图6中，operatornew函数是将返回一个指向不同的对象的指针（即对operator new的重复调用将返回不同的指针）；strncpy函数用于将指定长度的字符串复制到字符数组中；free函数则是C语言中释放内存空间的函数；这些函数均是在构建strobemer时所调用的函数。

其中，chop_X函数的耗时主要体现在：chop_X函数以序列作为输入，返回一个strobemer对象数组。具体的，首先chop_X函数计算每个k-mer和它们对应的hash值，然后以序列中所有有潜力能生成strobemer的位置为起始位置，窗口滑动n-1次，每次都按照规则选取一个strobe，并把这个strobe的序列用strncpy函数复制到strobemer对象中，形成一个三层嵌套循环。以minstrobes为例，chop_minstrobes函数的运算过程如图8所示，现有的strobemer构建方法包括外层循环、中层循环和内层循环，以此形成三层嵌套循环，该程序中循环的耗时与循环层数和循环次数都有关系，但循环次数取决于需要处理的序列长度，是不能改变的，因此，需要减少循环的嵌套次数来降低函数的耗时。

如图8所示，上述采用minstrobes构建方法构建strobemer的方法包括：（1）外层循环：选取strobemer起始位置为第一个strobe，判断该条碱基序列是否还有strobe以及其他strobe是否合法，若否，则构建完成所有strobemer并退出外层循环，反之则进行中层循环；（2）中层循环：取第一个strobe为初始序列，并判断strobe数量是否满足strobemer的要求，若是，则构建完一个strobemer并退出中层循环回到外层循环，反之，则初始化窗口的最小hash值，取最小hash值的strobe的位置并选取strobe的窗口上下限，进入内层循环；（3）内层循环：判断窗口内是否还有候选strobe，若是，则判断该候选strobe的hash值是否小于当前窗口的最小hash值，若是，则更新窗口的最小hash值并取最小hash值的strobe位置，重新判断窗口内是否还有候选strobe，若是，重复上述判断条件及操作，反之则退出内层循环回到中层循环。

另一个热点函数是strobemer的构建函数。因为chop_X以strobemer的对象数组作为返回值，在构建每一条序列的strobemer之前，需要动态申请一个对象数组作为返回结果。与静态声明不通，动态分配十分耗时，因为其需要在运行时才向操作系统申请，而且这样的分配方式还会使代码更复杂、更容易出错。因此，需要改变算法中的数据结构来降低耗时。

为此，本实施例提供了一种基于改进strobemer的基因序列生成方法，如图9所示，包括以下步骤：

读取多条碱基序列；

进一步的，在生成碱基序列中的前n-k+1个碱基字符所对应的strobe序列时，将筛选的合法strobe序列的hash值保存到数组（即uint64_t的数组）中，该数组的最大有效长度为n-k+1，便于后续直接调取该数组中的数据。

上述步骤中，对于一条碱基序列而言，它的strobemer集合是以每个碱基起始的strobe为第一个strobe，依次在窗口内选取strobes构建为strobemer的集合。现有的策略在每取到一个strobe后，都需要判断strobe的合法性并计算strobe的hash值，存在许多冗余的判断和计算，因为每个strobe除了可以作为第一个strobe以外，它还可以在窗口内成为候选strobe。本实施例所述改进后的方案中，在读取碱基序列之后，先对碱基序列包含的所有strobe进行一次合法判断和hash值计算，保存在计算结果，为后续每个基于strobe构建strobemer的过程节省计算时间，避免重复及重复计算所造成的时间浪费。

进一步的，在分别对每条碱基序列构建strobemer序列的过程中：

保存每个合法strobe序列在碱基序列中的起始位置；

进一步的，将每个合法strobe序列在碱基序列中的起始位置保存至数组（即unsigned int的数组）中，便于后续直接调取该数组中的数据。

现有的策略中，chop_x热点函数以一个三层嵌套循环为核心内容，循环内嵌套一个主要内容为分支判断的内层循环。分支判断的代价是昂贵的，因为无论是编译器的自动优化还是手动优化，分支判断都会使得编译过程不得不去考虑多种情况甚至是跳转，在考虑完这些情况后再优化，优化效果会大大降低；再加上分支预测并不是百分百都能预测中，如果预测失败则会更加耗时，因此，应该尽可能的减少无用的分支判断。为此，本实施例中，改进现有方案中hash值计算及strobe合法性判断的方式，避免每次都需要判断strobe是否有效并计算strobe的hash值的情况，直接使用一次循环遍历序列，采用一个uint64_t的数组来存储每个strobe的hash值，同时，采用一个unsigned int的数组来存储每个strobe在序列中的起始位置，对于不能生成strobe的位置不再保存。通过上述方案，虽然不能降低计算hash值这一步骤的复杂度，但是，却能在chop_x中消去一个复杂度为O(n)的循环，与改进前处理3GB数据需要550s以上的时间相比，改进后的单线程处理性能提升在50s以上，对于randstrobes这种构建方法，优化效果甚至超过80s。

而且，现有的strobemer构建方法使用C语言的strncpy函数复制strobes到结果对象中，复制函数总是耗时的，因为除了要申请存储空间，它还需要拷贝一份一样的数据到新的位置。因此，应该尽可能的减少不必要的数据拷贝。

考虑到一个strobemer是由若干个strobes组成的，而且这些strobes都是基因序列的一个子串，因此，要想知道一个strobemer到底包含哪些strobes，只需要记录下每个strobe在序列中的起始位置即可。故，本实施例中，在返回结果中只记录每个strobemer的起始位置的unsigned int数组，避免使用strncpy函数进行多余的拷贝，而且能够避免重复调用作为热点函数之一的strobemer构建函数来动态分配strobemer对象数组。通过本实施例上述优化返回结果的数据结构的方案，单线程处理性能提升在140s以上。

尽管优化返回结果的数据结构避免了strobemer对象数组的动态分配，但这使得strobes的hash值数组和index数组的动态分配。既然动态分配总是不可避免的，那么至少可以尽可能的减少动态分配的次数。因此，本实施例中，设置一个缓冲区，把所有需要在构建过程中反复使用和动态分配的数据结构都放在缓冲区中。即，构建缓冲区，在缓冲区中保存每个合法strobe序列的hash值，并保存每个合法strobe序列在碱基序列中的起始位置。

上述结构在缓冲区初始化时就分配好，每次使用前对其大小进行检查，当大小不够时重新分配更大的缓冲区。通常情况下，扩充数组的代价是巨大的，因为除了扩充数组以外，还需要把扩充前已经保存在数组中的内容复制过去。但在strobemer构建中，数组的扩充只会发生在计算某条碱基序列的strobemers之前，换言之，此时strobes的hash值数组和index数组保存的是上一条序列的相关值。因此，数组的扩充只需要动态分配一个新的数组，释放掉原有的数组就可以，不需要进行复制。比起计算每条序列的strobemer前都要动态分配一个数组，计算完后都要释放掉数组，扩充数组的代价显然更小。在使用缓冲池之后，单线程处理性能提升在150s以上，hybridstrobes这种构建方法性能提升甚至超过180s。

进一步的，采用多线程读取/处理多条碱基序列，每个线程负责处理一条碱基序列，不同线程同时处理不同的碱基序列。现有的策略采用单线程的方式处理碱基序列，只有一条碱基序列被处理完后才会开始下一条碱基序列的处理，这样的串行处理方式效率低下。考虑到碱基序列与碱基序列之间strobemer集合的构建起始是互不相关的，因此，本实施例中，采用多线程方式处理碱基序列，每个线程负责一条序列的处理，不同线程可以同时处理不同的序列，当一个线程处理完当前序列后即可获取一条新的序列再次处理，线程和线程间是互不干扰的，以此提高运行的效率。

进一步的，利用生产者线程，读取多条碱基序列，将读取的数据放入缓冲队列中；

现有的strobemer构建方法采用单线程串行实现，然而随着多核处理器的出现和发展传统的串行算法难以从中获益，不能最大程度地发挥硬件资源的能力，达到更优的计算速度，因此，在保证程序准确性的前提下，本实施例利用生产者—消费者模型对strobemer算法的设计和实现进行了优化改进，实现了程序的多线程并行。其中，生产者线程负责从fasta文件（fasta格式是一种基于文本用于表示核酸序列或多肽序列的格式）中读取数据放入缓冲队列，消费者线程负责从缓冲队列中拿取数据、解析成fasta数据并处理成strobemers。线程以块为单位处理数据，起始状态下可自定义缓冲池的大小，缓冲池用来限制生产者和消费者一次最多能够处理的数据大小，这使得数据不需要一次性全部读入占用大量的内存而是流式读入只占用规定的内存大小。同时也避免了反复对内存的申请和释放，节省了时间。相比于单线程串行处理3GB数据需要500s以上的耗时，多线程并行性能得到了大幅提高，程序耗时随线程数的增多而线性降低，在40个线程的情况下，处理3GB数据耗时低于5s。

此外，考虑到单指令流多数据流（Single Instruction Multiple Data，SIMD）作为Intel和AMD x86指令集的扩展，能够使多个功能单元同时作用于批量数据元素。由于单个指令可以同时处理多个数据元素，这种指令级并行形式提供了一种数据并行的新方法。除了支持64位、128位和256位寄存器外，还支持使用浮点和整数变量对短向量进行操作。strobemer算法的三种构建方法都需要从滑动窗口中选取hash值（或条件hash值）最小的strobe，这种算法很适合向量化。

现有程序中minstrobes构建方法的三层循环都不存在循环依赖，randstrobes和hybridstrobes构建方法虽然都要依赖前一个窗口选取的strobe的hash值来从当前窗口选取strobe，但最内层的循环，即从当前窗口中选取hash值（或条件hash值）最小的strobe却是能够展开的，因为窗口内strobes的hash值计算是相互不影响的，因此，可以使用手动向量化展开3种构建方法的最内层循环。通过调用intel的instrinsics函数能够简洁的实现手动向量化。

即，上述改进的chop_X函数是指在现有的strobemer构建方法的基础上，向量化展开strobemer构建方法中strobe选取过程的内层循环。

进一步的，上述内层循环向量化展开包括：

将获取的向量的最小值与上一循环得到的最小值进行比较，确定当前循环的窗口最小值。

优选的，若strobe的hash值的数量小于x时，则利用保存hash值的数组中的最大值填充strobe的剩余位置的hash值。

下述通过三种不同的构建方法进行详细的解释说明。minstrobes构建方法的最内层循环实现的是遍历窗口内所有的strobes的hash值，通过比较寻找最小hash值对应的strobe，instrinsics提供了规约向量元素中最小值的操作，因此，可以把mistrobes的最内层循环展开为每一次循环同时读入x个strobe的hash值（x为向量的大小），并将其初始化为一个向量，然后使用instrinsics函数规约该向量的最小值，最后，比较此次规约得到的最小值和此次循环之前得到的最小值，决定此次循环为止窗口内的最小值。当窗口内剩余的strobes数目小于x时，不再展开循环而是继续使用初始minstrobes的思路寻找最小值。通过上述方式，把每次循环读取一个strobe的hash值和当前最小hash值比较转换为每次循环找到x个strobe的最小hash值和当前最小hash值比较。例如图12所示，当窗口大小为9、向量大小为4时，初始程序的minstrobes需要9次循环，而向量化之后的minstrobes仅需要3次循环即可。在实际情况中，使用的是AVX 512指令集，数据宽度为512位，每个hash值的数据宽度是64，因此使用向量化能够同时处理8个数据，再加上最外层循环的大小约等于一条序列的长度，使用向量化在理论上能够很好的降低运算的次数。

randstrobes构建方法的最内层实现思路是利用上一个窗口的条件hash值与当前窗口的所有strobes的hash值做异或运算求条件hash值，然后找出条件hash值最小的strobe。instrinsics提供了向量的异或操作，类似于minstrobes的向量过程，randstrobes的向量化过程也一次读入8个strobes的hash值初始化成一个向量，同时把上一个窗口的最小条件hash值初始化成一个8个维度取值都相同的向量，调用instinsics函数的异或操作把该向量和strobe的hash值向量进行异或运算得到8个strobes的条件hash值向量，用规约最小值函数找出8个strobes的最小条件hash值，把规约结果和上一次循环得到的当前窗口最小条件hash值做比较，如果规约结果比上一次循环得到的当前窗口最小条件hash值小，那么规约结果将成为新的最小条件hash值，否则不做任何操作。在randstrobes中，相对于minstrobes，还多一步异或操作，除了能够达到上述minstrobes所实现的技术效果外，上述方式使得原本需要进行8次的异或操作，仅通过1次向量异或即可完成，进一步降低了计算成本。

hybridstrobes构建方法其实可以看做是更小规模的minstrobes，因为hybridstrobes把窗口再分成了x个子窗口，strobes从x个子窗口中的某一个窗口中用minstrobes的思想选取。对上一个窗口中取到的strobe对应的hash值取余，余数决定子窗口的选择。因为hybridstrobes分割窗口的思想使得其最内层循环不再是以为上下限，而是更小。例如，一个长度为9的窗口，在x为3的前提下，hybridstrobes方法只需要在长度为3的窗口内取最小值。尽管循环的次数变了但是循环内部找最小值的思想却没有变化，所以hybridstrobes向量化的思想和minstrobes是一样的，只不过hybridstrobes向量化的过程一次性读入的是子窗口内的8个strobes的hash值。与上述提到的两种构建方式相同，hybridstrobes会面临子窗口中strobes的数目小于8的情况。一种最坏的情况是窗口（或子窗口）内strobes的数目小于8，这会使向量化展开退化到向量化之前的状态，hybridstrobes比其他两种构建方式更容易遇到这种最坏的情况，这是由于当窗口strobes的数目为[1,8)时，前两种构建方法才会遇到上述最坏的情况，但因为hybridstrobes缩小了实际计算的窗口大小，所以其在窗口strobes的数目为[1,8x)时（x为向量的大小），都会遇到这种情况。

为了解决这种情况，使向量化能在最大的窗口范围都发挥应有的作用，本实施例采取了一种新的循环展开策略。前两种构建方法的循环展开策略是，只展开8的整数倍的strobes，对于窗口初始状态内strobes的数目就小于8和窗口已经按8展开过多次，剩下的strobes数目小于8的情况，采用原始循环（即未向量化前的循环）继续找最小值。这样的原始循环次数一定小于8，对于一条长度上万的序列来说，这样的计算是可以忍受的。对于hybridstrobes来说，面对这种情况本实施例采用一种新的策略，即当strobes数目小于8时，把strobes的hash值全部读入向量，不足8的剩余位置全部用UINT64_MAX填充，因为求取的是最小值，这样的操作并不会对结果造成影响，但是却能让不足8的strobes也能使用向量化。显然，这种策略同样也适合minstrobes和randsrobes，因此，最终选取这一策略来进行循环展开。

手动向量化前实现了自定义的数据对齐，帮助处理器更快地取到数据，为向量化做准备。向量化过程中使用了lowbit函数解决intel提供的向量化接口只能最小化规约得到二进制序列的问题。在向量化strobe的选取过程实现了数据并行之后，程序单线程条件下性能提升在50s以上，randstrobes构建方法提升超过120s。

为了进一步证明本实施例所述方案的优越性，对现有方法和本实施例上述方法进行了测试比较。首先，针对单线程处理优化的效果，利用3GB、格式为fasta的测试数据进行测试，测试结果如图10所示。初始条件下，执行现有的程序三种不同的构建方式耗时都在百秒级，hybridstrobes作为耗时最短的构建程序也需要450s左右的耗时。执行本实施例所述的优化后的方案，单线程程序的性能得到有效提升，耗时最长的minstrobes构建函数构建时间保持在100s左右，性能提升4倍，其余两种构建函数都能在50s以内完成strobemer的构建，性能提升在8倍以上。

其次，针对多线程处理优化的效果，作为衡量并行计算效果的核心指标，加速比定律准确地描述了算法程序并行化后的性能增益。对优化后的单线程处理进行并行化后，性能得到了进一步提升，如图11所示，随着线程数的增加，程序的加速比也线性增加，其中加速效果最好的minstrobes构建方法加速比在20以上，其余两种构建方法加速比也在10左右。在多线程的条件下，处理1GB的生物数据能保证在个位数的时间内完成。

实施例二

本实施例提供了一种基于改进strobemer的基因序列生成系统，包括：

碱基序列读取模块，用于读取多条碱基序列；

以上实施例二中涉及的各步骤与方法实施例一相对应，具体实施方式可参见实施例一的相关说明部分。

本领域技术人员应该明白，上述本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims

1.一种基于改进strobemer的基因序列生成方法，其特征是，包括：

读取多条碱基序列；

2.如权利要求1所述的基于改进strobemer的基因序列生成方法，其特征是，在对每条碱基序列构建strobemer序列的过程中，

保存每个合法strobe序列在碱基序列中的起始位置；

3.如权利要求2所述的基于改进strobemer的基因序列生成方法，其特征是，构建缓冲区，在缓冲区中保存每个合法strobe序列的hash值，并保存每个合法strobe序列在碱基序列中的起始位置。

4.如权利要求1所述的基于改进strobemer的基因序列生成方法，其特征是，采用多线程读取多条碱基序列，每个线程负责处理一条碱基序列，不同线程同时处理不同的碱基序列。

5.如权利要求1所述的基于改进strobemer的基因序列生成方法，其特征是，利用生产者线程，读取多条碱基序列，将读取的数据放入缓冲队列中；

6.如权利要求1所述的基于改进strobemer的基因序列生成方法，其特征是，所述改进的chop_X函数是指向量化展开strobemer构建方法中strobe选取过程的内层循环。

7.如权利要求6所述的基于改进strobemer的基因序列生成方法，其特征是，所述内层循环向量化展开包括：

8.一种基于改进strobemer的基因序列生成系统，其特征是，包括：

碱基序列读取模块，用于读取多条碱基序列；

9.如权利要求8所述的基于改进strobemer的基因序列生成系统，其特征是，在对每条碱基序列构建strobemer序列的过程中，

保存每个合法strobe序列在碱基序列中的起始位置；

10.如权利要求9所述的基于改进strobemer的基因序列生成系统，其特征是，构建缓冲区，在缓冲区中保存每个合法strobe序列的hash值，并保存每个合法strobe序列在碱基序列中的起始位置。