第三章 序列比较

第三章序列比较

序列比较是生物信息学中最基本、最重要的操作，通过序列比对可以发现生物序列中的功能、结构和进化的信息。序列比较的根本任务是：通过比较生物分子序列，发现它们的相似性，找出序列之间共同的区域，同时辨别序列之间的差异。在分子生物学中，DNA或蛋白质的相似性是多方面的，可能是核酸或氨基酸序列的相似，可能是结构的相似，也可能是功能的相似。一个普遍的规律是序列决定结构，结构决定功能。研究序列相似性的目的之一是，通过相似的序列得到相似的结构或相似的功能。这种方法在大多数情况下是成功的，当然，也存在着这样的情况，即两条序列几乎没有相似之处，但分子却折叠成相同的空间形状，并具有相同的功能。这里先不考虑空间结构或功能的相似性，仅研究序列的相似性。研究序列相似性的另一个目的是通过序列的相似性，判别序列之间的同源性，推测序列之间的进化关系。这里，将序列看成由基本字符组成的字符串，无论核酸序列还是蛋白质序列，都是特殊的字符串。本章着重介绍通用的序列比较方法。

序列的相似性

3.1

3.1序列的相似性

序列的相似性可以是定量的数值，也可以是定性的描述。相似度是一个数值，反映两条序列的相似程度。关于两条序列之间的关系，有许多名词，如相同、相似、同源、同功、直向同源、共生同源等。在进行序列比较时经常使用“同源”（homology）和“相似”（similarity）这两个概念，这是两个经常容易被混淆的不同概念。两条序列同源是指它们具有共同的祖先。在这个意义上，无所谓同源的程度，两条序列要么同源，要么不同源。而相似则是有程度的差别，如两条序列的相似程度达到30%或60%。一般来说，相似性很高的两条序列往往具有同源关系。但也有例外，即两条序列的相似性很高，但它们可能并不是同源序列，这两条序列的相似性可能是由随机因素所产生的，这在进化上称为“趋同”（convergence），这样一对序列可称为同功序列。直向同源（orthologous）序列是来自于不同的种属同源序列，而共生同源（paralogous）序列则是来自于同一种属的序列，它是由进化过程中的序列复制而产生的。

序列比较的基本操作是比对（align）。两条序列的比对（alignment）是指这两条序列中各个字符的一种一一对应关系，或字符对比排列。序列的比对是一种关于序列相似性的定性描述，它反映在什么部位两条序列相似，在什么部位两条序列存在差别。最优比对揭示两条序列的最大相似程度，指出序列之间的根本差异。

3.1.1字母表和序列

在生物分子信息处理过程中，将生物分子序列抽象为字符串，其中的字符取自特定的字母表。字母表是一组符号或字符，字母表中的元素组成序列。一些重要的字母表有：

（1）4字符DNA字母表{A,C,G,T}；

（2）扩展的遗传学字母表或IUPAC编码（见表2.3）；

（3）单字母氨基酸编码（见表2.1）；

（4）上述字母表形成的子集。

下面所讨论的内容独立于特定的字母表。首先规定一些特定的符号：

①Α—字母表；

②A*—由字母表A中字符所形成的一系列有限长度序列或字符串的集合；

③a、b、c—单独的字符；

④s、t、u、v、x—A*中的序列；

⑤|s|—序列s的长度。

为了说明序列s的子序列和s中单个字符，我们在s中各字符之间用数字标明分割边界。例如，设s=ACCACGTA，则s可表示为

0A1C2C3A4C5G6T7A8。

i:s:j指明第i位或第j位之间的子序列。当然，0￡i￡j￡|s|。子序列0:s:i称为前缀，即prefix(s,i)，而子序列i:s:|s|称为后缀suffix(s,|s|-i+1)。有两种特殊的情况，即i=j或i=j-1。

①i:s:i表示空序列

②(j-1):s:j表示s中的第j个字符，简记为sj。

一般认为，子序列与计算机算法中子串的概念相当。但是，严格地讲，子序列与子串的概念是有区别的：子串是子序列，而子序列不一定是子串。可以通过选取s中的某些字符（或删除s中的某些字符）而形成s的子序列，例如TTT是ATATAT的子序列。而s的子串则是由s中相继的字符所组成，例如TAC是AGTACA的子串，但不是TTGAC的子串。如果t是s的子串，则称s是t的超串。子串也可以称为连续子序列。

两条序列s和t的连接用s++t来表示，如：

ACC++CTA=ACCCTA

字符串操作除连接操作之外，另有一个k操作，即删除一个字符串两端的字符。其定义如下：

prefix(s,l)=sk|s|-l，

suffix(s,l)=k|s|-ls，

i:s:j=ki-1sk|s|-j。

序列比较可以分为四种基本情况，具体任务和应用说明如下：

（1）假设有两条长度相近的、来自同一个字母表的序列，它们之间非常相似，仅仅是有一些细微的差别，例如字符的插入、字符的删除和字符替换，要求找出这两条序列的差别。这种操作实际应用比较多，例如，有两个实验室同时测定某个基因的DNA序列，其结果可能不一样，需要通过序列比较来比较实验结果。

（2）假设有两条序列，要求判断是否有一条序列的前缀与另一条序列的后缀相似，如果是，则分别取出前缀和后缀。该操作常用于大规模DNA测序中序列片段的组装。

（3）假设有两条序列，要求判断其中的一条序列是否是另一条序列的子序列。这种操作常用于搜索特定的序列模式。

（4）假设有两条序列，要求判断这两条序列中是否有非常相似的子序列。这种操作可用于分析保守序列。

当然，进行序列比较时，往往还需要说明是采取全局比较，还是采取局部比较。全局比较是比较两条完整的序列，而局部比较是找出最大相似的子序列。

3.1.2编辑距离（Edit Distance）

观察这样两条DNA序列：GCATGACGAATCAG和TA TGACAAACAGC。一眼看上去，这两条序列并没有什么相似之处，然而如果将第二条序列错移一位，并对比排列起来以后，就可以发现它们的相似性。

如果进一步在第二条序列中加上一条短横线，就会发现原来这两条序列有更多的相似之处。

上面是两条序列相似性的一种定性表示方法，为了说明两条序列的相似程度，还需要定量计算。有两种方法可用于量化两条序列的相似程度：一为相似度，它是两条序列的函数，其值越大，表示两条序列越相似；与相似度对应的另一个概念是两条序列之间的距离，距离越大，则两条序列的相似度就越小。在大多数情况下，相似度和距离可以交互使用，并且距离越大，相似度越小，反之亦然。但一般而言，相似度使用得较多，并且灵活多变。

最简单的距离就是海明（Hamming）距离。对于两条长度相等的序列，海明距离等于对应位置字符不同的个数。例如，图3.1是3组序列海明距离的计算结果。

使用距离来计算不够灵活，这是因为序列可能具有不同的长度，两条序列中各位置上的字符并不一定是真正的对应关系。例如，在DNA复制的过程中，可能会发生像删除或插入一个碱基这样的错误，虽然两条序列的其他部分相同，但由于位置的移动导致海明距离的失真。就图3.1中例子最右边的情况，海明距离为6，简单地从海明距离来看，两条序列差别很大（整个序列的长度只有8bp），但是，如果从s中删除G，从t中删除T，则两条序列都成为ACACACA，这说明两条序列仅仅相差两个字符。实际上，在许多情况下，直接运用海明距离来衡量两条序列的相

似程度是不合理的。

为了解决字符插入和删除问题，引入字符“编辑操作”（Edit Operation）的概念，通过编辑操作将一个序列转化为一个新序列。用一个新的字符“-”代表空位（或空缺，Space），并定义下述字符编辑操作：

Match（a，a）—字符匹配；

Delete（a，-）—从第一条序列删除一个字符，或在第二条序列相应的位置插入空白字符；

Replace（a，b）—以第二条序列中的字符b替换第一条序列中的字符a，a1b；

Insert（-，b）—在第一条序列插入空位字符，或删除第二条序列中的对应字符b。

很显然，在比较两条序列s和t时，在s中的一个删除操作等价于在t中对应位置上的一个插入操作，反之亦然。需要注意的是，两个空位字符不能匹配，因为这样的操作没有意义。引入上述编辑操作后，重新计算两条序列的距离，就成为编辑距离。

以上的操作仅仅是关于序列的常用操作，在实际应用中还可以引入复杂的序列操作。下面是两条序列的一种比对：

上述比对不能反映两条序列的本质关系。但是，如果将第二条序列头尾倒置，可以发现两条序列惊人的相似：

再比如，下面两条序列有什么关系？如果将其中一条序列中的碱基替换为其互补碱基，就会发现其中的关系：

CTAGTCGAGGCAATCT

GAACAGCTTCGTTAGT

3.1.3通过点矩阵分析两条序列的相似之处

进行序列比较的一个简单的方法是“矩阵作图法”或“对角线作图”，这种方法是由Gibb首先提出的。将两条待比较的序列分别放在矩阵的两个轴上，一条在X轴上，从左到右，一条在Y 轴上，从下往上，如图3.2所示。当对应的行与列的序列字符匹配时，则在矩阵对应的位置作出“点”标记。逐个比较所有的字符对，最终形成点矩阵。

图3.2序列比较矩阵标记图

显然，如果两条序列完全相同，则在点矩阵主对角线的位置都有标记；如果两条序列存在相同的子串，则对于每一个相同的子串对，有一条与对角线平行的由标记点所组成的斜线，如图3.3中的斜线代表相同的子串“ATCC”；而对于两条互为反向的序列，则在反对角线方向上有标记点组成的斜线，如图3.4所示。

图3.3相同子串矩阵标记图

图3.4反向序列矩阵标记图

对于矩阵标记图中非重叠的与对角线平行斜线，可以组合起来，形成两条序列的一种比对。在两条子序列的中间可以插入符号“-”，表示插入空位字符。在这种对比之下分析两条序列的相似性，如图3.5所示。找两条序列的最佳比对（对应位置等同字符最多），实际上就是在矩阵标记图中找非重叠平行斜线最长的组合。

图3.5多个相同连续子序列矩阵标记图

除非已经知道待比较的序列非常相似，一般先用点矩阵方法比较，因为这种方法可以通过观察矩阵的对角线迅速发现可能的序列比对。

实例1：

实例2：

两条序列中有很多匹配的字符对，因而在点矩阵中会形成很多点标记。当对比较长的序列进行比较时，这样的点阵图很快会变得非常复杂和模糊。使用滑动窗口代替一次一个位点的比较是解决这个问题的有效方法。假设窗口大小为10，相似度阈值为8。首先，将X轴序列的第1-10个字符与Y轴序列的第1-10个字符进行比较。如果在第一次比较中，这10个字符中有8个或者8个以上相同，那么就在点阵空间（1,1）的位置画上点标记。然后窗口沿X轴向右移动一个字符

的位置，比较X轴序列的第2-11个字符与Y轴序列的第1-10个字符。不断重复这个过程，直到X 轴上所有长度为10的子串都与Y轴第1-10个字符组成的子串比较过为止。然后，将Y轴的窗口向上移动一个字符的位置，重复以上过程，直到两条序列中所有长度为10的子串都被两两比较过为止。基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且可以明确地指出两条序列间具有显著相似性的区域。

3.1.4序列的两两比对

序列的两两比对（Pairwise Sequence Alignment）就是对两条序列进行编辑操作，通过字符匹配和替换，或者插入和删除字符，使得两条序列达到一样的长度，并使两条序列中相同的字符尽可能地一一对应。设两条序列分别是s和t，在s或t中插入空位符号，使s和t达到一样的长度。图3.6是对序列AGCACACA和ACACACTA的两种比对结果以及对应的字符编辑操作。

下面就不同类型的编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表A中的任意字符a、b，定义：

这是一种简单的代价定义，在实际应用中还需使用更复杂的代价模型。一方面，可以改变各编辑操作的代价值，例如，在蛋白质序列比较时，用理化性质相近的氨基酸进行替换的代价应该比完全不同的氨基酸替换代价小；另一方面，也可以使用得分（score）函数来评价编辑操作。下面给出一种基本的得分函数：

在进行序列比对时，可根据实际情况选用代价函数或得分函数，即选用（3-1）式或（3-2）式。

下面给出在进行序列比对时常用的概念：

（1）两条序列s和t的比对的得分（或代价）等于将s转化为t所用的所有编辑操作的得分（或代价）总和；

（2）s和t的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；

（3）s和t的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。

使用前面代价函数w的定义，可以得到下列比对的代价：

s:AGCACAC?A

t:A?CACACTA

cost（s，t）=2

而使用得分函数p的定义，可以得到下列比对的得分：

s:AGCACAC?A

t:A?CACACTA

score(s，t)=5

进行序列比对的目的是寻找一个得分最高（或代价最小）的比对。

3.1.5用于序列相似性的打分矩阵（scoring matrix）

无论是3-1式还是3-2式，都是简单相似性评价模型，在计算比对的代价或得分时，对字符替换操作只进行统一的处理，没有考虑“同类字符”替换与“非同类字符”替换的差别。实际上，不同类型的字符替换，其代价或得分是不一样的，特别是对于蛋白质序列。某些氨基酸可以很容易地相互取代而不用改变它们的理化性质。例如，考虑这样两条蛋白质序列，其中一条在某一位置上是丙氨酸，如果该位点被替换成另一个较小且疏水的氨基酸，比如缬氨酸，那么对蛋白质功能的影响可能较小；如果被替换成较大且带电的残基，比如赖氨酸，那么对蛋白功能的影响可能就要比前者大。直观地讲，比较保守的替换比起较随机替换更可能维持蛋白质的功能，且更不容易被淘汰。因此，在为比对打分时，我们可能更倾向对丙氨酸与缬氨酸的比对位点多些奖励，而对于丙氨酸与那些大而带电氨基酸（比如赖氨酸）的比对位点则相反。理化性质相近的氨基酸残基之间替换的代价显然应该比理化性质相差甚远的氨基酸残基替换得分高，或者代价小。同样，保守的氨基酸替换得分应该高于非保守的氨基酸替换。这样的打分方法在比对非常相近的序列以及差异极大的序列时，会得出不同的分值。这就是提出打分矩阵（或者称为取代矩阵）的原由。在打分矩阵中，详细地列出各种字符替换的得分，从而使得计算序列之间的相似度更为合理。在比较蛋白质时，我们可以用打分矩阵来增强序列比对的敏感性。打分矩阵是序列比较的基础，选择不同的打分矩阵将得到不同的比较结果，而了解打分矩阵的理论依据将有助于在实际应用中选择合适的打分矩阵。以下介绍一些常用的打分矩阵或代价矩阵。

3.1.5.1核酸打分矩阵

设核酸序列所用的字母表为Α={A，C，G，T}。

（1）等价矩阵

等价矩阵（见表3.1）是最简单的一种打分矩阵，其中，相同核苷酸匹配的得分为“1”，而不同核苷酸的替换得分为“0”（没有得分）。

（2）BLAST矩阵

BLAST是目前最流行的核酸序列比较程序，表3.2是其打分矩阵。这也是一个非常简单的矩阵，如果被比的两个核苷酸相同，则得分为“+5”，反之得分为“-4”。

（3）转换-颠换矩阵

核酸的碱基按照环结构分为两类，一类是嘌呤（腺嘌呤A，鸟嘌呤G），它们有两个环；另一类是嘧啶（胞嘧啶C，胸腺嘧啶T），它们的碱基只有一个环。如果DNA碱基的变化（碱基替换）保持环数不变，则称为转换（transition），如A→G，C→T；如果环数发生变化，则称为颠换（transversion），如A→C，A→T等。在进化过程中，转换发生的频率远比颠换高，而表3.3所示的矩阵正好反映了这种情况，其中转换的得分为“-1”，而颠换的得分为“-5”。

3.1.5.2蛋白质打分矩阵

设蛋白质的字母表为表2.1（见第二章）。

（1）等价矩阵

（3-3）

其中，R ij代表打分矩阵元素，i、j分别代表字母表第i个和第j个字符。

（2）遗传密码矩阵GCM

GCM矩阵通过计算一个氨基酸残基转变到另一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。如果变化一个碱基，就可以使一个氨基酸的密码子改变为另一个氨基酸的密码子，则这两个氨基酸的替换代价为1；如果需要2个碱基的改变，则替换代价为2；以此类推（见表3.4）。注意，Met到Tyr的转变是仅有的密码子三个位置都发生变化的转换。在表3.4中，Glx代表Gly、Gln或Glu，而Asx则代表Asn或Asp，X代表任意氨基酸。GCM常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树，但是它在蛋白质序列比对尤其是相似程度很低的序列比对中很少被使用。

表3.4遗传密码矩阵

A S G L K V T P E D N I Q R F Y C H M W Z

B X Ala=A01122111112222222222222 Ser=S10112211221121111221222 Gly=G11022122112221221221222 Leu=L21202121222111122111222 Lys=K22220212121111222212122 Val=V12112022112122122212222 Thr=T11221201221121222212222

Pro=P11212210222211222122222 Glu-E12121122012212222222122 Asp=D12122122101222212122212 Asn=N21221212210122212122212 Ile=I21211112221021122212222 Gln=Q22211221122201222122122 Arg=R21111211222110221111222 Phe=F21212122222122011222222 Tyr=Y21222222211222101132212 Cys=C21122222222221110221222 His=H22212221211211212022212 Met=M22211112222121232202222 Trp=W21112222222221221220222 Glx=Z22221222122212222222122 Asx=B22222222211222212122212 22222222222222222222222 X

（3）疏水矩阵

该矩阵（见表3.5）是根据氨基酸残基替换前后疏水性的变化而得到得分矩阵。若一次氨基酸替换疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。

表3.5蛋白质疏水矩阵

R K D E B Z S N Q G X T H A C M P V L I Y F W Arg=R1010998866655555433333210 Lys=K1010998866655555433333210 Asp=D9910108876665555544433321 Glu=E9910108876665555544433321 Asx=B8888101088887777666555443 Glx=Z8888101088887777666555443 Ser=S667788101010109999887777664 Asn=N666688101010109999888777664 Gln=Q666688101010109999888777664 Gly=G556688101010109999888877665 ???=X555577999910101010998888775 Thr=T555577999910101010998888775 His=H555577999910101010999888775 Ala=A555577999910101010999888775 Cys=C4455668888999910109999885 Met=M334466888899991010101099887 Pro=P33446678888899910101099987 Val=V3344557778888891010101010987

Leu=L33335577778888999101010998 Ile=I33335577778888999101010998 Tyr=Y2233446666777788999910108 Phe=F1122446666777788889910109 Trp=W001133444555556777888910（4）PAM矩阵

为了得到打分矩阵，更常用的方法是统计自然界中各种氨基酸残基的相互替换率。如果两种特定的氨基酸之间替换发生得比较频繁，那么这一对氨基酸在打分矩阵中的互换得分就比较高。PAM矩阵就是这样一种打分矩阵。PAM矩阵是第一个广泛使用的最优矩阵，它是基于进化原理的，建立在进化的点接受突变模型PAM（Point Accepted Mutation）基础上，通过统计相似序列比对中的各种氨基酸替换发生率而得到该矩阵。Dayhoff和她的同事们研究了71个相关蛋白质家族的1572个突变，发现蛋白质家族中氨基酸的替换并不是随机的，由此，断言一些氨基酸的替换比其他替换更容易发生，其主要原因是这些替换不会对蛋白质的结构和功能产生太大的影响。如果氨基酸的替换是随机的，那么，每一种可能的取代频率仅仅取决于不同氨基酸出现的背景频率。然而，在相关蛋白中，存在取代频率大大地倾向于那些不影响蛋白质功能的取代，换句话说，这些点突变已经被进化所接受。这意味着，在进化历程上，相关的蛋白质在某些位置上可以出现不同的氨基酸。

一个PAM就是一个进化的变异单位，即1%的氨基酸改变。但是，这并不意味着经过100次PAM后，每个氨基酸都发生变化，因为其中一些位置可能会经过多次改变，甚至可能变回到原先的氨基酸。因此，另外一些氨基酸可能不发生改变。PAM有一系列的替换矩阵，每个矩阵用于比较具有特定进化距离的两条序列。例如，PAM?120矩阵用于比较相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值反映两条相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的概率。从理论上讲，PAM?0是一个单位矩阵，主对角线上的元素值为1，其它矩阵元素的值为0。其他PAM-N矩阵可以通过统计计算而得到。首先针对那些确信是相距一个PAM 单位的序列进行统计分析，得到PAM?1矩阵。PAM?1矩阵对角线上的元素值接近于1，而其它矩阵元素值接近于0。例如，可以按下述方法构建PAM?1矩阵。首先，构建一个序列间相似度很高（通常大于85％）的比对。接着，计算每个氨基酸j的相对突变率m j。相对突变率就是某种氨基酸被其它任意氨基酸替换的次数。比如，丙氨酸的相对突变率是通过计算丙氨酸与非丙氨酸残基比对的次数来得到。然后，针对每个氨基酸对i和j，计算氨基酸j被氨基酸i替换的次数。最后，将以上替换次数除以对应的相对替换率，利用每个氨基酸出现的频度对其进行标准化，并将以上计算结果取常用对数，于是得到了PAM-1矩阵中的元素PAM-1（i，j）。这种矩阵被称作对数几

率矩阵(log odds matrix)，因为其中的元素是根据每个氨基酸替换率的对数值来得到的。

将PAM-1自乘N次，可以得到矩阵PAM-N。虽然Dayhoff等人只发表了PAM?250，但潜在的突变数据可以外推至其他PAM值，产生一组矩阵。可以根据待比较序列的长度以及序列间的先验相似程度来选用特定的PAM矩阵，以发现最适合的序列比对。一般，在比较差异极大的序列时，通常在较高的PAM值处得到最佳结果，比如在PAM?200到PAM?250之间，而较低值的PAM矩阵一般用于高度相似的序列。实践中用得最多的且比较折衷的矩阵是PAM-250。

（5）BLOSUM矩阵

BLOSUM矩阵是由Henikoff首先提出的另一种氨基酸替换矩阵，它也是通过统计相似蛋白质序列的替换率而得到的。PAM矩阵是从蛋白质序列的全局比对结果推导出来的，而BLOSUM矩阵则是从蛋白质序列块（短序列）比对而推导出来的。但在评估氨基酸替换频率时，应用了不同的策略。基本数据来源于BLOCKS数据库，其中包括了局部多重比对（包含较远的相关序列，与在PAM中使用较近的相关序列相反）。虽然在这种情况下没有用进化模型，但它的优点在于可以通过直接观察而不是通过外推获得数据。同PAM模型一样，也有一系列的BLOSUM矩阵，可以根据亲缘关系的不同来选择不同的BLOSUM矩阵进行序列比较。然而，BLOSUM矩阵阶数的意义与PAM矩阵正好相反。低阶PAM矩阵适合用来比较亲缘较近的序列，而低阶BLOSUM矩阵更多是用来比较亲缘较远的序列。一般来说，BLOSUM-62矩阵适于用来比较大约具有62％相似度的序列，而BLOSUM-80矩阵更

适合于相似度为80%左右的序列。

第三章序列比较

两两比对算法

3.2

3.2两两比对算法

进行序列的两两比对最直接的方法就是生成两条序列所有可能的比对，分别计算得分（或代价）函数，然后挑选一个得分最高（或代价最小）的比对作为最终结果。但是，两条序列可能的比对数非常多，是序列长度的指数函数，随着序列长度的增长，计算量呈指数增长。从算法时间复杂性的角度来看，这种比对方法显然不合适。用上一节中所介绍的点矩阵分析方法，在寻找斜线及斜线组合时，仍然需要较大的运算量。因此，必须设计高效的算法以找出最优的比对。著名的Needleman-Wunsch算法，就是针对寻求最佳序列比对这一问题所设计的动态规划寻优策略。下面，首先介绍基本的序列比较算法，即动态规划算法（Dynamic Programming），该算法把一个问题分解成计算量合理的子问题，并使用这些子问题的结果来计算最终答案。这个算法是生物信息学的基本算法之一。

动态规划是一种常用的规划方法，往往用于在一个复杂的空间中寻找一条最优路径。对于

一个具体的问题，如果该问题可以被抽象为一个对应的图论问题，并且问题的解对应于图中从起点到终点的最短距离，那么就可以通过动态规划算法解决这个问题。在运用动态规划时，有以下几个要求：（1）首先，搜索问题能够划分成一系列相继的阶段；（2）起始阶段包含基本子问题的解；（3）在后续阶段中，能够按递归方式逐步计算前面阶段的每个局部解；（4）最后阶段包含全局解。

3.2.1序列两两比对基本算法

设序列s、t的长度分别为m和n。考虑两个前缀，0:s:i和0:t:j,i,j≥1。假如已知序列0:s:i和0:t:j 所有较短的连续子序列的最优比对，即已知：

（1）0:s:(i-1)和0:t:(j-1)的最优比对；

（2）0:s:(i-1)和0:t:j的最优比对；

（3）0:s:i和0:t:(j-1)的最优比对；

则0:s:i和0:t:j的最优比对一定是上述三种情况之一的扩展，即

（1）替换（s i，t j）或匹配（s i，t j)，这取决于s i是否等于t j；

（2）删除（s i，-）；

（3）插入（-，t j）。

令S(0:s:i，0:t:j)为序列0:s:i和与序列0:t:j比对的得分，可根据下列递归算式计算最大值：

（3-4）

其初值为：

（3-5）

按照这种方法，对于给定的打分函数p（s i,t j），两条序列所有前缀的比对得分值定义了一个(m+1)×(n+1)的得分矩阵

D=(d i,j)

其中，d i,j=S(0:s:i,0:t:j)。对于一个长度为n的序列，有n+1个前缀(包括一个空序列)，所以，得分矩阵的大小为(m+1)×(n+1)。其中，矩阵的纵轴方向自上而下对应于第一条序列（s），横轴方向从左到右对应于第二条序列（t）。矩阵横向移动表示在纵轴序列中加入一个空位，纵向的移动表示在横轴序列中加入一个空位，而斜对角向的移动表示两序列各自相应的字符进行比对。注意，各轴第一个元素的索引下标为0。

d i,j的计算公式如下：

（3-6）

d i,j最大值的三种选择决定了各矩阵元素之间的关系，用图3.7

表示。

矩阵右下角元素即为期望的结果：d m,n=S(0:s:m,0:t:n)=S(s,t)。

首先初始化得分矩阵D，然后计算D的其它元素。计算过程从d0,0开始，可以是按行计算，每行从左到右，也可以是按列计算，每列从上到下。当然，任何计算过程，只要满足在计算d i,j 时d i-1,j、d i-1,j-1和d i,j-1都已经被计算这个条件即可。在计算d i,j后，需要保存d i,j是从d i-1,j、d i-1,j-1或d i,j-1中的哪一个推进的，或保存计算的路径，以便于后续处理。上述计算过程到d m,n结束。

与计算过程相反，求最优路径或最优比对时，从d m,n开始，反向前推。假设在反推时到达d i,j，现在要根据保存的计算路径判断d i,j究竟是根据d i-1,j、d i-1,j-1和d i,j-1中的哪一个计算而得到的。找到这个点以后，再从此点出发，一直到d0,0为止。走过的这条路径就是最优路径（即得分最大路径），其对应于两条序列的最优比对。

假设我们采用公式（3-2）的打分函数，即序列字符匹配得分为+1，失配得分为0，空位罚分为-1，并且假设待比较的两条核苷酸序列分别是s=AGCACACA和t=ACACACTA。从公式（3-4）或公式（3-6）来看，动态规划算法的执行过程实际上是逐步计算得分矩阵D每一个元素的过程。首先按照公式（3-5）对矩阵D进行初始化，即首先计算矩阵第0行和第0列的元素值。其中，d0,0代表序列s和t两个空前缀比对的得分，其值显然为0。第0行的其它元素d0,j表示序列s空前缀与序列t前面连续j个字符组成的前缀比对的得分，相当于在序列s前面插入了j个空位，而每个空位的罚分为-1，所以d0,j=-j，即矩阵第0行元素值依次减1。同理，矩阵第0列元素值也应该依次减1。矩阵初始化的结果见图3.8。

t→

s ↓

A C A C A C T A

0-1-2-3-4-5-6-7-8 A-1

G-2

C-3

A-4

C-5

A-6

C-7

A-8

图3.8序列AGCACACA和ACACACTA比对的初始化得分矩阵初始化以后，假设逐行计算得分矩阵D的其它元素值。假设现在计算矩阵第1行、第1列的元素值，即计算矩阵（1,1）的元素值。可以通过三种途径到达该位置：

（1）从（0，0）出发，经过两条序列第一个字符的比对（A，A）；

（2）从（0，1）出发，在第二条序列加上的空位（A，-）；

（3）从（1，0）出发，在第一条序列加上的空位，（-，A）。

因此，矩阵（1,1）的元素值来自于下列三个值中最大的一个：

（1）左上方（0,0）位置的值加上匹配（A，A）的得分1，和为1；

（2）上方（0,1）位置的值加上空位罚分-1，和为-2；

（3）左边（1,0）位置的值加上空位罚分-1，和为-2。

最终，矩阵（1,1）的元素值等于1。当所有元素值计算完以后，形成图3.9所示的得分矩阵。

计算完得分矩阵D后，从元素（8，8）所在的位置反推最优路径。在图3.9中画出了一条穿越得分矩阵的路径，该路径表明如何通过合理的比对达到最大的得分。其中，斜线表示匹配或替换，垂直线表示删除，而水平线则表示插入。由该路径可以得到下面这种序列比对：

从图3.9可知，总的比对得分为5。值得注意的一点是，在有些情况下，最优比对并非唯一，亦

即存在几条最优路径。

以上计算是在打分函数的基础上进行的，得分值表示序列之间的相似程度。在实际应用中，也可以利用“代价”函数进行计算。两条序列比对的代价越低，序列就越相似；比对的代价越高，序列的差异就越大。因为计算方式刚好与得分函数相反，所以，具体计算时应求出最小代价所对应的路径。一般来讲，由于比对的得分可正可负，使用起来就更加灵活，所以大量的序列比对实用程序在计算时都采用得分的概念。

在计算序列相似程度时还应该考虑序列长度的影响。令S(s,t)表示两条长度分别为m和n 的序列的相似性得分，假设S(s,t)=99，两条序列有99个字符一致。如果m=n=100，则可以说这两条序列非常相似，几乎一样。但是，如果m=n=1000，则仅有10%的字符相同。所以，在实际序列比较时，使用相对长度的得分就更有意义，定义：

（3-7）

以sim（s，t）作为衡量序列相似性的指标。

从所使用的数据结构d i,j本身及其计算过程来看，序列两两比对基本算法的空间复杂度和时间复杂度都是O(mn)。如果两条序列的长度相近，空间和时间复杂度就为O(n2)。

3.2.2子序列与完整序列的比对

序列比较的基本动态规划算法找出两条序列的最佳比对，而不在乎它是否具有生物学意义。有些同源序列虽然全序列的相似性很小，但是存在高度相似的局部区域。如果在进行序列比对时注重序列的局部相似性，则可能会发现重要的比对。因此，不能仅仅只关注全局最佳的那一个。Smith和Waterman在Needleman-Wunsch算法的基础上进行改进，提出序列局部比对算法；后来其他人又进一步改进，形成改良Smith-Waterman算法，该算法将寻找多种最好的但不相互交叉的比对方式作为结果。下面分别介绍各种局部比对方法。

在有些情况下，我们需要将一个较短的序列（或探测序列，或模式序列）与一个较长的完整序列比较，试图找出局部的最优匹配。假设我们希望在较长序列ATGCAGCTGCTT中搜寻短序列AGCT。在所有可能的序列比对中，我们感兴趣的是：

这之所以是我们最感兴趣的比对，是因为它表明了，较短的序列完整地出现在较长的序列之中。我们有时希望避免对序列一端或两端出现的空位进行罚分，例如，在寻找一条短序列和整个基因组的最佳比对时就希望这样。