CN110310705A - 支持simd的序列比对方法及装置 - Google Patents

Abstract

本发明提供一种支持SIMD的序列比对方法，包括根据查询序列和参考序列的长度确定得分矩阵的行数和列数；根据指定的宽度值确定得分矩阵在行方向和列方向上的扩展区间；基于查询序列和参考序列的比对，确定得分矩阵在所述扩展区间内相应各元素的得分，同时将处于扩展区间之外的元素全部清零；然后根据所确定的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。该方法能充分利用SIMD来提高序列比对的计算效率。

Description

支持SIMD的序列比对方法及装置

技术领域

本发明涉及基因数据处理，尤其涉及对DNA序列比对进行加速的方法和系统。

背景技术

随着下一代测序技术(NGS)技术的快速发展，基因测序通量不断提升。在DNA测序中，DNA分子首先会被随机打断成很多个片段，然后将这些无序片段克隆，再经由测序仪进行测序，产生数百万很短的DNA reads(下文中也可称为测序片段)。之后将这些产生的DNAreads与参考基因组进行比对，找到这些DNA reads在参考基因组上的位置。BWA-MEM已成为当下最流行的也是最常用的基因序列匹配工具。提高和完善BWA-MEM工具的效率对于基因测序有着重要的意义。

发明内容

BWA-MEM采用经典的seed-and-extend范式，其主要耗时的部分包括两个：一是根据BWT(Burrow-Wheeler Transform，伯罗斯-惠勒变换)算法计算每个seed的位置，另一个是使用序列比对算法(Smith-Waterman,SW)，基于动态规划对seed进行扩展以找到两序列之间的最优本地比对。其中BWT计算复杂度为O(n)，可优化空间不大；而SW算法计算复杂度为O(n^2)，既复杂又耗时。

因此，本发明的目的在于提供一种可利用SIMD(Single Instruction MultipleData，单指令多数据)技术进行加速的SW算法的方法和系统。

本发明的目的是通过以下技术方案实现的：

一方面，本发明提供了一种支持SIMD的序列比对方法，包括：

s1)根据查询序列和参考序列的长度确定得分矩阵的行数和列数；

s2)根据指定的宽度值确定得分矩阵在行方向和列方向上的扩展区间；

s3)在查询序列和参考序列的比对过程中，确定得分矩阵在所述扩展区间内相应各元素的得分，同时将处于扩展区间之外的元素全部清零；

s4)根据所获得的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。

在上述方法中，还可包括根据下述公式判断是否执行步骤s1)-s4)：

其中a的值等于或略小于要用于执行所述方法的处理器中向量寄存器包含的元素个数；m为查询序列的长度，w表示所述指定的宽度值；以及当上述公式成立时确定执行步骤s1)-s4)。

在上述方法中，所述扩展区间可以为[max(i-w,0)，min(i+w,n)]，其中n为参考序列的长度；w表示所述指定的宽度值，i为[-1,m]之间的整数，m为查询序列的长度。

又一方面，本发明提供了一种计算机可读介质，其中存储有计算机程序，所述计算机程序在被执行时能够用于实现上文所述支持SIMD的序列比对方法。

又一方面，本发明提供了一种序列比对装置，包括存储器和支持SIMD指令集的处理器，所述存储器中存储有计算机程序，所述计算机程序在被所述处理器执行时能够用于实现上文所述支持SIMD的序列比对方法。

与现有技术相比，本发明的优点在于：

通过固定读取宽度和余位补零等方式解决了现有BWA-MEM工具中SW算法中难以在支持SIMD指令集的处理器上有效执行的问题，从而能充分利用SIMD加速序列比对的计算速度，从而改善了BWA-MEM基因匹配工具的效率。

附图说明

以下参照附图对本发明实施例作进一步说明，其中：

图1为得分矩阵的数据依赖关系示意图；

图2为解耦得分矩阵数据依赖关系的变换示意图；

图3为根据本发明一个实施例的支持SIMD的序列比对的方法的流程示意图；

图4为访存数据对比示意图。

具体实施方式

为了使本发明的目的，技术方案及优点更加清楚明白，以下结合附图通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

单指令流多数据(SIMD)是一种实现数据级并行的技术，其可以在1条单独的指令中同时执行多个运算操作，以增加处理器的吞吐量。支持SIMD的CPU有多个执行部件，但都在同一个指令部件的控制之下，中央控制器向各个执行部件发送指令，在指令译码后几个执行部件同时访问内存，一次性获得所有操作数进行运算。SIMD的典型代表是向量处理器(Vector Processor)和阵列处理器(Array Processor)。目前主流处理器支持的SIMD指令集包括MMX、SSE、SSE2等。MMX提供了8个64bit的寄存器进行SIMD操作，SSE提供了128bit的8个寄存器进行SIMD指令操作。SSE2使用了144个新增指令，扩展了MMX技术和SSE技术，支持多种数据类型的算术运算，支持灵活、动态范围更广的计算功能。

在BWA-MEM工具中使用的序列比对SW算法主要分两步：计算得分矩阵和寻找最佳相似片段对。首先基于两个序列获取得分矩阵，接着用动态规划回溯的方法找到局部最大相似片段对，也就是先找到得分矩阵中最大的元素，然后按照元素原路径一步一步往前回溯，直到回溯到0时停止。例如，假设两个序列，一个是查询序列，一个是参考序列。其中定义查询序列为Q＝q₁q₂..q_m、参考序列为D＝d₁d₂..d_n，得分矩阵为W(q_i,d_j)，且当q_i≠d_j时W(q_i,d_j)≤0，当q_i＝d_j时，W(q_i,d_j)＞0。对于比对过程中新出现的空位(gap)和延续的gap分别定义为G_init和G_ext，表示新出现一个gap和在已有gap基础上再添加一个空位的不同代价。SW算法过程主要依赖三个矩阵，矩阵E表示reads末端以gap结尾的得分矩阵，矩阵F表示参考序列以gap结尾的得分矩阵，H为整体比对的得分矩阵。这三个矩阵的递推公式如下:

比对得分矩阵，其中1≤i≤m和1≤j≤n

基于上述公式，矩阵E、F、H的数据依赖关系可以表示为图1中所示关系。从图1中，可以看出矩阵E的每个元素的计算依赖于前一列的元素；矩阵F的每个元素依赖于前一行的元素；而H矩阵元素的计算依赖于前一行、前一列和左上角元素，其中H_i,j依赖于F_i,j，而F_i,j依赖于H_i-1,j,进而H_i,j依赖于H_i-1,j，同理H_i,j依赖于H_i,j-1。这种数据的依赖关系使得无法使用向量化指令来完成上述计算，而SIMD要求采用向量化的指令。因此需要解除这种数据依赖关系才能利用SIMD加速SW算法。图2给出了常用的解耦得分矩阵数据依赖关系的变换示意图。假设向量化指令宽度为4个元素大小，则逐行以每4个元素为一个块，通过右移一个元素、与下一行元素对应、补零等对数据进行重排，使得下一行的计算依赖于上一行的某个向量块值，以及同行的上一个向量块的值，从而解耦了同一个向量化块中的元素依赖关系。因此现有的BWA-MEM工具已经可以利用上述方法对其SW算法进行SSE2指令加速。

然而，现有BWA-MEM在利用SW算法进行扩展时存在内存访问不连续的问题，这无疑会在一定程度上影响SIMD性能的有效利用。如上文介绍的，现有的BWA-MEM工具采用seed-and-extend方式：在seed阶段，BWA取read的碱基片段在参考序列上进行精确匹配，并选择满足一定匹配次数和长度要求read的碱基片段作为seed(种子)，即seed是read序列的一部分；在extend阶段，BWA利用SW算法将seed在read和reference上向两边延伸比对(且容忍gap)，进而找到整个read在reference上符合条件的全局匹配。在初始扩展时会提供一个w，表示可进行正向或逆向扩展的宽度(bandwidth)，也就是上文提到的1≤j≤n变成max(i-w,0)≤j≤min(i+w,n)，而且在每次计算完一行的信息(也即是计算完i表示的一行)，其会根据得分矩阵中相应E和H元素的值调整j的取值区域。这样，每次读入的数据宽度是变化的，从而导致内存空间不连续。以BWA-MEM工具中SW算法的扩展函数ksw_exend为例，其伪码如下：

从上述算法可以看出，由于beg和end随着计算区间大小不断变化，导致处理器在读取数据时内存访问不连续(例如图3左边图中阴影部分所示)，从而无法很好地利用向量化指令提高计算速度。

图4给出了根据本发明一个实施例的支持SIMD的序列比对方法，其主要包括步骤1)根据输入的查询序列和参考序列的长度确定得分矩阵的行数和列数；2)根据输入的宽度值确定得分矩阵在行方向和列方向上的扩展区间；步骤3)在查询序列和参考序列的比对过程中，仅确定得分矩阵在所述扩展区间内相应各元素的得分，同时将得分矩阵中处于扩展区间之外的元素全部清零；步骤4)根据所获得的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。可以看出，在该实施例中在每次迭代的过程中，将每次计算时得分矩阵的行列设置为固定值，也就是固定了i和j的区间；同时为兼容现有基于扩展宽度指定的方法，基于指定的扩展宽度确定在行方向和列方向上的扩展区间，而将该扩展区间之外的计算得分值全部设为零以对齐边界。这样既不影响最终比对结果，又使得内存访问得以连续(例如图3右边图中阴影部分所示)，从而更有利于SIMD向量化指令的执行效率。仍以扩展函数ksw_exend为例，在本发明的实施例中其伪码可如下表示：

其中对于各得分矩阵而言，其元素在行和列方向上的迭代区间是固定的，i的区间为[-1,m]，j的区间为[-1,n]，其中m为查询序列的长度，n为参考序列的长度；而每次迭代时，扩展区间是变化的，所收到的扩展宽度以w表示，从而可确定得分矩阵在行方向和列方向的扩展区间为[beg,end]，其中beg＝max(i-w,0)；end＝min(i+w,n)；同时在[beg,end]外的元素的得分值每次迭代都清0。

另外，采用根据本发明实施例的序列比对方法同时也引入冗余计算，有可能减少一定的向量化的受益。因此在优选的实施例中，还包括判断是否采用上述方法进行序列比对的步骤。假设向量化的加速比为a，指定的扩展宽度为w，查询序列的长度为m，则当下述公式(4)成立时，可采用根据本发明实施例的序列比对方法，否则可采用原有的BWA工具中的序列比对方法：

其中a的值等于或略小于要采用的处理器中向量寄存器包含的元素个数，例如对于Int32的系统，128位的寄存器包括4个元素。当上述公式成立时，采用根据本发明实施例的序列比对方法的受益大于冗余计算增加的计算开销。

虽然本发明已经通过优选实施例进行了描述，然而本发明并非局限于这里所描述的实施例，在不脱离本发明范围的情况下还包括所做出的各种改变以及变化。

Claims (5)

1.一种支持SIMD的序列比对方法包括：

s1)根据查询序列和参考序列的长度确定得分矩阵的行数和列数；

s2)根据指定的宽度值确定得分矩阵在行方向和列方向上的扩展区间；

s3)在查询序列和参考序列的比对过程中，确定得分矩阵在所述扩展区间内相应各元素的得分，同时将处于扩展区间之外的元素全部清零；

s4)根据所获得的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。

2.根据权利要求1所述的方法，还包括根据下述公式判断是否执行步骤s1)-s4)：

其中a的值等于或略小于要用于执行所述方法的处理器中向量寄存器包含的元素个数；m为查询序列的长度，w表示所述指定的宽度值；以及当上述公式成立时确定执行步骤s1)-s4)。

3.根据权利要求1所述的方法，其中所述扩展区间为[max(i-w,0)，min(i+w,n)]，其中n为参考序列的长度；w表示所述指定的宽度值，i为[-1,m]之间的整数，m为查询序列的长度。

4.一种计算机可读介质，其中存储有计算机程序，所述计算机程序在被执行时能够用于实现如权利要求1-3中任一项所述的方法。

5.一种序列比对装置，包括存储器和支持SIMD指令集的处理器，所述存储器中存储有计算机程序，所述计算机程序在被所述处理器执行时能够用于实现如权利要求1-3中任一项所述的方法。