蛋白质序列比对算法

蛋白质序列比对算法（精选六篇）

蛋白质序列比对算法 篇1

蛋白质是一种由20种L-型α-氨基酸通过α-碳原子上的取代基间形成的酰胺键连成的具有特定空间构象和生物功能的肽链构成的生物大分子。氨基酸是带有氨基的有机酸, 它由一个氨基、一个羧基、一个氢原子和一个氨基团组成。在蛋白质合成时, 一个氨基酸的α-羧基和另一个氨基酸的α-氨基连接起来形成肽键。

蛋白质一级结构是指多肽链中氨基酸的排列顺序以及二硫键的位置, 靠共价键维持多肽链的连接, 而不涉及空间排列, 是一个无空间概念的一维结构。测定蛋白质的一级结构具有重要意义, 蛋白质中氨基酸的序列分析时判断蛋白质三维结构的前提。许多遗传性疾病也都是由于蛋白质中个别氨基酸的突变造成的, 所以氨基酸序列分析对疾病的诊断和治疗有重要意义。

通过比较蛋白质的一级结构, 也就是将新测定的蛋白质一级结构与已知的蛋白质氨基酸序列进行比较, 可以预测蛋白质的功能, 并对分析研究这些蛋白质以及它们存在的生物体之间的进化关系有非常重要的关系。

2. 蛋白质双序列比对方法

2.1 序列比对实例

先观察两段蛋白质序列, 上为待检序列, 下位目标序列:

如果两序列中的两个氨基酸残基比对成功则计为1分, 则上面比对为6分。接着向右移动待检序列将移动后的序列与目标序列比对, 直至移动到最后一个位置。然后同理移动目标序列, 直至所有情况都比对完。

对于两个长度为n的序列, 按如上方法共有

这是一个指数级复杂度的计算问题。

我们再仔细观察以上比对, 发现如果在目标序列的第7位插入一空位, 那么则除了一位不匹配以外, 其他完全匹配, 其计分为9。

空位的引入有其生物学的依据, 它意味着两条序列有残基的插入或缺失。当然如果引入的空位的数目不加以限制, 即使联配得分再高那也可能没有生物学意义。因此必须对空位进行罚分, 罚分包括起始空位罚分和延伸空位罚分, 通常起始空位罚分较大, 而延伸空位罚分较小。通常采用线形或仿射罚分函数, 即r (g) =-gd或r (g) =-d- (g-1) e, 其中d为初始罚分, e为延伸罚分, g为引入空位数目。

2.2 计分矩阵

2.2.1 PAM计分矩阵

当双序列比对时, 我们将匹配的残基对计为1, 不匹配的计为0, 这种计分方式我们把它称为一致计分矩阵。这种计分矩阵非常简单, 但是它忽略了蛋白质的进化和结构信息。

Dayhoff等人对进化过程中氨基酸彼此间的替换做了详细的频率研究, 于是产生了一个在短的进化过程中氨基酸彼此间替换的相对频率表, 称为单点可接受突变 (point accepted mutation, PAM) 计分矩阵。它取1个蛋白质序列中的氨基酸突变1%作为进化距离的单位, 称之为一次PAM。Dayhoff等人手工比较了当时数目有限的同源蛋白质序列, 并用统计方法得到对应1次PAM矩阵的数据。

将1次PAM矩阵自乘多次, 可以得到进化距离较远的矩阵, 低号的PAM矩阵可用于进化进化距离较近的亲缘关系, 而亲缘关系较远的可以用更高号码的矩阵。常用的PAM250计分矩阵, 相当于两条序列仍然有20%的相同残基。

2.2.2 BLOSUM计分矩阵

不少情况下PAM计分矩阵可能失效, 因为其进化距离较远的矩阵是推算出来的而不是直接计算得到的, 其准确率受到一定限制, 这就需要使用新的计分矩阵。BLOSUM矩阵是由Henikoff首先提出的另一种氨基酸替换矩阵, PAM矩阵是从蛋白质序列的全局比对结果推导而来的, 而BLOSUM是从蛋白质序列块 (短序列) 比对而推导出来的。然而BLOSUM矩阵号数的意义与PAM相反, 低号BLOSUM矩阵适合用来比较亲缘关系较远的序列, 而高号BLOSUM用来比较亲缘关系较近的序列。常见的BLOSUM62矩阵可用来比较大约有62%相似度的序列。

2.3 Needleman-Wunsch算法

Needleman-Wunsch算法是基于动态规划的全局比对算法。对于长度分别为m和n的两个序列S和T, 构造矩阵, 构造mn矩阵, 假设序列S放在上面, 序列T放在矩阵左边, 并使用BLO-SUM50计分矩阵, 空位罚分取d=8的线形罚分函数r (g) =-gd。σ (i, j) 表示矩阵中的第i行, 第j列元素的最优匹配得分, 其中σ (00) =0, σ (i, 0) =-id, σ (0, j) =-jd。

用动态规划的方法建立计算σ (i, j) 的公式 (1) 如下:

公式 (1) 中的s (i, j) 这里代表两个氨基酸残基的比对得分, 可从BLOSUM50表中获取。由公式 (1) 得出表11的最优联配得分矩阵。 (如表1)

σ (m, n) =1是这两条氨基酸序列的总得分。接下来我们要做的就是从最右下角的单元回溯, 找到最优联配结果。通过得分矩阵进行动态规划回溯法获得相似性比较的算法如下:

其中insert (a, b, c) 函数表示在一个序列b的第c个位置插入一个字符a, 符号'-'代表空位。

最终的联配结果如下:

3. Smith-Waterman算法

Smith-Waterman算法是基于局部联配的动态规划方法, Needleman-Wunsch适用于整体水平上相似度较高的两个序列。如果两个序列的亲缘关系较远, 它们在整体上可能不具有相似性, 但在一些较小的区域上却可能存在局部相似性。Smith-Waterman算法是序列局部比对算法的基础。它的思想与Needleman-Wunsch基本相似, 区别有以下几点:

1) 不需设置负分, 因为两个零长字符串即为得分为0的局部比对, 这一事实表明在构建局部比对时, 不需要使用负分, 所以σ (0, 0) =0, σ (i, 0) =0, σ (0, j) =0。那么求解σ (i, j) 的动态规划算法的公式改为如下

2) 矩阵中每个单元均可以是比对结果序列片段的终点, 于是回溯的时候应该从得分最大的元素开始回溯。

最终的联配得分矩阵如表22所示

根据表2我们得到最优局部联配为

4. 算法改进--寻找所有最优比对

Needleman-Wunsch算法和Smith-Waterman算法在寻找最优比对的时候未考虑存在多种最优联配的情况。例如在回溯过程中有可能M[i][j]==M[i-1][j-1]+s (i, j) 同时M[i][j]==M[i-1][j]+d, 但算法中只取一种匹配情况。

为了遍历所有路径, 本节中提出了使用链队列结构的方法, 具体思路如下:

1) 从单元 (i, j) 出发开始回溯, 将单元 (i, j) 作为树的根结点

2) 判断单元 (i, j) 结点可以由 (i-1, j-1) , (i-1, j) , (i, j-1) 单元中的哪些生成, 如果只有一个单元可以到达 (i, j) , 那么将其作为 (i, j) 的孩子, 然后对 (i, j) 的孩子重复第二步;否则则将所有可以到达的单元 (i, j) 都作为 (i, j) 的孩子, 并将所有孩子入队列, 然后对队列中的单元依次重复步骤2) 。

当遍历完成后, 所有路径都已存储在树中。

5. 基于Needleman-Wunch算法或Smith-Waterman算法的双路并行算法

Needleman-Wunch算法和Smith-Waterman算法的大部分时间花销在于构造联配得分矩阵, 本节中提出了一种双路并行的算法可以将此算法提高一倍。

观察初始得分矩阵可知, 当单元 (i, i) 计算出来后, 那么它所在的一行, 一列可以同时计算出来, 于是我们可以设计两个并行块, 一块专门计算行, 一块专门计算列, 并利用多处理器或者多核处理器实现并行计算。

6. 结束语

蛋白质双序列比对的算法很多都基于基于NeedlemanWunch算法或Smith-Waterman算法所做的改进, 本文在执行效率和遍历方法上对算法进行了些改进。

摘要：Needleman-Wunsch算法和Smith-Waterman算法是蛋白质序列比对的两种重要方法, 根据实际需要选择不同的计分矩阵和算法可以达到较好的比对结果。但这两种算法有其缺陷, 本文在此两种算法的基础上提出以链队列的形式遍历所有最优匹配的算法, 并优化了得分矩阵的计算方法, 提出双路并行计算得分矩阵的方法。

关键词：PAM250,Needleman-Wunsch算法,Smith-Waterman算法,双路并行

参考文献

[1]蔡禄著.生物信息学教程[M].北京:化学工业出版社, 2007

多序列比对算法综述 篇2

关键词：生物信息学,序列比对,动态规划,算法

一、背景与意义

随着人类基因组计划的顺利实施和信息技术的迅速发展, Gene Bank、EMBL和DDBJ国际三大核酸序列数据库数据量呈指数增长。生物学家、数学家和计算机科学家都面临着一个相同的并且严峻的问题, 如何利用、表达这些数据进而分析与解释基因序列间的潜在关系, 从中发掘出对人类有利的信息。为迎接挑战, 一门涉及生物、数学、物理、化学、计算机科学等诸多科学的新型学科应运而生, 并且日益成为二十一世纪自然科学的核心领域之一。

生物信息学的主要研究对象是DNA序列和蛋白质序列, 主要通过分类、分析和检索核苷酸序列或氨基酸序列, 获取基因编码和调控、代谢途径、DNA和蛋白质结构功能及其相互关系等各方面的知识。所以在生命起源、生物进化以及细胞、器官和个体的出现、生长、病变、消亡等生命科学问题中, 生物信息学都起着非常重要的作用。生物信息学是发现生命科学问题中的基本规律和时空联系, 发掘生物序列数据中蕴含的生物学意义的交叉学科[1]。

在生物信息学中, 序列比对是最基本、最重要的操作。对于基因序列, 通过比对可以推测出哪个基因家族可能包含该序列, 并可以推测出该序列可能具有的生物学功能;对于蛋白质序列, 通过比对可以推测出该序列可能的功能和结构, 并可以找出与它同源的蛋白质序列。所以在生物信息学中, 序列比对具有非常重要的意义和实用价值。目前, 国际上提出了众多经典的比对算法, 也开发了众多的序列比对软件。但对于同一组序列, 不同的软件采用不同的序列比对算法, 其运算速度和比对结果都有很大差异。有些软件考虑了比对结果而运行时间较长, 而有些软件运算速度很快比对结果却不理想。一般情况下两者不能同时兼得。所以, 对于序列比对算法的研究还有待继续深入。

二、多序列比对的研究现状及发展趋势

多序列比对是双序列比对的一般性推广。由于核酸数据库容量的增长呈指数级, 序列比对的数量通常都会远远大于两个, 使用动态规划算法来解决比对问题已经是不可行的了, 这使得多序列比对成为一NP难题。为了解决这一难题, 人们提出了许多近似算法。

2.1动态规划算法

多序列比对最早采用的是动态规划算法来解决。动态规划算法中最为经典的是Needleman-Wunsch算法, 其解决思路是把整个问题分解成多个相互联系的子问题, 通过依次解决每个子问题, 从而解决整个问题。动态规划算法最初用于求解两个序列的比对, 当把动态规划的基本思想推广到多序列比对时, 3个很短序列的比对还可以顺利进行。比对序列的数量如果超过3个, 由于需要很大的存储空间和很长的运行时间, 比对根本无法进行下去。所以多序列比对问题不能采用动态规划算法来解决。Carrillo和Lipman等人对该算法进行了改进, 提出了Carrillo-Lipman算法, 通过减少存储空间将序列比对的数量提高到10。2004年, 唐玉荣等人对动态规划算法进行了优化[2], 与基本动态规划法敏感性相同, 但降低了算法的时间复杂度, 并在减少存储空间方面也有一定的效果。

2.2启发式算法

目前, 绝大多数算法属于启发式算法, 包括星比对算法, 渐进式比对算法, 迭代细化方法等。其中应用最早的是星比对算法, 而应用最广并且效果较好的是渐进式比对算法。Hogeweg和Hesper首先提出渐进式比对算法, 而后Feng和Taylor对其加以完善。与动态规划算法相比, 该算法在计算速度、存储空间和序列数目等方面都更加优良。并且, 渐进式比对算法能够直接用于构造进化树, 反映序列间的进化关系。2005年, 段敏等人提出了一种用减少序列比对过程中总评分的方法来达到局部优化目的的多序列比对算法[3]。启发式算法虽然在一定层度上减少了算法的运行时间和存储空间, 但都有一些不足之处。星比对算法中, 无论采用何种方法并不能保证找到的序列是最好的中心序列。渐进式比对算法中, 构造的指导树有时不一定真正反映系统的进化信息, 根据指导树渐进比对容易产生局部最优化问题。迭代细化算法中, 无法采用何种迭代策略得到的结果最优。

2.3随机算法

多序列比对中, 应用最多的随机算法有遗传算法、模拟退火算法和粒子群算法等。遗传算法是一种全局意义上的自适应随机搜索方法, 它借鉴生物进化规律, 模拟生物进化过程中的一系列事件, 包括突变、交配和选择, 最终得到一个优化解。模拟退火算法则是模拟物理中的退火过程并结合复杂系统中的组合优化之间的相似性来寻找最优解。2008年, 向昌盛等人提出了将遗传算法和模拟退火算法相结合的遗传退火进化思想[4], 设计了运用该思想进行多序列比对的算法过程, 实验结果表明该算法是行之有效的。2011年, 徐小俊等人针对粒子群优化易陷入局部最优、收敛速度慢的现象, 提出了一种分段取值惯性权重 (SW) 方法[5], 该方法在解决多序列比对问题时可以有效地避免算法早熟, 并提高解的精度。

2.4分治算法

分治算法是把一个大问题分解成若干个小问题来解决。Stoye提出了一种新的分治算法DCA, 将Carrillo-Lipman算法引入进来。在不影响特征表现的前提下, 把序列分割成完全满足Carrillo-Lipman算法长度要求的子序列, 使用Carrillo-Lipman算法进行序列比对。2000年Stoye又提出了一种OMA算法, 以达到减少存储空间的目的。2009年, 业宁等人设计了一个DCA-Clustal W算法来解决多序列比对问题[6], 从纵向和横向两个方面将复杂问题简单化, 并在Bali Base基准数据集上测试了算法的可行性。

2.5其他算法

2006年, 陈娟等人给出了多重序列比对的蚁群算法[7], 结果显示蚁群算法可以有效解决多重序列比对问题并具有自适应性、鲁棒性等特点。而文献[8, 9]针对蚁群算法易于陷入局部最优解、收敛速度慢等问题, 提出了改进的方法。

三、结论

生物信息学中最基本和最核心的问题之一就是多序列比对。由于多序列比对处理的数据越来越庞大和复杂, 所以其算法对计算精度和运算速度的要求也越来越高。如何能快速有效的获得比对的结果, 一直苦恼着众多的学者们。

参考文献

[1]孙啸, 陆祖宏, 谢建明.生物信息学基础.北京:清华大学出版社, 2005, 10-17

[2]唐玉荣, 一种优化的生物序列比对算法.计算机工程与设计, 2004, 25 (11) :1936-1945

[3]段敏, 许龙飞.生物DNA序列比对算法研究.佳木斯大学学报, 2005, 23 (2) :153-158

[4]向昌盛, 周建军, 周子英.模拟退火遗传算法在生物多序列比对中的应用.湖南农业科学, 2008, (4) :29-34

[5]徐小俊, 雷秀娟, 郭玲.基于SWGPSO算法的多序列比对.计算机工程, 2011, 37 (6) :184-186

[6]业宁, 张倩倩, 许翠云.一种多序列比对分值算法DCA-ClustalW.计算机与数学工程, 2010, 38 (11) :30-33

[7]陈娟, 陈.多重序列比对的蚁群算法.计算机应用, 2006, 26 (6) :124-128

[8]彭东海, 骆嘉伟, 袁辉勇.基于改进蚁群算法的多序列比对.计算机工程与应用, 2009, 45 (33) :114-119

多重序列比对的模型与算法 篇3

多重序列比对是生物信息学的一个重要问题, 它在序列装配、蛋白质结构预测、蛋白质结构功能分析、系万方数据统进化分析、数据库检索以及引物设计等问题的研究中被广泛使用。由于MSA的计算是一个十分困难的问题, 因而开发求解MSA的算法一直是一个很活跃的研究领域, 许多实用的算法已被提出.目前关于MSA的算法很多, 使用最为广泛的MSA算法是Clustal W。

二、模型建立与求解

1、衡量多重序列比对好坏的定量优化模型

从数学的观点讲, 多重序列比对是一个优化问题。有K=12个生物序列:

且满足下面个条件:

(1) 矩阵的行数等于序列的数目;

(2) 如果移去每行中的“-”字符, 将得到原来的序列;

(3) 任一列中不允许同时为“-”。

为了评价一个比对质量的优劣, 需要定义一个适当的目标函数一种常用的目标函数是基于SP (Sum of P air) 打分系统给出的, 其定义如下:

其中Zsore (bi, bj) 表示由比对B对序列Si和Sj导出的配对比对的计分。

本文只讨论核酸 (DNA) 序列的比对, 并定义Zsore (bi, bj) 如下:

这样, 具有生物学意义的多重序列的最佳比对问题就成为如下优化问题:

其中Γ表示序列S1, S2, L, SK的所有可能多重序列比对集合。

2、基于遗传算法的MSA

(1) 初始种群的产生

(2) 适应度计算

(3) 交叉操作

(4) 变异操作

(5) 终止条件

3、基于免疫遗传算法的MSA

本文采用的模式是在遗传算法的框架中增加一个免疫算子。其中免疫算子 (在遗传算法第四步和第五步加入) 操作分三步进行: (a) 提取疫苗; (b) 接种疫苗; (c) 免疫检验。

4、算法实现和分析

为了检验算法的有效性, 本文已将算法用C语言 (见附录一) 实现, 并将题中的12个序列实验, 给出使序列有最大相似程度的比对, 结果如图2。

计算时, 取群体规模N=100, 杂交率为50%, 变异率为20%, 疫苗长度t为个体长度的1/10, 接种疫苗的个体数为群体规模的1/5, 终止条件为群体中前n (=20) 最好个体的适应度值之和保持D (=500) 代不变。经计算所得的最好比对B的分值为795。相同的参数情况下采用遗传算法 (即去掉免疫遗传算法中的免疫算子) 所得结果见图3。

上述实验结果表明, 与遗传算法相比, 本文给出的免疫遗传算法不仅收敛的速度快, 而且还可以求得更优的比对。

三、结论

本文给出的免疫遗传算法是在保留遗传算法优良特性的前提下, 利用MSA问题最优解的特征来克服遗传算法收敛速度慢、局部搜索能力差和容易出现未成熟收敛的缺陷。若B是问题的最优比对 (个体) , 则B中任一片段也是最优的, 因此, 用最优个体中的一个片段 (疫苗) 去替换非最优个体中对应片段 (接种疫苗) , 一般会增强其个体的抗体 (适应值增高) 又因为每一代种群更新都要求种群中N个个体必须是互异的, 这样, 种群的多样性得到保证, 经免疫操作后, 种群的整体质量得到了提高。

由此可想到算法的一种修正方案:选出种群中的最佳个体之后, 截取其中哪个片段作为接种疫苗由选定的接受接种疫苗的个体来确定, 接种的个体不同, 其接种的疫苗可以不同, 使其达到最佳匹配。这样, 算法的效率将会得到进一步改进。

摘要：本文就基因的多重序列比对, 为设计合理的衡量比对好坏, 建立多重序列比对SP打分系统的的优化模型, 并对随机12个原始序列进行检验得分为161, 长度44。在分析国内外序列比对算法及序列比对的优化问题本质的基础上, 结合生物免疫系统的特点提出一个基于免疫遗传算法的MSA方法。

关键词：MSA,免疫遗传算法,SP打分系统

参考文献

[1]李素贞、莫忠息、张轩、陶玉敏:《基于免疫遗传算法的多重序列比对》, 《武汉大学学报》 (理学版) , 第50卷第5期2004.10。

[2]周康:《关于DNA序列比对算法的简述》, 《武汉工业学院学报》, 第24卷第2期, 2005.6。

[3]刘桂霞、常群、朱元贤、周文刚、周春光:《基于免疫遗传算法多序列比对的实现》, 《吉林大学计算机科学与技术学院》。

[4]刘帅、马志强、刘清雪、陆林英:《基于自适应免疫遗传算法的多序列比对》, 《东北师范大学计算机科学系》。

生物序列比对算法的研究现状 篇4

序列比对有双序列比对和多序列比对之分, 常见的双序列比对算法点阵图方法和动态规划算法, 而多序列比对算法主要有渐进比对和迭代比对两大类。

1、点阵图法

点阵图法[1]的基本思想是通过将一条序列排在上端, 另一条序列纵列在左端, 两个序列在任何位置上若出现相同残基, 就在两个序列对应位置上标注一个点, 做成一个图。排列成对角线的点列体现出两条序列间具有相同的字符串, 从而形象地表明序列间的相似性, 双序列点阵图示如图1所示。

点阵图法的主要优点在于可以找到序列间的所有可能的残基匹配, 但主要的局限是点阵计算机程序并不能显示真实的比对排列。

2、动态规划算法

动态规划算法主要有全局排列和局部排列两大类。

全局排列动态规划算法是由Needleman和Wunsch于1970年首先提出的[2], 算法的基本思想是:用比对的两条序列构建一个相似打分矩阵S, 矩阵中的元素可通过公式 (1.1) 获得。

这里, 是序列a在位置i和序列b在位置j的分值。是序列a位置i和序列b位置j上排列性状的分值。Wx是序列a中长度为x的间隔罚分, Wy是序列b中长度为y的间隔罚分。由于算法遍历矩阵的每个点, 具有最佳罚分。

例1:对序列a=GCTGATATAGCT, b=GGGTGATTAGCT, 选择参数s (a, a) =1, s (a, b) =-1, 插入删除单个字母的罚分为2, 计算相似打分矩阵S如图2。

根据上述矩阵可得最优比对如图3所示:

局部排列动态规划算法是由Temple Smith和Michael Waterman于1981年提出的, 同样算法也是通过比对的两条序列构建一个相似打分矩阵 (记为H) , 矩阵中的元素Hi, j可通过公式 (1.2) 获得。

这里, Hi, j是序列a在位置i和序列b在位置j的分值。w (ai, bj) 是序列a位置i和序列b位置j上排列性状的分值。Wx是序列a中长度为x的间隔罚分, Wy是序列b中长度为y的间隔罚分。

3、渐进比对算法

渐进比对算法属于启发式的多序列比对算法。最常见的渐进比对算法就是由Feng和Doolittle提出的Clustal算法, Clustal的基本思想是基于相似序列通常具有进化相关性这一假设。比对过程中, 先对所有的序列进行两两比对并计算它们的相似性分数值, 然后根据相似性分数值将它们分成若干组, 并在每组之间进行比对, 计算相似性分数值。根据相似性分数值继续分组比对, 直到得到最终比对结果。比对过程中, 相似程度较高的序列先进行比对, 而相似性较低的序列则添加在后面。Clustal算法的主要三个步骤如下:

(1) 两两比对:先将比对序列进行两两比对分别构建距离矩阵。

(2) 系统发生树构建:根据计算所获得的距离矩阵构建系统发生树。

(3) 进化式比对:对关系密切的序列进行加权, 然后从紧密的两条序列开始, 逐步引入临近的序列并不断重新构建比对, 直到所有序列都被加入为止。

4、迭代比对算法

迭代比对算法的基本思想是基于一个比对算法, 通过迭代方式精细多序列比对, 直到比对结果不再改变为止。根据迭代策略的不同, 迭代比对算法大致可分为Prrp法, 隐马尔科夫法, 模拟退火法和遗传算法等。

小结

进行序列比对的目的之一是让人们能够判断两个序列之间是否具有足够的相似性, 从而判定二者之间是否具有同源性。值得注意的是, 相似性和同源性虽然在某种程度上具有一致性, 但它们是完全不同的两个概念:相似性是指一种很直接的数量关系, 比如部分相同或相似的百分比或其它一些合适的度量;而同源性是指从一些数据中推断出的两个基因在进化上曾具有共同祖先的结论, 它是质的判断。基因之间要么同源, 要么不同源, 绝不像相似性那样具有多或少的数量关系。

摘要：本文综述了生物序列比对的基本思想和主要方法。通过序列比较可以发现生物序列中的功能、结构和进化的信息。进行序列比对的目的是让人们能够判断两个序列之间是否具有足够的相似性, 从而判定二者之间是否具有同源性。

关键词：生物序列,比对,算法

参考文献

[1]Gilbert DG.Dot plot sequence comparisonson Macintosh computers.Comput.Appl.Biosci., 1990, 6 (2) :117-117;

蛋白质序列比对算法 篇5

数据库序列比对是生物信息学研究领域中基本性的热点应用之一。研究人员经常需要面对将新发现的基因序列与已知的基因数据库进行搜索比对的课题状况。而通过基因序列比对,或者可以找出两个序列间的相似区域,判断该基因序列的新颖性,或者将两类生物体进行全基因组(full genome)序列比对,以判断两者间的同源性等。在已有的众多序列比对算法中,Smith-Waterman算法[1]最为著名。该算法基于Needleman-Wunsch算法[2]的改进,可用于局部序列比对,能够得到精确的序列比对结果。而Smith-Waterman算法在经过Gotoh[3]的进一步完善后,算法复杂度达到了O(M×N)。其中,M和N分别为两个比对序列的长度。但是,Smith-Waterman算法仍然表现了复杂度大,计算耗时的不足,而与此同时,随着各类生物序列数据库规模的逐年增大,其计算速度已经越来越难以满足实际应用要求。

近年来,出现了很多加速方法,按其原理可将其组织为两类:一是在算法层次进行改进,在降低算法复杂性后,提高速度;二是通过硬件加速平台,提升算法速度。

对于算法改进,BLAST[4]和FASTA[5]是算法改进后的典型代表。这些算法都基于启发式算法,可查找出两个比对序列中同时存在的某个特定长度的字或者种子,对这些字或种子所在的区域实现扩展,进行相似性比对。只有当种子在两个序列中都有命中的区域时,才进行拓展比对,因此大大缩短了运行时间,但是却无法得到与Smith-Waterman算法相同的精确结果。

对于通过特殊硬件平台加速,P-NAC[6]、BISP[7]以及SAMBA[8]等利用FPGA硬件实现核心算法,速度快,但是价格昂贵,缺乏灵活性,不易于扩展。在多处理器上并行Smith-Waterman算法也是常用的方式[9,10,11,12],但是商用多核处理器价位较高,相对低价的GPU逐渐成为时新的加速平台。Liu[13]利用GPGPU实现了基于Smith-Waterman的数据库序列比对,并相比于SSEARCH获得了多达5倍的性能提升,而Manavski[14]则第一次利用NVIDIA的计算统一架构(CUDA)平台付诸实现,相较Liu的研究结论又获得18倍的加速。之后,CUDA编程环境的GPU即成为研究人员聚焦关注的平台之一[15,16,17],采用CUDA编程模型[18],可方便地完成划分和并行任务。本文在基于CUDA平台的NVIDIA GeForce 8800GT上加速Smith-Waterman算法,获得了约3 000MCUPS(Million Cell Updates Per Second)的比对速度,其性能超过了基于启发式算法的BLAST在通用CPU上的速度。

本文组织如下:第1节和第2节分别简单介绍基于Smith-Waterman的序列比对算法,以及NVIDIA 8800GT显卡的结构特征和CUDA编程环境,第3节提出了基于CUDA编程环境的算法实现和优化,并对性能数据进行了分析,第4节概述本文的相关工作,最后一节总结全文工作。

1 序列比对算法

数据库序列比对是将查询序列与数据库中各个序列通过逐一运行序列比对算法,在完成比对后计算得到各次比对的结果分值,如果分值超过了指定阈值,就说明数据库中的该条待检序列与查询序列匹配成功,可以进入后续的其他操作流程。

序列比对算法的输入为两个序列文件,其输出为序列间的相似区域,下面是一个蛋白质序列比对的例子,具体细节如图1所示。

由图1可见,字符串代表由20个不同氨基酸组成的蛋白质序列,“|” 代表序列间这两个字符完全匹配,“:”代表序列间这两个字符近似匹配,可以根据指定的置换矩阵(substitution matrix)来判断哪些字符属于近似匹配,“-“表示插入了一个空位。

Smith-Waterman算法是进行序列比对的有效方法,通过利用动态规划方法来计算最优的局部序列匹配,而且多是采用Gotoh改进后的算法。算法具体过程如下:

假定序列A为查询序列,长度为|A|,B为数据库序列,长度为|B|,alpha为插入第一个空位的开销,beta为插入后续空位的开销,sbt_matrix为置换矩阵。此时,局部序列比对分值是在一个|A|×|B|的矩阵范围内应用动态规划的方法计算得到的各个点的最大值H(i, j),计算公式如下所示:

H(i, j)=max{E(i, j),H(i-1,j-1)+sbt_matrix (A[i],B[j])),F(i,j),0}

where (1)

E(i, j)=max{H(i, j-1)-alpha,E(i,j-1)-beta}

E(i, j)=max{H(i-1, j)-alpha,F(i-1, j)-beta}

其中,第一行和第一列的元素L(0, j),E(0, j),F(0, j),L(i,0),E(i,0),F(i,0)均初始化为0。H(i, j)则表示A序列前i个字符和B序列前j个字符进行比对时,当前所得到的最大分值。

由式(1)可以看出,矩阵中各个点的计算存在着严格的依赖关系,每个点的计算都取决于同一行左侧一列的点、同一列上一行的点以及左侧一列、上面一行的左上方点的数据值。因此,在实际计算时,可以按照行、列或者是反对角线的顺序依次计算。算法复杂度为O(|A||B|)。

2 GPU体系结构和CUDA编程环境

根据图像应用要求,GPU的体系结构能高效、可靠地支持数量庞大的细粒度的线程调度,而这一点在传统的通用并行计算领域也存在着广阔的应用空间。CUDA是NVIDIA于2006年推出的计算架构,凭此可发掘GPU强大的计算能力,给密集型应用计算带来了新的解决契机。

如图2所示,GeForce 8800 GT芯片由14个流多处理器(SM)组成,每个SM均包含有8个小的处理器核(processor),这8个小核共享一个指令单元,可通过SIMD的方式执行指令,运行频率为1.5GHz。同时,每个小核又包括一个32-bit的单精度浮点乘/加计算部件。此外,每个SM上还有2个特殊功能单元(SFU),用来完成复杂的浮点运算,如开方、求取正/余弦等。由此,GPU的整体运算能力已可达到378GFLOPS。

值得一提的是,8800GT的存储系统包括片上的寄存器、共享内存(shared memory)和cache、以及片外的device memory。其具体特性如表1所示。

CUDA是对标准C/C++的扩展,并提供了一些额外的数据类型、关键字以及函数调用方式等。在CUDA模式下,GPU就是一个专门用于加速的协处理器,程序中的串行部分由主机CPU来完成,而并行部分则包装成一个kernel函数,由主机调用GPU执行。CUDA上的线程组织方式分为grid-block-thread三级[19]。最顶层是grid,执行同一kernel函数的全部thread构成一个grid,Grid可包含多个block,每个block又包含多个thread。通常,一个SM上最多含有768个thread,而对block,则是最多容纳8个。执行时,位于同一block内的thread组织成warp形式,每个warp可包含32个连续thread。warp是GPU调度的基本单位,当SM上的小核执行到长延时指令时,系统将快速切换到另一个可以待命执行的warp,直到SM上没有可执行的warp时,则暂停,以此减少等待耗时。

3 算法实现及性能评估

对于访存性能提升策略,CUDA与通用CPU存在很大的不同。CPU主要是通过设置片内cache来加速访存,而CUDA遇到长延时访存指令时,将快速切换到其他线程执行,切换过程开销很小,因此,需要创建足够多的线程用来实现切换调度,从而充分利用访存操作的长延时。

在进行数据库比对时,查询序列和数据库序列的比对操作相互完全独立,直接的方法就是使每个线程都要完成一个查询序列与一个数据库序列的比对过程,这样,m个查询序列和n个数据库序列就可以创建m×n个线程。实际数据库中,n值常常较大,数量级在5以上,如此数目的线程数完全满足CUDA的调度需求。本文实验使用Swiss-Prot[19]蛋白质数据库(数据库共含有366 226个蛋白质,132 054 191个氨基酸),在数据库中提取查询序列,Smith-Waterman算法所使用的置换矩阵为BLOSUM50,首次空位插入开销为10,后续空位插入开销为2。

3.1 访存操作的优化

CUDA上访存操作对性能的影响主要表现在两个方面,一是访存延时大,需要超过200个时钟周期;二是对访存带宽的要求高,GPU的计算能力强达378GFLOPS,但是访存带宽却不到100GB/S,如果计算访存比低,将导致带宽形成短板。因此,除了需要创建大量调度线程来匹配访存延时外,还应该尽量减少程序中的访存操作。这样做目的主要有三个:一是降低对访存带宽的需求;二可降低对线程数的需求,避免出现活动的线程数不足引发的等待;三能保持程序的连贯性,减少线程执行中的切换。

由图3可知,算法中访存操作主要有以下几种:(1)读取查询序列A中的字符|A|次;(2)读取数据库序列B中的字符|B|次;(3)读取置换矩阵sbt_matrix中的相应项|A|×|B|次;(4)读写临时结果数组HE_in[]和HE_out[],读和写各|A|×|B|次。

在本文的实现中,借鉴了文献[11]和[16] 中采用的query-sequence profile方法,即根据查询序列和置换矩阵计算得到的每种字符对应的匹配分值,再按照查询序列的字符排列完成连续存放好。如图4所示,列序列是一个蛋白质查询序列示例,行字符代表氨基酸的种类(A～Z表示)。每一列记录了查询序列对应某种蛋白质的比对分值,矩阵中所有数据按列连续存放。当查询序列和数据库序列的某一个字符进行比对计算时,该字符对应的所有置换分值均可从内存中连续读到,无需再查找置换矩阵。query-sequence profile方法可以减少(1)、(3)种访存操作的总数量,同时,将随机访问转换成连续访问能够更好地利用cache或者数据读写合并,来实现加速访存。

profile信息表大小取决于查询序列和序列中字符种类的乘积。大多情况下,shared memory、constant memory以及texture memory都可以满足需要,但不同的数据存放方式决定了其性能也将有所不同。在CUDA模式下,可以显式指定数据的存放位置。图5给出了3种对应不同查询序列长度时的存储结果,可以看到,texture memory在这些长度下的速度都是最快的。下面将分析讨论最佳的存储方式。

使用shared memory时,要避免bank冲突,冲突发生时访问速度将会下降。而本文使用的查询序列长度却会使得warp内的并行thread访问同一bank,导致所有访问串行发生,势必影响了数据读取速度。此外,受到shared memory容量的限制,查询序列不可能太长;而且写shared memory时,所有线程都需要增加一个同步操作,同步操作会对单个block内的最佳thread数带来影响。

而对constant memory的访问,warp内的线程访问同一地址时,速度最快,否则,访问时间将随着访问的地址数呈线性增长。在本文中,形式各异的线程都在和数据库中的不同序列进行比对,因其访问不同种类字符的概率很高,决定了每个线程的访问地址都不一致。随着查询序列长度的增加,冲突越发明显,性能也逐渐变差。图5(b)显示,在查询序列长度达到127后,其比对速度即开始低于shared memory的性能,使用constant memory的扩展性也不好。

因此,选用texture memory来保存profile信息,一方面在容量上不受限制,另一方面速度又最快。借助profile,获得了39.3%的性能提升,如图6中的SW_prof所示。

第四种访存操作针对计算时临时生成的两列H值和E值(假设按列进行计算,按行计算则为H值和F值),需要写HE_out和读HE_in。由表1可知,支持读写操作的只有global memory、local memory或者shared memory。计算过程中,每个线程均需要保存长度等于查询序列的两列数据, 首先,shared memory空间限度无法满足该要求。若使用global memory,需要预先做好分配,当kernel中的线程数目较多时,同样会出现空间无法应对的问题,虽然可通过减少单个kernel线程数并增加kernel数来完成相同的任务,但线程数的减少却容易造成调度困难,从而影响整体效率提升。另外,多个kernel函数还将引入创建和退出的额外开销。实验统计表明,只有在查询序列较短时,global memory利用向量访存操作性能高于local memory。因此,文中采用local memory来保存H值和E值。但是local memory没有cache,速度很慢,因此要尽量减少访问次数。实际上,计算出来的分值都不大,用16位整数完全可以表示,将这种位数较少的读写操作进行合并,可以明显减少访存次数。目前,CUDA尚不支持对local memory向量读写,因此可将2个16位的访存进行合并,使用一条32位的读写指令即可完成。合并后的性能比SW_prof提高了66.3%,具体可参见图6中的SW_HEpk。

与对HE值的合并类似,第二种访存操作也可以合并。数据库序列都是8bit字符,存在global memory中,CUDA提供对其128位向量的取数操作,可一次提取16个字符到寄存器,这样就将原有的|B|次访存降为|B|/16次。结果性能如图6中的CUDA.dbpk项所示,相比SW_HEpk又获得了18.8%的提升。

3.2 其它优化

在CUDA模式下,一个warp内的所有thread采用SIMD的形式,执行相同的指令。为了充分利用计算资源,同一warp内thread的计算量应该尽量接近,使负载达到均衡。在执行算法时,计算量与两个序列的长度乘积成正比,且不同thread的查询序列相同,因此就需要保证相邻thread得到的数据库序列长度基本接近。为此,将对所有数据库序列进行预先降阶重排,降阶排列可以使得比对序列长、计算任务重的thread所在的block优先执行,留下的则是计算任务较轻的block,这样,GPU的多个SM就可以在更加接近的时间内全部完成,从而实现良好的负载均衡。其性能如图6中的SW_sort所示,这是在SW_init的基础上进行降排序得到的结果,排序后获得了1.82倍的速度提升。图6中的SW_sort右侧的各项结果都是在已排序的前提下而得到的。

在其他资源确定不变情况下,需要尽可能地提高GPU的占用率(occupancy),也就是要提高每个SM上正在执行的thread数(以warp为单位)。不同数据库序列和查询序列之间的比较完全独立,block内的thread并不存在同步的模式,所以可将所有thread组织在一个block内,使得单个SM上活动的thread数尽可能多,由此提高占用率。

最后,根据编译器生成的PTX文件,对kernel函数的内存循环代码进行手工的优化,重点在于减少寄存器的使用以及消除冗余的操作。此时,发现代码中存在大量的数据转换指令,可通过调整变量的类型,而将其消除。优化过后,每个线程使用的寄存器由24个减少到21个,而每个SM上的可活动线程数由320增加到了384个,由此获得了大约4%的性能提升,如图6中的SW_tune所示。图6中的其它结果却是在每个block包含320个thread的情况下得到的。

经过上述优化后,最终的计算访存比就达到了18:1。在这种计算访存比的意义作用下,可假设访问device memory的开销是300个时钟周期,那么每个SM上只需要6个warp,即192个thread就可以匹配访存延时,而目前却有384个thread,已经完全能够满足为匹配访存延迟所需进行的规定调度,因此,在访存上已经不会存在明显的性能瓶颈。

3.3 性能对比

在实验最后,将本文得到的结果和已有的加速方法进行了对比,对比结果如表2所示。表中,SW_tune为本文的结果实现,其配置为:使用了单块NVIDIA GeForce 8800GT型号的GPU,主机为Intel Core 2 6300的双核CPU,运行频率1.86GHz。Nvcc版本为release 1.1, V0.2.1221,且程序使用了-O3优化。SW_SAM为文献[16]中基于单块GeForce 8800GTX的CUDA实现。SW_Liu来源于文献[17],是基于GeForce 7800GTX,使用OpenGL接口的实现。Ssearch和Blast的结果数据均来自文献[16],基于3 GHz的Intel Pentium 4处理器平台。

由表2可以看到,本文的实验实现在性能上均优于其它几类的实现成果。SW_Liu是基于GPU的Smith-Waterman实现,本文在对应不同查询序列的长度下实现速度是前者的5.5～15.6倍。SW_SAM也是基于CUDA的实现,但其在降低访存开销、减少寄存器使用以及block/thread的组织等方面并未得到充分优化,与该方法相比,本文依然获得了约60%的性能提升。 SW_SAM实验中的GPU包含16个SM,文中的GPU只包含14个SM,如果在相同的平台背景下,本文方法的优势将更加明显。另外,即使相比于BLAST[6]这种快速的基于启发式的算法,本文也稳定地获得8.7%～106%的速度优势。

4 结束语

本文在基于NVIDIA GeForce 8800GT平台的CUDA模式下实现了Smith-Waterman算法的加速,使用了一系列方法优化程序的访存操作、线程划分和组织等。结论显示,算法在得到精确结果的同时,也保证了较快的速度,达到了约3 000MCUPS,现已超过了基于启发式的BLAST算法在主流通用CPU上的性能。

研究工作表明,CUDA平台对于生物信息学中的序列比对算法获得了良好的加速效果,因此,未来还可以对其它类型的序列比对算法实现加速。同时,也要继续将CUDA的加速效果和在当前主流使用的多核/众核平台上的实验结果进行对比研究,以分析不同平台的优、劣特性,为今后计算平台的选择提供全面、有益的参考依据。

蛋白质序列比对算法 篇6

序列比对是生物信息学研究的基础和前提,是生物信息学中计算的核心。它通过一定的算法对两个或多个核酸或氨基酸序列对齐进行比较,逐列比较其字符的异同,判断它们之间的相似程度和同源性,从而推测它们的分子结构、功能以及进化过程。这对于蛋白质或DNA的结构预测、进化分析和功能分析有很现实的意义。

由于多序列比对是一个NP完全问题,具有极高的计算复杂度,目前尚缺乏快速且准确的多序列比对算法。而启发式优化方法如遗传算法、模拟退火算法、蚁群算法等,在解决这类问题时往往能在合理的时间里取得相对较好的效果。因此,本文提出了一种将遗传算法和模拟退火算法相结合的优化方法用于DNA多序列比对问题。

2 多序列比对问题及数学描述

对于给定的一组DNA序列S=(S1,S2,…Sk), |Si|表示为某一条DNA序列Si(i=1,2,…,k)的长度。Si中的字符取自于字母表Σ,Σ={A,T,C,G},即表示DNA的四个核苷酸类型。对于DNA序列,给定包含N个序列的序列集S={S1,S2,…SN},N≥2,Sk=SK1=SK2…Sknk,k=1,2,…,N,其中,Sij, i=1,2,…,N, j=1,2,…, ni表示一个核苷酸。

我们用S1S2…SK表示为∑上的字符串,那么一个序列比对是字符表∑*的一个字符矩阵,且必须满足下列三个条件:

(1)矩阵共有k行;

(2)矩阵中的第i行去掉“-”后,即得到字符串 Si;

(3)矩阵中不包含字符全是空格的列。

对于x,y属于∑*,定义σ(x,y)为计分函数,表示x,y比较时的得分,表示每一种变异被赋予一定的分值即权值。一个序列比对通常用目标函数来衡量多序列比对结果好坏。多序列比对问题,就是通过在多个序列中插入空格来得到最优匹配,从而发现多个序列之间的相似程度和同源性。

3 遗传模拟退火思想的提出

遗传算法是模拟生物进化论的自然选择和遗传学机制,通过模拟自然进化过程搜索最优解的方法。遗传算法最大的特点是可以快速的搜索出所有解空间里的解,而不会陷入局部最优,同时,遗传算法从问题解的串集开始搜索,而不是从单个解开始,可以进行分布式计算,加快了求解速度。但是,遗传算法的明显缺陷是局部搜索能力较差,这容易导致过早收敛,即“早熟”现象,而且在进化后期搜索效率低。

模拟退火算法伴随温度参数的不断下降,能够从概率意义上跳出局部最优解从而在解空间中随机寻找目标函数的全局最优解。模拟退火算法解决优化问题时可以在局部以一定的概率接受恶化解,并且能够搜索复杂区域,有效地避免陷入局部极小。但是,模拟退火算法的不足在于它不能全面掌控整个搜索空间,可能进入不了适应度高的搜索区域,同时,求得一个高质量的近似最优解花费时间较长,这使得算法效率不高甚至丧失可行性。

遗传模拟退火算法(SAGA)是将遗传算法与模拟退火算法相结合的求解优化问题的算法,它是利用模拟退火算法针对遗传算子进行改进来提高算法的效率,以期待在DNA序列比对中收到良好的效果。一方面,将模拟退火算法引入到遗传算法群体更新的阶段,使交叉变异后的子代与父代一起竞争,这样即保证了群体多样性,又在后期逐步加快收敛速度,克服遗传算法早熟现象,得到全局最优解;另一方面,利用遗传算法求解是从初始解集出发,通过遗传机制解决问题。这两类算法相结合可以大大提高搜索效率。

4 基于遗传模拟退火的多序列比对算法描述

4.1 基本操作

4.1.1 选择适合的编码策略,确定个体空间

将k条序列比对结果存放在一张二维表Ak×l中,其中,K是参加比对的序列数,L=最长序列的长度×1.2,将这样一个比对结果作为算法中的一个个体。

4.1.2 定义适应度函数

我们用评价序列比对的目标函数作为一个比对结果的适应度函数,目标函数值越高说明适应度越强,选择进入下一代的概率就越高,反之,适应度越低,进入下一代的概率也越低。一个双序列比对结果的目标函数如下:

undefined(A)} 公式(4-1)

其中σ[s′[i],T′[i]]是序列S和序列T进行比对后,第i个位置比对的得分;σ[s′[i],T′[i]]-GP(A)是双序列比对Ai,j的计分值。式(4-1)中的GP(A)表示比对A的空位罚分,对于一个空位长度为k的罚值ωk,可用公式ωk=a+b*(k-1)来表示,其中,a表示空位设置罚分,b表示空位扩展罚分,k为空位长度。

对于多序列比对,常用的目标函数有SP目标函数和COFFEE函数。SP目标函数为多序列比对中的所有双序列比对的计分的总和,COFFEE函数反映了一个多序列比对和双序列比对数据库的一致性程度,它需要先创建一个双序列比对库,当序列之间的相似度较低的时候采用COFFEE函数优于其他目标函数。

4.1.3 设置遗传算法过程中涉及到的各个控制参数

包括群体规模N、交叉概率pc、变异概率pm、最大迭代次数Gmax,可接受的适应度无提升的最大迭代次数Gunimproved,空位设置罚分a,空位扩展罚分b。

4.1.4 产生初始种群

初始化时,在多个序列中随机插入空格来得到初始种群。对于每一条序列Si的随机选择x个位置插入空位“-”,x=L-Si,使得处理后每个序列的长度都为L。

4.1.5 新解的产生机制

遗传算法和模拟退火算法共同作用于新解的产生,假设父代和子代的目标函数值分别为fparent、fchild,接受child进入下一代群体的概率为Pchild,则:

当fchild>fparent时,Pchild=1;

当fchild≤fparent时,undefined。

其中,t为控制参数,在这里是退火温度,与迭代次数G有关:

undefined

4.1.6 遗传算子

4.1.6.1 选择退火算子

选择算子是依照适应度函数的评价从父代中挑选个体进入到新的群体。计算初始种群中每个个体的适应度值,然后进行交叉、变异操作得到新的个体,再计算新个体的适应度值。本文的算法将退火操作加入选择操作中,在得到子代及其适应度值后会依概率选择接受子代进入下一代个体。

4.1.6.2 交叉退火算子

交叉操作是将两个父代的基因进行基因重组,将优秀的基因遗传给下一代个体,并生成包含更复杂基因结构的新个体。交叉操作分以下几个步骤,本文采用一点交叉,具体操作如下:

(1)在群体中随机取出一对父代个体;

(2)根据序列长度L,随机选取[1,L-1]中的一个整数k作为交叉位置;

(3)根据交叉概率pc实施交叉操作,父代个体在交叉位置处互相交换各自部分内容,从而形成新的一对个体即一对子代个体。

图1表示的是一次交叉操作。

4.1.6.3 变异退火算子

在过去人们的研究中,往往将变异操作定义为序列中的某一个碱基变异为其他碱基或将其定义为随机选取序列中的一个或多个位置进行插入空格操作,有的甚至放弃变异操作,认为变异操作增加了序列比对的不确定性。本文将退火算法加入变异算子当中,依概率选择父代,进行变异操作,操作过程为随机选取个体中一条DNA,将其空格全部删除,然后再重新插入空位,得到新的个体。

4.1.7 终止条件

判断新的种群是否满足终止条件,是则结束,否则用新群体代替父群体,循环执行群体更新过程。

当出现以下情况时,算法终止:

(1)群体进化的代数超过最大代数值时;

(2)进化代数超过一定值而适应度值不再提高时。最优个体不再改变,最后一代中适应度值最高的个体即为最优解。

4.2 遗传模拟退火算法的描述

DNA多序列比对的遗传模拟退火算法描述如下:

(1)初始化种群,种群规模N,设定初始迭代G=0。

(2)删除种群中每个个体中全为空格的列。

(3)计算当前种群中每个个体的适应度函数值。

(4)按遗传操作策略开始迭代,运用选择退火算子、交叉退火算子、变异退火算子逐步作用于初始种群,根据当前解产生新解,形成新一代群体,计算新群体中个体的适应度函数值。

(5)计算适应度增量Δf,若Δf>0,接受子代个体替换父代个体;若Δf<0,计算接受适应度值较差的子代的概率Pchild,依概率Pchild接受子代为当前解。

(6)G=G+1;判断循环次数G,若G=Gmax,执行(7),否则返回到(3)。

(7)判断终止条件:进化代数超过一定值Gunimproved而适应度值不再提高,最优个体不再改变。若满足,执行(8),否则用新群体代替父群体,循环执行(3)～(6)。

(8)输出最优解及适应度函数值,程序结束。

5 实验

基于VS 2008.NET,在Intel(R) Core(TM) i3 CPU 2.13GHz 2.00G Mb内存的计算机上,对本文的算法进行了C#编程实现,并进行几组对比实验。实验数据来自数据库NCBI的DNA序列,为得到混合算法的比对效果,下面将分别与著名多序列比对程序ClustalX和经典遗传算法的比对结果进行比较。

经过多次试验,根据所取序列的长度,本文中试验的参数设置为:N=100,pc=0.8、pm1=0.001、pm2=0.01、Gmax=500,Gunimproved=50,a=5,b=1。实验1中取t0为100,实验2中取t0为200。本文采用基于SP的目标函数作为DNA序列比对的适应度函数。

5.1 基于SAGA和CLUSTALX的比对结果分析

实验1在DNA数据库NCBI中选取了5条DNA:GACCAGGATTCT,ACCAGCAATGT,GACCAGTCAATT,GACCTGTACTT,CACCATCATT,分别利用ClustalX和遗传模拟退火算法对上述DNA序列进行比对。图2和图3分别为利用ClustalX和遗传模拟退火算法的DNA比对结果,其中Score为比对的得分即适应度值。

由实验结果可以得出,利用ClustalX进行比对的得分为126,完全比对块个数为5;利用SAGA进行比对的得分为146,完全比对块个数为6。由实验结果可以得出,基于SAGA的比对方法比ClustalX的比对方法在得分上提高了15.87%,完全比对块个数也有所增加。该实验表明利用遗传模拟退火算法来解决DNA多序列比对问题是有效的,基于遗传模拟退火算法的DNA序列比对结果在比对得分和完全比对个数等几个方面均比ClustalX的比对结果优秀。

5.2 基于SAGA和GA的比对结果分析

实验2在DNA数据库NCBI中选取了5条长度约为30的DNA序列,分别利用经典遗传算法和遗传模拟退火算法对上述DNA序列进行比对。图4为利用经典遗传算法和遗传模拟退火算法进行序列比对时迭代次数和最大适应度值的关系。

从图4中可以看出基于经典遗传算法得到的最优适应度值为230,基于遗传模拟退火算法得到的最优适应度值为234,利用遗传模拟退火算法所得的结果比经典遗传算法所得的结果优秀。同时,遗传模拟退火算法在第96次达到了最优值,经典遗传算法在第107次达到了最优值,由此可以看出遗传模拟退火算法在一定程度上提高了算法的收敛速度,减少了计算时间。另一方面,由图中也可以看出,经典遗传算法在进行群体更新的时候毫无保留的选择适应度值高的个体进入下一代,容易陷入局部最优,而遗传模拟退火算法则可以依概率跳出局部最优,最终达到全局最优。

6 结束语

本文将遗传算法和模拟退火算法进行结合,针对遗传算法局部搜索能力较差,容易过早收敛的缺点,在选择、交叉、变异的过程中利用模拟退火算法的Metropolis准则来接受新解,从而大大提高搜索效率。本文对DNA多序列比对的应用进行了有意义的研究,下一步需要研究如何提高算法的高效性和稳定性。

参考文献

[1]Md.Sowgat Ibne Mahmud,Md.Akhter Hosen,Md.Saroer-E-Azam,M.A.Mottalib,Haw-lader Abdullah Al-Mamun.A Novel Two-Tier Multiple Sequence Alignment Algorithm.Proceedings of 13th International Conference onComputer and Information Technology,2010.

[2] S.B.Needleman, C.D Wunsch. A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology, 1970, 48:443～453.

[3] M.A. Larkin, G. Blackshields, N.P. Brown &others. ClustalW and ClustalX version 2.0 BIOINFORMATICS APPLICATIONS NOTE, Vol. 23 no. 21 2007, pages 2947—2948, doi: 10.1093/bioinformatics/btm404.

[4]MillerW.Comparison of Genomic DNA Se-quence:Solved and Unsolved Problems[J].Bioinformatics,2001,17(5):391-397.