生物信息学中一个优化的全局双序列比对算法

第24卷2004年6月

　

计算机应用

C omputer Applications

　

Vol. 24J une , 2004

文章编号:1001-9081(2004) 06Z -0307-02

生物信息学中一个优化的全局双序列比对算法

唐玉荣

(中国农业大学现代精细农业系统集成研究教育部重点实验室, 北京100083)

摘　要:最早的生物信息学中序列比对算法是基于动态规划思想的Needleman 2Wunsch 全局双序列比对算法, 由于其时间和空间复杂度巨大, 不适合实际的生物序列比对。态规划思想的全局双序列比对算法。实验结果表明, , 有效地降低了时间和空间复杂度。

关键词:算法; 双序列比对; ; 中图分类号:TP301. 61　引言

序列比对是生物信息学的核心研究内容之一。在生物学研究中, 为了判断两个序列是否具有足够的相似性, 从而判断两者是否具有同源性, 常常需要进行序列比对。序列比对根据同时进行比对的序列数目分为双序列比对和多序列比对。序列比对也可分为全局比对和局部比对, 全局比对考虑序列的全局相似性, 局部比对考虑序列片段之间的相似性。

最早的序列比对算法是Needleman 2Wunsch 全局双序列比对算法[1], 主要思想是利用动态规划的方法计算两条序列之间的一个最佳比对。该算法在其后的三十年中得到广泛的应用, 成为生物信息处理算法中的一个最基本的算法。随后涌现出了一大批优化的序列相似性比较算法, 它们虽有各自的定义, 但都是基于动态规划思想, 而且其时间和空间复杂度也都很巨大, 仅适用于字符数量不大的序列, 无法在实际中使用。本文提出了OG P 全局双序列比对算法, 该算法基于动态规划思想, 能有效地降低时间和空间复杂度。

过行i 的水平方向, 空位罚分的值取决于插入空格的个数。若两个序列为a =a 1a 2…a n 和b =b 1b 2…b n , 则s i ,j =s (a 1a 2…a n , b 1b 2…b n ) ,s i ,j 是到达序列a 中第i 位字符与序列b 中第j 位字符的比对得分值,s (a i ,b j ) 是字符a i 与b j 的计分值,w x 是在序列a 上空格长度为x 的空位罚分值,w y 是在序列b 上空格长度为y 的空位罚分值。在此注意,s i ,j 是得分矩阵中从三个方向上得到的一条最佳比对路径得分值。当得分矩阵的所有元素s i ,j 被计算出来后, 最佳路径的终点是在最后一行最后一列的位置。从这一点开始根据上面公式在得分矩阵中回溯寻找得到的路径就是一条最优路径

。

图1　动态规划算法描述

2　动态规划算法思想

[2]

生物序列包括DNA 序列和蛋白质序列,DNA 序列由4种碱基组成, 则DNA 序列可看作是由4种字符组成的字符序列; 蛋白质序列由20种残基组成, 则蛋白质序列可看作由20种字符组成的字符序列。序列比对算法实际上就是通过对两个或多个字符串序列插入“-”表示插入或删除, 来获得序列之间的最佳比对结果。例如两条DNA 序列ATG C 和AG TC , 希望通过对每条字符序列插入空格, 得到使两条序列的匹配字符数最大的最佳比对。经过观察, 最佳比对应该为:

ATG -C ; 　A -GTC

3　优化的全局双序列比对算法UHF

Needleman 2Wunsch 算法在计算得分矩阵D 过程中存储所

对于上述问题, 动态规划算法的解决方案基本可描述为:计算得分矩阵; 在得分矩阵中回溯寻找最优比对序列。

S i ,j =max{S i -1,j -1+s (a i , b j ) ,max (S i -x ,j -w x ) ,

x Ε1

max (S i ,j -y -w y ) }

y Ε1

具体如图1所示, 从三个方向可以到达矩阵元素(i , j ) :对角线方向元素、同一行或同一列的元素。在得分矩阵中, 到达位置为i , j 的某一个元素有三种可能的路径:通过位置i -1, j -1的对角方向, 没有空位罚分; 通过列j 的垂直方向, 通

有的元素, 时间和空间复杂度均为O (m ×n ) , 见图2(a ) 。

[3]

Hirschberg 算法能够解决基本动态规划算法空间复杂度太大的问题, 它在得分矩阵的计算过程中不存储所有元素, 只存储得分矩阵的当前行(或列) 和前一行(列) 。该算法的描述过程为:在矩阵D 的中间列(或行) 处将Bs 序列分割成两部分, 从上和下两个方向计算得分矩阵, 以寻找在这一列或行上最优路径所经过的点checkpoint ; 递归计算由checkpoint 点组成的子矩阵; 所有的checkpoint 点组成的路径就是最佳比对路径。此算法的空间复杂度为O (m +n ) , 计算时间却是基本动态规划算法的两倍。Ukkonen [4]是一种广泛使用的基于动态规划的快速序列比对算法, 在得分矩阵D 的基础上计算替换矩阵U , 矩阵元素为U [ab ,d], 其中ab =a -b 对应D 矩阵的对角线, d 对应D 矩阵中的得分值, U [ab ,d ]的值为a , 该值

ab +1和ab -1上值为d -1的元素所决定。由D 矩阵中ab 、

示例见图2(b ) 。此算法的时间复杂度为O (n +d 2) (在此n 是

　　收稿日期:2003-09-10; 修订日期:2003-12-23　　基金项目:北京市科技计划项目(H[1**********]021)

　　作者简介:唐玉荣(1974-) , 女(蒙古族) , 内蒙古准旗人, 博士研究生, 主要研究方向:计算机应用、生物信息学.

　　　　308计算机应用2004年

序列长度,d 是两条序列的得分值) , 但空间复杂度仍然为O (n 2) 。为了同时降低时空复杂度, 之后提出的将Hirschberg 算法和Ukkonen 算法进

行了结合应用, 时间和空间复杂度介于Hirschberg 和Ukkonen 之间。

在得分矩阵D 的计算过程中, 由于第i 行(或列) 的元素值都来源于第i -1行(或列) , 因此F A 算法[6]提出在计算得分矩阵时同时获得相邻两行元素的计算来源关系, 这样不需要回溯就可以获得checkpoint 点。并且为了降低时间复杂度, 在计算过程中增加了checkpoint 的数目, 减少了对矩阵的重复计算。实现过程见图2(d ) 。

从Ukkonen 算法中对U 、ab +列值为d 的元素值由列值为d -1行ab -1的元素所决定, 列计算得到, 见图2(c ) 2and 2C onquer 算法计算U 矩阵的过程中借鉴checkpoint 数目并记录元素来源关系的算法思想, 从而降低时间复杂度

。

Divide 2and 2Conquer 算法

[5]

{Outerloop , iterated until U[|As |-|Bs |,d]=|As |}U[ab,d]=max (U[ab+1,d -insertC ost ],

U[ab, d -mismatchC ost ]+1, U[ab-1,d -deleteC ost ]+1)

{InnerLoop , extends diagonal on a run of matches}

while (As[U[ab,d]+1]=Bs[U[ab,d]-ab +1]) U[ab,d]+=1

4　实验结果

G 000、5000和10000的序列, 2Wunsch 、Divide 2and 2C onquer 、和UHF 算法进行了对比实验, 以验证UHF 。实验是在一台Pentium4,CPU 主频为1. 6GHz , RAM 为640M B , 操作系统为Windows 2000的PC 机上进行。实验结果见表1和表2。表1结果表明, 序列长度为1000时,UHF 算法的时间复杂度明显小于Needleman 2Wunsch 、Divide 2and 2C onquer 和Hirschberg , 但略大于Ukkonen 算法; 序列长度为5000、10000时,UHF 算法的时间复杂度小于其他四种算法。表2结果表明, 在序列长度为1000、5000和10000时, UHF 算法的空间复杂度均明显小于Needleman 2Wunsch 和Ukkonen , 但略大于Divide 2and 2C onquer 和Hirschberg 算法。由于实际的生物序列长度很大, 而且UHF 算法在存储空间的需求增加在计算机可用存储空间的范围之内, 所以UHF 算法具有一定的实用性。

表1　不同序列长度下五种算法的运行时间比较

算法

Needleman 2Wunsch

Hirschberg Divide 2and 2C onquer

Ukkonen UHF

运行时间(s )

10009. 2218. 445. 994. 534. 98

500036. 5473. 0824. 3519. 9818. 73

1000066. 75133. 5040. 4135. 6232. 52

表2　不同序列长度下五种算法的空间需求比较

算法

图2　UHF 算法示例

Needleman 2Wunsch

Hirschberg Divide 2and 2C onquer

Ukkonen UHF

空间需求

100019M B 20K B 22K B 19M B 40K B

500054M B 49K B 52K B 54M B 141K B

10000101M B 87K B 95K B 101M B 265K B

UHF 算法在U 矩阵的计算过程中只存储前一列和当前

列, 寻找多个chekpoint 点, 从而达到节省存储空间的目的; 且

在U 矩阵计算过程中同时记录元素的来源关系, 最佳比对路径的获得不需要回溯。

UHF 算法的实现步骤如下:

1) 计算U 矩阵, 计算过程中只存储相邻两列的元素, 同时存储元素的来源关系;

2) 对U 矩阵的列进行k 等分, 由U 矩阵计算过程中保存的元素来源关系获得k 个checkpoint 。对k 个checkpoint 点进行排序, 得到各个小矩阵的开始点和结束点;

3) 对各个小矩阵递归调用UHF 算法, 直到小矩阵的开始点和结束点之间的行差和列差小于等于分割数k , 此时直接调用Hirschberg 算法计算;

4) 所有的checkpoint 点构成一条最佳比对路径。对U 矩阵进行计算的Ukkonen 算法描述如下:

{U[ab,d]=max a s. t. D[a,b ]=d where ab =a -b

=-infinity if no such a exists}

U[0,0]=max a s. t. As[1.. a ]=Bx[1.. a ]U[ab,d]=-infinity , if |ab |>d

参考文献

[1]　Needleman S , Wunsch C. A general method applicable to the search

for similarities in the amino acid sequences of tw o proteins [J].Journal of M olecular Biology , 1970,48. 443-453.

[2]　M ount DW. Bioinformatics :sequence and genome analysis [M].

USA :C old S pring Harbor Laboratory Press , 2002. 53-72. [3]　Hirschberg D. A linear space algorithm for com puting maximal

comm on subsequences [J].C omm ACM , 1975,18(6) :341-343. [4]　Ukkonen E. On approximate string matching [A].Proceedings Int

C orf F ound C om put Theory , 1983,158. 487-495

[5]　P owell DR , Allison L , Dix TI. A versatile divide and conquer

technique for optimal string alignment [J].Information Processing Letters , 1999,70(3) :127-139.

[6]　李昭, 杨琪, 祝明发. 存储约束条件下的序列联配算法[J].微电

子学与计算机, 2002,19(6) :1-5.