对基于随机上下文无关文法的RNA二级结构预测进行加速的方法

摘要

本发明公开了一种对基于随机上下文无关文法的RNA二级结构预测进行加速的方法，目的是加快使用SCFG进行RNA二级结构预测的速度。技术方案是先构建由主机和可重构算法加速器组成的异构计算机系统，接着主机将格式化后的CM模型和编码后的RNA序列发送至可重构算法加速器，可重构算法加速器的PE阵列执行无回溯CYK/inside算法计算，计算中采用“按区域分割”和“逐层按列并行处理”的任务划分策略实现细粒度并行计算，且n个PE采用SPMD方式同时计算位于矩阵不同列的n个数据，但计算时根据不同的状态类型采用不同的计算顺序；采用本发明实现了对基于SCFG模型的RNA序列二级结构预测应用加速，加速比高，成本低。

权利要求

1.一种对基于随机上下文无关文法的RNA二级结构预测进行加速的方法，其特征在于包括以下步骤：

第一步，构建由主机和可重构算法加速器组成的异构计算机系统，主机是一台个人计算机，可重构算法加速器与主机相连，主机上安装有多序列比对程序、CM模型构建程序、CM模型解析程序、可重构算法加速器逻辑下载程序和数据通讯程序，其中CM模型是Eddy和Durbin提出的一种进行RNA二级结构分析的概率模型，它利用随机上下文无关文法SCFG从一组同源RNA序列的多序列比对结果中发现一致性的结构，以刻画RNA家族共同的结构信息；主机通过JTAG编程电缆对可重构算法加速器进行配置，将RNA序列和CM模型通过PCI-E接口加载至可重构算法加速器，并启动可重构算法加速器；可重构算法加速器完成无回溯的CYK/inside算法的计算，将比对结果返回给主机；可重构算法加速器由算法FPGA、动态存储器DRAM和PCI-E接口组成：DRAM与算法FPGA相连，用于存储原始数据-RNA序列和CM模型及中间状态的计算结果，PCI-E接口连接主机与算法FPGA；算法FPGA的JTAG配置端口通过JTAG编程电缆与主机相连，主机通过JTAG编程电缆对算法FPGA进行逻辑配置；

算法FPGA与DRAM和PCI-E接口相连，算法FPGA从PCI-E接口接收数据和操作命令，并对操作命令进行解析：如果是数据写命令，将从PCI-E接口接收的数据存储在DRAM中；

如果是启动命令，执行无回溯的CYK/inside计算，将计算结果存储在DRAM中并向主机返回计算完成信号；如果是数据读命令，将计算结果从DRAM中读出并通过PCI-E接口发送给主机；

算法FPGA由IO接口控制器、DRAM存储控制器、PE阵列控制器、PE阵列、PE同步与写回控制器组成：IO接口控制器对外与PCI-E接口相连，对内与DRAM存储控制器、PE阵列控制器相连，IO接口控制器从PCI-E接口接收原始数据即CM模型和RNA序列，并向DRAM存储控制器发出数据写请求，将原始数据通过DRAM存储控制器写入DRAM；数据存储完成后，IO接口控制器将数据准备好信号发送给PE阵列控制器；IO接口控制器接收主机发出的启动命令，并转发给PE控制模块，启动PE阵列；计算完成后，IO接口控制器接收主机发出的的数据读命令，通过DRAM存储控制器将序列比对得分从DRAM中读出，通过PCI-E接口送回主机；

DRAM存储控制器对外与DRAM相连，对内与IO接口控制器、PE阵列控制器、PE阵列、PE同步与写回控制器相连，它接收来自IO接口控制器的数据写请求，将原始数据写入DRAM；

接收来自PE阵列控制器的数据读请求，将RNA序列从DRAM读出并发送至PE阵列；接收来自PE阵列的CM模型和数据读请求，将CM模型当前状态信

说明书

对基于随机上下文无关文法的RNA二级结构预测进行加速

的方法

技术领域

[0001] 本发明涉及对基于随机上下文无关文法SCFG(Stochastic Context-free Grammars)模型的RNA二级结构预测进行加速的方法。

背景技术

[0002] 随着人类及其他各物种基因组测序项目的完成，生物信息学的研究已步入后基因组时代。序列比对和搜索已成为分子生物学领域最重要的基础性工作。目前在DNA和蛋白质序列分析领域已有许多经典的研究工具，如BLAST、FASTA、HMMER、ClustalW等。而自20世纪80年代中期具有催化性质的RNA被发现以来，RNA所起的各种重要生物化学功能逐渐引起了人们的关注。与DNA序列的双螺旋结构不同的是，虽然RNA序列本身是一条单链，但它能通过碱基互补配对形成空间二维乃至三维结构，RNA的各种重要功能与其结构直接相关，随着近年来非编码RNA在生物体内所起的重要作用被发现，针对RNA结构的研究受到人们的日益关注与重视。

[0003] 由于RNA序列结构上的保守性大于序列的保守性，有时序列本身相似度很低的两个RNA分子有可能具有很相似的结构，从而具有相近的功能，所以仅仅采用传统的序列一维结构分析工具无法满足对RNA结构特性研究的需求。另一方面，由于结构信息的加入大大增加了RNA序列分析的复杂性，直接针对三维结构的理论预测方法还很不成熟。因此研究RNA二级结构是揭示整个RNA奥秘的根本途径。关于RNA序列分析模型及二级结构预测方法，一直是近年来生物信息学研究的热点和难点问题。

[0004] 目前，最直接的RNA结构测定方法是采用X射线衍射和核磁共振，这种方法虽然结果精确可靠但是只有在拥有相关设备的环境下才能进行，这种方法所用设备非常昂贵，且非常耗时。因此采用计算机和数学模型预测RNA序列二级结构的方法被广泛采用，目前主要有三类RNA二级结构预测方法：基于热力学模型的Zuker最小自由能方法，基于比较序列分析模型的多序列比对方法和基于SCFG的结构预测方法。其中，基于SCFG理论模型的标准算法为Coche-Younger-Kasami，简称CYK算法，是目前最好的RNA二级结构预测方法之一。

[0005] 无回溯的CYK算法用于实现单条序列与单个RNA家族的共变模型CM(CovarianceModel)间的比对，从而判断该RNA序列是否属于该家族。带回溯的CYK算法还可进一步得到该序列的二级结构。CYK算法是一种三维动态规划算法，根据矩阵填充方向不同，又可分为inside和outside两种，但本质上并无不同之处。无回溯CYK/inside算法的输入为一条长度为L的RNA序列x＝x1...xL和一个CM模型。CM模型是Eddy和Durbin提出的一种进行RNA二级结构分析的概率模型，它利用SCFG从一组同源RNA序列的多序列比对结果中发现一致性的结构，以刻画RNA家族共同的结构信息。CM模型由K个不同的状态以及对应的字符生成概率ek和状态转移概率tk组成，状态信息中包括当前状态的类型、编号、父/子状态的数量以及编号。在循环变量i、j、k的控制下，无回溯CYK/inside算法实际上是不断地迭代计算一个由K个三角矩阵叠加构成的三维矩阵(每一个三角矩阵对应CM模型中的一个状态)，迭代公式如下(1≤i≤j+1，0≤j≤L，1≤k≤K)：

[0006]

[0007] 上述公式中的M是CM模型对应的三维矩阵，M(i，j，k)代表当前计算的矩阵元素，i，j，k是元素在三维矩阵中的坐标。S(k)表示当前状态的类型，其中BIF表示分支状态，其他状态如D、S、E、MP、ML、IL、MR、IR统称为非分支状态，分别代表删除、开始、结束、匹配、左匹配、左插入、右匹配和右插入。C(k)表示当前状态的子状态集合，γ代表子状态集合中的元素。xi和xj分别表示处于i、j位置的序列字符，ek和tk分别代表当前状态下的字符生成概率和子状态转移概率，变量d＝j-i+1。

[0008] 无回溯的CYK/inside算法的计算过程是从三维矩阵的最下层(每一层对应CM模型中的一个状态)边缘的元素开始，沿对角线按照波前(wavefront)顺序进行计算。每层元素的计算启动前首先判断CM模型中当前状态的类型，根据当前状态类型选择公式中的一个分支执行。当下层三角矩阵所有元素计算完毕后再计算上一层元素。计算过程由下至上逐层推进，直到CM模型的根节点(对应三维矩阵最上层)中的所有元素算完。此时元素M(1，L，1)的值就是最终的比对分值，它代表该序列与当前RNA家族的相似度。2 3 2

[0009] 无回溯的CYK/inside算法时间复杂度为O(K·L+B·L)，空间复杂度为O(LK)，其中L是RNA序列长度，B是分支状态个数，K是CM模型中的状态个数。由于巨大的时空复杂度，该算法不适合处理较大规模的RNA序列。以LSU RNA序列为例，在单个Alpha处理器上，执行无回溯的CYK/inside算法需要17391秒(约4.8小时)，带回溯的CYK/inside算法则需要27798秒(约7.7小时)，因此对该算法进行加速成为一个极具挑战性的问题。

[0010] 目前人们主要采用软件并行方案基于多处理器平台如cluster系统对CYK算法实现加速。Jane C.Hill等人于1991年提出了一种按照对角线方向并行的无回溯CYK/inside算法，在7个处理器上可以获得3.5倍的加速效果。2005年，中科院GM.Tan等人提出了对CYK/inside算法中的三维动态规划矩阵采用层内并行处理策略，对矩阵进行等面积的动态任务划分，在多个处理器上并行执行。但由于CYK算法中每个元素的计算量与其在矩阵中所处的位置相关，这种按面积等分的任务划分策略会导致处理器负载不平衡，并且随着处理器个数的增加，计算过程中的通信开销随之增长。算法复杂的数据依赖关系导致处理器间的数据通信粒度小，通信次数频繁，因而并行处理的效率不高。在由8个节点(每个节点包括4个2.4GHz AMD Opteron CPU和8GB主存，一共32个处理器)构成的Cluster上执行长度为1542的LSU RNA序列的比对只能获得19倍的加速效果，并行处理效率不到60％。新加坡南洋理工大学T.Liu等人提出在PC集群上按状态层方向流水处理的层间并行处理方案，在设计时根据元素位置与计算量之间的关系，在进行计算区域划分时采用了近似等负载的任务划分策略，即单个元素计算量较小的区域包含的列数较多，单个元素计算量较大的区域包含的列数较少，尽量让每个处理器负责计算的元素的计算量之和相等而不是元素个数相等。研究表明，新算法的并行处理效率显著提高，在包含20个Inter-Xeon 2.0GHz CPU的多处理器平台和包含48个1.0GHz Alpha EV68处理器的Cluster上执行LSU RNA序列的比对，可以达到16倍和36倍的加速效果。但是大规模并行计算机系统硬件设备的购置、使用、日常维护的成本高昂，包含20个Inter-Xeon CPU的PC Cluster系统硬件成本约为30万元，系统功耗超过4KW，而由48个Alpha处理器构成的Cluster成本超过40万元，系统功耗超过6KW，只有少数机构才有实力采购使用，非一般用户能用得上。

[0011] 近年来，FPGA(Field Programmable Gate Array)器件在容量、时钟频率和功能方面快速发展，且高端FPGA芯片的硬件成本不足2万元，最大功耗不足30W。目前已有许多基于FPGA器件对经典的序列比对算法实现硬件加速的报道，例如针对Smith-Waterman算法进行加速的有Paracel的GeneMatcher，Compugen的Bioccelerator以及TimeLogic的DeCypher等等，这些专用加速器可以达到每秒处理2亿个矩阵单元的峰值计算能力；针对BLAST算法实现硬件加速的有Mercury BLAST、RC-BLAST、Tree-BLAST、BEE2、CLC Cube以及Mitrion的BLAST Accelerator等；针对RNA二级结构预测方法实现硬件加速的有并行Zuker算法(Tan G.2006)和并行Nussinov算法(Arpith Jacob 2008)。此外，有一些硬件实现的Smith-Waterman算法已经申请了专利，例如专利号分别为5553272，5632041，5706498，5964860和6112288的美国专利等等。但是现有的序列比对、搜索以及RNA结构预测领域的硬件加速方法都仅限于对二维动态规划矩阵的并行化，未涉及三维动态规划矩阵的硬件加速，也未见有使用FPGA器件实现基于SCFG的RNA二级结构预测加速方法的相关报道。

发明内容

[0012] 本发明要解决的技术问题在于：在保证计算结果正确性，兼顾较低成本的前提下，提供一种基于FPGA平台的对基于随机上下文无关文法的RNA二级结构预测进行加速的方法，加快使用SCFG进行RNA二级结构预测的速度，使得在RNA序列数据量急剧增长的背景下，加速后的无回溯CYK/inside算法能够满足生物和医学领域研究的日常需求。

[0013] 本发明的技术方案为：利用FPGA芯片实现细粒度的资源分配与任务划分，采用多个处理单元PE(Processing Element)对无回溯CYK/inside算法执行过程进行加速。

[0014] 具体技术方案为：

[0015] 1.构建由主机和可重构算法加速器组成的异构计算机系统。主机是一台个人计算机，可重构算法加速器插入主机主板PCI-E插槽与主机相连。主机上安装有多序列比对程序、CM模型构建程序、CM模型解析程序、可重构算法加速器逻辑下载程序和数据通讯程序。主机通过JTAG(Joint Test Action Group)编程电缆对可重构算法加速器进行配置，将输入的RNA序列和CM模型通过PCI-E接口加载至可重构算法加速器，并启动可重构算法加速器；可重构算法加速器完成无回溯的CYK/inside算法的计算，将比对结果返回给主机。

[0016] 可重 构 算 法加 速 器由 算 法FPGA、动态 存 储器 DRAM(Dynamic Random AccesssMemory)和PCI-E接口组成。

[0017] DRAM与算法FPGA相连，用于存储原始数据(包括RNA序列和CM模型)和中间状态的计算结果，其存储容量要求大于1GB。PCI-E接口连接主机与算法FPGA，实现主机与算法FPGA间的数据传递。

[0018] 算法FPGA的JTAG配置端口通过JTAG编程电缆与主机相连，主机通过JTAG编程电缆对算法FPGA进行逻辑配置。

[0019] 算法FPGA与DRAM和PCI-E接口相连。算法FPGA从PCI-E接口接收数据和操作命令，并对操作命令进行解析：如果是数据写命令，将从PCI-E接口接收的数据存储在DRAM中；如果是启动命令，执行无回溯的CYK/inside计算，将计算结果存储在DRAM中并向主机返回计算完成信号；如果是数据读命令，将计算结果从DRAM中读出并通过PCI-E接口发送给主机。

[0020] 算法FPGA由IO接口控制器、DRAM存储控制器、PE阵列控制器、PE阵列、PE同步与写回控制器组成。IO接口控制器对外与PCI-E接口相连，对内与DRAM存储控制器、PE阵列控制器相连。IO接口控制器从PCI-E接口接收原始数据(CM模型和RNA序列)，并向DRAM存储控制器发出数据写请求，将原始数据通过DRAM存储控制器写入DRAM；数据存储完成后，IO接口控制器将数据准备好信号发送给PE阵列控制器；IO接口控制器接收主机发出的启动命令，并转发给PE控制模块，启动PE阵列；计算完成后，IO接口控制器接收主机发出的的数据读命令，通过DRAM存储控制器将序列比对得分从DRAM中读出，通过PCI-E接口送回主机。

[0021] DRAM存储控制器对外与DRAM相连，对内与IO接口控制器、PE阵列控制器、PE阵列、PE同步与写回控制器相连，它接收来自IO接口控制器的数据写请求，将原始数据写入DRAM；接收来自PE阵列控制器的数据读请求，将RNA序列从DRAM读出并发送至PE阵列；接收来自PE阵列的CM模型和数据读请求，将CM模型当前状态信息和计算所需的操作数从DRAM读出发送给PE阵列；接收来自PE同步与写回控制器的写请求，将由PE阵列产生的数据写回DRAM。

[0022] PE阵列控制器与IO接口控制器、DRAM存储控制器、PE同步与写回控制器和PE阵列相连。它从IO接口控制器接收数据准备好信号后启动PE阵列。PE阵列控制器向DRAM存储控制器发数据读请求，使PE阵列从DRAM获得RNA序列。PE阵列控制器通过控制当前计算元素的下标i，j，k的值实现PE阵列的初始化并对三维矩阵M实施按列循环分配和对列元素实现由上至下或者由下至上的计算。PE阵列控制器在收到PE同步与写回控制器发出的PE同步与写回完成信号后产生PE控制信号，为PE阵列分配新的计算任务，更新PE计算元素的i，j，k，启动下一列元素的计算。

[0023] PE同步与写回控制器与PE阵列和DRAM存储控制器相连，它根据当前元素的下标i，j，k判断PE的计算是否终止。由于每一列元素的计算量并不完全相等，当某个PE完成本列元素的计算后，PE同步与写回控制器控制该PE进入同步等待状态(此时阵列中的其他PE可能仍在计算)，并接收该PE发出的数据写回请求，按照先请求先服务的原则将PE阵列存储在数据缓存区中的计算结果读出，通过DRAM存储控制器写回DRAM。写回操作完成后，PE同步与写回控制器向PE阵列控制器发出同步与写回完成信号。

[0024] PE阵列负责三维动态规划矩阵的并行计算，它由n个PE、n个数据传递寄存器组和n个数据缓存区构成(n为正整数，其大小受FPGA器件片内逻辑和存储容量的限制，且在阵列规模和时钟频率之间进行折中)。每个PE都有一片存放计算结果的数据缓存区，每个PE都与DRAM存储控制器、PE阵列控制器、PE同步与写回控制器相连，相邻PE之间通过一个数据传递寄存器组相连。n个PE串联构成线性阵列，第一个PE为主PE，其他PE都是从PE，称为第1从PE，第2从PE，…，第n-1从PE。只有主PE能够向DRAM存储控制器发送数据读请求，从PE都只能被动地从DRAM存储控制器接收数据。每个数据传递寄存器组与左右相邻的PE相连(最后一个数据传递寄存器组连接第n-1从PE和主PE)，使用FPGA片内分布式存储模块实现，采用双端口设计，可以同时被本PE和下一个PE访问，用于缓存从前一个PE的数据缓存区中读出的数据，实现数据在PE阵列间的传递。所有PE的数据缓存区都与PE同步与写回控制器相连。数据缓存区使用FPGA片内集中式存储模块实现，采用双端口设计，可以同时被本PE和下一个PE访问。每个数据缓存区都由多个独立的存储区域构成，每个存储区域存储三角矩阵的一列元素，所有存储区域中存放的元素的列编号都相同，但这些列元素都位于三维矩阵M的不同层。每个存储区域都由三个寄存器标明该区域的使用状态，其中位置寄存器deck_id记录所存储的数据所处的状态层编号，父状态寄存器p_cnt记录还未计算出的本状态层对应的父状态数量(即本状态层元素还要用到的次数)，区域数据有效标志寄存器BRAM_valid表示对应区域存储的数据是否有效。

[0025] PE由PE控制模块、分值计算模块、RNA序列存储器、CM模型存储器组成。

[0026] 分值计算模块对外与该PE的数据缓存区相连，从数据缓存区中读出当前状态的子状态数据，按照无回溯CYK/inside递归公式实现对三维动态规划矩阵元素的计算，并将计算结果写入数据缓存区。分值计算模块对内与CM模型存储器和PE控制模块相连，它从CM模型存储器中读出CM模型当前状态信息，从PE控制模块或者从数据缓存区中读出计算所需的子状态数据，写入分值计算模块内部的状态寄存器中。

[0027] 分值计算模块由结束状态计算模块、非分支状态计算模块和分支状态计算模块三个子模块构成，分别实现结束状态、非分支状态(不包括结束状态)和分支状态层的计算。三个子模块彼此独立，但都与数据缓存区、CM模型存储器和PE控制模块相连。三个子模块均由状态寄存器、加法器、比较器和数据输出寄存器组成。状态寄存器与加法器相连，保存CM模型的当前状态信息。加法器与状态寄存器、PE控制模块、比较器相连，非分支状态计算模块和分支状态计算模块中的加法器还与数据缓存区相连(结束状态计算模块与数据缓存区的数据通路连接是单向的，只写数据缓存区而不从数据缓存区中读回数据)，加法器从PE控制模块和数据缓存区中载入操作数，从状态寄存器中读取当前状态下的字符生成概率ek和状态转移概率tk，执行加法运算，并将结果传递给比较器。比较器与加法器和输出寄存器相连，它从加法器获得当前的计算结果，从输出寄存器获得当前最大值，在PE控制模块的控制下将计算结果和当前最大值二者中的较大者写入输出寄存器。输出寄存器与比较器和数据缓存区相连，它保存当前状态的所有子状态下运算结果的最大值，并在计算完成后将输出寄存器中的值即元素M(i，j，k)写入数据缓存区。

[0028] RNA序列存储器与PE控制模块和CM模型存储器相连，采用FPGA片内分布式存储模块实现，用于存储经过二进制编码后的RNA序列(使用00、01、10、11表示RNA序列中的四种碱基A、C、G、U)。RNA序列存储器的每个字节可存放4个序列字符，因此RNA序列存储器容量应当大于L/4(字节)，L为按字符计数的RNA序列长度。PE控制模块将从DRAM中读出的RNA序列存入RNA序列存储器(由主PE从DRAM读出，并广播至所有PE)。

[0029] 对于FPGA的片内存储容量而言，整个CM模型的存储需求过大(超过512KB)，不宜将其整体存储在FPGA片内存储器中，因此为每一个PE设置一个CM模型存储器，它与RNA序列存储器和分支计算模块相连，容量大于2Kb，只存储CM模型当前状态下的相关信息，包括当前状态的类型Stype、编号State_id、最后一个父状态的类型P_last、父状态数量P_num、第一个子状态的类型C_first、子状态数量C_num以及RNA字符的生成概率ek和转移概率tk。

[0030] PE控制模块对外与IO接口控制器、DRAM存储控制器和数据传递寄存器组相连。它接收IO接口控制器发出的启动命令，启动PE的计算；计算时通过DRAM存储控制器访问DRAM，读回CM模型、RNA序列及中间结果；如果当前状态为匹配、右插入和右匹配，PE控制模块还将访问前一个PE的数据缓存区，从数据传递寄存器组读回返回的数据。

[0031] PE控制模块对内与RNA序列存储器、分值计算模块和PE同步与写回控制器相连。PE控制模块是一个控制状态机，负责控制PE的初始化、数据载入、分值计算、数据同步以及数据写回和缓存区数据替换，控制状态机由开始状态、PE初始化状态、读CM模型状态、状态类型判断状态、结束状态元素计算状态、非分支状态元素计算状态、分支状态元素计算状态、列元素计算完成状态、写回状态、同步等待状态、缓存区数据替换状态、结束条件判断状态、结束状态这13个状态构成。

[0032] 主机使用EDA工具对算法FPGA的逻辑进行综合，生成二进制编码的配置文件，并通过JTAG编程电缆对算法FPGA进行配置。配置过程结束后，算法FPGA自动复位，可重构算法加速器进入正常工作状态，等待主机发出操作指令。

[0033] 第二步，主机使用多序列比对程序和CM模型构建程序将输入的一组同源RNA序列转换为CM模型，然后使用模型解析程序对CM模型进行格式化：将CM模型文件中的数据按照状态编号的升序逐行顺序地存放，为每个状态分配存储单元，并对不足一行(指的是CM模型文件中的每个状态对应的数据量的最大长度)的状态数据信息进行填充补0，以解决由于CM模型中不同状态行包含的信息量不等的问题，实现CM模型在可重构算法加速器DRAM中的规整存储，从而简化地址换算。格式化完成后，主机通过PCI-E接口将CM模型发送至可重构算法加速器的DRAM。

[0034] 第三步，主机对输入的RNA序列进行二进制编码，使用2位二进制数00、01、10、11对RNA序列中的A、C、G、U字符进行编码，并将编码后的序列发送至可重构算法加速器的DRAM。

[0035] 第四步，主机向可重构算法加速器发送启动命令，启动算法FPGA，主机进入等待状态。启动命令由操作码、RNA序列长度、RNA序列在DRAM中存放的基地址、CM模型状态数量、CM模型在DRAM中存放的基地址、中间结果和最终计算结果在DRAM中存放的基地址组成。算法FPGA接收到启动命令后启动PE阵列执行无回溯CYK/inside算法计算，计算过程包括以下步骤：

[0036] 4.1 PE初始化，算法FPGA在收到主机发出的启动命令后，立即启动PE阵列的初始化过程，具体步骤如下：

[0037] 4.1.1所有PE在各自的数据缓存区中选择编号最小的一个空闲区域保存计算结果(每个PE每次负责计算三角矩阵的一列)；

[0038] 4.1.2主PE通过DRAM存储控制器从DRAM中依次读出RNA序列和CM模型第K行数据，并发送至所有的从PE，从PE接收RNA序列和CM模型第K行数据分别存储在各自的RNA序列存储器和CM模型存储器中；

[0039] 4.1.3所有PE根据各自在PE阵列中的编号为将要计算的元素下标赋初始值：i＝P_id，j＝P_id，k＝K。i，j，k为元素M(i，j，k)在三维矩阵中的坐标，i，j为元素在矩阵中的行/列坐标，k为状态编号，P_id为PE的编号，其中主PE的P_id＝1，从PE的P_id为2～n之间的整数，K为CM模型中的状态数量同时也是矩阵最底层的编号，n为PE的数量。

[0040] 4.2数据载入：

[0041] 各PE根据当前状态的类型将操作数载入分值计算模块，一共分为五种情况：

[0042] ●如果当前状态类型为S、D、IL、ML，则PE分值计算模块从该PE的数据缓存区中载入当前状态的所有子状态数据层中的第j列元素；

[0043] ●如果当前状态类型为分支(BIF)状态，且P_id＝1，则PE从以下两个位置获得计算所需的操作数：1)从本地数据缓存区中载入分支状态右子状态层的第j列元素；2)向DRAM存储控制器发读请求，从DRAM中载入分支状态左子状态层的第i行元素，并将其广播至所有PE；

[0044] ●如果当前状态类型为分支(BIF)状态，且2≤P_id≤n，则PE从以下两个位置获得计算所需的操作数：1)从PE数据缓存区中载入右子状态数据层的第j列元素；2)等待接收主PE从DRAM中载入的左子状态数据层的第i行元素；

[0045] ●如果当前状态类型为IR、MR或者MP(分别对应CM模型中的右插入、右匹配和匹配状态)，且P_id＝1，则PE向DRAM存储控制器发读请求，从DRAM载入当前状态的所有子状态数据层的第j-1列元素；

[0046] ●如果当前状态类型为IR、MR或者MP，且2≤P_id≤n，则从前一个PE的数据缓存区中载入当前状态的所有子状态数据层的j-1列元素；

[0047] 4.3计算比对得分

[0048] 串行无回溯CYK/inside算法的计算过程是从三维矩阵最下层三角矩阵的主对角线元素开始，沿对角线按照波前顺序进行计算。当下层三角矩阵元素全部填充完后再计算上一层的元素，但这种方法在FPGA平台上并不适用。首先串行无回溯CYK/inside算法多维度、非一致性的数据依赖关系导致计算过程中任务划分困难，难以实现PE之间的负载平衡；其次，串行无回溯CYK/inside程序的空间局部性差(访存粒度小(数据量以KB为单位)，每次读写操作的数据量不等且在空间上不连续)，导致访存调度困难；第三，FPGA片3
内存储资源不能满足矩阵O(N)的存储需求，导致片外存储器访问延时将成为提高并行处理效率的瓶颈。计算过程中三维矩阵每一层元素的存储需求约为32MB，而目前容量最大的FPGA器件片内存储容量不足2MB，因此有限的片内存储容量和相对较低的存储带宽是基于FPGA平台对无回溯CYK/inside算法进行加速面临的挑战。

[0049] 为此，本发明采用“按区域分割”和“逐层按列并行处理”的任务划分策略实现对三维动态规划矩阵的细粒度并行计算，方法是：

[0050] 4.3.1首先对三维矩阵按状态层方向进行纵切，将其垂直分割为s个区域( L为RNA序列的长度)，每一个区域仍是一个三维矩阵，每个区域都包含K层，每一层包含n列元素；然后对每一个区域逐层计算：从第1个区域的第K层出发，当第K层计算完成后再算第K-1层，直到第1个区域的第1层算完，然后再计算第2个区域的第K层，第K-1层，…，第1层；之后再计算第3个区域…，直到最后一个区域的第1层算完。对每一层的计算而言，每个PE负责计算当前状态层中的一列，PE计算的元素的列号与PE在阵列中的序号P_id对应。主PE的编号P_id＝1，负责计算本区域第1列元素，从PE计算本区域的第2～n列元素。

[0051] 4.3.2阵列中的n个PE采用单程序多数据流(Single Program Multiple Data)方式同时计算位于矩阵不同列的n个数据，以达到对计算过程进行加速的效果。PE都使用无回溯CYK/inside算法的迭代公式对矩阵元素进行计算，但计算时根据不同的状态类型采用不同的计算顺序，方法是：

[0052] 4.3.2.1如果当前状态k是分支状态，则每一个PE都从各自列的顶部位置(i＝1)开始按照由上至下的顺序计算。当前状态层的计算启动时，主PE计算元素(1，j)，第1从PE计算元素(1，j+1)，…，第n-1从PE计算元素(1，j+n-1)，j为元素在矩阵中的列坐标，1≤j≤L。当第1行的n个元素计算完之后，所有PE同时更新元素的横坐标，计算下一行(i＝2)的n个元素。每一个PE由上至下串行计算本列元素，整个PE阵列每次完成本状态层中的一行元素的计算。当一个PE计算的元素坐标i＝j时(即到达各自列的底部位置)，PE计算暂停，进入等待状态。由于三角矩阵每一列元素的个数不等，PE按编号顺序依次进入同步等待状态，直到收到最后一个从PE的计算完成信号。特别地，当n＝1时，此时没有从PE，主PE按照从左至右的顺序对矩阵元素进行逐列计算，在每一列内部按照由上至下的顺序对每一个元素进行计算。

[0053] 4.3.2.2如果当前状态k是非分支状态，则每一个PE都从各自列的底部位置(i＝j)开始，按照由下至上的顺序计算。当前状态层的计算启动时，所有PE计算的元素都位于三角矩阵的主对角线上：主PE计算元素(i，j)，第1从PE计算元素(i+1，j+1)，…，第n-1从PE计算元素(i+n-1，j+n-1)，1≤i≤L，i≤j≤L。位于同一条对角线上的元素的计算量都相等，因此所有PE都可以同时完成当前元素的计算，当前对角线上的n个元素计算完之后，所有PE同时向上移动一个计算位置，主PE计算元素(i-1，j)，第1从PE计算元素(i，j+1)，…，第n-1从PE计算元素(i+n-2，j+n-1)。当PE计算的元素行坐标i＝1时，PE计算暂停，进入等待状态(如果计算结果需要写回DRAM，PE将在等待状态下发出写回请求)。同样，由于三角矩阵每一列元素的个数不等，PE将按编号顺序依次进入等待状态。主PE算完元素(1，j)后进入等待状态，此时第1从PE计算元素(1，j+1)，…，第n-1从PE计算元素(n-1，j+n-1)。第1从PE计算完(1，j+1)后进入等待状态，此时第2从PE计算元素(1，j+2)，…，第n-1从PE计算元素(n-2，j+n-1)，第2从PE计算完元素(1，j+2)后也进入等待状态，…，直到第n-1从PE的算完元素(1，j+n-1)后本区域当前状态层的计算结束。特别地，当n＝1时，此时没有从PE，主PE按照从左至右的顺序对矩阵元素进行逐列计算，在每一列内部按照由下至上的顺序对每一个元素进行计算。由于只有一个PE，此时也就不存在同步等待，当一列元素计算完成后如果需要写回，则直接发出写回请求。

[0054] 4.4所有PE的列元素计算结果都保存在PE各自的数据缓存区中，PE根据状态类型确定是否将数据写回DRAM。

[0055] 4.4.1如果当前状态是分支状态的左子状态，则所有PE都将计算结果写回DRAM(因为分支的左子状态数据在未来很长一段时间内不会被使用)。方法是：PE完成本列元素的计算后向PE同步与写回控制器发数据写回请求，随后PE同步与写回控制器将数据从PE的数据缓存区中读出，通过DRAM存储控制器将数据写回DRAM。如果阵列中同时有多个PE发数据写回请求，PE同步与写回控制器依次将各PE的数据写回DRAM。

[0056] 4.4.2如果当前状态是匹配、右插入或者右匹配状态，则只有最后一个从PE将计算结果写回DRAM(因为该PE的计算结果位于本区域当前状态层的的最右侧，在下一个区域同一状态层计算时将会被用到)，其它PE的计算结果都保存在各自数据缓存区中，以便在后续的计算中被本PE使用或通过数据传递寄存器组传递给下一个PE。

[0057] 4.5对各PE进行同步。由于三角矩阵每一列的元素个数和计算量不等，导致PE间的负载并不完全平衡，当阵列中的某个PE完成本列元素的计算时，阵列中该PE之后的处理单元仍然处于计算状态。为了保持数据依赖关系的一致性，该PE在完成本列元素的计算之后进入同步等待状态。PE在同步等待状态下收到第n-1从PE发出的当前状态层计算完成信号后进入缓存区数据替换状态。PE在缓存区数据替换状态下判断是否存在某个存储区域的父状态寄存器p_cnt数值为0，若为0则释放该状态占用的数据缓存区。

[0058] 4.6加载新任务。当每个PE都收到第n-1从PE(P_id＝n)的计算完成信号时，PE阵列控制器将本区域中的下一个状态层加载至PE阵列，为每一个PE分配新状态层中的一列元素，即为变量i，j，k赋值：如果k＞0，则i＝P_id，j＝P_id，k＝k-1，载入CM模型的第k-1行数据；否则，i＝P_id，j＝j+n，k＝K，载入CM模型的第K行数据，随后PE控制状态机进入结束条件判断状态。如果j＞L，则进入结束状态(此时当前PE处于空闲状态，直到阵列中的所有PE计算完成)，否则在数据缓存区寻找一个新的编号最小的空闲区域，然后转4.2，对下一个状态层进行计算；

[0059] 4.7当所有PE都进入结束状态后，主PE将最终的计算结果(元素M(1，L，1)的值，它代表最终的比对得分)写回DRAM，算法FPGA向主机发计算完成信号，PE计算过程结束。PE控制模块将所有寄存器清0，回到开始状态，等待主机发新的启动命令。

[0060] 第五步，主机接收到可重构算法加速器返回的计算完成信号后通过PCI-E接口向可重构算法加速器发数据读命令，算法FPGA接收数据读命令将计算结果从DRAM中读出，通过PCI-E接口转发给主机，RNA序列比对过程结束。

[0061] 与现有技术相比，采用本发明可达到以下技术效果：

[0062] 1.本发明通过对无回溯CYK/inside算法中三维动态规划矩阵分值计算进行区域分割、区域内逐层计算、层内按矩阵列顺序并行处理，对计算过程进行细粒度并行化，实现了对基于SCFG模型的RNA序列二级结构预测应用的加速。通过与运行在Intel Pentium42.66GH CPU、1.5GB主存计算机上的比对软件Infernal-0.55比较，采用本发明在RNA序列长度为379，CM模型状态数为1176时，可以获得接近14倍的加速效果，而当测试序列长度为959，CM模型状态数为3145时，可获得14.5倍的加速效果。

[0063] 2.就基于SCFG模型的RNA二级结构预测应用而言，本发明以很小的硬件代价(1片FPGA和一台通用计算机)取得了与包含20个Inter-Xeon CPU的并行计算机相当的加速效果。包含20个Inter-Xeon CPU的并行计算机系统功耗超过3KW，功耗分析结果表明，算法FPGA芯片的最大功耗小于30W，整个加速器系统平台(包括主机和可重构算法加速器)的总功耗不超过300W，仅为前者系统功耗的10％左右。

[0064] 3.本发明将可编程的FPGA芯片作为通用计算机的协处理器，具有成本低、计算部件密集、性价比高的优点；同时FPGA算法加速器具备的重构能力能够灵活地适应算法的改变，克服了传统专用硬件造价高、普及率低的缺点。

附图说明

[0065] 图1是本发明第一步构建的异构计算机系统总体结构图。

[0066] 图2是本发明异构计算机系统中可重构算法加速器中算法FPGA逻辑结构图。

[0067] 图3是本发明可重构算法加速器中算法FPGA中PE的逻辑结构图。

[0068] 图4是本发明各PE的PE控制模块状态转换图。

[0069] 图5是本发明异构计算机系统进行RNA序列比对的流程图。

[0070] 图6是图5第6步中PE阵列执行无回溯CYK/inside计算的流程图。

[0071] 图7是图6中PE阵列计算比对得分步骤的三维矩阵填充过程示意图。

具体实施方式

[0072] 图1是本发明第一步构建的异构计算机系统总体结构图。异构计算机系统由主机和可重构算法加速器组成。可重构算法加速器插入主机主板PCI-E插槽与主机相连。主机上安装有CM模型格式化程序、RNA序列编码程序可重构算法加速器逻辑下载程序和数据通讯程序。主机通过JTAG编程电缆对可重构算法加速器中的算法FPGA进行配置，将RNA序列和CM模型通过PCI-E接口加载至可重构算法加速器，并启动可重构算法加速器；可重构算法加速器完成无回溯的CYK/inside算法的计算，然后将比对结果返回给主机。

[0073] 可重构算法加速器由算法FPGA、DRAM和PCI-E接口组成。DRAM与算法FPGA相连，存储原始数据(包括RNA序列和CM模型)和中间状态的计算结果，其存储容量大于1GB。PCI-E接口连接主机与算法FPGA，实现主机与算法FPGA间的数据传递。JTAG编程电缆连接主机USB接口与算法FPGA的JTAG配置端口。

[0074] 算法FPGA与DRAM和PCI-E接口相连。算法FPGA从PCI-E接口接收数据和操作命令，并对命令进行解析：如果是数据写命令，则从将PCI-E接口接收的数据存储在DRAM中；如果是启动命令，则执行无回溯的CYK/inside计算过程，将计算结果存储在DRAM中并向主机返回计算完成信号；如果是数据读命令，则将计算结果从DRAM读回通过PCI-E接口发送给主机。

[0075] 图2是本发明异构计算机系统中可重构算法加速器中算法FPGA逻辑结构图。算法FPGA由IO接口控制器、DRAM存储控制器、PE阵列控制器、PE阵列、PE同步与写回控制器组成。IO接口控制器对外与PCI-E接口相连，对内与DRAM存储控制器、PE阵列控制器相连。

[0076] IO接口控制器从PCI-E接口接收原始数据，并向DRAM存储控制器发数据写请求，将原始数据写入DRAM；数据存储完成后，IO接口控制器将数据准备好信号发送给PE阵列控制器；计算完成后，IO接口控制器接收主机发出的的数据读命令，通过DRAM存储控制器将序列比对得分从DRAM读出，通过PCI-E接口送回主机。

[0077] DRAM存储控制器对外与DRAM相连，对内与IO接口控制器、PE阵列控制器、PE阵列、PE同步与写回控制器相连，负责实现对DRAM的访问。

[0078] PE阵列控制器与IO接口控制器、DRAM存储控制器、PE同步与写回控制器和PE阵列中的所有PE相连。它从IO接口控制器接收数据准备好信号后启动PE阵列，它向DRAM存储控制器发数据读请求，使所有的PE处理单元从DRAM获得RNA序列。在计算比对得分阶段，PE阵列控制器通过控制当前计算元素的下标i，j，k对三维矩阵实施按列循环分配并通过对i，j，k的赋值控制PE的初始计算位置和元素的计算顺序；在PE中间结果写回阶段，收到PE同步与写回控制器发出的PE同步与写回完成信号后为PE分配新的计算任务；在加载新任务阶段，通过对i，j，k的判断决定PE是否终止计算。

[0079] PE同步与写回控制器与PE阵列中每一个PE以及PE的数据缓存区相连，可以直接访问数据缓存区。PE同步与写回控制器根据i，j，k监控每个PE当前计算的元素所处的位置。

[0080] PE阵列由n个PE、n个数据传递寄存器组和n个数据缓存区构成。每个PE都与DRAM存储控制器、PE阵列控制器、PE同步与写回控制器、数据缓存区以及相邻的数据传递寄存器组相连。n个PE串联构成线性阵列，其中阵列左侧第一个PE为主PE，其他PE都是从PE，称为第1从PE，第2从PE，…，第n-1从PE。只有主PE能够向DRAM存储控制器发送数据读请求，从PE都被动地从DRAM存储控制器接收数据。每个数据传递寄存器组与左右相邻的PE相连(最后一个数据传递寄存器组连接第n-1从PE和主PE)，可以同时被本PE和下一个PE访问，用于缓存从前一个PE的数据缓存区中读出的数据，实现数据在PE阵列间的传递。每个PE都有一片存放计算结果的数据缓存区，所有PE的数据缓存区都与PE同步与写回控制器相连，在中间结果写回阶段可被PE同步与写回控制器直接访问。数据缓存区可以同时被本PE和下一个PE访问，用于存储本PE的计算结果。每个PE的数据缓存区都由多个独立的存储区域构成，每个存储区域存储三角矩阵的一列元素，所有存储区域存放的元素的列编号都相同，但这些列元素位于三维矩阵M的不同层。

[0081] 图3是本发明可重构算法加速器中算法FPGA中PE的逻辑结构图。PE由PE控制模块、分值计算模块、RNA序列存储器、CM模型存储器组成。

[0082] 分值计算模块对外与数据缓存区相连，从数据缓存区中读出当前状态的子状态数据，将计算结果写入数据缓存区；对内与CM模型存储器和PE控制模块相连，它从CM模型存储器中读出CM模型当前状态信息，写入分值计算模块内部的状态寄存器中，从PE控制模块或数据缓存区中载入计算所需的子状态数据，输入分值计算模块内部的加法器。

[0083] 分值计算模块由结束状态计算模块、非分支状态计算模块和分支状态计算模块三个子模块构成，分别实现结束状态、非分支状态(不包括结束状态)和分支状态层的计算。三个子模块彼此独立，但都与数据缓存区、CM模型存储器和PE控制模块相连。三个子模块均由状态寄存器、加法器、比较器和输出寄存器组成。状态寄存器与加法器相连，保存CM模型的当前状态信息。加法器与状态寄存器、PE控制模块、比较器相连，非分支状态计算模块和分支状态计算模块中的加法器还与数据缓存区相连(结束状态计算模块与数据缓存区的数据通路连接是单向的，只写数据缓存区而不从数据缓存区中读回数据)，非分支状态计算模块和分支状态计算模块的加法器从PE控制模块和数据缓存区中载入操作数，结束状态计算模块的加法器从PE控制模块载入操作数，从状态寄存器中读取当前状态下的字符生成概率ek和状态转移概率tk，执行加法运算，并将结果传递给比较器。比较器与加法器和输出寄存器相连，它从加法器获得当前的计算结果，从输出寄存器获得当前最大值，在PE控制模块的控制下将计算结果和当前最大值二者中的较大者写入输出寄存器。输出寄存器与比较器和数据缓存区相连，它保存当前状态的所有子状态下运算结果的最大值，并在计算完成后将输出寄存器中的值即元素M(i，j，k)写入数据缓存区。

[0084] RNA序列存储器与PE控制模块和CM模型存储器相连，采用FPGA片内分布式存储模块实现，用于存储经过二进制编码后的RNA序列。PE控制模块将从DRAM中读出的RNA序列存入RNA序列存储器(由主PE从DRAM读出，并广播至所有PE)。

[0085] CM模型存储器与RNA序列存储器和分值计算模块相连，存储CM模型中当前状态k的相关信息。

[0086] PE控制模块对外与IO接口控制器、DRAM存储控制器和数据传递寄存器组相连。它接收IO接口控制器发出的启动命令，启动PE。在PE初始化阶段，接收从DRAM中载入的CM模型状态数据，并将其缓存在CM模型存储器中。计算时，根据当前状态的类型访问DRAM和前一个PE的数据缓存区获取计算所需的数据。PE控制模块对内与RNA序列存储器、分值计算模块和PE同步与写回控制器相连。当前状态层的计算完成后，PE控制模块根据当前状态类型判断是否需要将计算结果写回DRAM。如果需要写回，则向PE同步与写回控制器发出写回请求，并等待接收PE同步与写回控制器发出的数据写回完成信号。写回操作完成后，PE控制模块发出当前状态层计算完成信号，并将存放子状态数据的存储区域对应的父状态寄存器数值减1。如果某个存储区域的父状态寄存器数值减为0，则将区域数据有效标志寄存器清0，并将该区域编号加入空闲存储区域集合。当PE控制模块接收到同步等待完成(即最后一个从PE发出的当前状态层计算完成)信号后，更新元素下标i，j，k的值，启动下一个状态层的计算。

[0087] 图4是本发明各PE的PE控制模块状态转换图。PE控制模块是一个控制状态机，由13个状态构成，状态定义和转换过程如下：

[0088] 开始状态：算法FPGA配置结束并复位之后，控制状态机自动进入开始状态，收到IO接口控制器发出的启动命令后进入PE初始化状态。

[0089] PE初始化状态：在数据缓存区中选择编号最小的一个空闲区域保存计算结果；向DRAM存储控制器发读请求，从DRAM载入RNA序列，并写入RNA序列存储器；为PE计算的元素下标赋初始值：i＝P_id，j＝P_id，k＝K(默认为从三维矩阵M最底层的边缘位置开始计算)。初始化操作完成后进入读CM模型状态。

[0090] 读CM模型状态：该状态下，PE控制模块向DRAM存储控制器发出读请求，从DRAM载入CM模型当前状态k的信息，并写入CM模型存储器，随后进入状态类型判断状态。

[0091] 状态类型判断状态：将当前状态编号k填入在PE初始化状态选择的空闲区域对应的位置寄存器deck_id中；并判断当前状态的类型，如果是结束状态，则进入结束状态元素计算状态；如果是非分支状态，则进入非分支状态元素计算状态；如果是分支状态，则进入分支状态元素计算状态，并给PE计算的元素下标赋新值：i＝1，j＝P_id，k＝K。

[0092] 结束状态元素计算状态：在本状态下启动结束状态计算子模块，按照CYK/inside算法迭代公式中的结束状态计算分支对矩阵元素M(i，j，k)进行计算；结束状态元素的计算只与元素位置相关，不需要从DRAM或者数据缓存区载入数据，产生的计算结果保存在数据缓存区，然后根据元素下标判断本列元素的计算是否完成。如果i＝1，则进入列元素计算完成状态；否则继续保持本状态，i＝i-1，计算下一个元素。

[0093] 非分支状态元素计算状态：在本状态下启动非分支状态计算子模块，按照CYK/inside算法迭代公式中的非分支状态(不包括结束状态)计算分支对矩阵元素M(i，j，k)进行计算。如果当前状态为开始状态、删除状态或左插入状态，则从本PE的数据缓存区载入操作数；如果当前状态为匹配或右插入状态，主PE从DRAM载入操作数，从PE从前一个PE的数据缓存区载入载入操作数。计算结果保存在数据缓存区，然后根据元素下标判断本列元素的计算是否完成。如果i＝1，则进入列元素计算完成状态；否则继续保持本状态，更新元素下标i＝i-1，计算下一个元素。

[0094] 分支状态元素计算状态：在本状态下启动分支状态计算子模块，按照CYK/inside算法迭代公式中的分支状态计算分支对矩阵元素M(i，j，k)进行计算。分支状态计算模块从DRAM中载入左子状态的一行，从PE的数据缓存区中载入右子状态的一列数据，计算结果保存在数据缓存区，然后根据元素下标判断本列元素的计算是否完成。如果i＝P_id，则进入列元素计算完成状态；否则继续保持本状态，元素下标i＝i+1，继续计算下一个元素。

[0095] 列元素计算完成状态：将存放子状态数据的存储区域对应的父状态寄存器数值减1；并根据当前计算的CM模型状态类型判断是否需要将保存在数据缓存区中的元素写回DRAM。如果当前状态是分支状态的左子状态，则所有PE的PE控制模块都向PE同步与写回控制器发出写回请求，随后进入写回状态；如果当前状态是匹配、右插入或者右匹配状态，则最后一个从PE的PE控制模块向PE同步与写回控制器发写回请求，随后最后一个从PE的PE控制模块进入写回状态，其它PE控制状态机进入同步等待状态；如果当前状态是左插入、左匹配、删除或者结束状态，则不需要将数据缓存区中的元素写回DRAM，所有PE的PE控制状态机直接进入同步等待状态。

[0096] 写回状态：PE在该状态下保持数据写回请求，直到PE同步与写回控制器响应该请求，随后PE同步与写回控制器将数据从PE数据缓存区中读出，通过DRAM存储控制器将数据写回DRAM。在收到PE同步与写回控制器返回的写回完成信号后向阵列中的其他PE发出当前状态层计算完成信号，进入缓存区数据替换状态。

[0097] 同步等待状态：该状态下PE处于空闲状态，在收到最后一个PE发出的当前状态层计算完成信号(表明整个PE阵列计算和写回操作都已完成)后进入缓存区数据替换状态。

[0098] 缓存区数据替换状态：PE控制状态机判断当前使用的存储区域的父状态寄存器p_cnt的值是否为0，如果为0(说明该存储区域所存放状态的所有父状态都已算完，该区域内的数据将不再被使用)，则将该存储区域的数据有效标志寄存器清0，并将该存储区域编号加入空闲存储区域集合；更新元素下标i，j，k的值：如果k＞0，则i＝P_id，j＝P_id，k＝k-1；否则i＝P_id，j＝j+n，k＝K，随后进入结束条件判断状态。

[0099] 结束条件判断状态：PE控制状态机判断分值计算模块终止条件：如果j＞L，则进入结束状态，否则在数据缓存区中选择编号最小的一个空闲区域，然后进入读CM模型状态，准备对下一个状态层中的列元素进行计算。

[0100] 结束状态：PE控制状态机将所有寄存器清0，回到开始状态，等待接收IO接口控制器转发主机发出的新的启动命令，开始下一轮计算。

[0101] 图5是本发明异构计算机系统进行RNA序列比对的流程图：

[0102] 第一步，构建由主机和可重构算法加速器组成的异构计算机系统。

[0103] 序列比对流程对应本发明的第二、三、四、五步，为了清楚地表明主机端和可重构算法加速器的工作，则可分为8个小步骤。

[0104] 第二步，包括1、2小步：

[0105] 1.主机使用多序列比对和CM模型构建程序将输入的一组同源RNA序列转换为CM模型，然后使用模型解析程序格式化CM模型，并将其通过PCI-E接口发送至可重构算法加速器。

[0106] 2.可重构算法加速器算法FPGA接收主机发送的CM模型，将其存储在DRAM中，然后向主机返回数据接收完成信号。

[0107] 第三步，包括3、4小步：

[0108] 3.主机对输入的RNA序列进行二进制编码，并将经过编码的RNA序列发送至可重构算法加速器。

[0109] 4.可重构算法加速器算法FPGA接收主机发送的RNA序列，将其存储在DRAM中，并向主机返回数据接收完成信号。

[0110] 第四步，包括5、6小步：

[0111] 5.主机收到可重构算法加速器返回的RNA序列接收完成信号后，向可重构算法加速器发送启动命令，启动算法FPGA执行计算，主机进入等待状态。

[0112] 6.算法FPGA接收到主机发送的启动命令随即启动PE阵列执行无回溯CYK/inside计算。算法FPGA计算完成后，将结果写回DRAM，随后向主机发出计算完成信号。

[0113] 第五步，包括7、8小步：

[0114] 7.主机在等待状态下接收可重构算法加速器发出的计算完成信号后通过PCI-E接口向算法加速器发出数据读命令。

[0115] 8.可重构算法加速器算法FPGA接收到主机发出的数据读命令，将计算结果从DRAM中读出并传给主机，比对过程结束。

[0116] 图6是图5第6步中PE阵列执行无回溯CYK/inside计算的流程图。图中的符号s表示分割后得到的较小的三维矩阵(即计算区域)的数量( )；m表示当前正在计算的区域的编号(1≤m≤s)；K表示CM模型的状态数量，即三维矩阵的层数；k表示当前正在计算的状态层的编号(1≤k≤K)；d为当前最小的空闲存储区域编号(1≤d≤n)；n为PE的数量；S(k)表示当前状态的类型；P(k)表示当前状态k的父状态；L为RNA序列的长度；D、S、E、MP、IL、ML、IR、MR和BIF表示具体的状态类型，分别代表CM模型中的删除、开始、结束、匹配、左插入、左匹配、右插入、右匹配和分支状态，其中除了分支状态BIF外，其余的状态统称为非分支状态。PE阵列的操作流程包括初始化、数据载入、计算比对得分、写回中间结果、PE阵列同步和加载新任务六个步骤：

[0117] 1.初始化。算法FPGA在收到主机发出的操作指令后，立即启动PE阵列的初始化：

[0118] 1.1所有PE在各自的数据缓存区中选择编号最小的空闲区域d用于保存计算结果；

[0119] 1.2主PE向DRAM存储控制器发数据读请求，DRAM存储控制器从DRAM依次读出RNA序列和CM模型第K行数据，并发送至所有的从PE；从PE接收RNA序列和CM模型第K行数据并分别存储在各自的RNA序列存储器和CM模型存储器中；

[0120] 1.3确定将要计算的元素的初始位置。所有PE根据各自在PE阵列中的编号P_id为将要计算的元素下标赋初始值：i＝P_id，j＝P_id，k＝K。

[0121] 2.数据载入。

[0122] 各PE根据当前状态的类型将操作数载入分值计算模块，一共分为五种情况：

[0123] ●如果当前状态类型为S、D、IL、ML，则PE分值计算模块从该PE的数据缓存区中载入当前状态的所有子状态数据层中的第j列元素；

[0124] ●如果当前状态类型为分支(BIF)状态，且P_id＝1，则PE从以下两个位置获得计算所需的操作数：1)从数据缓存区中载入分支状态右子状态层的第j列元素；2)向DRAM存储控制器发读请求，从DRAM中载入分支状态左子状态层的第i行元素，并将其广播至所有PE；

[0125] ●如果当前状态类型为分支(BIF)状态，且2≤P_id≤n，则PE从以下两个位置获得计算所需的操作数：1)从PE数据缓存区中载入右子状态数据层的第j列元素；

[0126] 2)等待接收主PE从DRAM中载入的左子状态数据层的第i行元素；

[0127] ●如果当前状态类型为IR、MR或者MP(分别对应CM模型中的右插入、右匹配和匹配状态)，且P_id＝1，则PE向DRAM存储控制器发读请求，从DRAM载入当前状态的所有子状态数据层的第j-1列元素；

[0128] ●如果当前状态类型为IR、MR或者MP，且2≤P_id≤n，则从前一个PE的数据缓存区中载入当前状态的所有子状态数据层的j-1列元素；

[0129] 3.计算比对得分

[0130] 3.1将三维矩阵按状态层方向垂直分割为s个区域( )，按照编号依次进行计算。对每一个区域m(1≤m≤s)采用逐层计算的策略：从区域m的第K层出发，当第K层计算完成后再算第K-1层，…，直到区域m的第1层算完。

[0131] 3.2对每一层的计算而言，每个PE负责计算当前状态层中的一列，PE计算的元素的列号与PE在阵列中的序号P_id对应。主PE的编号P_id＝1，负责计算本区域第1列元素，从PE计算本区域的第2～n列元素。所有PE都同时按照无回溯CYK/inside算法的迭代公式对矩阵元素进行计算，但要根据不同的状态类型采用不同的计算顺序：

[0132] 3.2.1如果当前状态k是分支状态，则每一个PE都从各自列的顶部位置(i＝1)开始，按照由上至下的顺序计算。整个PE阵列每次完成本状态层中一行元素的计算。当一个PE计算的元素坐标i＝j时(即到达各自列的底部位置)，PE计算暂停，进入等待状态。由于三角矩阵每一列元素的个数不等，PE按编号顺序依次进入同步等待状态，直到收到最后一个从PE的计算完成信号。

[0133] 3.2.2如果当前状态k是非分支状态，则每一个PE都从各自列的底部位置(i＝j)开始，按照由下至上的顺序计算。区域m中的第k个状态层的计算启动时，所有PE计算的元素都位于三角矩阵的主对角线上。当主对角线上的n个元素计算完之后，所有PE同时向上移动一个计算位置。当PE计算的元素行坐标i＝1时，PE计算暂停，进入等待状态(如果计算结果需要写回DRAM，PE将在等待状态下发出写回请求)。由于三角矩阵每一列元素的个数不等，PE将按编号顺序依次进入同步等待状态。直到第n-1从PE完成本列最上方元素的计算(i＝1)并写回DRAM(如果需要写回)后，本区域当前状态层k的计算结束。

[0134] 4.写回中间结果。所有PE的列元素计算结果都保存在PE各自的数据缓存区中，PE根据状态类型确定是否将数据写回FPGA片外的DRAM。

[0135] 5.对PE阵列中的PE进行同步，直到所有PE完成本状态层的计算。PE进行同步等待的同时判断当前使用的存储区域的父状态寄存器p_cnt的值是否为0，如果为0，则释放该缓存区。

[0136] 6.加载新任务。当收到最后一个PE(PE[n])的计算完成信号后(说明此时本状态层的计算任务已完成)，PE阵列控制模块将下一个状态层加载至PE阵列，为每一个PE分配新状态层中的一列元素，并对变量i，j，k进行赋值：如果k＞0，则i＝P_id，j＝P_id，k＝k-1，载入CM模型的第k-1行数据；否则，i＝P_id，j＝j+n，k＝K，载入CM模型的第K行数据，随后PE控制状态机进入结束条件判断状态。如果j＞L，则当前PE的计算终止，否则在数据缓存区寻找当前编号最小的一个空闲存储区域，回到第2步(数据载入)，准备对下一个状态层进行计算。

[0137] PE重复上述计算过程直到新指派的列元素编号j超出输入的RNA序列的长度(此时当前PE处于空闲状态，直到阵列中的所有PE计算完成)。当所有PE的计算都终止后，主PE将最终的计算结果(即元素M(1，L，1)的值)写回DRAM，算法FPGA向主机发计算完成信号，PE计算过程结束。PE控制模块状态机将所有寄存器清0，回到开始状态。

[0138] 图7是图6中PE阵列计算比对得分步骤的三维矩阵填充过程示意图。图中的符号K表示CM模型的状态数量，即三维矩阵的层数；k表示当前正在计算的状态层的编号(1≤k≤K)；i，j分别表示PE当前计算的元素在三角矩阵中的行/列坐标；L为RNA序列的长度；图中的虚线箭头代表三维矩阵的计算顺序。

[0139] PE阵列采用“按区域分割”和“逐层按列并行处理”的任务划分策略实现对三维动态规划矩阵的细粒度并行计算。如图7a所示，首先对三维矩阵按状态层方向进行纵切，将其垂直分割为s个区域( )，每个区域仍是一个三维矩阵，每个区域都包含K层，每一层包含n列元素；然后按照编号对分割后的区域进行逐一计算。对每一个区域m(1≤m≤s)，采用逐层计算的策略：从第1个区域的第K层出发(图7a区域1的第K层中标出的粗线代表初始计算位置)，当第K层计算完成后再算第K-1层，直到第1个区域的第1层算完，然后再计算第2个区域的第K层，第K-1层，…，第1层；之后再计算第3个区域…，直到最后一个区域s的第1层算完，此时区域s的第1层右上角元素M(1，L，1)的值便是最终的序列比对分值。对每一层的计算而言，每个PE负责计算当前状态层中的一列，PE计算的元素的列号与PE在阵列中的序号P_id对应。主PE的编号P_id＝1，负责计算本区域的第1列元素，从PE计算本区域的第2～n列元素。

[0140] 阵列中的n个PE采用SPMD方式同时计算位于矩阵不同列的n个数据，以达到对计算过程进行加速的效果。所有PE都使用无回溯CYK/inside算法的迭代公式对矩阵元素进行计算，但计算时根据不同的状态类型采用不同的计算顺序。

[0141] ●如果当前状态k是分支状态，则每个PE都从各自列的顶部位置(i＝1)开始计算，按照由上至下的顺序执行(如图7b所示)。当前状态层的计算启动时，主PE计算元素(1，j)，第1从PE计算元素(1，j+1)，…，第n-1从PE计算元素(1，j+n-1)。当第1行的n个元素计算完之后，所有PE同时更新元素的横坐标，计算下一行(i＝2)的n个元素。当一个PE计算的元素坐标i＝j时(即到达各自列的底部位置)，PE计算暂停，进入等待状态。由于三角矩阵每一列元素的个数不等，PE按编号顺序依次进入同步等待状态，直到收到最后一个从PE的计算完成信号。特别地，当n＝1时，此时没有从PE，主PE按照从左至右的顺序对矩阵元素进行逐列计算，在每一列内部按照由上至下的顺序对每一个元素进行计算。

[0142] ●如果当前状态k是非分支状态，则每个PE都从各自列的底部位置开始，按照由下至上的顺序计算(如图7c所示)。当前状态层的计算启动时，所有PE计算的元素都位于三角矩阵的主对角线上(i＝j)：主PE计算元素(i，j)，第1从PE计算元素(i+1，j+1)，…，第n-1从PE计算元素(i+n-1，j+n-1)。位于同一条对角线上的元素的计算量都相等，因此所有PE都可以同时完成当前元素的计算，当前对角线上的n个元素计算完之后，所有PE同时向上移动一个计算位置，主PE计算元素(i-1，j)，第1从PE计算元素(i，j+1)，…，第n-1从PE计算元素(i+n-2，j+n-1)。当PE计算的元素行坐标i＝1时，PE计算暂停，进入等待状态(如果计算结果需要写回DRAM，PE将在等待状态下发出写回请求)。同样，由于三角矩阵每一列元素的个数不等，PE将按编号顺序依次进入同步等待状态。主PE算完元素(1，j)后进入等待状态，此时第1从PE计算元素(1，j+1)，…，第n-1从PE计算元素(n-1，j+n-1)。第1从PE计算完(1，j+1)后进入等待状态，此时第2从PE计算元素(1，j+2)，…，第n-1从PE计算元素(n-2，j+n-1)，第2从PE计算完元素(1，j+2)后也进入等待状态，…，直到第n-1从PE的算完元素(1，j+n-1)。特别地，当n＝1时，此时没有从PE，主PE按照从左至右的顺序对矩阵元素进行逐列计算，在每一列内部按照由下至上的顺序对每一个元素进行计算。由于只有一个PE，此时也就不存在同步等待，当一列元素计算完成后如果需要写回，则直接发出写回请求。

[0143] 当PE控制模块收到最后一个从PE(P_id＝n-1)的当前状态层计算完成信号后(说明此时本区域当前状态层的计算任务已完成)，随后PE控制状态机进入结束条件判断状态。如果j＞L，则PE的计算过程结束，否则在数据缓存区寻找当前编号最小的一个空闲存储区域，回到数据载入步骤，准备对下一个状态层进行计算。当区域s中的最后一个状态层(k＝1)计算完成后，计算过程结束，此时矩阵右上角元素M(1，L，1)的值便是最终的序列比对得分。