支持SIMD的序列比对方法及装置制造方法及图纸

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

Sequence alignment method and device supporting SIMD

【技术实现步骤摘要】
支持SIMD的序列比对方法及装置
本专利技术涉及基因数据处理，尤其涉及对DNA序列比对进行加速的方法和系统。
技术介绍
随着下一代测序技术(NGS)技术的快速发展，基因测序通量不断提升。在DNA测序中，DNA分子首先会被随机打断成很多个片段，然后将这些无序片段克隆，再经由测序仪进行测序，产生数百万很短的DNAreads(下文中也可称为测序片段)。之后将这些产生的DNAreads与参考基因组进行比对，找到这些DNAreads在参考基因组上的位置。BWA-MEM已成为当下最流行的也是最常用的基因序列匹配工具。提高和完善BWA-MEM工具的效率对于基因测序有着重要的意义。
技术实现思路
BWA-MEM采用经典的seed-and-extend范式，其主要耗时的部分包括两个：一是根据BWT(Burrow-WheelerTransform，伯罗斯-惠勒变换)算法计算每个seed的位置，另一个是使用序列比对算法(Smith-Waterman,SW)，基于动态规划对seed进行扩展以找到两序列之间的最优本地比对。其中BWT计算复杂度为O(n)，可优化空间不大；而SW算法计算复杂度为O(n^2)，既复杂又耗时。因此，本专利技术的目的在于提供一种可利用SIMD(SingleInstructionMultipleData，单指令多数据)技术进行加速的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的典型代表是向量处理器(VectorProcessor)和阵列处理器(ArrayProcessor)。目前主流处理器支持的SIMD指令集包括MMX、SSE、SSE2等。MMX提供了8个64bit的寄存器进行SIMD操作，SSE提供了128bit的8个寄存器进行SIMD指令操作。SSE2使用了144个新增指令，扩展了MMX技术和SSE技术，支持多种数据类型的算术运算，支持灵活、动态范围更广的计算功能。在BWA-MEM工具中使用的序列比对SW算法主要分两步：计算得分矩阵和寻找最佳相似片段对。首先基于两个序列获取得分矩阵，接着用动态规划回溯的方法找到局部最大相似片段对，也就是先找到得分矩阵中最大的元素，然后按照元素原路径一步一步往前回溯，直到回溯到0时停止。例如，假设两个序列，一个是查询序列，一个是参考序列。其中定义查询序列为Q＝q1q2..qm、参考序列为D＝d1d2..dn，得分矩阵为W(qi,dj)，且当qi≠dj时W(qi,dj)≤0，当qi＝dj时，W(qi,dj)＞0。对于比对过程中新出现的空位(gap)和延续的gap分别定义为Ginit和Gext，表示新出现一个gap和在已有gap基础上再添加一个空位的不同代价。SW算法过程主要依赖三个矩阵，矩阵E表示reads末端以gap结尾的得分矩阵，矩阵F表示参考序列以gap结尾的得分矩阵，H为整体比对的得分矩阵。这三个矩阵的递推公式如下:比对得分矩阵，其中1≤i≤m和1≤j≤n基于上述公式，矩阵E、F、H的数据依赖关系可以表示为图1中所示关系。从图1中，可以看出矩阵E的每个元素的计算依赖于前一列的元素；矩阵F的每个元素依赖于前一行的元素；而H矩阵元素的计算依赖于前一行、前一列和左上角元素，其中Hi,j依赖于Fi,j，而Fi,j依赖于Hi-1,j,进而Hi,j依赖于Hi-1,j，同理Hi,j依赖于Hi,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.一种支持SIMD的序列比对方法包括：s1)根据查询序列和参考序列的长度确定得分矩阵的行数和列数；s2)根据指定的宽度值确定得分矩阵在行方向和列方向上的扩展区间；s3)在查询序列和参考序列的比对过程中，确定得分矩阵在所述扩展区间内相应各元素的得分，同时将处于扩展区间之外的元素全部清零；s4)根据所获得的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。

【技术特征摘要】
1.一种支持SIMD的序列比对方法包括：s1)根据查询序列和参考序列的长度确定得分矩阵的行数和列数；s2)根据指定的宽度值确定得分矩阵在行方向和列方向上的扩展区间；s3)在查询序列和参考序列的比对过程中，确定得分矩阵在所述扩展区间内相应各元素的得分，同时将处于扩展区间之外的元素全部清零；s4)根据所获得的得分矩阵获取查询序列相对于参考序列的局部最大相似片段。2.根据权利要求1所述的方法，还包括根据下述公式判断是否执行步骤s1)-s4)：其中a的值等于或略小于要用于执行所述方法的处理器中向量寄存器包含的元素个数；m为查询序列的长度，w表示所述...

【专利技术属性】
技术研发人员：谭光明，张中海，牛钢，曾平，
申请(专利权)人：北京哲源科技有限责任公司，
类型：发明
国别省市：北京,11