CN101887493B

CN101887493B - 一种基于布尔逻辑的基因比对处理方法

Info

Publication number: CN101887493B
Application number: CN 201010246011
Authority: CN
Inventors: 康继昌; 徐烽涛
Original assignee: HUZHOU RUIWANSI INFORMATION TECHNOLOGY Co Ltd
Current assignee: Ruiwansi (Hangzhou) Information Technology Co.,Ltd.
Priority date: 2010-08-05
Filing date: 2010-08-05
Publication date: 2013-02-20
Anticipated expiration: 2030-08-05
Also published as: CN101887493A

Abstract

本发明提供了一种基于布尔逻辑基因比对处理方法，该方法包括两条基因序列，其比对步骤为：将需要比对的基因序列输入到FPGA器件中后用单比特表示碱基，用布尔逻辑表示基因的突变规则，用FPGA器件实现比对计算，用位置计数器和突变计数器记录突变的类型和数量，结束时将记录发送给上层软件，上层软件计算处理后将结果显示在屏幕上。利用本发明提供的方法，比对逻辑与“读碱基”同时工作，比对逻辑无需占用额外的时间，比对计算的时间与序列的长度成线性关系，显著提高了比对计算速度。

Description

一种基于布尔逻辑的基因比对处理方法

技术领域

本发明涉及基因比对方法，具体讲涉及一种基于布尔逻辑的基因比对处理方法。

背景技术

1条DNA序列是由A，T，C，G等4种碱基组成。简单地说，作亲子鉴定，就要对“亲”、“子”2条DNA序列进行比对(alignment)。比对时，从2条DNA序列各取第1个碱基开始进行比对，然后各取第2个碱基进行比对，以此类推，从头至尾全部逐拍进行比对。当然，对亲子鉴定，绝大多数比对都是相同的(A-A，T-T，C-C，G-G)。如果有不同(A-T，A-C，A-G，...)，就说明有突变。双序列比对(Pairwise alignment)通常有3种算法：动态规划法(Dynamicprogramming)、点矩阵法(dot matrix method)和词法(Word method)。

现有的基因比对计算，大都采用普通的计算机或集群计算机系统，按照动态规划算法编程，然后通过程序进行计算。其计算时间随基因长度非线性增长，计算速度较慢。

动态规划法主要有3种算法：Smith-Waterman算法，Needleman-Wunsch算法和Gotoh算法。其中，Smith-Waterman算法是Temple Smith和Michael Waterman在1981年提出的，得到广泛应用。

近年国内外对于基于硬件的序列比对的研究主要集中在Smith-Waterman(简称S-W)算法的实现技术。与本发明最接近的成果，如国防科学技术大学计算机学院张阳等在“计算机科学与探索”2008年05期上发表的“生物信息学双序列比对算法加速器设计与实现”，该文简介：“在FPGA上实现大规模脉动式阵列对双序列比对算法进行加速能够大幅度提高比对的效率。以Smith-Waterman算法为例进行了实现验证，结果表明本设计在性能上优于传统设计。与Pentium42.60GHz通用微处理器计算机相比，使用加速器对长度为65536的序列进行比对可获得1555倍的加速比”。可见，其特点主要是用大规模脉动式阵列实现了加速。

[1]Zhang Yang etc.，pair-wise genetic alignment has been implemented by large-scalesystolic array on a FPGA chip，Computer Science and Exploration，2008，No.5，

发明内容

本发明的目的是提供一种基于布尔逻辑的基因比对处理方法，该方法包括两条基因序列，其改进之处在于：包括如下步骤：

将需要比对的基因序列输入到FPGA器件中后用单比特表示碱基；用布尔逻辑表示基因的突变规则，用包括门电路和寄存器的FPGA器件实现比对计算；用位置计数器和突变计数器记录突变的类型和数量，结束时将记录发送给上层软件；上层软件计算处理后将结果显示在屏幕上。

本发明提供的优选基因比对处理方法中，所说单比特表示碱基为4位表示，碱基A、T、C和G分别用0001、0010、0100和1000表示，4位中只有1位是“1”，称为单比特表示法。

本发明提供的再一优选基因比对处理方法中，所说FPGA器件通过2条3级流水线来实现比对计算，所说2条流水线指正常流水线和突变流水线，所说3级流水线指用以确定突变类型的三个寄存器比对过程。

本发明提供的再一优选基因比对处理方法中，所说用布尔逻辑表示基因的突变规则包含步骤：将输入的基因序列进行逻辑与运算；若与运算后输出为1，则说明碱基序列匹配，进行下一碱基比对运算；若与运算后输出结果为0，则说明碱基序列不匹配，比对流入突变流水线，通过3级寄存器流水线判断基因突变类型。

本发明提供的再一优选基因比对处理方法中，所说流水线中均包含两个用于存放碱基的寄存器。

本发明提供的再一优选基因比对处理方法中，所说位置计数器有1个，用来保存突变基因的位置和数量；突变计数器有4个，分别为转换计数器、颠换计数器、插入计数器和缺失计数器，用来记录4种突变类型的位置和数量。

本发明提供的再一优选基因比对处理方法中，所说的4种突变类型为插入I、删除D、转换TRS和颠换TRV。

本发明中的DNA序列比对过程中，碱基可能的互配方式只有4种情况：(1)匹配(即：A-A，T-T，C-C，G-G)。(2)插入I(Insertion)在序列X的某个位置插入了一个碱基。(3)缺失D(Deletion)在序列X的某个位置缺失了一个碱基。(4)错配(即：A-C，C-A，A-T，T-A，A-G，G-A，C-T，T-C，C-G，G-C，T-G，G-T)；错配又分为2种情况：A-G，C-T，G-A，T-C称为转换(TRanSition)；其余，A-C，A-T，G-C，G-T，C-A，T-A，C-G，T-G称为颠换(TRansVersion)。前者简称TRS，后者简称TRV。

匹配表示：比对的序列位点没有发生突变。不匹配的都表示有了突变。由上可见，有I、D、TRS、TRV等4种突变。

本发明中单比特化表示方法为：用一组2进制数表示基因的碱基，其中，只有一位是‘1’。基因碱基的这种表示方法在本发明中称为单比特化表示法。有A，T，C，G等4种基因碱基，在计算机中用ASCII码表示，每个碱基用1个字节，8位，来表示。而本发明提出用4位来表示：A用“0001”表示，T用“0010”表示，C用“0100”表示，G用“1000”表示。4位中只有1位是‘1’，故称为单比特化表示法。下面可以看到，它对简化比对逻辑有重要作用。

单比特化可在FPGA中实现，下面是用VHDL硬件描述语言完成的一个设计：

Process(RST，clk)

begin

if RST＝′0′then

SEQ_X＜＝″0000″；

elsif clk′event and clk＝′1′then

ifSTART＝′1′then

if SEQ_X＝A then

QX＜＝″0001″；

elsif SEQ_X＝T then

QX＜＝″0010″；

elsif SEQ_X＝C then

QX＜＝″0100″；

elsif SEQ_X＝G then

QX＜＝″1000″；

end if；

end process；

其中，START是启动信号；SEQ_X是以ASCII码表示的碱基变量；QX是以“单比特化”表示的碱基变量。每个clk(时钟周期)可完成1个碱基的转换。

本发明提出了基于硬件比对算法的基因比对处理方法。虽然也是基于FPGA实现，但采用的比对算法不同于传统的比对算法。本发明的硬件比对算法，是直接采用布尔逻辑精确地表示各种基因突变规则，来识别各种基因突变，与现有的基因比对算法有本质区别。而且在本发明中，硬件电路不是单纯用于比对，也用于加快比对运算速度。

具体实施方式

本发明的具体结构是一块PCB板，除了有一块ALTERA公司的FPGA器件之外，还有电源，晶振，配置器件，USB接口，安装在一个机箱内。在机箱面板上有开关、RST按钮和指示灯等。USB接口用于和微机(或PC机)相连。在微机上有人机界面，也就是用户的软件操作平台。用户通过“动态联编库DLL(dynamic linking library)”提供的函数来实现本“基因比对处理方法”。

本发明提供的处理方法的比对过程为：第一步，将要比对的基因序列输入FPGA器件；第二步，FPGA器件收到基因序列之后，用单比特表示碱基；用布尔逻辑表示基因的突变规则，用包括门电路和寄存器的FPGA器件实现比对计算；第三步，用位置计数器和突变计数器记录突变的类型和数量，结束时将记录发送给上层软件；上层软件计算处理后将结果显示在屏幕上。另外，还可显示彩色的基因图谱等信息。第二步比对算法具体描述为：

1)基因比对算法逻辑设计

2条经过“单比特化”转换后的基因序列QX和QY，分别保存在2个FIFO(先进先出)存储器中。当比对时，每个时钟周期(我们的实验环境时钟频率为100MHz，时钟周期为10ns)，从2个FIFO存储器中各读出一个碱基，又逐个时钟地分别而同时地流入2条n级流水线寄存器中。大多数的时间里，2个碱基匹配，1个时钟即可完成识别。如不匹配就说明基因有了突变，需要上述的2条流水线寄存器来处理。按照基因突变规则设计的布尔逻辑可以识别基因突变的类型。有了流水线，识别基因突变类型的操作不会中断读入碱基的工作。因此，基因比对和读入碱基同时在工作，一旦序列从头到尾读完碱基后，比对也就同时完成了。计算时间与序列长度成线性关系。

(a)“匹配”的比对算法

设有2条基因序列依次为SEQ_X和SEQ_Y，首先要转换为单比特化表示，分别为QX和QY，由4位表示。例如，QX的4位依次为QX(0)，QX(1)，QX(2)，QX(3)。4位中，只有1位是1。如果QX(0)为1，QX就表示A；如果QX(1)为1，QX就表示T；如果QX(2)为1，QX就表示C；如果QX(3)为1，QX就表示G。比对“匹配”逻辑可以用1个“与门”来实现。“匹配”的英文词为(Match)，如果QX和QY都是A，可表示如下：

M0-AA＝QX(0)AND QY(0)；

其实，同是T，同是C，同是G都可表示匹配。故有：

M0-TT＝QX(1)AND QY(1)；

M0-CC＝QX(2)AND QY(2)；

M0-GG＝QX(3)AND QY(3)；

因此，QX和QY完整的比对匹配算法是

M0＝M0-AA OR M0-TT OR M0-CC OR M0-GG；

以3级流水线为例，流水线可表示为：

REG_X1、REG_X2、REG_X3，以及REG_Y1、REG_Y2、REG_Y3，

同理，REG_X1andREG_Y1完整的比对匹配算法是

M1＝M1-AA OR M1-TT OR M1-CC OR M1-GG；

同理，REG_X2 and REG_Y2完整的比对匹配算法是

M2＝M2-AA OR M2-TT OR M2-CC OR M2-GG；

同理，REG_X3and REG_Y3完整的比对匹配算法是

M3＝M3-AA OR M3-TT OR M3-CC OR M3-GG；

(b)“突变”的比对算法

①插入(Insertion)如下表：

2条基因序列基本相同，只是上面一条QX的第4位点插入了1个G，后面的各位依次后移。2条基因序列逐位比对时，前3位相同，没有“突变”。第4位开始，其QX和QY不同，就是遇到了“突变”：

QX：G A T G C A C G T

| | | / / / / /

QY：G A T C A C G T

可见，QX的第4位的G与QY的第4位C不同，但是，其第5位C与QY的第4位的C却相同。而且，以后的各位也如此。这也是识别“插入”必需的条件。可是，QX和QY经常在变化，故需要流水寄存器来帮助。在REG_X1、REG_X2、REG_X3，以及REG_Y1、REG_Y2、REG_Y3中，保留着不同时刻的碱基状态。有了这些碱基状态，才能确定是遇到了什么类型的“突变”。

以“插入”为例，表示如下：

时刻	1	2	3	4	5	6	7	8	9	10	11	12
													QX	G	A	T	G	C	A	C	G	T
QY	G	A	T	C	A	C	G	T
													M0	1	1	1	0

REG_X1

G

A

T

G

C

A

C

G

T

REG_Y1

G

A

T

C

A

C

G

T

M1

1

0

REG_X2

G

A

T

G

C

A

C

G

T

REG_Y2

G

A

T

C

A

C

G

T

M2

1

0

REG_X3

G

A

T

G

C

A

C

G

T

REG_Y3

G

A

T

C

A

C

G

T

M3

1

0

根据上述“匹配”的比对算法M0＝M0-AA OR M0-TT OR M0-CC OR M0-GG，其不同时刻的逻辑运算值可由QX与QY的值生成M0。如表中第4行，在时刻1，2，3，由于QX与QY相同，M0＝‘1’。在时刻4，由于QX与QY不同，M0＝‘0’。(M0＝‘0’)就成为突变的时间标志。

为了确定这个突变是“插入”，还要看QX在时刻4与QY在时刻5是否相同，这要等下一拍才行。到了下一拍，QX在时刻4的C已经移到REG_Y1。而且，(M0＝‘0’)也移了一拍，突变的时间标志变成了(M1＝‘0’)。故确定这个突变是“插入”的条件有2个。(M1＝‘0’)上面已介绍；QX和REG_Y1都是C的条件是：QX(2)and REG_Y1(2)。同理，还应考虑都是A，T和G的情况。最后，在(M1＝0)的时刻，判别“插入”的比对算法为：

M01＜＝(QX(0)AND REG_Y1(0))OR(QX(1)AND REG_Y1(1))OR(QX(2)ANDREG_Y1(2))OR(QX(3)AND REG_Y1(3))；

以此类推，逐拍还可生成突变的时间标志(M2＝0)和突变的时间标志(M3＝0)。

根据插入的特点是Y序列的碱基和X序列的后1位相同，故其比对算法还可以有：在(M2＝0)的时刻，用：

M12＜＝(REG_X1(0)AND REG_Y2(0))OR(REG_X1(1)AND REG_Y2(1))OR(REG_X1(2)AND REG_Y2(2))OR(REG_X1(3)AND REG_Y2(3))；

在(M3＝0)的时刻，用：

M23＜＝(REG_X2(0)AND REG_Y3(0))OR(REG_X2(1)AND REG_Y3(1))OR(REG_X2(2)AND REG_Y3(2))OR(REG_X2(3)AND REG_Y3(3))；

以上可见，有3种可选的方案。由于出现突变的位置具有随机性，在具体实现中可灵活性地使用这3种方案。流水线级数可变，供复杂情况的应用。关于复杂情况的处理，是属于进一步优化的问题，尚待以后研究解决。

②缺失(Deletion)如下表：

2条基因序列基本相同，只是原来QX的第4位C缺失了，后面的依次前移1位。

QX：G A T A C G T

| | | \ \ \ \

QY：G A T C A C G T

根据缺失的特点是Y序列的碱基和X序列的前1位相同，故比对算法为：

在(M1＝0)的时刻，

M10＜＝(REG_X1(0)AND QY(0))OR(REG_X1(1)AND QY(1))OR(REG_X1(2)ANDQY(2))OR(REG_X1(3)AND QY(3))；

在(M2＝0)的时刻，

M21＜＝(REG_X2(0)AND REG_Y1(0))OR(REG_X2(1)AND REG_Y1(1))OR(REG_X2(2)AND REG_Y1(2))OR(REG_X2(3)AND REG_Y1(3))；

在(M3＝0)的时刻，

M32＜＝(REG_X3(0)AND REG_Y2(0))OR(REG_X3(1)AND REG_Y2(1))OR(REG_X3(2)AND REG_Y2(2))OR(REG_X3(3)AND REG_Y2(3))；

值得注意的是，虽然经过2-3级流水寄存器才能完成了比对，但这并不干扰“读碱基”的工作。实际上，只有遇到基因有突变时，流水寄存器才工作，不会影响总的比对时间。

③转换(TRS)如下表：

两条序列基本相同，只有第4位不同，C转换成T了。这与上述的I或D一样。故只判第4位还不能识别是什么类型的突变。要等下一拍，如果没有识别到I或D，才能确定其为“转换”。如上所说，还是在(M1＝0)的时刻来识别：

判C和T均为‘1’的算法是，

M-CT＜＝REG_X1(2)AND REG_Y1(1)；

类似地，还有TC，AG，GA等TRS的3种情况：

QX：G A T T A C G T C G A T A A C G T C G A T G A C G T C

QY：G A T C A C G T C G A T G A C G T C G A T A A C G T C

它们对应的算法是，

M-TC＜＝REG_X1(1)AND REG_Y1(2)；

M-AG＜＝REG_X1(0)AND REG_Y1(3)；

M-GA＜＝REG_X1(3)AND REG_Y1(0)；

4项组合起来，

M-TRS＜＝M-CT OR M-TC OR M-AG OR M-GA。

④颠换(TRV)

颠换与转换的不同之处就在于组合不同。其8项的组合(即A-C，A-T，G-C，G-T，C-A，T-A，C-G，T-G)可按同样的方法设计，故从简。

3)硬件实现

以上逻辑有的用寄存器表示，有的用门电路表示，都可用VHDL语言描述，便于用FPGA器件实现。用布尔逻辑精确表示各种序列比对规则，既包括匹配也包括突变，总的执行过程也是一种流水形式。使得1对基因序列从头到尾依次经过该比对逻辑，即可完成比对。在比对过程中，有一个位置计数器，每比对1对碱基时，计数器加1，用于给出位置信息。还有4个计数器，分别记录突变类型的信息。当发生突变时，信息都要记录下来。比对结束时，记录发回上层软件。上层软件进行适当的处理后，把处理结果在屏幕上显示，以满足用户的需要。

技术人员从本发明提供的实施例可以，确认匹配的基因对1个时钟周期即可给出结果。突变的基因对的“转换(TRS)”和“颠换(TRV)”需要3个时钟周期给出结果。上述事实说明，由于比对逻辑不影响“读碱基”的过程，所以也就不影响总的比对速度。“插入”和“缺失”由于要调整2条基因恢复对齐，较快的序列要停止“读碱基”1个时钟周期。实际上，突变的基因对数量是很少的，对总体影响不大。比对时间主要是“读碱基”的时间，用于比对额外的时间只占总时间很少的一部分。人类基因长度若按10亿计，“读碱基”时间需10亿拍。每拍10ns，总计比对时间为：100亿ns.＝1000万微秒＝1万毫秒＝10秒。就是说，人类基因比对的基本时间是10秒。用于比对额外的时间决定于“突变”的数量。即使有10％的“突变”，实际比对时间也只有11秒。

本发明提供的方法中全部比对电路占用FPGA的硬件资源较少，只占用Cyclone器件总资源的5％，故一块FPGA器件足够支持多对基因序列的并行处理。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案，技术人员阅读本发明提供的技术方案后可以参照本发明实施例进行种种变更或修改，但这些变更或修改均在申请待批的权利要求保护范围之内。

Claims

1.一种基于布尔逻辑基因比对处理方法，包括两条基因序列，其特征在于：比对步骤如下：

A、将需要比对的基因序列输入到FPGA器件中；

B、FPGA器件收到基因序列之后依次执行如下步骤：

B1、用单比特表示碱基；

B2、用布尔逻辑表示基因的突变规则，用包括门电路和寄存器的FPGA器件实现比对计算；

C、位置计数器和突变计数器记录突变的类型和数量，结束时将记录发送给上层软件；

D、上层软件计算处理后将结果显示在屏幕上；

所说单比特表示碱基为4位表示，碱基A、T、C和G分别用0001、0010、0100和1000表示；

所说用布尔逻辑表示基因的突变规则包含如下步骤：

B21、将输入的基因序列进行逻辑与运算；

B22、若与运算后输出为1，则说明碱基序列匹配，进行下一碱基比对运算；

B23、若与运算后输出结果为0，则说明碱基序列不匹配，比对流入突变流水线，通过3级寄存器流水线判断基因突变类型。

2.按照权利要求1所述的基因比对处理方法，其特征在于：所说FPGA器件通过2条3级流水线来实现比对计算，所说2条流水线指正常流水线和突变流水线，所说3级流水线指用以确定突变类型的三个寄存器比对过程。

3.按照权利要求1所述的基因比对处理方法，其特征在于：所说流水线中均包含两个用于存放碱基的寄存器。

4.按照权利要求1所述的基因比对处理方法，其特征在于：所说位置计数器有1个，用来保存突变基因的位置；突变计数器有4个，分别为转换计数器、颠换计数器、插入计数器和缺失计数器，用来记录4种突变类型的数量。

5.按照权利要求4所述的基因比对处理方法，其特征在于：所说的4种突变类型为插入I、删除D、转换TRS和颠换TRV。