布尔逻辑DNA序列比对系统的软件平台设计

2023年12月7日发(作者：)

布尔逻辑DNA序列比对系统的软件平台设计

陈亚东;冯萍;王进科

【摘 要】针对目前基于动态规划的DNA序列全局比对算法时间复杂度较高,设计了一个DNA序列全局比对系统.该系统用FPGA进行序列的比对,并配备一个软件平台存储数据、发送命令以及发送和接收数据.测试数据表明,该系统的DNA序列比对时间在序列相似度较低情况下,为Needleman的42%;在序列相似度较高的情况下,为Needleman的6%.%Present DNA sequence alignment algorithms

based on dynamic programming have a high time complexity.A global

alignment system is introduced to address this system aligns

sequences in FPGA and store sequences and aligning results in PC

system also send commands to experiment results

show compared to Needleman algorithm, the aligning time decreases by

58 when the sequences have a low similarity, and by 94 when the

sequences have a high similarity.

【期刊名称】《科学技术与工程》

【年(卷),期】2011(011)007

【总页数】7页(P1468-1473,1479)

【关键词】FPGA;DNA序列;软件平台;动态规划

【作 者】陈亚东;冯萍;王进科

【作者单位】西北工业大学计算机学院,西安,710072;西北工业大学计算机学院,西安,710072;西北工业大学计算机学院,西安,710072 【正文语种】中 文

【中图分类】TP311.1

一条DNA序列可以看成A(腺嘌呤)、C(胞嘧啶)、G(鸟嘌呤)、T(胸腺嘧啶)四种碱基字符的线性排列。在生物进化过程中，由于突变，DNA序列会发生改变。如果两条序列是由共同的祖先序列进化而来，就称这两条序列为同源序列。序列比对通过碱基字符的一一配对，使两条序列最大程度地反映它们的共同性。共同性的多少反映了序列与其祖先序列在亲缘关系上的远近程度。序列比对是生物信息处理中最基本、最重要的操作［1］。

Needleman算法［2］是经典的全局动态规化比对算法，Smith-Waterman算法［3］是对 Needleman算法进行改进而形成的局部动态规划比对算法，适用范围很广。BLAST算法［4］是启发式的局部相似性搜索算法，采用搜索种子-扩展种子的策略，可以快速给出数据库中与查询序列有较高相似度的序列集合，但该算法的序列比对精度稍差。BLAT算法［5］与BLAST相似，不同之处在于BLAT对数据库建立索引而非对查询序列建立索引;BLAT将高得分片段对接合在一起作为比对结果返回，而BLAST将每个高得分片段对都返回。国防科大的张阳等用FPGA构成的脉动阵列对Smith-Waterman算法进行了优化［6］，取得了不错的进展。

现场可编程门阵列(FPGA)内部支持程序并行执行，时钟周期在ns级，很适合用来进行序列比对。遗憾的是，FPGA的存储容量一般不大。目前，Altera公司的新一代FPGA产品Arria II GX，其Memory大小不到1 MB。作者所在小组设计的布尔逻辑DNA序列比对系统采用PC内存存储－FPGA比对的策略，弥补FPGA容量不足的问题。这种策略使得PC和FPGA之间的数据传输成为比对系统的关键。USB 2．0最高速度达480 Mb/s，可以满足FPGA对数据的处理速率，作为比对系统的数据传输方式在理论上是可行的。

本文论述了布尔逻辑DNA序列比对系统软件平台的设计和实现。

1 布尔逻辑比对系统组成

布尔逻辑比对系统由硬件(协处理器)和软件(BooleanAlign)两部分组成。BooleanAlign配合协处理器完成整个系统的功能，二者之间通过USB系统互连。协处理器是一个外接USB设备，通过主机上的USB接口与主机相连。系统的组成框图如图1所示。

当开始比对的时候，BooleanAlign首先向协处理器发送ASCII码表示的DNA序列。由于DNA序列往往很长，BooleanAlign与协处理器协商好分段发

图1 布尔逻辑DNA序列比对系统框图

送，每次发送512 B。每发送一段数据，协处理器就处理一段数据，然后返回处理结果。接着，Boolean-Align继续向协处理器发送下一段数据。这样循环，直到序列全部发送完毕。序列发送完毕的同时，比对结果也全部返回到BooleanAlign。最后，Boolean-Align依据这些数据构建两条序列的最佳匹配和系统发生树。比对系统的工作流程图如图2所示。

图2 布尔逻辑比对系统工作流程图

2 BooleanAlign平台设计

BooleanAlign是布尔逻辑比对系统的总控平台。DNA序列的输入、比对的开启、数据的发送与接收、控制命令的发送以及最佳匹配的构建和系统发生树［7］的重建，都由BooleanAlign控制和完成。BooleanAlign的运行环境为 Windows，采用 Microsoft公司的高效编程语言C++加以实现。BooleanAlign的程序调用关系如图3所示。

图3 BooleanAlign程序调用关系图

应用程序层(D2XX API Application)通过调用动态链接库(FTD2XX．DLL)函数实现与 USB设备(FTD2XX．SYS)的通信。FTD2XX．SYS和 FTD2XX．DLL分别是FTDI(Future Technology Devices International Ltd)公司为其USB接口转换芯片FT245BM提供的驱动程序和编程接口函数库。BooleanAlign平台的设计包括人机界面设计、数据通信、比对结果处理和基于距离的系统发生树重建。

2.1 人机界面设计

利用MFC(Microsoft Foundation Class)可以快速构建Windows应用程序框架［8］。VB也可以快速构建应用程序框架，但执行效率不如VC++．Java源代码经编译后变成字节码(．class)文件，程序需要边解释边执行，不适应BooleanAlign程序特性。基于以上情况，BooleanAlign选择MFC作为人机界面的实现工具。BooleanAlign人机界面采用MFC多文档结构，以满足用户同时打开多个窗口查看比对数据的要求。采用“多视窗”技术，创建动态分割窗口，将界面分割成独立的多个功能区域。菜单命令的响应函数，用来实现用户的各种操作。实现过程中的难点在于:视图窗格(pane)中显示的四种碱基字符(A、T、C、G)的字符宽度不相等，导致鼠标移动到相应字符上时，字符位置(posi)显示错误。解决方法是:调用CDC的文本绘图函数ExtTextOut()，将每一个字符的实际显示宽度调整成一致。这样，程序就可以直接利用鼠标的逻辑坐标值来计算字符在碱基序列中的位置了。字符位置计算公式如式(1)。

式(1)中，posi是字符在碱基序列中的位置，x_LPPoint是鼠标的逻辑坐标，widthCharacterRect是每个字符的实际显示宽度。人机界面如图4所示。

图4 BooleanAlign主要界面

2.2 数据通信

2．2．1 数据通信格式

上层应用程序与底层硬件的数据交换必须按照一定的协议进行。BooleanAlign发送给协处理器的数据和协处理器返回的数据必须具有一定的格式，双方才能识别所接收到的数据的意义。

BooleanAlign发送给协处理器的数据单元的格式如图5。

图5 发送数据单元格式

“数据块”是长度为k个字节的DNA序列，即k个碱基。k在一般情况下都取512，只有在最后一个数据单元k的值可能会小于512。由于插入和缺失是一对相对的概念，为了区别这两个概念是相对于哪条DNA序列而言，有必要对发送的数据加个标识。“标识”给数据块加上了一个标志，用来说明该数据块来自两条序列中的哪一条，占一个字节，取值范围是三个 ASCII字符“E”、“N”、“D”，如表1 所示。当标识是“D”时，表示参与比对的两条序列已全部发送到协处理器，即数据发送完毕。

表1 标识说明标 识 说 明序列N数据块属于B序列D E 数据块属于A数据块属于B序列，是最后一个数据块

协处理器返回给BooleanAlign的数据单元的格式如图6所示。

图6 接收数据单元格式

cnt代表发生变异的碱基位点(碱基字符在DNA序列中的位置)个数，占一个字节，最大可表示256个变异位点。每个变异位点信息用3个字节描述，position描述变异发生在DNA序列的什么位置，占两个字节;type描述变异是四种类型中的哪一种，占一个字节，分别用 00000001B表示颠换，00000010B表示插入，00000100B表示缺失，00001000B表示转换。2．2．2 数据通信过程

(1)过程描述

在BooleanAlign与协处理器通信过程中，数据要分段多次发送接收，每次通信过程包括一次发送一次接收。要保证整个通信过程不出差错，需要一种机制来控制通信过程的进行。BooleanAlign采用MFC的线程机制，通过控制线程的阻塞与唤醒，实现数据通信过程中发送与接收的同步。

如图7(a)所示，首次发送数据时，send线程将firstWr置为1，然后进入

WriteDev状态。在 WriteDev状态，send线程截取DNA序列的一段并发送到协处理器。发送完毕且发送正确后，send线程等待50 ms(等待协处理器将本次序列比对完毕)，然后将信号 eventRd置为有效，eventWr置为无效，send线程进入Block(阻塞)状态;否则，writeDevError有效，endRd被置为1，send线程进入 writeDev Error状态，接着立即终止。在Block状态receive线程将接收并存储协处理器返回的数据，这将在下面叙述。当eventWr被receive线程置为有效时，send线程又进入WriteDev状态。此时，如果还有DNA序列要发送，send线程继续截取DNA序列的一段并发送到协处理器;如果DNA序列已经全部发送完毕，那么send线程进入终止态。

图7 数据发送、接收状态

如图7(b)所示，当eventRd被send线程设置为有效后，receive进程被启动并进入ReadDev状态。在ReadDev状态，receive线程接收并存储协处理器返回的数据。这里有三种可能。第一，如果数据接收正确，则eventWr置为有效，eventRd置为无效，receive线程进入Block(阻塞)状态。第二，如果这是最后一段数据(lastWr=1)，则 receive线程进入ReadDev End状态。第三，如果接收数据出错，则readDevError有效，endWr置为 1，receive线程进入ReadDev

Error状态，接着立即终止。在Block状态send线程发送下一段 DNA序列。当

eventRd被send线程置为有效时，receive线程又进入ReadDev状态。

(2)通信的实现细节

协处理器是一个USB2．0设备，FTDI公司提供了访问USB2．0设备的编程接口函数(如表2)，使应用程序可以直接调用接口函数实现与协处理器的通信。Windows操作系统本身带有USB设备的驱动程序，不过这是个通用的USB驱动程序，不能实现针对协处理器的最大传输率。FTD2XX．SYS可以针对协处理器的特点，对USB通信过程进行简化，从而实现通信速率的提高。如图3所示，在上层应用程序，WriteDev()函数调用FTD2XX．DLL里的FT_Write()库函数，实现向协处理器发送数据;ReadDev()函数调用 FTD2XX．DLL里的 FT_Read()库函数，实现接收数据。BooleanAlign使用隐式连接调用动态连接库文件中的函数。程序在运行时动态调用所需要的代码和数据，节省了内存占用。FTD2XX．SYS向上提供上层软件访问协处理器的接口，向下直接操纵协处理器，允许动态链接库函数通过驱动程序直接访问协处理器。FTD2XX．SYS采用批量传输方式，有一个接口，两个传输管道，每个管道支持最大64字节的数据包。

表2 FT_Read与FT_Write函数参数说明表函数 参 数 说 明ftHandle FT_Read设备句柄lpBuffer 指向读缓冲区的指针，用来存储要从协处理器读取的数据dwBytesToRead 32位整数，表示要从协处理器读取的数据字节个数lpdwBytesReturned指向32位整型数据的指针，表示实际上从协处理器读取到的字节个数ftHandle FT_Write设备句柄lpBuffer 指向写缓冲区的指针dwBytesToWrite 32位整数，表示写缓冲区中数据有多少字节lpdwBytesWritten 指向32位整型数据的指针，表示实际写到协处理器的字节个数

2.3 比对结果的处理

主机接收到的比对结果包含position和type，分别表示突变发生的位点和突变类型。应用程序根据这两种数据决定如何处理该位点的这对碱基。如果是转换和颠换，将该位点碱基以黑色显示;如果是插入，A序列上该位点的碱基与B序列上的“_”配对，B序列该位点以后的碱基依次后移一位;如果是缺失，A序列以“_”与B序列该位点的碱基配对，A序列该位点以后的碱基依次后移一位。假设，有序列S:TACAGTTGGATC，序列T::TTTGGA，其比对结果的集合R={＜2，插入＞，＜3，插入＞，＜4，插入 ＞，＜5，插入 ＞，＜11，插入 ＞，＜12，插入＞}。经应用程序处理以后，序列呈现为:

T A C A G T T G G A T C

T _ _ _ _ T T G G A _ _

2.4 基于距离的系统发生树重建

非加权组平均法(UPGMA)［9］是一种简单的基于距离的系统发生树够建方法。Kimura模型［10］是常用的一种距离计算模型，其距离计算公式为:

d1表示比对排列后转换的位点数目与序列总长度的比值;d2表示比对排列后颠换的位点数目与序列总长度的比值，K表示两条 DNA序列的距离［11］。构建基于距离的系统发生树，首先要建立一个距离矩阵，用来选择一对距离最小的DNA序列。程序中用一个二维数组(GDVector)存储距离矩阵，矩阵元素是一个结构体，名为Element，如图8所示。

图8 距离矩阵元素

Element存储了行序列与列序列的距离、行序列和列序列各自的名称以及本节点的名称。例如，图9(a)是a、b、c、d四条DNA序列的距离矩阵。a列b行处的矩阵元素(命名为E)，其结构体内容为:E．GD=0．000 161，E．up=“a”，E．down=“b”，E．name=“ab”。

图9 距离矩阵删减过程

聚合节点是BooleanAlign绘制系统发生树的基础，它代表了一个如图10所示的形状。每个聚合节点有四个顶点，图10示意了聚合节点u和v的8个顶点。BooleanAlign借助距离矩阵找到构建系统发生树所需的全部聚合节点，然后用MFC画图函数LineTo()绘制系统发生树。

图10 聚合节点代表的形状 以来自GeneBank加德纳菌属的四条序列a(GI号:297532249)、b(GI号:297532250)、c(GI号:297532251)、d(GI号:297532252)为例，构建系统发生树的步骤为:

Step 1:建立初始化距离矩阵(GDVector)。以a、b为例，经协处理器比对后得到的转换位点个数为12，颠换位点个数为26，序列总长度为195。根据式(2)，得K=0．098 008，将K填入距离矩阵的a列 b行。依次计算序列 ac、ad、bc、bd、cd 之间的距离，得到初始距离矩阵如图9(a)所示。

Step 2:如果 GDVector．size()=0，转到 Step 6;否则，找出距离矩阵中GD最小的元素minElement。在本例中为元素ac，它的GD为0．049 519。

Step 3:计算聚合节点坐标。a、c首先聚合出一个节点u。根据BooleanAlign给定的绘图初始点P(x，y)和 u 的分支长度 du=minElement．GD/2，计算 u 的四个顶点 point1(x，y)、point2(x，y)、point3(x，y)、point4(x，y)坐标。

Step 4:删减距离矩阵。已经参加过聚合的序列不再参加下次聚合，要从距离矩阵中删掉。在本例中，删掉序列a、c。然后，计算u与其它序列之间的距离。u与序列m的距离dum=(dam+dcm)/2，m=b，d。删减后的距离矩阵如图9(b)所示。

Step 5:重复 Step2。

Step 6:根据聚合节点，绘制系统发生树。

按照上述的方法，来自狼狗的五条DNA序列的系统发生树如图11所示。

图11 五条狼狗DNA序列的系统发生树

分支 a，b，c，d，e，f的长度如表3 所示。

表3 五条狼狗DNA序列的系统发生树分支长度分支a 分支b 分支c 分支d 分支e 分支f 0．000 08 0．000 578 0．000 498 0．000 121 0．000 8410．000

384 3 测试结果和分析

在Windows xp操作系统，1．61 GHz处理器主频，512 MB内存条件下，进行了两项测试。一项测试BooleanAlign在序列相似度较高情况下的性能，另一项测试BooleanAlign在序列相似度较低情况下的性能。每项测试包括五组测试数据，序列长度从300到16 195不等，每组测试都是双序全局列比对。测试数据如表4、表5所示，给出了在不同长度情况下，BooleanAlign和Needleman各自的性能。

表4 BooleanAlign与Needleman的性能对比1注:maxLen表示序列最大长度Time/ms Similarity/%SourceBoolean-Align needle Boolean-Align needle

maxLen/(bp)470 5 900 98．2 98．1 300大肠杆菌 563 6 820 99．8 99．8 1

533小鼠 593 12 832 99．7 99．7 3 932小鼠 609 35 900 99．9 99．9 9 054狼狗果蝇5 922 159 900 94．2 93．4 16 195

表5 BooleanAlign与Needleman的性能对比2Time/ms Similarity/%Source

Boolean-Align needle Boolean-Align needle maxLen/(bp)437 7 900 49．2

37．0 302大肠杆菌 2 390 8 902 48．0 41．8 1 335小鼠 7 172 9 910 46．9

40．7 4 026小鼠 11 422 11 911 47．7 40．9 5 688狼狗果蝇12 250 19 906

47．3 41．5 6 395

令Tb表示BooleanAlign比对耗费时间，Tn表示Needleman算法比对耗费时间。在第一项高相似度测试中，五组数据的Tb/Tn的均值为0．06。在第二项低相似度测试中，五组数据的Tb/Tn的均值为0．42。

4 结论

测试结果显示，在两条序列相似度很高的情况下，BooleanAlign的比对耗费时间是Needleman算法的6%;在两条序列相似度较低情况下，BooleanAlign的比对耗费时间也仅是Needleman算法的42%。可以得出，BooleanAlign结合FPGA的速度优势和PC内存容量的优势，提高了双序列全局比对的速率。 参考文献

【相关文献】

1 孙 啸，陆祖宏，谢建明．生物信息学基础．北京:清华大学出版社，2005:81

2 Needlman S B，Wunsch C D．A genera method application to the search for similarities

in the amino acid sequence of two proteins．Journal of Molecular Biology，1970，48(3):443—453

3 Smith T F，Waterman M S．Identification of common molecular sequence．Journall of

Molecular Biology，1981;147:195—197

4 Altschul SF，Gish W，Miller W，et al．Basic local alignment search tool．Journall of

Molecular Biology，1990;215(3):403—410

5 Kent W J．BLAT—the BLAST-like alignment tool．Genome Res，2002;10:656—664

6 张 阳，窦 勇，夏 飞．生物信息学双序列比对算法加速器设计与实现．计算机科学与探索，2008;2(5):519—528

7 郭 静，王 超，陈 崚．基于遗传算法的系统发生树构建方法．计算机工程与应用，2009;45(16):72—76

8 ［美］Shephend G，Wingo S．深入解析MFC．赵剑云，卿 瑾，译．北京:中国电力出版社，2003:5—17

9 Sokal R R，Michener CD．A statistical method for evaluating systematic

relationships．University of Kansas Scientific Bulletin，1958;28:1409—1438

10 Kimura M．A simple method for estimating evolutionary rates of base substitutions

through comparative studies of nucleotide sequences．Journal of Molecular Evolution，1980;16(2):111—120

11 许忠能．生物信息学．北京:清华大学出版社，2008:235—240