（计算机软件与理论专业论文）生物序列比对算法并行性的研究.pdf

摘要 d n a 和氨基酸序列比对是现代分子生物学的一个重要工具。大多数序列比对工具使 用的算法是n e e d l e m a n 和w u n s c h 在1 9 7 0 年提出的动态规划算法。这一算法的时间复杂 度是o ( m ”1 ，m ，? 是两个序列的长度。 本文通过对现有的生物序列比对算法进行研究分析，提出了一种基于前缀计算的序 列比对并行算法。这种算法采用数据划分技术，利用高性能计算系统，例如并行计算机 或分布工作站网络，有效的解决了超长序列比较引起的空间和时间问题。 在此基础上对几种经典的串行算法进行了改进，包括线性空位罚分和仿射空位罚 分模型、全局比对和局部比对算法等。本文提出的并行算法在o ( m n p ) 时间内找到最优 比对，其中p 是处理器数。传统的最优的串行算法使用d ( 什h ) 空间，o ( m 竹) 时间。本文 提出的并行存储有效的算法保持时间最优，使用0 ( r e + n p ) 空间。最后还对算法实现过程 进行了介绍。 关键词：序列比对并行算法前缀计算数据划分 a b s t r a c t o n eo ft h ei m p o r t a n tt o o l sf o rm o d e mm o l e c u l eb i o l o g yi ss e q u e n c ea l i g n m e n tf o rd n a a n da m i n oa c i d t h eb a s i ca l g o r i t h mf o rm o s to fs e q u e n c ea l i g n m e n tt o o l si sc a l l e dd y n a m i c p r o g r a m m i n ga l g o r i t h m , p r e s e n t e db yn e e d l e m a na n dw u n s c hi n1 9 7 0 t h et i m ec o m p l e x i t y i so ( m n ) ，w h e r ema n dna l et h el e n g t h so f t h et w os e q u e n c e s b ym e a n so fr e s e a r c h i n ga n da n 蜘n gt h ep r e s e n tb i o s e q u e n c ea l i g n m e n ta l g o r i t h m s ，a n e wp a r a l l e ls e q u e n c ea l i g n m e n ta l g o r i t h mi sp r e s e n t e d ，b a s e do np r e f i xc o m p u t a t i o n u s i n g d a t ad i v i s i o nt e c h n i q u e s ，t h em e t h o ds o l v e st h ep r o b l e mo ft i m ea n ds p a c et h a ta r i s ei n p a i r w i s ec o m p a r i s o no fe s p e c i a l l yl o n gs e q u e n c e b a s e do nt h i sa l g o r i t h m ，w em a k ea ni m p r o v e m e n to nt h ec l a s s i c a ls e q u e n t i a la l g o r i t h m s ， s u c ha sl i n e a rg a pp e n a l t ya n da f f i n eg a pp e n a l t ym o d e l ，g l o b a ls e q u e n c ea l i g n m e n ta n dl o c a l s e q u e n c ea l i g n m e n t a l lt h ea l g o r i t h m sp r e s e n t e di nt h i sp a p e rf i n dt h eo p t i m a la l i g n m e n ti n o ( m n p ) ，w h e r epi st h en u m b e ro fp r o c e s s o r s t h es e q u e n t i a lo p t i m a la l g o r i t h ms o l v et h e p r o b l e mi nd ( 州+ ms p a c ea n do ( m n ) t i m e w ep r o p o s e das p a c e - s a v i n ga l g o r i t h mt h a tu s e o ( m + n p ) s p a c ea n d r u n si no p t i m a lt i m e k e y w o r d s ：s e q u e n c ea l i g n m e n t p a r a l l e la l g o r i t h mp r e f i xc o m p u t a t i o n d a t a d 押i s i o n 创新性声明 本人声明所呈交的论文是我个人在导师的指导下进行的研究工作及取得的研究成 果。尽我所知，除了文中特别加以标注和致谢中所罗列的内容以外，论文中不包含他人 已经发表或撰写过的研究成果；也不包含为获得西安电子科技大学或其它教育机构的学 位或证书而使用过的材料。与我一同工作的同志对本研究所儆任何贡献均己在论文中 做了明确的说明并表示了谢意。 申请学位论文与资料若有不实之处，本人承担一切相关责任。 本人签名：z 刍垒兰笙 日期 关于论文使用授权的说明 本人完全了解西安电子科技大学有关保留和使用学位论文的规定，即：研究生在校 攻读学位期间论文工作的知识产权单位属西安电子科技大学。本人保证毕业离校后，发 表论文或使用论文工作成果时署名单位仍然为西安电子科技大学。学校有权保留送交论 文的复印件，允许查阅和借阅论文：学校可以公布论文的全部或部分内容，可以允许采 用影印、缩印或其他复制手段保存论文。( 保密的论文在解密后遵守此规定) 本人签名：f 刍垒垒 导师签名：塑i ! 兰 日期! ! 堕j j ! ： 日期 第一章绪论 第一章绪论 1 1 生物信息学概述 近年来随着快速序列测定、基因重组、基因芯片、多维核磁共振等技术的应用， 生物学实验数据呈爆炸趋势增长，迫切需要一种强有力的工具来组织数据，以利 于对已知生物学知识的储存和进一步加工利用。与此同时以数据处理和分析为本 质的计算机科学技术和网络技术的发展使对大规模数据的贮存、处理和传输成为 可能，- f l 崭新的拥有巨大发展潜力的科学生物信息学悄然兴起并蓬勃发展 起来。 生物信息学( b i o i n f o r m a t i c s ) 是生物学与计算机科学以及应用数学等学科相互 交叉而形成的- - f 新兴学科。它通过对生物学实验数据的获取、加工、存储、检 索与分析，进而达到揭示数据所蕴含的生物学意义的目的。由于当前生物信息学 发展的主要推动力来自分子生物学，生物信息学的研究主要集中于核苷酸和氨基 酸序列的存储、分类、检索和分析等方面，所以目前生物信息学可以狭义地定义 为：将计算机科学和数学应用于生物大分子信息的获取、加工、存储、分类、检 索。 国外一直非常重视生物信息学的发展，各种专业研究机构和生物科技公司、制 药公司等与日俱增。而事实上，欧美等发达国家在生物信息方面已经有较长时间 的积累。早在1 9 7 9 年，美国洛斯阿拉莫斯国家实验室就建立起g e n b a n k 数据库， 欧洲分子生物学实验室1 9 8 2 年就已经提供核酸序列数据库e m b l 的服务，日本也 与1 9 8 4 年着手建立国家级的核酸序列数据库d d b j 并于1 9 8 7 年开始提供服务。 欧美各国及日本相继成立了生物信息数据中心，如美国1 9 8 8 年成立的国家生 物技术信息中心( n c b i ) 、国家基因组资源中心、英国的欧洲生物信息研究所 ( e b i ) 、日本的国家遗传学研究所等。以西欧各国为主的欧洲分子生物学网络组 织( e m b n e t ) ，是目前国际最大的分子生物信息研究、开发和服务机构，通过计 算机网络使英、德、法、瑞士等国的生物信息资源实现共享。 目前，绝大部分的核酸和蛋白质数据库由美国、欧洲和日本的3 家数据库系统 产生；他们共同组成了g e n b a n k e m b l d d b j 国际核酸序列数据库，每天交换数 据，同步更新。其它一些国荔：，如德国、法国、意大利、瑞士、澳大乖j 亚、丹麦 和以色列等，在分享网络共享资源的同时，也分别建有自己的生物信息学机构、 二级或更高级的具有各自特色的专业数据库以及自己的分析技术，服务于本国生 物( 医学) 研究和开发，有些服务也向全世界开放。 生物序列比对算法并行性的研究 我国在1 9 9 3 年启动了基因组相关的研究项目，近两年又在上海和北京相继成 立了国家人类基因组南、北两个中心。1 9 9 9 年7 月，我国在国际人类基因组注册， 承担了其中1 的铡序任务，此举标志着我国己掌握生命科学领域中最前沿的大片 段基因组测序技术，在结构基因组学中占了一席之地。在水稻基因组研究领域， 北京华大基因研究中心在今年四月份宣布完成了水稻基因组工作框架图，并于日 前完成了改进型框架图的绘制工作。关于水稻单条染色体序列测定，中国已经完 成了第四条染色体的精确测序工作，在这条染色体上共发现了4 6 5 8 个基因，精确 率达到9 9 9 9 。此外，还参与了部分疾病基因的研究。 2 0 0 0 年7 月3 日，由中国科学院上海生命科学研究院生物信息中心和国家人 类基因组南方研究中心承担开发的、由国家“8 6 3 ”计划生物技术领域资助的生物 信息学项目我国首家自主开发的d n a 序列公共数据库( b i o s i n od a t a b a s e ) 正式 启用，同时开始接受我国核酸序列的注矾登记。 此外，清华大学在基因调控及基因功能分析、蛋白质二级结构预测方面、天津 大学物理系和中科院理论物理所在相关算法方面、中科院生物物理所在基因组大 规模测序数据的组装和标识方面、北京大学化学学院物理化学研究所在蛋白质分 子设计方面、华大基因组研究中心( 中科院遗传所人类基因组研究中心) 在大规 模测序数据处理自动化流程体系及数据库系统建立方面均已展开相关研究。北京 大学己建立了e m b l 中国镜像数据库，将该数据库移植到中国本地，并提供部分 的检索服务；复旦大学遗传学研究所为克隆新基因而建立的一整套生物信息系统 也已初具规模：中科院上海生化所、生物物理所等单位在结构生物学和基因预测 研究方面也有相当的基础。但从全国总体上来看与国际水平差距很大。 1 2 生物序列比对算法的发展 生物的遗传信息和功能信息存储在d n a ，r n a 和蛋白质中。这些生物大分子 非常复杂，由成千上万甚至上百万的原子按化学方法结合在一起而组成，有复杂 的3 维原子结构。但是，它们都是由一种易于理解的化学物的确定的字母表装配 而成的聚合体链。d n a 由四种脱氧核糖核酸( 腺嘌呤，胸腺嘧啶，乌嘌呤和胞核 嘧啶，或者a 、t 、c 、g ) 。r n a 由四种核糖核酸( a 、u 、c 、g ) 组成。组成蛋 白质的氯基酸有2 0 种，可以分别用2 0 个字母表示，包括a 、r 、n 、d 、c 、q 、 e 、g 、h 、i 、l 、k 、m 、f 、p 、s 、t 、w 、y 、v 。因为它们是组成成分的线性 链，所以它们可以用符号序列表示。 这些大分子序列包含大量有关它们生物功能的信息。在一个d n a 链中，四种 脱氧核糖核酸可能以任意顺序出现它们出现的顺序决定了d n a 的生物功能。在 一个蛋白质中，氨基酸序列决定了它的三维原子结构和功能。这些大分子的序列 第一章绪论 分析在现代分子生物学中起重要作用。 七十年代以来，d n a 测序方法的飞速发展，极大地引发了序列信息量的扩增， 从而使可供比较的序列数量呈现爆炸式增长。分子生物学家应该意识到，将未知 序列同整个数据库中的已知序列进行比较分析已经成为他们手中一个强有力的研 究手段。生物学家对基因和蛋白质序列进行比较，就是从核酸以及氨基酸的层次 去分析序列的相同点和不同点，以期能够推测它们的结构、功能以及进化上的联 系。序列比较包括同一个序列内不同片段的比较以及两个或多个序列的比对。在 过去的三十年里，序列比对的各种算法已经发展得越来越迅速，也越来越成熟。 下面简要介绍序列比对算法的发展过程。 1 9 7 0 年，g i b b s 和m c i n t y r e 发表了单序列分析方法矩阵打点作图法，用于 寻找单条序列中的重复片断，从而推测其功能【1j 。 1 9 7 0 年，n e e d l e m a n 和w u n s c h 发表了广受重视的两个序列比对算法1 2 】。 1 9 7 5 年，h i r s c h b e r g 提出一种基于分治法的存储有效的序列比对算、法【3 1 。 1 9 7 8 年，d a y h o f f 等研究氨基酸序列的进化，建立p a m 替换矩阵【4 l 。 1 9 8 1 年，s m i t h 和w a t e r m a n 基于仿射空位罚分模型开发局部序列比对算法口】。 l9 8 5 年，p e a r s o n 和l i p m a n 开发快速比对软件包f a s t a ，提供d n a 在d n a 数据库、蛋白质在蛋白质数据库，以及翻译序列在d n a 和蛋白质数据库中的查找 【毛7 1 。 1 9 8 7 年，g r i b s k o v ，m c l a c h l a n 和e i s e n b e r g 提出用p r o f i l e 检测一种蛋白质是 否属于某个己知的蛋白质家族 1 9 9 0 年，n c b i 第一次发行了启发式序列比对软件包b l a s t ，提供d n a 在 d n a 数据库、蛋白质在蛋白质数据库，以及翻译序列在d n a 和蛋白质数据库中 的查找 9 1 。 1 9 9 1 年，h u a n g 和m i l l e r 开发出s i m 程序，找出两个序列( 均为d n a 或均为 蛋白质) 之间的k 个最好的互不相交的局部比对f l o 】。 1 9 9 2 年，h e n i k o f t s 和h e n i k o f t j g 在序列比对算法中引入之后被广泛使用 的b l o s u m 替换矩阵【。 1 9 9 5 年，e d d y 提出一种新方法用隐马尔可夫模型( h i d d e nm a r k o vm o d e l ， h m m ) 进行多序列比对【1 2 , 1 3 1 。 1 9 9 7 年，z h a n g ，p e a r s o n 和m i l l e r 提出了将d n a 序列与蛋白质序列进行比对 的算法f r a m e 1 4 】。 1 9 9 7 年，n c b i 对b l a s t 进行改进，可以处理带空位的g a p p e d b l a s t 和位 置特异性p s i b l a s t t l 5 】。 2 0 0 0 年，t o r b j o mr o g n e s 提出了一种新的灵敏性更高的并行启发式算法 p a r a l i g n ，用于查找同源序列i 】“。 生物序列比对算法并行性的研究 此外，还有很多改进的算法，它们综合各种原始算法的优点，解决包括序列比 对和数据库搜索等在内的各种问题。最近几年，人们开始研究全基因组之间的比 对算法旧9 1 。 1 _ 3 生物序列比对的现状及存在问题 目前常用的序列比对软件包有： 基于“近似字符串匹配”算法的c l e a n u p1 8 ，它能够确定从核苷酸序列数据库 中指定的任何一对序列之间的整体同源性，并自动从冗余数据库中生成一组纯化 的无冗余的核苷酸序列集合。无冗余序列有助于进行统计学分析和加快广泛性检 索核苷酸序列数据库的速度。 b i o f a c e t ( 即大规模序列比较软件包l a s s a p ) 是一种跨越多种u n i x 平台的新 颖而全面的序列比对软件包。它使用了目前所有主要的序列比对算法：b l a s t 、 f a s t a 、s m i t h - w a t e r m a n 算法、n e e d l e m a l l w u n s c h 算法、k - b e s t 比对方法、字符 串匹配算法( 主要针对冗余问题) 、模式匹配算法( 譬如搜索p r o s i t e 特征模式) 等。b i o f a c e t 中的所有算法都是基于两个序列比对，不同算法之间的优势共享。它 的优势包括：数据库内或数据库间比较；直接计算；序列翻译；结构化的计算结 果和强大的再分析能力( 支持3 种输出格式) ；并行计算和利用特殊硬件设备而使 性能加强；同时它提供的应用编程接口( a p i ) 允许用户植入任何其它基于两个序 列比对的算法。 蛋白质多序列编辑器p r o m s e d 2 是运行于w i n d o w s 平台的应用程序，它能自 动或手动完成d n a 和蛋白质的序列比对、编辑和分析。它能读入几种常见格式 ( n b r f p i r ，f a s t a ，m s f ，e m b l s w i s s - p r o t ，i n t e l l i g e n e t i c s 和c l u s t a l 等) 的 序列数据，用户界面友好。它手动分析时能用不同的颜色组表示氨基酸序列在突 变等性质上相似的位点，是一套方便小巧的工具程序。 现在全世界数目众多的测序仪每天都在产生海量的生物信息数据，截止2 0 0 2 年1 2 月，e m b l 核酸序列数据库中已经有将近2 8 0 亿条序列。要对这些数据进行 处理，必须利用高性能计算机以及计算机群集方式来支持运算，针对生物计算的 特点设计特殊的体系结构和并行算法。同时还需要改进现有的理论分析方法，发 展研究基因组完整信息结构和信息网络的研究方法，开展大计算量算法和软件的 并行化工作。就国内而言，还需要设计全新的、有自主知识产权的新算法，以及 相应的软件包。 第一章绪论 1 4 本文的主要工作 d n a 和氨基酸序列比对是现代分子生物学的一个重要工具。大多数序列比对 工具使用的算法是n e e d l e m a n 和w u n s c h 在1 9 7 0 年提出的动态规划算法1 2 l 。这一 算法的时间复杂度是o ( n 、 是序列长度。 本文中用并行前缀计算作为基本模块，为在并行计算机或网络工作站上解决各 种序列比较问题开辟一条新思路，包括如何用恒定空位罚分函数和仿射空位罚分 函数进行全局比较和局部比较。这种算法提供线性加速使用o ( n d o g n ) 个处理器 ( ms n ) ，其中一是较长的序列的长度。改进的并行的存储有效算法保持时间最优， 使用0 ( m + n p ) 空间，其中p 是使用的处理器数。全文中，并行算法的最优时间意 味着相对于著名的串行算法有线性加速。 论文章节是这样安排的： 第一章简要介绍生物信息学和序列比对的国内外发展现状，以及作者的主要研 究工作。 第二章详细分析了生物序列比对的基本算法以及一些改进的算法。 第三章提出基于并行前缀计算的序列比对算法。 第四章是并行前缀序列比对算法的实现过程中需要考虑的问题。 第五章是本文的工作总结，在简要回顾论文工作的基础上，对序列比对算法的 设计和改进提出进一步的工作。 6生物序列比对算法并行性的研究 第二章生物序列比对 2 1 生物学背景 一旦将基因组转化成序列，而且成功的预测出它的( 部分) 基因内容，就可 以开始分析这个基因组中的基因之间以及与别的基因组的基因之间的进化关系。 这两种情况都需要两个基因的核酸或蛋白质的比较方法。 在几百万年中基因通过突变逐渐发生变化。在不同物种的基因组中同一个 祖先的两个直向同源基因会独立进化。因此在它们的d n a 序列中相似之处会很 少。它们编码的蛋白质的氨基酸也会不同。在一个物种的基因组中，可以通过称 之为基因复制的过程中的变换或重组产生基因的新备份。在基因复制之后，两个 共生同源基因2 分别经过多次突变趋于不同。 如果将基因的核苷酸序列写作一个字符串，那么突变引起的变化有下面几种 情况： 插入在序列中插入一个或几个字符 删除从序列中删除一个或多个字符 替代用另一个字符替代某个字符 插入和删除互为逆过程：己知两个序列，如果在一个序列中插入一个( 或多 个) 字符得到另一个，那么同样从后一个序列中删除这个字符也会得到前一个。 核酸在序列中的删除或插入常常用一个空位( 一) 表示。一般很难确定序列比 对中的空位是由插入还是删除造成的。因此这种空位经常被简称为插删( i n d e l s ) 。 在基因组序列化并测出基因之后，主要任务就是查找与这些基因相似的其它 生物基因。进行序列比较的最广泛采用的方法就是序列比对( a l i g n m e n t ) ，这种方 法尽可能确切的反映出序列之间的相似和相异。 序列比对在计算生物学中使用非常广泛。 假设在模型生物体中发现了一个基因，它在细胞新陈代谢中起的作用是已知 的。人类基因组数据库包含数百万的序列段。很有可能在人体中存在类似的突变 基因。驭一个功能己知的基因序列，将它与人类基因组数据库中所有序列进行比 对【2 0 1 。如果一个序列与已知序列比对的非常好，那就可能有相似的功能，有可能 1 直向同源( o r t h o l o g o u s ) ：共同祖先的直接后代( 没有发生基因复制事件) 之间的同源基因称为直向同源a ：共生i 司源( p a r a l o g o u s ) ：两个物种a 和b 的i _ d 源基挑l ，分别是共同姐丸璀州组中由复制事件l 叮产生的不h 拷- l 的后代，这棱称为共生同源基闻。 第二章生物序列比对 7 同源( h o m o l o g i e s ) 。序列相似和序列同源是不同的概念，序列之间的相似程度是可 以量化的参数，而序列是否同源需要有进化事实的验证。一旦在几种不同物种中 发现很多同源基因，在所有物种之间的最优比对会有助于展示随着时间的过去物 种从共同祖先分离时基因的进化史。 序列比对也用于找到遗传疾病的原因。遗传学家可以将健康人的序列与患者 的序列进行比对。如果所有患者都有一个健康人没有的变异，很明显这个变异就 是疾病的起因。 最常见的比对是蛋白质序列之间或核酸序列之间的两两比对，通过比较两个 序列之间的相似区域和保守性位点，寻找二者可能的分子进化关系。进一步的比 对是将多个蛋白质或核酸同时进行比较，寻找这些有进化关系的序列之间共同的 保守区域、位点和p r o f i l e 从而探索导致它们产生共同功能的序列模式。此外， 还可以把蛋白质序列与核酸序列相比来探索核酸序列可能的表达框架：把蛋白质 序列与具有三维结构信息的蛋白质相比，从而获得蛋白质折叠类型的信息。 比对还是数据库搜索算法的基础，将查询序列与整个数据库的所有序列进行 比对从数据库中获得与其最相似序列的已有的数据，能最快速的获得有关查询 序列的大量有价值的参考信息，对于进一步分析其结构和功能都会有很大的帮助。 近年来随着生物信息学数据大量积累和生物学知识的整理，通过比对方法可以有 效地分析和预测一些新发现基因的功能。 2 2 序列比对的基本概念 度量两个序列相似之处有两种方法：相似性( s i m i l a r i t y ) 测量和距离( d i s t a n c e ) 测量【2 1 1 。两个字符串的相似性的概念来自d n a 由共同祖先进化而来，而距离的 概念则侧重d n a 在进化过程中发生的变异。 定义2 1如果给每个匹配赋权值，那么两个字符串的相似性就是这些权值 之和的最大值。 定义2 2如果给基因的每次变异都赋一个权值，那么两个字符串的距离是 将一个字符串变成另一个需要进行的那些变异的权值之和。 定义2 3一条生物序列s 中的字符由某个有限字符集合确定，s e 。 定义2 4 如果s 是一个序列，那么l 司表示s 中的字符个数，研f 】表示序列的 第i 个字符。 定义2 5 对于t y u ，盯0 ，y ) 称为计分函数，表示工和y 在进行比 较时的匹配得分。公式( 2 1 ) 中的计分函数是最简单的一种。 生物序列比对算法并行性的研究 盯g ，y ) = 2 ， x = y ，x ，y e 一1 ， 工x ，y ( 2 1 ) 2 x = 一o r y = 一 定义2 6 两个序列s 和r 的比对( 也称为联配，a l i g n m e n t ) 可以用f 和 r 中字符一一对齐表示，其中，r u ) ，且i s7 卜：p l ，将f 和r 中的。 字符去掉后所得的序列分别与s 和r 相同。 比对( a ) 得分可以用下式表示： f ，i s c o e ( 彳) = 艺仃 【j 1 丁【f d ( 2 2 ) j j l 如果使用上面的计分函数( 2 1 ) ，序列a t g a c c 和a g a a t c 的一种比对如下，得 分为3 。 a g aatc 22 2 2 21 2 定义2 7两个序列s 和r 的最优比对是指在这两个序列的所有比对中得分 最大的一个。 这种计分方法与l e v e n s h t e i n 在1 9 6 6 年提出的两个序列之间的编辑距离的概 念密切相关。从序列s 到序列t 的l e v e n s h t e i n 编辑距离指将s 变成丁要做的最少 的编辑操作，这里编辑操作是替代、插入和删除。这里不匹配对应替代，空位对 应插入或删除。这两个都是负分值。为了使比对的总分最大，比对中的不匹配和 空位应该最少。在生物学上，计分函数表示序列中发生的各种突变( 替代、插入 和删除) 的生物学概率。最优比对就是两个序列匹配字符最多的情况。 常用的表示方法如下： ( 日，口) 表示匹配( 从s 到丁没有变化) ， ( 口，) 表示删去字符( 在s 中) ， ( 口，6 ) 表示( t 中的) b 替代( s 中的) a ，这里口b ， ( 一，6 ) 表示插入字符b ( 在t 中) 。 s 与r 中的问题是对称的，因此在s 中的删除也可以看作是在丁中的插入， 反之亦然。 2 3 动态规划算法 为找出两个序列的最优比对，一种最简单的方法就是一一列举所有可能的比 对，然后找出得分最大的情况。但是这种方法显然是不可行的。晟初使用的是“点 第二章生物序列比对9 矩阵分析”法( d o t - m a t r i xa n a l y s i s ) ，也称之为点图法( d o t - p l o t ) ，是两个序列相 似性的一种图形化的直观表示【i j 2 。现在生物学家使用的太部分双序列比对工具 的基本算法是动态规划算法，而进行数据库查找时最常用的是启发式方法。 动态规划是解决多阶段决策最优化问题的一种方法。能够用动态规划解决的 问题，往往是最优化问题，且问题的最优解( 或特定解) 的局部往往是局部问题在 相应条件下的最优解，而且问题的晟优解与其子问题的最优解要有一定的关联， 要能建立递推关系。 动态规划的思想是这样的：如果一条路径终止于最优路径上的一点，那么这 条路径本身就是起点到这个中间点的最优路径，也就是说，任何一个终止于最优 路径上的一点的子路径必然就是终止于这一点的最优路径本身。这样，最优路径 就可以通过把各个最优的子路径连接而成。 假设有两条长度分别为所和一的序列s 1 小】和n 1 堋】。不失一般性，本文中 假设m h 。s 1 一司和r f l 胡分别称为序列s 和丁的前缀。构建一个大小为+ 1 ) x ( ，升1 ) 的矩阵 以其中m f ，j 】( o f m ，0 孑竹) 记录了比对s 1 一f 】和玎1 刀 的最优得分。那么，脚m ，川是s 和r 的最优比对得分。将一个串与空串进行比对 没有什么价值，我们只能在空串中插入空格。因此， 村h o 】- 0 ， ( 2 3 ) 枷q 力= o r c j ，钆) ( 2 4 ) k = l 聊，o 】= 出。，_ ，) ，。、 k = ll z ) ， = 埘o j + 盯k ，u ) 当比对非空串研1 力和玎1 月时，有三种可能性：1 ) s 【1 1 与一个空格匹配，2 ) z 们与一个空格匹配，3 ) 研司与丁们匹配。得到下面的递归公式： m p 一1 ，卅+ 盯0 【j l 一) ，d ，- 】= 所a x j p ，一1 】+ 盯c 一r l ，d 彳【f 一1 ，_ ，一1 】+ 仃( s f l 丁【，d ( 2 6 ) 这个算法计算( f + 1 ) ( ，+ 1 ) 的得分矩阵从在公式中的三路选择导致矩阵元素 之间的如下依赖关系： 黜域一1m k ，一1 】一 m 月 计算矩阵的顺序可以是按行( 每行从左到右) 、按列( 每列从上到下) ，或者 l o 生物序列比对算法并行性的研究 按对角线顺序填充。这里对角线表示表中的所有使得i h = ，是常数的数据项【f ，1 。在 矩阵的右下角是要求的最优比对得分。 这种方法可以用编辑图( e d i t g r a p h ) 来表示，要求的比对就是图中左上角起 点到右下角终点之间的最长路径。在图中有三种边，水平线、竖直线和对角线( 粗 细) 分别表示插入、删除和匹配不匹配每条边按照相应的得分赋一个权值， 图2 1 所示是两个序列“a t c t g a t ”和“t g c a t a ”比对的编辑图。 2 4 全局比对与局部比对 2 4 1 全局比对 定义2 8 对于两条序列s 和n 全局比对是指将整个s 与丁进行比对，而 局部比对是指将s 的子序列与丁的子序列进行比对。 n e e d l e m a n 和w u n s c h 在1 9 7 0 年提出的n e e d l e m a n - w u n s c h 算法中 2 】，最优比 对必然对每个序列都由始至终进行比较，就是说从搜索空间( 或得分矩阵) 的左 上角直至右下角。换句话说，这是一种全局比对，最优比对中包括了全部的最短 匹配序列。 实际上n e e d l e m a n - w u n s c h 算法由两步构成： 按照前面给出的初始条件( 2 - 3 ) 、( 2 - 4 ) 、( 2 - 5 ) 和递归公式( 2 - 6 ) ，计 算出两个矩阵的相似分值，存于一个矩阵m 中： 根据此矩阵，按照回溯的方法寻找最优的比对序列。 n e e d l e m a n w u n s c h 算法也称为全矩阵( f u l l - m a t r i x ) 算法。一般，沿着一个矩 阵的两维放置两个序列，从左上角开始，从左到右、从上到下计算每个矩阵元素。 两个序列的比对与穿过矩阵的路径( 从左上角到右下角) 相对应。在矩阵中水平 或垂直移动，表示在相应的序列中插入一个空位。沿对角线移动表示对应的两个 a t c t g a t 第二章生物序列比对 字符相匹配。从它的三个相邻元素中选择向上、左和左上的移动最好的一步，计 算每个元素的值。 当矩阵计算完之后，右下角的值是最优得分。如果从m t 力向在公式( 2 - 6 ) 中得到最大值的m i ，一l 】、m i i - 1 ，力或者m f - l ，j - 1 】之- - ：e l j 连接线，最优比对就可以 表示成矩阵m 中的一条路径，它从m 小，行 出发，在m o ，o 】终止。这样的一条路 径被称为最优路径。穿过从m m ，”】到m o ，0 】的一条最优路径，找回最优比对称为 回溯( t r a c e b a c k ) 过程。图2 2 是两个序列“c g g t t a ”和“c a g g a ”比对的得分 矩阵和回溯的最优路径的示例，采用的计分函数是公式( 2 1 ) 。 。 c gg tt8 蓦 c a g g a 图2 2 n e e d l e m a n w u n s e h 算法得分矩阵以及回溯路径不意 得到的最优比对是： ( 1 ) c-gg tta或者 ( 2 ) cgg tta ca gg 一 aca g g a 从比对结果看( 1 ) 的匹配字符有4 个，空位有3 个；( 2 ) 的匹配字符有3 个， 替换( 不匹配) 1 个，空位3 个。因此( 1 ) 是最优比对，而( 2 ) 不是。这是因为利用 得分矩阵计算最大分值得到的最优比对很大程度上与计分函数有关。 根据公式( 2 6 ) 矩阵的每个数据元素至少需要l o 次运算，一共有( 肌+ 1 ) ( 肿1 ) 个。因此，计算整个矩阵的时间复杂度为o ( m n ) 。要找到穿过矩阵的路径需要0 ( ”) 次运算，因此整体时间复杂度是o ( n 2 ) ，这比完全列举所有比对要好得多。 如果按行计算的话，计算某一行的时候只需要保存前一行的数据和同一行前 一列的数据，即d + 卅) 空间。为了找回比对，需要设置回溯的指针，因此，空 间复杂度是o ( m n l 。 2 4 2 局部比对 生物序列中有重要功能的片断往往比较保守，即变异的概率很低。序列的 其它部分可能有较高的变异概率在进化过程中会变得面目全非。如果片面强调 全局比对，可能会漏掉真正的同源序列。良好的局部比对往往会更有效的揭示同 源关系。设计局部比对的另外一个很明显的原因就是在比较一个拼接后的m r n a 1 2 生物序列比对算法并行性的研究 和它的基因序列时，每个外显予都应该进行局部比对。 由于局部比对是在a 的予序列和口的子序列之间的一个得分最高的比对，所 以常常又称为子序列匹配。与此对应，前面介绍的全局比对又常常称为序列匹配。 s m i t h w a t e r m a n 是一种局部比对算法，它是在n e e d l e m a n w u n s c h 算法的基础上 发展而来的【s 】。这种比对的路径不需要到达搜索图的尽头，只需要在内部开始和 终结。 己知两个序列s 1 m 】平口丌1 n 口0 ，y ) 表示两个符号x 和y 匹配得分，要找两 个序列的高分相似性片断对h s p s ( h i g h s c o r ep a i r s ) 。用矩阵m 表示两个序列的 比对得分，令m f ，o 】- 虹0 ，刀= o ，其中0 f i d 訇锄。设m i ，月是分别以研f 和7 们结尾的两个子串的最大相似得分，s 1 f 】和玎1 卅的最优比对只能是下面的 情况之一： 如果研f 与7 们匹配，则得分 虹f - l ，产l 】+ 岱【f 】，7 u ) ： 如果s 【f 与一个空格匹配，得分m i - l ，月+ 口( 研f 】，o ) ； 如果丁们与一个空格匹配，得分川f ，- l 】+ “- ，丁们) ； 最后，为了防止计算得分为负，用0 表示在研f 和丁们之前没有任何匹配。 那么研l 一司和1 刃( h ，j 对锄) 的最优比对分数是下面的最大值 m i ，月= m a x m i 一1 ，- ，一1 】+ 盯$ 【f 1 丁【，d m i - ! + 盯附：! ( 2 _ 7 ) m i ，_ ，一1 】+ 盯卜丁【，d 0 要找最大分值片断对，首先要找到矩阵m 中的最大值确定得到最大分值的 i 和，然后从此向后回溯，直到m 的某个元素为0 。这一过程找出最大分值片断， 同时产生相应的比对。要找到下一个最大分值片断，只要找到m 的( 与第一次回 溯无关的) 下一个最大值，然后回溯。 一 ct a taatccc 一： c t g t a t c 图2 3s m i t h w a l e n n a n 算法得分矩阵和回溯的最优路径示例 第二章生物序列比对 图2 3 显示了序列“c t a t a a t c c c ”和“c t g t a t c ”比对的得分矩阵和回溯的最优路 径。采用的计分函数是公式( 2 一1 ) 。 图2 3 所示的矩阵得到的最优局部比对是 etata atcc c或ctata atc c c ct g tatcct g tatc 一般说来，对于4 个碱基具有相同频率的随机长序列，郧【f ，丁们) 值的平均值 为零。与空位匹配的值应大于或等于匹配与不匹配权重的差值。 在使用恒定空位罚分的情况下，s m i t h w a t e r m a n 局部比对算法与前面介绍的 n e e d l e m a n w u n s c h 全局比对算法的时间复杂度和空间复杂度相同，即时间d ( r n n ) ， 空间o ( m n ) 。 应该意识到，寻优方法总是把最佳的比对方法表达出来，而不在意它是否具 有生物学意义，另一方面寻求局部比对时可能会发现若干个重要的比对，因此， 不能仅仅注意最佳的一个。一个标准的s m i t h w a t e r m a n 算法只会报告出最好的一 个比对，改进的算法会报告出第二和第三的比对方式，从而显示出功能结构域。 2 5 空位罚分和替换矩阵 2 5 1 空位罚分 空位处理是针对序列进化过程中可能发生的插入和缺失而设计的。插入和缺 失可能只涉及1 个或2 个残基，也可能是整个功能域。空位的引入不仅考虑到空 位数而且考虑到连续的空格的数量。一般，连续的插删操作的总罚分不同于它们 的罚分之和。 定义2 9 空位( g a p ) 是在给定比对的一个字符串中的最大连续的空格。空 位长度是指空位中的插删操作数，即空格数。 最简单的罚分方法是恒定空位罚分( c o n s t a n t g a p p e n a l t i e s ) ，在这种情况下， 单个空格是不罚分的，每个空位罚分为，与空位中空格数无关。如果用# g a p s 表示空位数，用# s p a c e 表示空位中的空格数盯o ，力( 墨y 是字符或空格) 是计分 函数，那么口o ，一) = “，) = o ，要找的比对是 r1 州刮出i ，日：) + # g a p s ， l j 其中，a7 和b 表示在a 和口中插入空位后的序列。 将恒定空位罚分模型的进行推广，对空位中的每个空格增加一个权值瞩。在 这种情况下，表示开始一个空位的成本，w ；表示将空位扩展一个空格的成本。 1 4生物序列比对算法并行性的研究 这就得到仿射空位罚分模型( a 卵n e g a p p e n a l t i e s ) 。它的得名是因为长度为k 的空 位的罚分是由仿射函数+ 七巩决定的。因此，要找的比对就是 r 卅础 盯( 一? ，口7 ) + # g a p s + 睨# s p a c e 。 l ，j 恒定空位罚分模型就是仿射模型瞰= 0 的情况。当玎仁0 时，就得到另一种常 用的罚分模型线性空位罚分。这时= 女畎，空位罚分与空位长度成正比，要 找的比对是 r、 烈 c r ( a b i ) + w # s p a c e l j 有人建议，一些生物现象用下面的空位罚分函数模型更好：空位中的每个附 加的空格比它前面的空格罚分更少。换句话说，空位罚分是一个凸函数( c o n v e x ) ， 而不是它的长度的仿射函数。一个例子是函数+ l o g k ，其中k 是空位长度。 最后，可能考虑的最一般的空位罚分是任意空位罚分( a r b i t r a r yg a pc o s t ) ， 在这种情况下，空位罚分是长度女的任意函数w ( 曲。恒定、仿射和凸函数是任意 罚分模型的受限隋况。 使用仿射空位罚分的序列比对算法 前面介绍的算法使用的都是线性空位罚分，实际上，更常用的是仿射空位罚 分。它用两个参数设定空位罚分，一个与空位开放( g a po p e n i n g ) 有关，另一个与 空位扩展( g a pe x t e n s i o n ) 有关。任一空位的出现均处以空位设置罚分，而任一空位 的扩展必须处以空位扩展罚分。对于一个空位长度为七的罚分可用下式表示： w , = h + g k ，其中h 是空位设置罚分，g 为空位扩展罚分， g o 。 表2 1空位设置和空位扩展罚分对比对的影响 空位设置罚分空位扩展罚分说明 大大极少插入或缺失：适用于非常相关蛋白质之间比对 大小少量大块插入：用于整个功能域可能插入的情况 小小大量小块插入：用于亲缘关系较远蛋白质同源性分析 空位罚分的取值还与匹配得分和替代得分有关。如果空位罚分太小( 与替代 得分相比) 那么很可能在比对结果中根本不出现空位。如果空位罚分太大，那么 比对结果中可能会有很多空位从而产生高分比对。实际上，空位罚分的选择有 重要的经验因素在内，比如下面两个程序使用不同的罚分方法： p r o g r a m m a t c h m i s m a t c h g a p o p e ng a p - e x t e n d a l i g n5- 41 64 g c gb e s t f i t1 095 0 3 第二章生物序列比对 大多数比对程序均对特定的替换矩阵设定了空位罚分值的缺省值，如果用户 希望使用不同的替换矩阵，则原来的空位罚分值设定不一定合适。如何设定罚值 并无明确的理论可循，一般用大的空位设置罚分值配以很小的空位扩展罚分值被 普遍证实是最佳的设定思路。 两端空位不罚分的比对 另外还有一种隋况，就是两端空位不罚分的比对。在这种情况下，不管其它 空位的罚分是多少，在比对末尾或开头的插删操作罚分值是0 。例如， s = 一cacdb d