生物信息学讲义第四章双序列比对

1、第四章双序列比对本章介绍序列比对时常用的几个概念,包括同一性(identity、相似性(similarity同源性(homology等。这些概念,在蛋白质和核酸序列比对时经常使用。双序列比对(pairwise alignement是指通过一定的算法对两个DNA或蛋白质序列进行比对分析,从而找出两者之间最大相似性匹配。双序列比对是序列分析常用方法之一,是多序列比对和数据库搜索的基础。双序列比对基于一定的算法,而这些算法是基于序列本身的属性而不是关于该序列的注释信息。双序列比对基本上可以分为两类,一类是基于序列的局部相似性,而另一类则基于序列的全局相似性。两者之间无论是生物学基础,还是所用的算法,

2、都有很大区别。数据库检索(database query和数据库搜索(database search是分子生物信息学中两个不同的概念,有时常被混淆,有必要把这两个术语作简单说明。所谓数据库检索,是指对序列、结构以及各种二次数据库中的注释信息进行关键词匹配查找。例如,对核酸或蛋白质数据库输入关键词“human adrenergic recptors”(人的肾上腺素受体,即可找出数据库中所有和人的肾上腺素受体有关的序列条目。数据库查询有时也称数据库检索,它和因特网上通过搜索引擎查找需要的信息是一个概念。而数据库搜索是指通过特定的序列相似性比对算法,找出核酸或蛋白质序列数据库中与代查序列(query

3、sequence具有一定程度相似性的序列。例如,给定人的肾上腺素受体序列,通过数据库搜索,在核酸或蛋白质序列数据库中找出与该代查序列具有一定相似性的序列。显然,数据库检索和数据库搜索是在生物信息学中是两个完全不同的概念,它们所要解决的问题、所采用的方法和得到的结果均不相同。基于文字注释信息的数据库检索是核酸和蛋白质序列分析的一个重要组成部分,这一点不应予以忽略。然而,在生物信息学应用研究中,我们经常面临的问题是已经测得一个核酸序列片段,或一个由核酸序列翻译得到的蛋白质序列片段,而尚无注释信息。此时,数据库搜索就成了确定对该序列进一步分析研究的有效方法。本章将重点介绍序列比对的基本概念和常用算法

4、,特别是对代查序列和目标(subject序列之间的相似性关系,做一些定量的分析。为了识别一个新测定的序列和一个已知基因家族之间的进化关系,确定它们是否具有同源性,通常需要通过序列比对,找出它们之间核苷酸碱基或氨基残基的最大匹配,从而定量给出其相似性程度。如果两者的相似性程度很低,则很难确定它们是否具有同源性,除非使用亲源性分析(Phylogeny Analysis等其它分析方法,或有实验结果加以证实。通过对基因或蛋白质之间家族关系的分析,可以从浩繁的基因组信息中找出一些线索,从而对该基因或蛋白质家族的完整性(completeness进行预测。例如,假定某个基因家族的10个成员在鼠中已知,而在人

5、中只找到7个,那么很有可能还有3个成员在人的基因组中有待发现。这些分析结果还可用于药理学和分子生物学研究,用来解释某些受体对某种药物的特殊反应,尽管这些受体的序列在人基因组中还没有测得。它可以为分子生物学家以鼠的序列数据为模板克隆人的该基因的受体提供可能性。为了能够更有效地分析数据库搜索结果,有必要对序列比对的基本原理和数据库搜索的常用算法和有一个比较详细的介绍。4.3 字母表和复杂度核酸和蛋白质序列可以看成是由字母表中选出的字母组成的。一个字母表的复杂度定义为它所包含的不同字母的数目。例如英语字母表的复杂度是26;DNA序列的字母表复杂度是4。假定序列中含有不确定的碱基N(例如表达序列标签E

6、ST,则其复杂度为5。蛋白质序列由20种氨基酸残基组成,其复杂度为20。若某个残基未定,则用X表示。有时用B 表示天冬酰胺或天冬氨酸,若用三字符表示,则写成Asx;用Z表示谷氨酰胺或谷氨酸,若用三字符表示,则写成Glx。有些序列比对程序对这些特殊残基进行预处理。本章中将不涉及这些特殊残基,只讨论由20中普通氨基酸组成的蛋白质序列。在开始介绍序列比对的基本原理前,有必要分清算法和程序之间的区别。所谓算法,是指按照一定的方式描述计算过程或处理某个问题的一系列步骤,而程序则是算法的具体实现,也就是用某种计算机语言编写的实现某个算法的一在组指令集合。一个算法可能会有多种实现的方法。如果算法的描述或定义

7、明确,那么这些不同的实现方法,即不同的程序应给出同样的结果。然而,对某个算法可能有不同的理解,在具体实现时,可能会有一定的区别。首先,我们用一个简单的例子说明序列比对的基本原理。图6.1所示是对两个蛋白质序列片段进行比对的一般方法,基本思想是将两个序列上下排列,若上下对应的残基相同,则用竖线表示。可以通过插入空位(gap使上下两个序列具有最好的匹配,即两个序列之间对用所对应的相同残基最多。插入空位前序列1 (代查序列 AGGVLIIQVE|序列2 (目标序列 AGGVLIQVG插入空位后序列1 (代查序列 AGGVLIIQVE| |序列2 (目标序列 AGGVLI-QVG图6.1 利用插入空位

8、的方法获得最佳序列匹配图中竖线表示相同残基之间的匹配,插入空位前共有6个相同残基,插入空位后共有9个相同残基下一步,则可以计算相同残基个数,并用分数给出定量指标。如图6.1中未经比对以前的得分为6,而比对后的得分为9。显然,从这个例子中可以看出,匹配对准的相同残基数越多,两个序列之间相似性比对的得分就越高。当然,这只是一个用来说明比对原理的简单例子,序列很短,只有10来个残基,而大多数蛋白质序列的长度为200到500个残基,甚至更长。其次,这两个序列的长度几乎相等,而在实际情况下代查序列和目标序列的长度往往差别很大。此外,这两个序列的大部分残基相同,没有其它可选择的匹配方式。另一方面,序列比对

9、结果也可以根据引入空位的数目和非匹配残基的数目来度量。由此而引出距离矩阵的概念,即可以用距离矩阵的方式表示两个序列之间的相似性距离。序列比对所用的距离矩阵可能不止一个,同一算法的不同实现所用的距离矩阵可能会有所不同。为了进一步说明序列比对的基本原理,下面我们用一对比较接近实际的序列。假定序列A有400个残基,序列B有650个残基。如果序列A与序列B的一部分相同,则可以认为A是B的子序列。此时,只要在适当部位插入空位,就可以使序列A和序列B完全匹配。假定序列A中有两个片段分别与序列B中的两个片段相同。通过序列比对,可以找出这些相同片段,并在序列A中插入空位,使序列A与序列B有最大的匹配,如图6.

10、2(b所示。此时,可以找出序列A与序列B之间具有最好匹配的子序列。利用简单的渐进算法,可以直接实现上述比对,因为两个序列之间的相似片段比较明显。上述序列比对实例只是简要说明了如何找出两个序列之间完全相同的片段。实际操作中可以利用计算机程序实现上述序列比对的基本算法。然而,序列比对不仅需要考虑子序列之间的匹配,而且需要对整个序列进行比较。也就是说,必须考虑两个序列中所有残基的匹配。这就意味着,不可能使所有残基都能严格匹配。在这种情况下,比对过程中确定空位的过程变得十分复杂。最简单的办法使通过不加限制地插入空位的办法获得相同残基的最大匹配数。我们知道,空位的引入,意味着两个序列之间残基的插入或删除

11、。如果对引入空位不加限制,所得比对结果即使分值较高,也缺乏生物学依据。因此,必须有一种机制,对空位的引入加以限制。常用的方法就是空位罚分,即每插入一空位就在总分值中罚去一定分值,即加上一负分值,包括起始空位罚分和延伸空位罚分。所谓起始空位,是指序列比对时,在一个序列中插入一个空位,使两个序列之间有更好的匹配;所谓延伸空位,是指在引入一个或几个空位后,继续引入下一个连续的空位,使两个序列之间有更好的匹配。延伸空位罚分值可以与起始空位罚分值相同,也可以比起始空位罚分值小。因此,序列比对最终结果的分数值是两个序列之间匹配残基的总分值与空位罚分的总和。上述序列比对过程中,只考虑了残基的同一性,即两个序

12、列之间完全相同的匹配残基数目。可以把这种只考虑残基同一性的矩阵理解为一个分数值为1和0的分数矩阵(见表6.1,即相同残基的分数值为1,不同残基的分数值为0。这种矩阵通常称为稀疏矩阵,因为矩阵大多数单元的值为0。显然,这种单一的相似性分数矩阵具有很大局限性。改进分数矩阵的表征性能,找出那些潜在的具有生物学意义的最佳匹配,提高数据库搜索的灵敏度,而又不至于降低信噪比,是序列比对算法的核心。相似性分数矩阵就是为解决上述问题而产生的。相似性分数矩阵的构建,是基于远距离进化过程中观察到的残基替换率,并用不同的分数值表征不同残基之间相似性程度。恰当选择相似性分数矩阵,可以提高序列比对的敏感度,特别是两个序

13、列之间完全相同的残基数比较少的情况下。必须说明,相似性分数矩阵有其固有的噪声,因为它们在对两个具有一定相似性的不同残基赋予某个相似性分值时的同时,也引进了比对过程的噪声。这就意味着随着微弱信号的增强,随机匹配的可能性也会增大。本书不准备深入讨论有关相似性分数矩阵的问题,而只对两个常用的相似性分数矩阵作简单介绍,即突变数据矩阵和残基片段替换矩阵。突变数据矩阵(Mutation Data Matrix,简称MD,Dayhoff等,1978是基于单点可接受突变的概念,即Point Accepted Mutation,简称PAM。1个PAM的进化距离表示在100个残基中发生一个可以接受的残基突变的概率

14、。对应于一个更大进化距离间隔的突变概率矩阵,可以通过对原始矩阵进行一定的数学处理获得。例如,PAM250相似性分数矩阵相当于在两个序列之间具有20%的残基匹配。在序列比对中,通常希望使用能够反映一个氨基酸发生改变的概率与两个氨基酸随机出现的概率的比值的矩阵。这些比值可以用相关几率(relatedness odds矩阵表示。在序列比对过程中,两个序列从头到尾逐个残基进行比对,所得几率值的乘积就是整个比对的分值。在实际使用时,通常取几率值的对数以简化运算。因此,常用的突变数据矩阵PAM250实际上是几率值的对数矩阵(表6.2。矩阵中值大于0的元素所对应的两个残基之间发生突变的可能性较大,值小于0的

15、元素所对应的两个残基之间发生突变的可能性较小。表6.2 突变数据相似性分数矩阵PAM250C 12S 0 2T -2 1 3P -3 1 0 6A -2 1 1 1 2G -3 1 0 -1 1 5N -4 1 0 -1 0 0 2D -5 0 0 -1 0 1 2 4E -5 0 0 -1 0 0 1 3 4Q -5 -1 -1 0 0 -1 1 2 2 4H -3 -1 -1 0 -1 -2 2 1 1 3 6R -4 0 -1 0 -2 -3 0 -1 -1 1 2 6K -5 0 0 -1 -1 -2 1 0 0 1 0 3 5M -5 -2 -1 -2 -1 -3 -2 -3 -2

16、-1 -2 0 0 6I -2 -1 0 -2 -1 -3 -2 -2 -2 -2 -2 -2 -2 2 5L -6 -3 -2 -3 -2 -4 -3 -4 -3 -2 -2 -3 -3 4 2 6V -2 -1 0 -1 0 -1 -2 -2 -2 -2 -2 -2 -2 2 4 2 4F -4 -3 -3 -5 -4 -5 -2 -6 -5 -5 -2 -4 -5 0 1 2 -1 9Y 0 -3 -3 -5 -3 -7 -2 -4 -4 -4 0 -4 -4 -2 -1 -1 -2 7 10W -8 -2 -5 -6 -6 -7 -4 -7 -7 -5 -3 2 -3 -4 -5 -2

17、 -6 0 0 17C S T P A G NDE Q H R K M I L VF Y W* 表中把理化性质相似的氨基酸按组排列在一起,正值表示进化上的保守替代,值越大,保守性越大。序列分析的难点是要确定那些仅有20%相似性的序列之间是否具有同源关系。PAM250突变数据矩阵因此成为很多序列分析软件的缺省矩阵,因为它在20%的水平上反映出两个序列之间的相似性。按理说,使用与比对序列的实际进化距离更接近的相似性矩阵更为有效,但在实际使用中却无法实现,因为这意味着需要事先知道两个序列之间的进化距离,而导致先入为主的错误。因此,在实际进行序列比对时,应该选择各种不同的相似性分数矩阵进行多次比对,并

18、对比对结果进行分析比较,才能得到比较理想的结果。4.7.2 BLOSUM矩阵突变数据矩阵的产生基于相似性较高(通常为85%以上的序列比对,那些进化距离较远的矩阵(如PAM250是从初始模型中推算出来而不是直接计算得到的,其准确率受到一定限制。而序列分析的关键是检测进化距离较远的序列之间是否具有同源性,因此突变数据矩阵在实际使用时存在着一定的局限性。为了克服上述弊病,Henikoff夫妇(Henikoff和Henikoff,1992从蛋白质模块数据库BLOCKS中找出一组替换矩阵,用于解决序列的远距离相关。在构建矩阵过程中,通过设置最小相同残基数百分比将序列片段整合在一起,以避免由于同一个残基对

19、被重复计数而引入的任何潜在的偏差。在每一片段中,计算出每个残基位置的平均贡献,使得整个片段可以有效地被看作为单一序列。通过设置不同的百分比,产生了不同矩阵。由此,例如高于或等于80%相同的序列组成的串可用于产生BLOSUM80矩阵(BlOcks SUbstitution Matrix 发音为blossom;那些有62%或以上相同的串用于产生BLOSUM62矩阵,依此类推。序列比对实际上是根据特定的数学模型找出两个序列之间的最大匹配残基数。而序列比对的数学模型一般用来描述两个序列中每一个子字符串之间匹配的情况。通过改变某些参数可以得到不同的比对结果,例如空位罚分值大小。此外,序列长度差异和字母表

20、复杂度也会比对结果产生影响。合理地调节参数,会减少空位数目,得到较好的结果,而放宽对空位罚分的限制,理论上可以对任意两个序列进行比对而得到某个结果。因此,序列比对的结果并不能作为两者之间一定存在同源关系的依据。常用序列比对程序通常给出一些统计值,用来表示结果的可信度。BLAST程序中使用的统计值有概率p和期望值E。p值表示比对结果得到的分数值的可信度。一般说来,p值越接近于零,则比对结果的可信度越大;相反,p值越大,则比对结果来自随机匹配的可能性越大。期望值E描述的是搜索某一特定数据库时,随机出现的匹配序列数目。例如,E值为1可以解释为当前搜索中,由随机产生的相同分值的匹配的可能性为1。而E值

21、为0则表明搜索结果不大可能是随机产生的。点阵图(dotplot是用图示方法进行双序列比对的最基本方法。假定由两个序列A和B,它们的长度不一定相同,但最好相差不大。把序列A的残基序列沿X轴排列,序列B的残基序列沿Y轴排列,并构建一个矩阵。该矩阵所有元素的初始值均为0。对每个矩阵元素XiYj,赋予一个相似性值,表示该矩阵元素对应的两个残基之间的相似性程度,其中i的值为1到序列A的长度,j的值为1到和序列B的长度。如果只考虑同一性而不考虑相似性,则可以简单地将序列A和序列B相同残基所对应的矩阵元素的值置1,不同残基所对应的矩阵元素的值置0。若序列不是很长,上述矩阵很容易用可视化的图形方式表示,例如用

22、某个字符表示值为1的元素,如表6.3中的字符X。若序列较长时,则可以使用适当的图形程序(图6.6。若序列不是很长,上述矩阵很容易用可视化的图形方式表示,例如用某个字符表示值为1的元素,如表6.3中的字符X。若序列较长时,则可以使用适当的图形程序(图6.6。* M T F R D L L S V S F E G P R P D S S A G G S S A G G M *T *F * *R * *D * *L * *L * *S * * * *V *S * * * *F * *E *G * * *P *R * *P *D * *S * * * *S * * * *A *G * * * G *

23、* *若某一位置所对应的两个残基相同,则在该位置用“*”号表示。点阵图的最基本特征是有一些由噪音组成的随机点和具有信号特征的主对角线。其中主对角线由连续的点组成,它代表两条序列中具有相似性的区域。为清楚起见,表6.3中用X 表示。点阵图为从噪声背景中提取信号提供了一个图形表示法。尤其是在对于数据库搜索的结果的解释方面。点阵图中两条完全相同的序列被表示成一条连续的对角线,如图6.6(a所示。为了能够在一张图中表示出两个较长序列的比较结果,必须提高点阵图的分辨率。因此,表6.3中的X被替换成了点,连续的点连接在一起构成线。点阵图的名称由此而来。两条不完全相同却具有一定相似性的序列在点阵图上表示为一

24、些间断的对角线。其中不相连的区域表示那些不匹配的区域,如图6.6(b所示。那些相似性程度较低的序列所构成的点阵图具有较大的噪声。而一定长度的相似序列片段则在点阵图上以成组的、较短的对角线出现。它们平行于主对角线,其距离则反映了插入空位的多少。以上,我们讨论了如何利用单元矩阵来构建点阵图。更加复杂的点阵图可基于不同的打分规则而构建,这些打分规则规定了不同残基之间相似性程度的分值。例如,可以根据不同残基之间在进化、结构、理化性质等方面的相似性来规定它们之间的相似性分数值。在这种情况下,由于点阵图不只是简单的稀疏矩阵,那些偏离主对角线的点也可能具有相当重要的作用,所以过滤掉噪音就变得十分重要了。常用

25、的方法是引入滑动窗口算法作为平滑函数提高点阵图的信噪比。从上面的介绍中可以看出,序列比对是一个简单的数学模型,模型的参数可以加以调节。不同的模型所反映的生物学性质不同。例如,可以根据分子间的结构、功能和进化等方面的相关性来进行构建。必须指出,比对的结果是没有正确和错误之分的,其区别是由于模型所反映的生物学性质的不同。总体来说,比对模型可以分为两类,一类是考察两个序列之间的整体相似性,称全局性比对;另一类则着眼于序列中的某些特殊片段,比较这些片段之间的相似性,即局部性比对。区分这两类相似性和这两种不同的比对方法,对于正确选择比对方法是十分重要的。应该指出,在实际应用中,用整体比对方法企图找出只有

26、局部相似性的两个序列之间的关系,显然是徒劳的;而用局部比对得到的结果也不能说明这两个序列的三维结构或折叠方式一定相同。目前常用的BLAST和FastA等数据库搜索程序均采用局部相似性比对的方法,具有较快的运行速度,采用某些优化算法可进一步提高速度。局部相似性搜索主要用于找出序列中的功能位点,如酶的催化位点等。它们通常只有一个或几个残基,具有较高的保守性,并且不受序列中其它部分的插入和突变的影响。从这个意义上说,局部相似性搜索比整体相似性比对更加灵敏,也更具有生物学意义。需要指出的是,这并不能说明那些具有一定相似性的序列片段一定具有相同的三维结构。4.10 整体比对算法在对上述基本概念有所了解后

27、,我们开始讨论整体比对的Needleman和Wunsch算法(Needleman 和Wunsch,1970。从本质上讲,这一算法和已经广为使用的点阵图方法类似。整体比对方法中,两条蛋白质序列具有最多匹配残基定义为最佳匹配,其中允许进行必要的插入或缺失。为控制无限制的空位插入,我们引入罚分(penalty的概念。与点阵图类似,整体比对基于一个二维矩阵,并通过某种算法找出最佳匹配路径。矩阵的最基本形式是:将两序列中匹配残基所对应的单元的值置为1,不匹配的值置为0。然后对矩阵中的每个单元进行连续求和,即把能够到达该位置的所有单元中的最大值与该位置的值相加。若当前位置为第i行、第j列,那么能够达到它的

28、单元为:(a第i+1行中的第j 个单元之后的所有单元;(b第j+1列中的第i个单元之后的所有单元。对矩阵的所有单元都重复这一操作,直到全部结束为止。这样,可以构建一条最大匹配路径,它由N末端具有最大值的单元格开始,按照取最大值的原则一直到C末端,即从序列的起始开始到最后一个残基为止。不在主对角线上的单元格表示需要在此插入空位。在允许空位插入的情况下,可以籍此来寻求最大比对。假如不允许空位插入,则只能找一条分值较低的路径。下面我们通过两条短序列“ADLA VFALCDRYFQ”和“ADLGRTQNCDRYYQ”的比对详细讨论这一算法。根据算法,首先构建一个二维矩阵,用来表示两个序列的匹配状况。第

29、一个序列沿水平方向,即X轴;第二个序列沿垂直方向,即Y轴(表6.4。A D L G A V F A L C D R Y F QR000000000001000T000000000000000Q000000000000001N000000000000000C000000000100000R000000000001000Y000000000000100Y000000000000100Q000000000000001表6.4 Needleman和Wunsch算法初始矩阵接下来开始对矩阵单元连续求和。此处,从矩阵中最后一个单元开始。具体操作时,从最后一行、最后一列开始,先按列,后按行,逐步进行。本例中

30、最后一列为谷氨酰胺Q,而Y轴方向的序列有两个Q,按照相同匹配的分值为l的原则,该列的两个单元的分值为1,其它均为0。下一步统计倒数第2列,从该列的最后一行开始。该矩阵单元对应的残基不同,X轴方向为F,而Y轴方向为Q,因此,该单元的分值为0。从该列往上,下一个单元对应的残基为F和Y,也是一对非匹配残基,其本身的分值为0。这一单元的总体分值为1。因为根据连续求和原则,能够达到该单元的唯一的单元的值为1。从该单元往上,这一列中其它单元的分值均为1,这是因为该列残基F与Y轴方向序列的残基均不匹配,单元本身的分值均为0,而且可达到它们的单元中最大值(Q-Q单元提供为1。A D L G A V F A L

31、 C D R Y F QR000000005433110T000000005432110Q000000005432111N000000005432110C000000004532110R000000002223110Y000000002222210Y000000001111210Q000000000000001表6.5 Needleman和Wunsch算法中间矩阵接下来统计倒数第3列。该列最后一行为不同残基对Y和Q,分值为0。往上移动一个单元,到达倒数第二行,此时为一个相同匹配Y-Y,本身的分值为1。再加上能够到达该单元的最大值(Q-Q位置提供,所以该单元的值为2。如表6.5所示(表中标出了该

32、单元和能够达到该单元的所有单元。这样,算法逐列的向矩阵深入,并将子路径中最大值加到当前单元中。现在我们来看第7列,考察那个明显的L-L位置,其本身值为1,再加上两条子路径中的最大值5;这样,这个位置的值就是6(见表6.6。依次类推,直到矩阵中所有单元的分值全部计算完毕(表6.6。A D L G A V F A L C D R Y F QA976676675432110D786666665442110L667555556432110G555655555432110R555555555433110T555555555432110Q555555555432111N555555555432110C44

33、4444444532110D343333333342110R222222222223110Y222222222222210Y111111*Q000000000000001表4.6 Needleman和Wunsch算法最终矩阵在完成所有矩阵单元的分值计算后,接下来就是从最高分值单元开始找出最大分值路径,也就是找出最佳匹配。根据上述求和过程的特性,最大分值单元一定是在序列的N端,也就是矩阵的左上角。从这一起始单元回溯,找出具有最大分值的路径,即最佳路径。所谓回溯,就是由算法结束时的单元开始,反向向前逐一找到达到这一单元所经过的所有单元。本例中,最佳路径中间有一个间隔,可以通过在Y轴方向序列中的N和

34、C之间插入一个空位来实现。最终的比对结果如图6.7所示。ADLGA VFALCDRYFQ| | |ADLGRTQN-CDRYYQ图4.7 Needleman和Wunsch算法示例的比对结果可以看出,对用1和0表示匹配和非匹配的初试分数矩阵,上述连续求和得到的最大单元分值,即本例中矩阵起始单元的最大匹配值9,实际上就是最佳匹配路径中相同匹配残基的数目。从上述算法过程,可以看出Needleman和Wunsch算法考虑了两个序列中所有残基的贡献。其最佳路径的回溯一定从N端开始,而每个单元的分值计算则是从C端开始。因此,这种方法称整体性序列比对,其结果一定是两个序列的整体比对。Needleman-Wu

35、nsch算法适用于整体水平上相似性程度较高的两个序列。如果两个序列的亲缘关系较远,它们在整体上似乎不具有相似性,但在一些较小的区域上却可能存在局部相似性。1981年,Smith和Waterman提出了一种用来寻找并比较这些具有局部相似性区域的方法,即我们所熟知的Smith-Waterman算法。与Needleman-Wunsch算法类似,它也是一种基于矩阵的方法,而且也同样是运用回溯法建立允许空位插入的比对。多年来,Smith-Waterman算法一直是序列局部比对相关算法的基础,许多其它算法都是在这一算法的基础上开发和改进的。它也经常作为比较不同比对方法的标准。Smith和Waterman在

36、识别局部相似性时,确实具有很高的灵敏度,但使用时要注意,它只是寻找序列中一些小的,具有局部相似性的片段,而不是序列的整体相似性。*x A D L G A V F A L C D R Y F Qx0000000000000000A0D0L0G0R0T0Q0N0C0D0R0Y0Y0Q0表4.7 Smith-Waterman算法的起始矩阵Smith-Waterman算法一个重要特性是矩阵中每个单元均可以是比对结果序列片段的终点,该片段的相似性程度由该单元中的分数值表示。该算法的具体过程如下。首先,在矩阵的最上面一行和最左边一列的前面插入一个边界行和边界列,图中用字符“x”表示,我们把它称为第0行和第

37、0列。该边界行和边界列的所有单元的分值均为0.0(表6.7。可以把这些单元理解为长度为0的片段的起始端,它们的相似性分数值自然为0。至于为什么用小数而不是用整数,本质上没有区别。接下来,就是计算矩阵中每个单元的分数值。与Needleman-Wunsch算法不同,该算法在计算矩阵单元分值时,从左往右、从上到下,并沿对角线从左上角到右下角。即用三个函数分别计算由三条路径到达该单元的分数值并找出其中的最大值,假如数值小于0,则用0替代。这三个函数分别计算:(a当前单元对角线方向前一格的数值与相似性数值的和,匹配时值为1.0,不匹配时值为-0.333;(b当前行前面的各数值与相应的空位罚分值之和的最大

38、值;所用求空位罚分值的函数为W k=1.0+0.333k,k表示连续的第k个空位。(c当前列前面的各数值与相应的空位罚分值之和的最大值。如果出现负值就用0代替,表示到当前位置没有继续进行相似性比对的可能(表6.8。*x A D L G A V F A L C D R Y F Q表4.8 Smith-Waterman算法初对角矩阵一旦矩阵中所有单元的分值计算完毕,就可以找出具有最高分值的单元,也就是代表两个序列间高分匹配的终点。到达这个单元的其它矩阵元素可以通过回溯方法确定(表6.9。然后就根据回溯的路径求得一个片段的比对。如果需要,还可以找出在上述回溯范围以外的其它具有较高分值的矩阵单元,再进

39、行回溯,即找出多个具有较高分值的相似性片段。* x A D L G A V F A L C D R Y F Q表4.9 Smith-Waterman算法得到的最终矩阵从上面的介绍可以看出,Smith-Waterman算法中的分数矩阵中的最高分值可能不在N 端,它表示一个具有最高相似性的序列片段的终点,而序列中其它部分的相似性程度都比它低。从这一算法的本身可以看出,它是一个局部相似性比对而不是整体相似性比对的方法。4.12 动态规划算法以上我们介绍的算法属于动态规划算法的范畴。动态规划的基本思路是将复杂问题分解成与之相似的子问题,再通过求解各个子问题来得到整个问题的解决方案。正如上面所述,对于差

40、异程度不很大的两条序列,比对结果往往不是一种,而是可能有几种。利用动态规划方法求解这一问题,需要采用回溯(backtracking技术,并根据空位罚分等参数比较达到高分比对的不同路径。动态规划的核心,是要将各个最佳子路径连接起来,最终构成具有最优结果的路径,而得到比对结果。一个序列与整个数据库中所有序列比对的过程可以看作双序列比对的扩展。但是,要快速实现数据库的搜索并非易事。随着数据库规模的不断扩大,人们力图探索提高搜索效率的方法。无论是Needleman-Wunsch的算法,还是Smith-Waterman算法,对于数量不大的序列来说,其运行时间上可接受。对于大规模的数据库搜索,它们都非常耗

41、时,有时甚至无法实现。Smith-Waterman算法已经在一些专用计算机系统上实现,例如具有大规模并行运算功能的巨型机上运行的MPSrch程序。但是这些系统造价相当昂贵,并且随着计算机硬件的迅速发展和更新而被淘汰。显然,运行速度无疑是数据库搜索程序需要考虑的一个问题。就像上面两算法所提到的,数据库搜索程序的运行速度与代查序列长度及数据库容量密切相关。FastA和BLAST程序是目前最常用的基于局部相似性的数据库搜索程序,它们都基于查找完全匹配的短小序列片段,并将它们延伸得到较长的相似性匹配。它们的优势在于可以在普通的计算机系统上运行，而不必依赖计算机硬件系统而解决运行速度问题。 4.13.1

42、 FastA FastA 算法是由 Lipman 和 Pearson 于 1985 年发表的（Lipman 和 Pearson，1985） 。FastA 的基本思路是识别与代查序列相匹配的很短的序列片段，称为 k-tuple。蛋白质序列数据库 搜索时，短片段的长度一般是 1-2 个残基长；DNA 序列数据库搜索时，通常采用稍大点的 值，最多为 6 个碱基。通过比较两个序列中的短片段及其相对位置，可以构成一个动态规划 矩阵的对角线方向上的一些匹配片段。FastA 程序采用渐进（heuristic approach）算法将位于 同一对角线上相互接近的短片段连接起来。 也就是说， 通过不匹配的残基将

43、这些匹配残基片 段连接起来，以便得到较长的相似性片段。这就意味着，FastA 输出结果中允许出现不匹配 残基。这和下面将要介绍的 BLAST 程序中的成对片段类似。如果匹配区域很多，FastA 利 用动态规划算法在这些匹配区域间插入空位。 由 FastA 搜索产生的典型输出结果如图 6.8 所示。输出结果的第一行列出程序名称和版 本号，以及该程序发表的杂志。接下来列出所提交的序列，如本例中的人肾腺能-1A 受体， 数据库名称和版本号，此处为 SWISS-PROT 库第 35 版。然后是所用参数和运行时间，本例 为 12.42 秒。 紧跟这些一般信息的是数据库搜索结果。 首先列出搜索得到的目标序

44、列简单说明， 其数 目可由用户定义。为节省篇幅，这里只列出 10 个目标序列，实际运行时，可以列出更多的 序列，例如前 50 个。所列出的目标序列的信息包括：序列所在数据库名称的缩写，如这里 的 SW 表示 SWISS-PROT；目标序列的标识码、序列号和序列名等部分信息。括号中标明 匹配部分的残基数。 紧接着是由程序计算得到的初始化和优化后的分数值。 最后一列是期望 值即 E 值（详见 6.7.3 节） ，用来判断比对结果的置信度。接近于 0 的 E 值表明两序列的匹 配不大可能是由随机因素造成的。在这个例子中，前 10 个序列的 E 值都很低，也就是说它 们的置信度非常高。 在对搜索结果给

45、出概要信息后， 输出结果中详细列出代查序列与目标序列间的详细的比 对结果，其数目也可以由用户选择。为减少输出篇幅，有时在第一次运行时只输出概要信息 和部分详细比对结果。 比对结果给出代查序列和目标序列之间匹配残基的排列情况， 完全匹 配用冒号“:”表示，相似性匹配用句号“.”表示。残基的相似性定义取决与所用的相似性 分数矩阵，也可以由用户选择。本例中使用 BLOSUM 替换矩阵，在输出结果的概要信息部 分列出。 比对结果上方的数字表示代查序列的残基位置， 下方的数字表示目标序列的残基位 置，此处均不计插入的空位。 4.13.2 BLAST BLAST 是基本局部比对搜索工具（Basic Loc

46、al Alignment Search Tool）的缩写，由 Altschul 等于 1990 年提出。这一算法在许多计算机系统上的高效运行速度，特别是很早就 实现了 UNIX 系统下的并行化， 使其成为数据库搜索的流行软件， 为基于公共服务器的数据 库搜索带来了革命性的转变。 BLAST 算法本身很简单，它的基本要点是序列片段对（segment pair）的概念。所谓序 列片段对是指两个给定序列中的一对子序列， 它们的长度相等， 且可以形成无空位的完全匹 配。 BLAST 算法首先找出代查序列和目标序列间所有匹配程度超过一定阈值的序列片段对， 然后对具有一定长度的片段对根据给定的相似性阈值延

47、伸， 得到一定长度的相似性片段， 称 高分值片段对（high-scoring pairs, HSPs） 。这就是无空位的 BLAST 比对算法的基础，也是 BLAST 输出结果的特征。在图 6.9 中所显示的 BLAST 输出结果中共有 4 个高分值片段对。 与 FastA 类似，BLAST 的输出结果列出程序名称和版本号，以及文献出处。接下来列 出代查序列名称，本例中仍是人肾腺能-1A 受体；以及数据库名称，这里是处理过的非冗 余 SWISS-PROT 序列数据库。 紧接着就是搜索结果。与 FastA 一样，先列出目标序列，其数目可由用定义，这里仍只 列出 10 个。 目标序列信息包括数据库

48、缩写， 如这里的 sp 表示 SWISS-PROT； 目标序列编号、 标识码和序列名等部分信息； 接下来是目标序列中高分值片段对的最高得分， 高分值片段对 的数目由参数 N 确定。最后一列是概率值或称 p 值（详见 6.7.3 节） ，可用来判断比对的置 信度， 一个接近于 0 的 p 值表明两序列的匹配不大可能是由随机因素造成的。 在这个例子中， 前 10 个序列的 p 值都相当低，也就是说它们的置信度很高。 在搜索结果概要后， 就是依照用户要求所列出的代查序列和目标序列间无空位高分值片 段对的比对结果，和 FastA 一样，用户可以选择列出的序列多于或者少于概要中的序列数， 或者不列出比对

49、详细结果。 对于每一个高分值片段对， 都列出它们的起始和终止位点在序列 中的绝对位置。代查序列和目标序列之间的相同匹配用相应的氨基酸代码表示。此外，本例 中使用了 SEG 程序，用于屏蔽低复杂度的区域，即那些具有较高频度的特定残基或残基组 合的区域，如 GAPGAPGAPGAP 这样的重复序列。它们通常是胶原蛋白那样具有特殊结构 的序列，如不加以屏蔽，会导致假的高分值匹配的出现。被屏蔽的区域在搜索结果中用字符 X 表示。 允许空位的 BLAST 如上所述，最初的 BLAST 程序只能用于无空位的比对。经验表明比对结果通常会出现 一些无空位但不连续的区域，如图 6.9 中前 10 个目标序列的高

50、分值片段对都不止一个。不 难想象， 这些高分值片段对可以通过一些相似性较低且有空位的片段连接起来， 组成了一些 更长的或许更具实际生物学意义的比对。 基于上述思路，BLAST 算法经过改进允许空位插入（Altshul 等，1997） 。为缩短对数 据库初始搜索的时间， 新的算法只找出一个最好的高分值片段， 并以此为基础运用动态规划 方法将这一片段向两端延伸， 最终产生的比对结果可能有空位插入。 由于免去了查找所有高 分值片段对的步骤，新的算法比原算法快 3 倍。对 BLAST 算法的进一步扩充，可以考虑双 序列比对和多序列比对的有效结合，下一章将详细讨论具体算法。 4.14 本章小结 按关键词

51、进行数据库查询和序列相似性搜索是数据库应用的两个主要方 面。代查序列和目标序列的相似性程度，是识别未知序列和已知基因家族之间进化 关系的一个常用方法。 一组经过严格定义的解决某个给定问题的计算过程称为算法， 程序则是算 法的具体实现。一个算法可以有多种不同的实现方法，其结果理论上应该完全相同， 但实际上却可能有所差别。 双序列比对最基本的方法是通过插入空位将相同残基上下对齐， 并统计匹 配残基的个数而得到相似性分数。 对于长度不等、相似性程度较低的序列之间的比对应该采用分数矩阵。分 数矩阵不但包括相同残基之间的相似性分数值，而且包括不同残基之间的相似性分 数值。为尽可能减少空位插入并达到最佳比

52、对结果，需要引入对空位的罚分。 只考虑残基完全匹配的矩阵称为稀疏矩阵，其大部分矩阵元素的值为 0， 因此识别能力较差。根据进化距离较远的残基替换率可以得到权重矩阵，这类矩阵 对不同残基之间的匹配赋予相似性分数值，但其在增加微弱信号的同时，也增加了 随机噪声。如何从高强度的噪声背景中提取较弱的信号，是序列分析中的一大难题。 Dayhoff 突变数据矩阵（PAM）基于可接受点突变概念。PAM250 矩阵进 化距离相当于两个序列之间的相似性为 20%， 即序列比对中的模糊区， 因此 PAM250 矩阵是序列比对时常用的缺省矩阵。 BLOSUM 矩阵是根据 BLOCKS 数据库中经过比对的序列模块的替换概率 构建的。与 Dayhoff 突变数据矩阵相比，BLOSUM 矩阵对识别进化距离较远的相似 性序列更加可靠，因为前者的替换率是从高相似性序列中推导得来的。 序列比对结果并不能用来证明两个序列之间一定相关。 序列比对结果通常 附有可信度统计值，如双序列比对时采用的概率值 p 或期望值 E。