（计算机软件与理论专业论文）多序列比对优化方法研究.pdf

摘要 序列比对是生物信息学中一项重要的基础性研究课题，它的最基本任务之一是进行 多序列比对，目前还没有一个最佳的多序列比对算法。本文提出使用遗传算法和粒子群 优化算法来解决多序列比对问题。 首先分析了空位罚分、替换矩阵和目标函数对序列比对的影响，具体实现了s p 和 c o f f e e 目标函数。然后研究并实现了双序列比对的精确算法一动态规划算法，并对基 于渐进方法构建的多序列比对算法c l u s t a d w 进行了深入分析。接着通过对多序列比对 算法的现状的研究以及对遗传算法和粒子群优化算法特点等的分析，提出基于遗传算法 的多序列比对算法m s a g a ( ag e n e t i ca l g o r i t h m d e d i c a t e df o r m u l t i p l es e q u e n c e a l i g n m e n t s ) 和基于粒子群优化算法的多序列比对算法m s a p s o ( ap a r t i c l e s w a r m o p t i m i z a t i o n d e d i c a t e df o rm u l t i p l e s e q u e n c ea i i g n m e n t s ) ，并分别实现了基于s p 和 c o f f e e 目标函数的m s a g a 和m s a p s o ，两种算法复杂度都只与进化代数和种群大 小有关。最后用基准多序列比对库b a l i b a s e 中的用例对算法进行测试，结果表明 m s a g a 和m s a p s o 算法在解决基因序列比对问题上是有效的。 关键词：生物信息学多序列比对遗传算法粒子群优化算法目标函数 a b s t r a c t s e q u e n c ea l i g n m e n to f b i o i n f o r m a t i c si sa ni m p o r t a n tf u n d a m e n t a ls u b j e c ti nb i o i n f o r m a t i e sr e s e a r c h ， o n eo fi t sm o s tb a s i ct a s ki st om u l t i p l es e q u e n c ea l i g n m e n t s s t i l lt h e r ei sn o ta l lo p t i m a la l g o r i t h mo f m u l t i p l es e q u e n c ea l i g n m e n t s t h i sp a p e rp r o p o s e sam e t h o d ，w h i c hu s eg e n e t i ca l g o r i t h ma n dp a r t i c l e s w a r mo p t i m i z a t i o na l g o r i t h mo fe v o l u t i o n a r ya l g o r i t h m ，t os o l v et h ep r o b l e m so fm u l t i p l es e q u e n c e a l i g n m e n t s f i r s t l y ，t h ee f f e c to ns e q u e n c ea l i g n m e n t c a u s e db yt h eg a pp e n a l t y ，s u b s t i t u t i o nm a t r i x a n do b j e c t i v ef u n c t i o na l ea n a l y z e d ，a n dt h es po b j e c t i v ef u n c t i o na n dc o f f e eo b j e c t i v e f u n c t i o na r ei m p l e m e n t e d t h e nt h ed y n a m i cp r o g r a m m i n ga l g o r i t h mt h a ti so n eo fa c c u r a t e a l g o r i t h m so fp a i r w i s ea l i g n m e n t i ss t u d i e da n di m p l e m e n t e d a n dt h ec l u s t a lwt h a ti so n e o ft h em u l t i p l es e q u e n c ea l i g n m e n t sa l g o r i t h mb a s e do nt h eh e u r i s t i cm e t h o di sa n a l y z e di n d e t a i l l a t e rt h r o u i 曲t h es t u d yo nc u r r e n ts i t u a t i o no ft h em u l t i p l es e q u e n c ea l i g n m e n t s a l g o r i t h ma n dt h ea n a l y s i st o t h ep r i n c i p l ea n dc h a r a c t e r i s t i co ft h eg e n e t i ca l g o r i t h ma n d p a r t i c l es w a r mo p t i m i z a t i o na l g o r i t h ma n ds oo n ，p r e s e n tt h em u l t i p l es e q u e n c ea l i g n m e n t s a l g o r i t h mb a s e do ng e n e t i ca l g o r i t h mm s a g a ( a g e n e t i ca l g o r i t h md e d i c a t e df o rm u l t i p l e s e q u e n c ea l i g n m e n t ) a n dt h em u l t i p l es e q u e n c ea l i g n m e n t sa l g o r i t h mb a s e do np a r t i c l e s w a r m o p t i m i z a t i o na l g o r i t h mm s a p s o ( ap a r t i c l e s w a r m o p t i m i z a t i o n d e d i c a t e df o r m u l t i p l es e q u e n c ea l i g n m e n t ) a l s ot h et w oa l g o r i t h m sb a s e d o ns po b j e c t i v ef u n c t i o na n d c o f f e eo b j e c t i v ef u n c t i o na r ei m p l e m e n t e d ，a n db o t hc o m p l e x i t i e so ft w oa l g o r i t h m sa l e o n l yr e l a t e dt oe v o l u t i o ng e n e r a t i o n a ln u m b e ra n dt h es i z eo fp o p u l a t i o n f h a a l l y ，t h et w o a l g o r i t h m sa l eu s e dt o l e s tb e n c h m a r km u l t i p l es e q u e n c ea l i g n m e n t sd a t a b a s eb a l i b a s e t h er e s u l t ss h o wt h a tt h ep r o p o s e dm s a g a a n dm s a p s o a l g o r i t h m a r ef e a s i b l et os o l v et h e p r o b l e m o f s e q u e n c ea l i g n m e n t k e y w o r d s ：b i o i n f o r m a t i c sm u l t i p l es e q u e n c ea l i g n m e n t s g e n e t i ca l g o r i t h m p a r t i c l es w a r mo p t i m i z a t i o n0 b j e e t i v ef u n c t i o n 创新性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知，除了文中特别加以标注和致谢中所罗列的内容以外，论文中 不包含其他人已经发表或撰写过的研究成果：也不包含为获得西安电子科技大学 或其它教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所 做的任何贡献均已在论文中做了明确的说明并表示了谢意。 申请学位论文与资料若有不实之处，本人承担一切相关责任。 本人签名日期2 叫几如 关于论文使用授权的说明 本人完全了解西安电子科技大学有关保留和使用学位论文的规定，即：研究 生在校攻读学位期间论文_ l ：作的知识产权单位属西安电子科技大学。本人保证毕 业离校后，发表论文或使用论文工作成果时署名单位仍然为西安电子科技大学。 学校有权保留送交论文的复印件，允许查阅和借阅论文：学校可以公布沦文的全 部或部分内容，可以允许采用影印、缩印或其它复制手段保存论义。( 保密的论 文在解密后遵守此规定) 本学位论文属于保密在一年解密后适用本授权书。 本人签名： 导师签名 邈丛立 重i ! 垄 日期趔盛! i 翌j 1 f | 湖加心，f f 扩 第一章绪论 第一章绪论 1 1 背景及意义 随着人类基因组以及其它物种基因组计划的启动和实旌，生命的奥秘不断地 被揭示，有关生物体和核酸、蛋白质的序列和结构数据的分子数据大量产生并呈 指数增长，分子生物学家正在依靠现代计算机和网络技术，实现对这些海量数据 的管理、使用及进一步挖掘( d a t am i n i n g ) ，以期从中认识新的生命现象、寻找新的 生命规律。2 0 世纪8 0 年代末，由分子生物学、计算机科学以及信息技术等学科 的交叉和结合产生了生物信息学( b i o i n f o r m a t i c s ) ，并迅速发展起来。它的形成和发 展对目前在全球范围内广泛开展的各物种的基因组学、后基因学、信息科学、计 算机与网络本身以及生物医药的新药开发等多个领域将产生重大影响，并将成为 2 l 世纪生命科学的基石。 生物信息学这一术语，广义上说可指利用信息技术管理和分析生物学数据。 这就意味着生物信息学所涉及的范围相当广泛，从人工智能、机器人一直到基因 组( g e n o m e ) 数据分析。就基因组数据分析这一角度来看，生物信息学主要是指核 酸和蛋白质序列数据的计算机处理和分析。 人类基因组计划( h u m a ng e n o m ep r o j e c t ，h g p ) 是由美国科学家r e n a t o d u l b e c c o 于1 9 8 6 年率先提出的，1 9 9 0 年由美国政府正式启动。它是人类认识自 身、揭开人体奥秘、奠定2 l 世纪医学和生物学发展基础的一项庞大工程，目前已 经发展成为全球性的重大国际合作项目。 随着人类基因组计划的实施，通过基因组测序，蛋白质序列测定结构解析等 实验，分子生物学家提供了大量的有关生物分子的原始数据，发展情况如表1 1 和图1 1 所示： 表l _ 1g e n b a n k 核酸序列数据库增长情况 年份碱基数序列年份碱基数序列 1 9 8 26 8 0 3 3 86 0 61 9 9 21 0 1 0 0 8 4 8 67 8 6 0 8 1 9 8 32 2 7 4 0 2 92 4 2 71 9 9 31 5 7 1 5 2 4 4 21 4 3 4 9 2 1 9 8 43 3 6 8 7 6 54 1 7 5 1 9 9 42 1 7 1 0 2 4 6 22 1 5 2 7 3 1 9 8 55 2 0 4 4 2 0 5 7 0 01 9 9 53 8 4 9 3 9 4 8 55 5 5 6 9 4 1 9 8 69 6 1 5 3 7 19 9 7 81 9 9 66 5 1 9 7 2 9 8 41 0 2 1 2 1 1 1 9 8 71 5 5 1 4 7 7 61 4 5 8 41 9 9 71 1 6 0 3 0 0 6 8 71 7 6 5 8 4 7 2 多序列比对优化方法研究 19 8 82 3 8 0 0 0 0 02 0 5 7 919 9 82 0 0 8 7 6l7 8 42 8 3 7 8 9 7 l9 8 93 4 7 6 2 5 8 52 8 7 9l1 9 9 93 8 4 1 1 6 3 0 1l4 8 6 4 5 7 0 1 9 9 04 9 1 7 9 2 8 53 9 5 3 32 0 0 0l l l 0 1 0 6 6 2 8 81 0 1 0 6 0 2 3 1 9 9 17 1 9 4 7 4 2 65 5 6 2 72 0 0 l1 4 3 9 6 8 8 3 0 6 41 3 6 0 2 2 6 2 引自2 0 0 1 年1 2 月发布的g e n b a n k 第1 2 7 版说明，到目前为止，g e n b a n k 中的d n a 碱基数目已经超过了1 5 8 亿5 千万，d n a 序列数目已经超过了1 5 0 0 万。 # b a s ep a i r si nb i l l i o n s1 9 8 2 - 2 0 0 0 图1 11 9 8 2 - 2 0 0 0 年g e n b a n k 数据的增长情况 我们需要利用现代计算技术对上面的原始数据进行收集、整理、管理以便于 检索使用。为了解释和理解这些数据，需要对数据进行比对、分析，建立计算模 型，进行仿真、预测与验证，因而出现了生物信息学。生物信息学是在生命科学 的研究中，以计算机为工具对生物信息进行储存、检索和分析的科学。它是当今 生命科学和自然科学的重大前沿领域之一，同时也将是2 l 世纪自然科学的核心领 域之一。 1 2 研究现状 目前，被称为“生物学史上划时代工程”的人类基因图谱测序已经基本完成， 就人类基因组来说，得到序列仅仅是第一步，下一步我们将进入更为艰巨和复杂 的“后基因时代( p o s t - g e n o m ee r a ) ”，即收集、整理、检索和分析序列中表达的蛋 白质结构与功能的信息。因此研究新的计算方法，从序列数据库中提取有用的生 物信息，已经成为当务之急。 寻找序列相似性的过程称为序列比对，序列比对是序列分析的基本手段。通 过序列比较，检测新测定序列与数据库已知结构和功能的序列之间的相似性关系， 从而以足够的可信度确定新序列的结构和功能信息。这一过程往往需要通过数据 第一章绪论 库搜索，找出具有相似性的同源序列，对于d n a 序列，可推测该未知序列可能属 于哪个基因家族，具有哪些生物学功能，而对蛋白质序列来说，有可能找到已知 三维结构的同源蛋白质，而推测其可能的空间结构与功能。 然而，就目前的数学和计算机科学的能力而言，对数据容量达到上十亿字节 的数据库进行生物计算仍然是一项很艰巨的任务。虽然最简单的序列比对可以被 简化成字符串匹配的算法，但是扩展的和多重的序列比对仍处于实验摸索中。 近期，生物信息学的目标有以下几个方面： ( 1 ) 大规模基因组测序中的信息分析：大规模测序是基因组研究的最基本任 务，它的每一个环节都与信息分析紧密相关。 ( 2 ) 新基因和新的单核苷酸多态性( s n p s ) 的发现与鉴定：使用基因组信息学的 方法通过超大规模计算是发现新基因的重要手段，可以说大部分新基因是靠理论 方法预测出来的。构建s n p s 及其相关数据库是基因组研究走向应用的重要步骤。 ( 3 ) 完整基因组的比较和分子进化研究：通过多种同元比较和分析方法，可以 预测出一个基因可能具有的功能。 ( 4 ) 大规模基因功能表达谱的分析：基因组测序工作完成以后，虽然弄清了核 酸序列，但不知道它们的功能如何，或者说它们是如何按照特定的时空进行表达 的，表达量有多少等。 ( 5 ) 生物大分子的结构模拟与药物设计：基因的鉴定仅是阐明了其一级结构， 然而作为生命本质的生物信息流所包含功能的实现，必然要经过空间重构( 蛋白质 的折叠) ，才能表现出生命的功能。反之，从己知功能的蛋白质结构出发，研究这 些蛋白质功能的分子基础及其变化对蛋白质的三维重构和功能的影响，从而为基 因疗法设计相应的蛋白质受体药物，这些是摆在生物医学科学家面前的紧迫任务。 1 3 论文主要工作及安排 本文主要对生物信息学的多序列比对算法进行研究。首先，在遗传算法和粒 子群优化算法基本理论基础上，提出了基于遗传算法的多序歹0 比对算法和基于粒 子群优化的多序列比对优化算法。然后对两种算法进行了具体实现，并通过实验 结果对算法进行评价。最后对两种算法进行比较并作具体的分析。 论文内容具体安排如下： 第一章主要介绍课题的背景、意义及现状，提出了作者的工作设想，并简要 说明了作者所做的主要工作。 第二章首先研究了进行序列分析所涉及到的一些问题、概念以及对序列比对 的影响，包括：空位罚分、替换矩阵和目标函数。然后详细阐述s p 目标函数和 4 多序列比对优化方法研究 c o f f e e 目标函数的原理及特点，并对其进行具体的实现。最后深入分析了序列 比对的双序列比对和多序列比对问题。 第三章主要研究序列比对的经典算法和近年来比较新的比对算法及其发展 动态。 第四章本章主要分析了遗传算法的原理及特点，接着提出使用遗传算法来求 解多序列比对问题，并具体实现了该算法。然后使用基准多序列比对库b a l i b a s e 中的用例对所提出的算法进行测试，最后通过实验结果对算法进行分析。 第五章本章介绍了粒子群优化算法的背景，分析了其原理以及特点。接着提 出基于粒子群优化算法的多序列比对算法，并具体实现了该算法。然后使用基准 多序列比对库b a l i b a s e 中的用例对所提出的算法进行测试，最后通过实验结果 对算法进行分析。 第六章主要对所提出的两种多序列比对优化算法进行比较和分析。 最后对全文进行总结并指出不足之处和今后的研究方向。 第二章序列分析基础 第二章序列分析基础 2 1概论 随着基因组计划的实施和被称为“生物学史上划时代工程”的人类基因图谱 测序的基本完成，大量序列数据伴随而来，从而使序列分析成为了计算机在生物 学中应用的热点。 进行讨论之前，先给出同源性( h o m o l o g y ) 和相似性的明确定义，所谓同源序列， 简单地说，是指从菜一共同祖先经趋异进化而形成的不同序列。而相似性的概念 ，有两层含意：( 1 ) 指那些折叠方式相似却没有明显的序列相似性的蛋白质。( 2 ) 蛋 白质中一组具有相同催化活性和空间构象( 折叠方式) 的氨基酸残基，但分子间 整体上的序列和结构却不具有相似性。是指序列比对过程中用来描述序列之间相 同相似d n a 碱基或氨基酸残基序列所占比例的高低。 对于不同的相似性程度，我们采用不同的序列分析方法，如图2 1 所示。 般来说，序列间的相似程度越低，序列分析结果所得的可靠性就越差，当相似性 低于某一界限时，就很难得出明确的结论。这一界限通常被称为序列相似性界限， 即序列相似性朦胧区( t w i l i g h tz o i l e ) 。通常这一界限为2 0 左右。 双序列比对 多序列比对图 图2 1 各种相似性程度下所用不同序列分析方法 在生物信息学中，对序列数据进行相似性比较即序列比对，是一种基本的信 息处理方法。它对于发现生物序列中的功能、结构和进化的信息具有非常重要的 意义2 1 。序列比对就是运用某种特定的数学模型或算法，找出两个或多个序列之间 的最大匹配碱基或残基数，比对结果体现了算法在多大程度上反映了序列之间的 丫 区 区同引引剥倒觥暗 6 多序列比对优化方法研究 相似性关系以及它们的生物学特征p 】。这个过程往往需要通过数据库搜索，找出具 有相似性的同源序列，对于d n a 序列，可推测该未知序列可能属于哪个基因家族， 具有哪些生物学功能。而对蛋白质序列来说，有可能找到已知三维结构的同源蛋 白质，而推测其可能的空间结构与功能。 d n a 序列由四种核苷酸组成，用“爿”，“丁”，“c ”，“g ”代表四种碱基，其复 杂度为4 。例如：“c c 4 粥c 翻g “r ”可代表一个简单的d n a 序列，这里没有考虑 不确定的碱基，我们常用用来表示不确定的碱基。蛋白质序列由2 0 种氨基酸组成， 由埘口c d e f g h 删p q r s t v w x y z 代表不同的残基，其中“x ”，“b ”，“z ”表 示3 个特殊的残基。“z ”表示某个不确定的残基。”b ”表示天冬胺或天冬氨酸，用 三字符表示为“a s x ”。“z ”表示谷氨酰胺或谷氨酸，用三字符表示为“g l x ”。其 复杂度为2 3 。例如：“b e g s s t t n m a b n n m ”代表一个简单的蛋白质序列。因此生 物序列比对可以看作字符串的比对。 下面各节对序列比对中常用的空位罚分与相似性计分矩阵机制以及对双序列 比对与多序列比对作了详细的介绍。 2 2 空位罚分与相似性计分矩阵 在进行序列比对时，我们需要通过插入空位( g a p ) 使其具有最佳匹配，即序列 之间所对应的相同残基最多。序列比对中确定空位的过程十分复杂。最简单的办 法是通过不加限制地插入空位的办法获得相同残基的最大匹配数。我们知道，空 位的引入，意味着两个序列之间残基的插入或删除。如果引入空位不加限制，所 得的比对结果即使分值较高，也缺乏生物学依据。因此，必须有一种机制，对空 位的引入加以限制。常用的方法就是空位罚分。 通过空位罚分机制实现序列比对过程中，只考虑了残基的同一性，即序列之 间完全相同的匹配残基数目。可以把这种只考虑残基同一性的矩阵理解为单一的 相似性矩阵，即相同残基的计分值为1 ，不同残基的计分值为0 。显然，这种单一 的相似性矩阵具有很大局限性。改进计分矩阵的表征特性，找出那些潜在的具有 生物学意义的最佳匹配，提高数据库搜索的灵敏度，而又不至于降低信噪比，是 序列比对算法的核心。相似性计分矩阵就是为解决以上问题而产生的。 2 2 1 空位罚分 为了能得到最佳的匹配，在序列比对过程中往往需要插入空位( g a p ) ，确定空 位的过程是非常复杂的。空位处理是针对序列进化过程中可能发生的插入和缺失 第二章序列分析基础 而设计的。我们要对引入空位加以限制，常用的方法就是空位罚分。插入和缺失 可能只涉及1 个或2 个残基，也可能是整个功能域( d o m a i n ) ，所以，在进行空位罚分 设计时必须反映这些情况。 引入空位罚分是为了补偿插入和缺失对序列相似性的影响。由于没有什么合 适的理论模型能很好地描述空位问题，因此空位罚分缺乏理论依据而更多的带有 主观特色。一般的处理方法是用两个罚分值，一个与空位设置( g a po p e n i n g ) 有关， 它是指序列比对时，在一个序列中插入一个空位，使序列之间有更好的匹配。另 一个与空位扩展( g a pe x t e n s i o n ) 有关，它是指在引入一个或几个空位，继续引入下 一个连续的空位，使序列之间有更好的匹配。任一空位的出现均处以空位设置罚 分，而任一空位的扩大必须处以空位扩展罚分。对于一个空位长度为k 的罚值w 。可 用如下公式( 2 1 ) 表示： 毗= 口+ b k式( 2 一1 ) 其中a 表示空位设置罚分，6 为空位扩展罚分，k 为空位长度。这两个参数值设 置的变化对匹配结果产生的影响如表2 1 所示。 衰2 1 空位设置罚分和空位扩展罚分对匹配结果的影响 空位设置罚分空位扩展罚分说明 极少插入或缺失：适用于非常相 大大 关蛋白质问的匹配： 少量大块插入；适用于整个功能 大小 插入的情况 大量小块插入：适用于亲缘关系 小小 较远的蛋白质同源性分析 还有一种更为理想的空位罚分【4 l 对上面的公式进行改进，改进的公式如式( 2 2 ) 所示： = a l o g ( k + 1 )式( 2 - 2 ) 经过多年的实验，已经确定了一个合适的空位罚分值。大多数比对程序对特 定的替换矩阵设定了空位罚分的默认值( d e f a u l t ) ，如果使用者希望使用不同的替换 矩阵，则原来的空位罚分设置不一定合适。如何设置罚分并无明确的理论可循， 但大的空位罚分配以很小的空位扩展罚分，这样的设置被普遍证实是最佳的设置 思路。 在实际的情况中，通过改变某些参数可以得到不同比对的结果，比如空位罚 分值和选择不同的替换矩阵，因此合理的调节参数，能够减少空位数目，得到较 好的比对的结果。 多序列比对优化方法研究 2 2 2 替换矩阵 序列比对算法中我们常常要使用替换矩阵( s u b s t i t u t i o nm a t r i c e s ) 。它的原理是 找到一个可以估计任何比对的某一统计数，使生物学关系匹配最显著的比对统计 数最大。替换矩阵包括了在比对中各种匹配方式如何打分的信息，故替换矩阵又 常被称为计分矩阵( s c o r i n gm a t r i c e s ) 。 为了找出潜在的具有生物学意义的最佳匹配，提高数据库搜索的灵敏度，而 又不至于降低信噪比，我们常使用的计分矩阵为相似性计分矩阵，它的构建是基 于远距离进化过程中观察到的残基替换率，并用不同的计分值表征不同残基之间 的相似程度。恰当的选择相似性矩阵，可以提高序列比对的敏感度，特别是序列 间完全相同的残基数比较少的情况下。没有哪个计分矩阵适用于任何情况，因此， 在实际进行序列比对时，应该选择各种不同的相似性矩阵进行多次比对，并对比 对结果进行分析比较，才能得到比较合理的结果。 虽然针对不同的研究目标和对象应该构建适宜的取代矩阵，但国际上常用两 种相似性计分矩阵：( 1 ) 突变数据矩阵( m u t a t i o nd a t a ，m d ) ，它的记分方法是基于蛋 白序列中单点可接受( p o i n ta c c e p t e dm u t a t i o n ，p a m ) 的概念，通过对蛋白质进化模 式的研究而建立的【5 】。1 个p a m 进化距离对应于残基发生突变的概率以每1 0 0 个残 基中有1 个可接受单点突变而取。有p a m 2 0 ，p a m 6 0 ，p a m l 2 0 ，p a m 2 5 0 等，比如 p a m 2 5 0 相似性分值相当于序列间仍然有2 0 相同残基，数字表示进化的距离，数 字越小，距离越近，反之距离越远。( 2 ) 残基片断替换矩阵。突变数据矩阵是基于 相似性较高( 通常为8 5 以上) 的序列比对产生的，对于那些进化距离较远的矩阵如 ( p a m2 5 0 ) 是从初始模型推算出来的而不是直接计算得到的，因此准确率受到一定 限制。由此产生了残基片断替换矩阵b l o s u m ( b l o c k ss u b s t i t u t i o nm a u i x ) 拍j 。它 是从蛋白质模块数据库b l o c l s 中找出一组替换矩阵，直接利用多序列比对分析亲 缘关系较远的蛋白质，而不是用相近的序列，用于解决序列的远距离相关。比如， 高于或等于8 0 相同残基的序列组成序列模块可用于产生b l o s u m s 0 矩阵等等。 大量的实验表明，b l o s u m 矩阵总体比p a m 矩阵更适合于生物学关系的分析和局 部相似性搜索。 2 3 1背景介绍 2 3 目标函数 双序列比对计分很容易，我们可以根据相似矩阵和空位罚分机制来得到，并 第二章序列分析基础 且能很好的反映在生物学上的意义。具体的方法是把所得的双序列比对的每一列 根据所选的相似性矩阵计分求和，然后根据空位罚分减去由于增加空位而带来的 代价。对于多序列比对不仅在计算上是困难的，而且在选择如何计分即目标函数 也是很困难的。 首先介绍目标函数( o f - o b j e c t i v e f u n c t i o n ) 的概念，目标函数是来考查多序列比 对结果好坏的一种度量标准。因为所有多序列方法都依赖于一个目标函数来说明 所得的多序列比对结果的好坏从而反映此方法的精确度和有效性。理想的情况是 希望目标函数所得的分值有更多的生物学的意义。 在进化算法中，目标函数不仅是评价一个结果( 个体) 的好坏( 适应度f i t n e s s ) 的标准，而且目标函数的值还反映这个比对的生物的信息并且预示序列的结构或 者比对序列之间存在的进化关系。理论上，如果被比对的每一列碱基有相同的进 化或在r n a 和蛋自质的三维折迭结构中起到相似的作用，这样的多序列比对结果 就是正确的。 目前，还没有一个公认的目标函数来评价多序列比对的结果【”，不过有一些常 用的目标函数，下面介绍两个广泛使用的目标函数：基于s p ( s u m o f - p a i r ) 计分的 目标函数和c o f f e e ( c o n s i s t e n c yb a s e d o b j e c t i v e f u n c t i o nf o r a l i g n m e n t e v a l u a t i o n ) 目标函数。 2 3 2s p 目标函数 在说明具体问题之前，我们需要定义k 个字符串的多序列的子集。 定义：假设一是k 个字符串多序列s ly 5 2 ，i ( 1 ，k ) ，然后我们定义4 ，是 从4 的所有行中，删掉所用列都是空位而得出的。我们称a ，是4 到，的投影，我 们用简单的符号4 砧代替4 “圳。 例如：s l = 一gctgataaggtcct s 2 = ggg tgtt - tggt ct s 3 = 一gct t 一一tggtcct s 4 = gggtgtt tggt ct 那么彳，为： s 2 = gggtgtttggt ct s 3 = 一gct t tggtcct 怎么来评价多序列比对结果的质量，或者说怎么给多序列计分。假设 c ：j _ r 是所有占的可能比对集合到一个实数的对应关系。如果一个最优比对是 n ，那么我们的目标是发现一个最优多序列比对爿+ ，c ( a + ) = a r g n i n 。c ( a ) 或者 1 0多序列比对优化方法研究 c ( a 。) = a r g m a x c ( 爿) ，依赖于是最大相似度还是最小距离。 多序列比对把认为是相似的碱基排列在同一列中，类似双序列比对，我们用 一种计分方式表示多序列比对，一个简单的方法，也是实用很广泛的一种的方法， 多序列比对中包含所有双序列得出的计分，即：多序列比对中的所有双序列比对 的计分总和，我们把它称为s p 目标函数引，定义如式( 2 3 ) 所示： k - ik c ( 一) = e ( a 。) j t tj ，i + l 式( 2 3 ) c ( a 。u ) 是双序列比对4 。，的计分值。 如果我们假设比对的每一列的分数都是独立的，一般经常是指双序列比对， 我们能使用一个式( 2 4 ) 改写式( 2 - 3 ) 的s p 目标函数： k - ik c 似) = d ( a i ，h i ，4 ，明) 式( 2 - 4 ) h = li - i 州“ 例如：为了简单说明问题，假设d 口，倒= - 2 ，a b ，d 似，0 = d 6 = - 1 ， a “”，其它情况d = 3 。 s l = 一gctg atata act s 2 = ggg tgattag ct s 3 = ag cggaac acc t s c o r o f c o l u m n ：- 49 - 1 一l 9911 19 699 = 4 3 根据习惯使用式( 2 3 ) 或者式( 2 4 ) 。实际解决问题的时候我们使用相似性计分 矩阵来得出相似性函数d ，并且增加空位罚分机制，具体程序实现的目标函数形 式如式( 2 5 ) 和式( 2 6 ) 所示： ， c o s t ( s ，0 ) = d ( s ，) + g n u m f 式( 2 - 5 ) - 1 月i - i c o s t ( a ) = c o s t ( s 。，s ，) 式( 2 - 6 ) 1 = 2j = l 其中d ( s 。，5 j k ) 为序列s 。和5 j 对应残基的替代分值，g 为空位罚分，n u m q 为序 列s ，和s j 中的自然空位数。如果s 的一个比对a 满足：c o s t ( a ) = m a x 。( c o s t ( a ) ) ， 则称一是一个最优比对。 程序实现如下： 1 f o r 出l = 1 ；s l = s e q u e n c en u m b e r ；s i + 叫，s e q u e n c e _ n u m b e r 为序列的个数+ 2 f o r ( s 2 = 1 ；s 2 s e q u n c e _ s 2 _ l e n g t h ) ) * m s a 用来存储比对结果+ 5 r e t u r no ： 第二章序列分析基础 6 s c o r e + = m a t r i x m s a i s i i m s a s 2 i * m a t r i x 为相似性矩阵 7 s c o r e 一= g a p o p e n 4g a p n u m b e r ； p 计算空位罚分+ ， 8 r e t u r ns c o r e ； 2 3 3c o f f e e 目标函数 基于s p 的目标函数是基于一个相似性替换矩阵，它的应用最广泛，对于每一 个可能的氨基酸替代需要一个替换矩阵，对于删除插入【9 】的情况使用空位罚分给 出其代价和给序列加权【10 1 。在这种情况下，一个最佳多序列比对被定义为：对于 替代和插删有最小的代价。这个计分方法的主要限制是依赖于通用的相似性替代 矩阵。通常基于大量比对的统计分析来建立。因此这不可能适合于每个序列。 为了弥补这个缺陷，出现了基于其它的计分方式的目标函数。 下面我们介绍c o f f e e 目标函数，这个目标函数反映了一个多序列比对和这 个序列对应的双序列比对的数据库的一致性的程度。c o f f e e 主要的优点是它的灵 活性。使用c o f f e e ，任何适合于操作双序列比对的方法都能延伸去操作多序列 比对。它能提高多序列比对的精确性，当序列之间相似度比较低的时候c o f f e e 优于其它的目标函数。 c o f f e e 目标函数有两部分组成： ( 1 ) 一组相关的双序列比对库，这个数据库由双序列比对产生，但至少要包括 定义这个多序列足够的信息。例如一个多序列比对包括个序列，这个 双序列比对数据库中要包括所有可能的( 2 朋2 个双序列比对结果，双 序列比对的方法我们可以任意选择。 ( 2 ) 计算多序列比对和库中双序列的一致性的c o f f e e 计分函数，假设一个 满意的双序列比对数据库已经建好了。基于c o f f e e 计分函数的目标函 数定义如式( 2 7 ) 所示： r 。1n r 一1n c o f f e e s c o r e = i s c o r e ( a i d ) l i 彬l e n ( a , j ) i 式( 2 7 ) l ，一l j i + l jl i - i - j + lj 一个多序列比对包括个序列s t s n ，a l ，为多序列中的双序列的投影， l e n ( a t 为双序列比对的序列的长度，s c o r e 4 是多序列比对中a q 和数据库中 4 。双序列比对结果的一致性分数。即：s c o r e 口为多序列比对中a ，和数据库 中4 u 双序列比对结果的匹配残基的数目。，为双序列比对相关的权值，最简单 的阱，等于两个比对序列的相似度的百分比。从公式( 2 7 ) 我们可以看出c o f f e e s c o r e 是在0 - 1 之间的数。 1 2多序列比对优化方法研究 c o f f e e 目标函数和s p 目标函数的主要的不同是：在c o f f e e 中没有额外的 空位罚分，因为这些信息都已经包含在双序列比对库中。 c o f f e e 目标函数的计分过程如图2 2 所示： 我们可以计算出来c o l o n ns c o r e = 8 0 + 5 0 + 2 5 “0 ，8 0 + 3 0 + 2 5 + 5 0 + 6 0 + 6 0 = 0 7 所有的列都如此计算然后计算得出a l i g n m e n t s c o r e 就是c o f f e e 分值。 图2 2c o f f e e 计分过程 第二章序列分析基础 程序实现如下： 1 f o r ( 2 l ；i s e q u e n c en u m b e r ；i + + )s e q u e n c e n u m b e r 为序列个数 2 f b r 仃= f + 1 ；， “2 ，u 2 - g3 ，z k 1 - ) t k 的一条路径。 如果把每一条边一 v 都赋一个值盯( 砧一 ，那么如此路径的分值就是 ：盯( “。斗d i + i ) ，即序列比对的分值。 我们考虑两个序列a = 口l 口2 a m 和b = b i b 2 一占，a 包含吖个碱基，b 包含 个碱基。a 和b 序列比对的有向图被记为 o 。o 口的顶点是( f 力碱基对， i 0 ，m 】，- ，【o ，n 】，实际g a ，b 包括卅1 行和+ 1 列个顶点，4 和口比对的路径就 是从( 0 ，0 ) 到( t 帅的条路径。 ，口的边由如下的边组成： ( 1 ) ( f 一1 ，d 一( ) ，f 1 ，m 】，j 0 ， ，表示为 a i ，】 ( 2 ) ( f ，一1 ) 斗( f ) ，f 【0 ，m 】，j 【1 ，n 】，表示为【- ，b 。 ( 3 ) a i ，一1 ) - 9 g 力，i i ，圳， i ，j ，表示为( 口。，6 ，】 峨寸 眦寸 p 叶叶o ，嘲【- ，g 】 - , m j 图3 1 序列a = g m 和口= m g m 比对有向图 ，b 比对的过程就是寻找一条从( o ，0 ) ( 图4 1 的左上角) 到( 2 ，3 ) ( 图4 1 的右下角) 的路径，最终可以通过这个路径得到如下的比对结果： 一gm m gm 我们从图3 1 可以看出从( o ，0 ) 至t j ( m ，j 】v ) 的路径有许多条，每一条这样的路径 决定序列4 和b 的个结果，并且序列4 和曰比对结果都被一条唯一这样的路径 决定，到底怎样从这些路径中选择到最优的比对结果呢，这就要借助计分函数和 空位罚分得到。 对于每一个结点的分值，有它的直接前驱节点和输入边的分值取得，如图3 2 所示。 p 2 0 多序列比对优化方法研究 ( f 乒1 )仃( 【- 6 ，】)( f ，力 图3 2 节点( j ，j 的输入边及其分值 d ( f j ) = m a x ( 盯( f - l ) + 盯( 【q ，- 】) ) ，( 盯( f - 1 j - 1 ) + 盯( 【a ，b j 】) ) ，( o - ( ，d + 盯( ，b ： ) ) ) ，计算完所有节点分值后，从节点( ，| ) 开始回溯就可得到比对结果 序列。 我们通过实际比对两条短序列来研究动态规划算法是如何工作的，给定两个 序列a = “a d a v l c d 阳”和b = “a d r t c d r o ”。由于我们可以看出它们有高的相 似度，所以选择p a m 2 5 0 替换矩阵( 如图3 3 ) ，我们把它存在一个m a t r i x 数组中。 盯( 【a ，b ，】) 由替换矩阵取得。空位罚分为l o ，即盯( 一，6 l 】产- 1 0 ，盯( 口，】) = 一1 0 。 s h o r tm a t r i x = p a m 2 5 0 hik l m np 图3 3p a m 2 5 0 替换矩阵 每个节点的分值计算程序实现如下： 1 s o ，o 】_ 0 2 f o r j = 1t on d o 3 s 0 j 】= s 0 1 】+ 盯( - ，b 。】) 4