CN106778077A - 一种将Smith‑Waterman算法在FPGA平台实现的新方法 - Google Patents

Abstract

该发明公开了一种将Smith‑Waterman(以下简称S‑W)算法在FPGA平台实现的新方法。本发明是围绕S‑W算法的打分与回溯两方面进行优化。S‑W算法打分方面，本发明提出优化的删减策略和优化的动态规划(Dynamic Programming)分块方案来减少FPGA实现的时间消耗和资源消耗；S‑W算法回溯方面，本发明采用简化的回溯方向，更加便于FPGA实现，本发明采用两组回溯模块，与打分系统构成流水线处理，大大提高打分系统硬件利用效率，减少序列比对的整体时间。

Description

一种将Smith-Waterman算法在FPGA平台实现的新方法

技术领域

本发明涉及生物信息计算领域，特别涉及一种将Smith-Waterman算法在FPGA平台实现的新方法。

背景技术

生物信息学是生物学、医学、数学、计算机科学等多个学科的前沿交叉学科，具有重要的理论和应用价值。生物信息序列通常指的是通过基因测序技术得到的DNA、RNA序列，这些序列包含了该物种的遗传信息，通过对生物信息序列的分析处理得到的信息在医学疾病研究、生物进化研究等领域都具有重要的意义。

生物信息序列的处理流程中最为重要的是序列比对，同时也是耗时较久的一个流程，这是因为通常待比对物种基因碱基数量庞大，比如人类全基因组的总碱基个数大约是30亿，采用存储介质存储后，目标序列和参考序列所占用空间将达到数百GB，另一方面对这数量级别的数据进行计算处理将会耗费大量的计算资源和耗费大量的时间。

生物信息序列比对最基本的是双序列比对，两条由A，C，G，T四种基本碱基构成的DNA序列通过比对算法得到待比对序列相对于参考序列的比对结果。双序列比对算法一般有3种：动态规划算法(Dynamic Programming)、和词法和点矩阵法，本发明使用动态规划算法。

用于生物序列比对的动态规划算法主要为两种：适用于局部比对的Smith-Waterman算法和适用于全局比对的Needleman-Wunsch算法。本发明基于S-W算法提出一种针对短序列比对在FPGA平台实现的新方法。

当前国内外对于利用硬件平台实现序列比对主要集中在S-W算法的实现技术。主要的硬件计算平台有ASIC(Application Specific Integrated Circuits)，FPGA(FieldProgrammable Gate Arrays)，GPU(Graphics Processors Units)，和CPU(CentralProcessing Unit)等。相比于ASIC的最佳性能，制约ASIC的更是它投入研发的高额费用和它研发周期过长带来的时间成本，而市场竞争激烈需求瞬息万变，ASIC较难立住脚跟。GPU平台虽有高吞吐量，但是其存储资源受限。CPU较为常规，并且现有的生物序列比对多利用计算机集群搭配业界序列比对软件来处理，实现起来最为简单，但是由于其串行执行命令效率也最低。相比于前面的平台，虽然FPGA时钟频率一般在兆赫兹级别，但是FPGA拥有大量的逻辑资源，非常适合并行计算，并且集成了FPGA芯片的开发板外围设备齐全，特别是存储外设。这些特点使得FPGA非常适合S-W算法在其实现。

基于S-W算法的特点，现有在FPGA上实现的方法都是基于一种脉动式阵列，但规模有大有小，实现的方法也各有差异，但在各自支持的目标序列都得到了加速效果。本发明在公认的脉动式阵列基础上采用了优化的删减策略和优化的分块动态规划策略。

发明内容

本发明提供一种将S-W算法在FPGA平台上实现的新方法。该方法改进之处在于：包括以下步骤：

上位机将需要比对的目标序列、参考序列以及目标序列的相关信息从存储介质输入到FPGA器件并数字量化表示单个碱基；根据目标序列相关信息计算出所需参数，同时对目标序列和参考序列添加头部辅助碱基；准备就绪后PE阵列控制器的全局计数器开始计数并为PE阵列提供启动计算/停止计算信号，同时目标/参考序列补给模块开始协同为PE阵列提供碱基；在基于优化的S-W动态规划删减策略和优化的动态规划分块策略下，PE脉动阵列开始对目标序列和参考序列进行打分，打分的同时记录每个相似矩阵元素的回溯方向和记录相似矩阵最后一列的最大分值坐标；一组矩阵的回溯方向交替存进2个RAM中的一个，打分结束后基于优化的回溯策略进行回溯得到比对结果，并通过PCIE接口传到上层上位机。

本发明提供的新型S-W实现方法中，所说的目标序列、参考序列指的是生物信息学中需要比对的短序列(read，100-150bp)和参考基因序列，序列中每个碱基使用3位表示，碱基A、C、G和T分别用000、001、010、011表示，无效碱基使用111表示，序列相关信息包括目标序列ID、待比对目标序列长度、是否为精确匹配、在参考基因序列的地址，这些信息由上位机通过PCIE接口传到FPGA板卡的DDR存储器件中，待需要进行动态规划时信息提供给本发明的处理核心。

本发明提供的新型S-W实现方法中，需要根据输入中需要比对的序列长度，计算出实际需要动态规划的序列长度和PE阵列控制器全局计数器的计数范围。实际进行动态规划的序列长度和输入序列长度的差值作为头部辅助碱基的个数。

本发明提供的新型S-W实现方法中，PE阵列控制器提供全局计数值，为PE阵列提供启动计算/停止计算信号，这两个信号是决定PE开始和停止动态规划的控制信号。

本发明提供的新型S-W实现方法中，目标/参考序列补给模块负责给20个PE输送目标序列碱基和参考序列碱基，每个时钟周期输送一个碱基。

本发明提供的新型S-W实现方法中，基于优化的S-W动态规划删减策略，在计算相似矩阵时只计算限定范围内的元素。

本发明提供的新型S-W实现方法中，采用20个PE组成PE阵列，采用优化的分块策略进行动态规划计算，这20个PE遵循本发明所述的运作规则对长度大于80bp的目标序列进行动态规划。打分的同时记录每个元素的回溯方向和相似矩阵最后一列的最大分值以及最大分值坐标。本发明提供的新型S-W实现方法中，优化的回溯策略指的是将原本7种的回溯方向简化为最后3种方向处理，即左、上、左上3种回溯方向，即分别对应插入(Insert)、删除(Deletion)和匹配/失配(Match/Missmatch)。回溯结果效仿SAM格式中的CIGAR值，采用有限个数的CIGAR计数器来统计CIGAR值中‘M’即Match/Missmatch，‘I’即Insert，‘D’即Deletion的个数，经过整理后通过PCIE接口传到上层上位机。

附图说明

图1为本发明方法的示意框图；

图2为本发明所依赖的脉动阵列示意图；

图3为优化的动态规划分块方案示意图；

图4为传统动态规划分块方案示意图；

图5为PE模块核心电路；

具体实施方式

为了使本发明更便于理解，以下结合附图和具体实施方式来对本发明做进一步说明。

在介绍本发明之前需要介绍S-W算法及其在FPGA硬件实现的特点。S-W算法本质是给一个矩阵中的元素打分，如附图2所示，水平方向为目标序列方向，竖直方向为参考序列方向，S-W算法依赖以下一套仿射罚分公式：

其中pen_match、pen_missmatch分别是碱基匹配和失配的罚分，gapo、gape分别是一个空位打开和延展的罚分，罚分数值参考的是业界普遍采用的BWA软件提供的罚分模型，也可以自定义。

打分过程有以下特征：

1.在一个比对任务内部存在统一的数据依赖关系，即H(i,j)的计算依赖于相邻的三个元素：E(i,j)，F(i,j)和H(i-1,j-1)。只有当这三个元素结果有效，H(i,j)的计算才能够启动。因此，矩阵既可以按行也可以按列顺序计算，但是由于存在数据相关，不管采用那种方式都只能串行执行。

2.处于同一反对角线上的元素不存在数据依赖关系。

基于以上特征，因此，可以从矩阵的左上角开始，沿对角线方向，采用波前计算方法对位于反对角线上的元素实现并行计算，矩阵上方的数字是对角线的编号，同时也代表了计算的先后顺序，每一列对应一个PE模块处理，由以上特征可以看出在非边界区域H值和F值依赖于上一个PE，E值只依赖于本PE。本发明对传统的S-W算法做了优化的删减策略和优化的动态规划分块策略以及优化的回溯策略，以下将进行说明。

本发明提出的将S-W算法在FPGA平台实现的一种新方法整体框架图如附图1所示。整体结构包括PE阵列控制器1-1；序列整理模块1-2；序列补给模块1-3；20个PE组成的PE阵列1-4；参数生成模块1-5；两个RAM组成的存储模块1-6；两个回溯模块组成的1-7；以及模块之间的多路选择器。

以下说明本发明的优化删减策略。如附图3所示，目标序列长度为18bp，参考序列长度为20bp，如按照传统的结构在FPGA上实现，不分块则需要18个PE模块组成PE阵列，需要计算20*18＝360个元素的值。本发明在对足够多的DNA参考序列和目标序列进行动态规划观察之后发现，并不是全部的元素都需要被计算，回溯路径只会在对角线附近波动，故引入带宽概念，如附图元素3-4所示，本发明定义在主对角线元素上下左右四个方向一定的距离为带宽，比如此附图的带宽即为2，带宽在基因比对上可以理解为该长度的目标序列下，最大允许存在的插入或者删除的个数。在基于带宽这一概念下，沿主对角线3-1可以确定另外两条平行主对角线的线3-2与线3-3，本发明的优化删减策略规定S-W算法只对3-3右上方区域与3-2左下方区域相交元素进行动态规划，采用此策略后需要计算的元素个数为87，只为传统计算方法的24％，当然这与带宽大小和目标序列长度有关系。本发明采用优化删减策略后大大减小了计算压力。在基于本发明所述删减策略下，3-2、3-3所包含元素中除去第一行与第一列的元素外，处在3-2线上的元素所依赖的E值本来自上方，处在3-3线上的元素所依赖的F值本来自左方，因为这些元素不进行计算，采用优化的删减策略后这两种情形下所需要的E值和F值均指定为0，且经过测试，采用此方法处理后回溯结果不受影响。

以下解释本发明动态规划分块策略。如附图3与附图4，附图4是传统动态规划分块方法，附图3与附图4采用一样的目标序列和参考序列长度，对于18*20的打分矩阵，传统方法将其分为3块区域：4-4、4-5、4-6，PE数量采用6，整个打分矩阵分三次进行，只有当4-4右下角的元素完成计算后才可进行4-5区域的计算，即分块区域的衔接存在计算空白区。附图3是基于优化的删减策略下的分块方案。3-5、3-6、3-7将整个计算区域划分为4各部分，在3-5计算周期时，PE6首次进入计算，下一个计算周期PE阵列进入新的一轮计算，头部PE如PE1、PE2、PE3开始从尾部PE PE4、PE5、PE6接收所需要的H值，E值，F值等。其中深灰色区域计算依赖于中灰色区域。采用本发明所述分块方案，在计算过程中不会出现传统分块方案中的衔接空白区，整个计算过程流畅，同时采用少量PE实现长序列的动态规划，大大提高了PE的复用率和计算的效率，同时也减少了FPGA的资源消耗。

附图5为本发明PE的核心电路。采取一个时钟周期处理一组反对角线上的元素。输入信号有clk、Qry_i、Ref_i、fstH_i、H_i、H_joint_i、F_i、F_joint_i和E_joint_i，输出信号有H_o、F_o和E_o。中间信号有Mij、Hij、H_sel、F_sel、E_sel、F*、和E*。

Qry_i和Ref_i分别为目标碱基和参考碱基信号；

fstH_i是PE模块第一次进入计算所使用的H(i-1,j-1)值；

H_joint_i是衔接处来自尾部PE的H(i-1,j-1)值；

H_i是来自上一个PE的H(i-1,j-1)值；

F_i是来自上一个PE的F(i,j)值；

F_joint_i是衔接处来自尾部PE的F(i,j)值；

E_joint_i是衔接处来自尾部PE的E(i,j)值；

H_sel是当前计算周期所需要的H(i-1.j-1)值；

F_sel是当前计算周期所需要的F(i.j)值；

E_sel是当前计算周期所需要的E(i.j)值；

Mij是S-W罚分公式中的H(i-1,j-1)+S(i,j)值；

Hij是当前计算周期下的矩阵元素的H(i,j)值，由H_sel、E_sel和F_sel比较得到；

F*是为下一个计算周期准备的F(i,j)值；

E*是为下一个计算周期准备的E(i,j)值；

H_o是Hij经过两个时钟周期延迟的输出信号；

F_o是F*经过一个时钟周期延迟的输出信号；

E_o是E*经过一个时钟周期延迟的输出信号。

S-W打分主要依赖三个值H、E、F。按照上述删减策略和分块策略本发明采用20个PE组成PE阵列，PE1-PE20，带宽采用10，PE阵列在连接时PE10的H_o连到PE1的H_jiont_i，PE10的F_o接到PE1的F_jiont_i，PE11的E_o连接到PE1的E_joint_i，同理依此类推到PE2-PE10，采用相类似的模块间连接方法，将PE11-PE20的信号连接到PE2-PE10。

S-W算法包含打分和回溯两个部分，打分部分已述，以下解释本发明采用的回溯策略。如下表1所示

表格1

对于打分矩阵中的每个元素，传统回溯方向有表中的7种，本发明采取的做法是将传统方向中有岔路方向的化简为单路方向，这么做的目的是为了简化硬件电路在FPGA上的实现难度。

在存储回溯方向时，以一列即2*WIDTH+1个元素组成的数据存进RAM，等待全部的回溯方向存储完毕后开始回溯，本发明采用每两个时钟周期处理一次回溯方向，根据得到的回溯方向对相应的统计计数器进行计数。同时又根据得到的方向对取RAM数据的地址维持或者减1调整，整个目标序列回溯完成之后，仿照SAM格式文件中的CIGAR值采用108M1I31M这种表示方法来给出最终的比对结果。

Claims (6)

1.一种将S-W算法在FPGA上实现的新方法。该方法包括优化的动态规划删减策略、优化的动态规划分块策略和优化的S-W回溯方法。

2.如权利要求1所述的一种将S-W算法在FPGA上实现的新方法，其特征在于动态规划元素如何进行删减。

S-W算法本质是给相似矩阵中元素进行打分，实际上并不需要对矩阵的全部元素进行动态规划，本发明基于带宽(WIDTH)的概念，以矩阵主对角线的延长直线为主轴基准，向上移动WIDTH个元素单位确定上边界，同理向下移动同样多的元素单位可以确定下边界。上边界覆盖的元素所采用E值指定为0，下边界覆盖的元素所采用F值指定为0。被发明只对上下边界包含的元素进行动态规划，其他元素被删减。

3.如权利要求1所述的一种将S-W算法在FPGA上实现的新方法，其特征在于S-W算法相似矩阵如何分块。

4.如权利要求2所述的一种将S-W算法在FPGA上实现的新方法，在基于优化的删减策略基础上，采用较少数量的PE对100-150bp长度的目标序列进行动态规划。

本发明PE数量采用20，带宽采用10，分块界定标志为PE20开始进入计算，当PE20完成计算后，PE1-PE10开始新一轮动态规划，如此类推知道完成整个需要动态规划的区域。

5.如权利要求1所述的一种将S-W算法在FPGA上实现的新方法，其特征在于优化的S-W回溯方法。

6.本发明对传统的7种回溯方向采用简化处理，只保留3种基本的回溯方向。在打分回溯两个流程中采用流水线处理，提高硬件电路资源利用率。