第三章 序列比较
3.3 序列多重比对
与序列两两比对不一样，序列多重比对（Multiple Alignment ）的目标是发现多条序列的共性。如果说序列两两比对主要用于建立两条序列的同源关系和推测它们的结构、功能，那么，同时比对一组序列对于研究分子结构、功能及进化关系更为有用。例如，某些在生物学上有重要意义的相似性只能通过将多个序列对比排列起来才能识别。同样，只有在多序列比对之后，才能发现与结构域或功能相关的保守序列片段。对于一系列同源蛋白质，人们希望研究隐含在蛋白质序列中的系统发育的关系，以便更好地理解这些蛋白质的进化。在实际研究中，生物学家并不是仅仅分析单个蛋白质，而是更着重于研究蛋白质之间的关系，研究一个家族中的相关蛋白质，研究相关蛋白质序列中的保守区域，进而分析蛋白质的结构和功能。序列两两比对往往不能满足这样的需要，难以发现多个序列的共性，必须同时比对多条同源序列。
图3.14是从多条免疫球蛋白序列中提取的8个片段的多重比对。这8个片段的多重比对揭示了保守的残基（一个是来自于二硫桥的半胱氨酸，另一个是色氨酸）、保守区域（特别是前4个片段末端的Q-PG ）和其他更复杂的模式，如1位和3位的疏水残基。实际上，多重序列比对在蛋白质结构的预测中非常有用。
多重比对也能用来推测各个序列的进化历史。从图3.14可以看出，前4条序列与后4条序列可能是从两个不同祖先演化而来，而这两个祖先又是由一个最原始的祖先演化得到。实际上，其中的4个片段是从免疫球蛋白的可变区域取出的，而另4个片段则从免疫球蛋白的恒定区域取出。当然，如果要详细研究进化关系，还必须取更长的序列进行比对分析。
对于多重序列比对的定义，实际上是两个序列的推广。设有k 个序列s 1, s 2, ... ,sk ，每个序列由同一个字母表中的字符组成，k 大于2；通过插入操作，使得各序列s 1, s2, ... ,sk 的长度一样，从而形成这些序列的多重比对。如果将各序列在垂直方向排列起来，则可以根据每一列观察各序列中字符的对应关系，如图3.14。
通过序列的多重比对，可以得到一个序列家族的序列特征。当给定一个新序列时，根据序列特征，可以判断这个序列是否属于该家族。对于多序列比对，现有的大多数算法都基于渐进比对的思想，在序列两两比对的基础上逐步优化多序列比对的结果。进行多序列比对后，可以对比对结果进行进一步处理，例如构建序列的特征模式，将序列聚类，构建分子进化树等。
3.3.1 SP模型
SP 模型（Sum-of-Pairs ，逐对加和）是一种多重序列比对的评价模型。在多重比对中，首先要对所得到的比对进行评价，以确定其优劣。例如，对图3.14中的8条序列进行比对，可以得到另外两种结果，如图3.15所示。那么，这样的三个多重比对，哪一个更好呢？这就需要有一种方法来评价一个多重比对。
评价一个多重序列比对比评价序列两两比对结果更复杂。这里，我们假设得分（代价）函数具有加和性，即多重比对的得分是各列得分总和。那么，我们首先考虑如何给比对的每一列打分，然后将各列的和加起来，成为一个总得分。
在处理每一列时，自然的处理方式是寻找一个具有k 个变量的打分函数（k 是参与多重比对的序列的个数），而每一个变量或者是一个来自特定字母表中的字符，或者是一个空位。我们很难得到这样一种具有k 个变量的表达式函数。另一方面，这种隐式函数不具有统一的形式，随着k 的变化，函数的表现形式也发生变化，不利于计算机处理。可以考虑使用显式函数，在具体实现显式函数时，用一个k 维数组来表示该显式函数（类似于打分矩阵），指定对应于k 个变量各种组合的函数值。这带来一个新问题，即所需的数组空间很大，而且随着k 的变化，数据结构也要随之动态变化。
我们所期望的函数在形式上应该简单，具有统一的形式，不随序列的个数而发生形式变化。根据得分函数的意义，函数值应独立于各参数的顺序，即与待比较的序列先后次序无关。另外，对相同的或相似字符的比对，奖励的得分值高，而对于不相关的字符比对或空白，则进行惩罚（得分为负值）。满足上述条件的一个函数就是常用的逐对加和SP 函数。SP 函数定义为一列中所有字符对得分之和：
(3-34)
其中，c1,c2,„,ck是一列中的k 个字符，p 是关于一对字符相似性的打分函数。对于p 可采用不同的定义。下面一个是SP 函数计算例子。
总得分根据字符两两比较得分计算演化而来，即逐对计算p(1,2)，p (1,3)，... ，p (1,8)，p (2,3)，p(2,4)，... ，p (2,8)，... ，p (7,8) 的所有得分，再加和得到结果。在上述计算中，我们假定：如果两个对比的字符相同，则得分为0，否则，得分为-1。所以，上述一列的得分为：（-7-6-5-4-3-2-1）+2= -26。
将一个多重比对所有列的得分全部加起来，其和即为该多重序列比对的得分。
注意：在进行多重序列比对时，可能会出现两个空位字符的比对，因此需要扩充函数p 的定义域，即增加p(-，-) 的定义。通常的定义是：
p(- , -) =
0 （3-35）
这似乎有点意外，因为一般只要是遇到空位字符，其得分就应该是负的，所以两个空位字符的比对应得到更多的负分。但是，这样处理是有道理的。在序列多重比对中，我们往往在得到整体比对的基础上进一步分析两条序列的对应关系。例如，根据图3.14的比对结果，取出最后两条比对的序列，见图3.16（a ）。这里存在空位字符对比的列，相当于这两条序列都进行了插入操作。但是由于插入位置相同[如图3.16（a ）箭头所指位置]，这两条序列本身在此位置上是完全相同的，所以，此位置上的编辑代价为零，或者得分为0。因而在分析这两条序列时，可以去掉这些空位字符对比的列，得到由多重序列比对结果推演的序列两两比对，如图3.16（b ）所示。其结果又称为多重序列比对在两条特定序列上的投影（projection ）。
若先处理每一个序列对，而在处理序列对时，逐个计算字符对，最后加和，则SP 得分模型的计算公式如下：
(3-36)
其中，α是一个多重比对，αi,j 是由α推演出来的序列s i 和s j的两两比对。
具体计算SP 时，我们可以先对多重序列比对a 的每一列进行计算，然后将每一列的得分值相加；也可以先计算每一对推演出来的两两序列比对的得分，然后再加和。但是注意，上述两种计算过程仅在p(- , -)=0这一条件成立下才等价。
3.3.2 多重比对的动态规划算法
多重序列比对的最终目标是通过处理得到一个得分最高（或代价最小）的序列对比排列，从而分析各序列之间的相似性和差异。如同处理序列两两比对一样，依然可以用动态规划算法。 回想前一节中穿越得分矩阵的路径（见图3.9）。得分矩阵相当于二维平面，而对于三条序列，每一种可能的比对可类似地用三维晶格中的一条路径表示，而每一维对应于一条序列，如图
3.17所示。路径的起点为晶格的左下角，路径的终点为晶格的右上后角。图3.17中的路径记为 (0,0,0)，(1,0,0)，(2,1,0)，(3,2,0)，(3,3,1)，(4,3,2)。如果存在多条序列，则形成的空间是超晶格（Hyperlattice ）。
假设在三维晶格的顶端、前面、后面各有一平行光源，这些光源将代表多重序列比对的路径投影到相对的侧面，即投影到一对序列所在的平面（图3.17仅画出了右边的光源）。将多重比对投影到两个序列所在的平面在加速优化计算中将起着重要的作用。
一个多重序列比对投影的结果可能比原来要短，例如，下面的多重序列比对
G---SNS
GN----S
GNAVSNS
如果在前两条序列方向进行投影，则投影结果为：
G-SNS
GN--S
如果路径的某一段沿投影方向进行，那么该段路径就不可能产生投影。图3.17中的第一段沿着第一条序列向前进，推进方向垂直于其他两条序列，则平行光源对这一段不产生直线投影，而仅仅产生一个投影点。
图3.17三条序列所对应的三维晶格
序列两两比对的动态规划算法经改进后可直接用于序列的多重比对。就二重序列而言，将动态规划算法的计算过程看成是在二维平面上按一定顺序访问每一个节点，访问节点的先后顺序取决于节点之间的关系，如图3.7所示。二维平面（或二维晶格）与前面所介绍的超晶格相似，上一节得分矩阵中的位置与每一个晶格节点相对应。在计算过程中，对每个节点逐一计算其得分（或代价）。每个节点的得分代表两个序列前缀的最优比对的得分，而晶格最后一个节点（右下角节点）的得分即为两条完整序列的比对得分。
在超晶格中，序列的放置与上一节有所不同，计算是从左下角进行的，而不是得分矩阵中的左上角。如图3.18所示，计算过程从超晶格的坐标点（0,0,„,0）开始，按节点之间的依赖关系向右上后方推进，直到计算完最后一个节点。实际计算时，可以采用类似二维的计算方式，即先低维后高维（对应于先行后列）。在图3.18中，当前点的得分计算取决于与它相邻的7条边，分别对应于匹配、替换或引入空位等三种编辑操作。计算各操作的得分（包括前趋节点的得分），选择一个得分最大的操作，并将得分和存放于该节点。在三维或超晶格中，计算公式与上一节相似，但计算一个节点依赖于更多的前趋节点。在三维情况下要考虑7个前趋节点，在k 维情况下要考虑2-1个前趋节点。
k
假设以k 维数组A 存放超晶格，则计算过程如下：
① a[ 0, 0, „ ,0 ] = 0
② a[ i ] = max {a[ i - b ] + SP-score[Column(s , i , b )]}
这里i 是一个向量，代表当前点，b 是具有k 个元素的非零二进向量，代表i 与前一个点的相对位置差（例如，在二维的情况下，b =（1，1）、（1，0）或（0，1）），s 代表待比对的序列集合，而
其中，s j [ij ]表示第j 条序列在第i j 位的字符，SP-score(Column(s , i , b )) 代表SP 模型的得分值。计算过程是一个递推的过程，在计算每个晶格节点得分的时候，将其各前趋节点的值分别加上从前趋节点到当前点的SP 得分，然后，取最大值作为当前节点的值。
随着待比对的序列数目增加，计算量和所要求的计算空间猛增。 对于k 条序列的比对，动态规划算法需要处理k 维空间里的每一个节点，计算量自然与晶格中的节点数成正比，而节点数等于各序列长度的乘积。另外，计算每个节点依赖于其前趋节点的个数。在二维情况下，前趋节点个数等于3，三维等于7，四维等于15。对于k 维超晶格，前趋节点的个数等于2- 1。因此，用动态规划方法计算多重比对的时间复杂度为O(2∏i=1,...,k ∣s i ∣) 。这个计算量是巨大的，而动态规划算法对计算空间的要求也是很大的。 k k
称一个问题在多项式时间内可解，当且仅当存在一个时间复杂度为O(n) 的求解算法，其中c 是常数，n 是问题的输入规模（如序列的长度、序列的个数）。例如前面介绍的序列两两比较动态规划算法的时间复杂度为O(n) 。假设有k 个长度n 的序列，则利用SP 模型寻找最优多重序列比对算法的时间复杂度为O(2n ) ，这里k 和n 都是问题的规模。 k k 2c
可以证明，利用SP 模型寻找最优多重序列比对是一个NP-完全问题。1971年，Cook 给出NP-完全问题的定义。P 类问题为多项式界的问题；而NP-完全问题是这样一类问题，如果其中的某个问题存在多项式界的算法，则NP 类中的每一个问题都存在一个多项式界算法。NP-完全问题通常被认为是一些人们难以在有限的时间、空间内对问题求出最佳解的问题，几乎所有专家都认为不可能在多项式时间内准确求解NP-完全问题。到目前为止，已经证明属于NP-完全问题的达两千多个，这些问题同旅行售货员问题、有向图哈密顿回路等问题是一样难的，换句话说，要么大家都有多项式时间算法，要么都没有。目前，人们对NP-完全问题通常采用以下二种方法进行
求解：舍去寻找最优解的要求，寻找对一般问题比较接近最优解的近似解；利用非常规的求解技术求解，例如，利用神经网络、遗传算法等方法进行问题求解。
对于生物信息学的问题，有许多是NP-完全问题。在实际应用中，有一些变通的方法可以用来求解NP-完全问题。概括如下：（1）只求解规模比较小的问题；（2）利用动态规划、分支约束等技术减小搜索空间，提高求解问题的效率；（3）针对具体问题的特点，根据实际输入情况，设计实用求解算法，这样的算法虽然在最坏的情况下其时间复杂度是非多项式的，但是算法执行的平均效率和复杂度与多项式的算法相当；（4）采用近似算法或者启发式方法，如局部搜索、模拟退火、遗传算法等。
对基于SP 模型寻找最优多重序列比对这样一个问题，可以用近似的方法求解，其算法的时间复杂度可用多项式表示。下面分别介绍几种实用的多重序列比对算法。
3.3.3 优化计算方法
如果待比对的序列很多，多重序列比对所形成的超晶格空间将会非常大，算法不可能访问所有的节点，而且也没有必要这样做。正如序列两两比对一样，可以想象代表最优的路径处于主对角线附近，至少不会在超晶格中的某个平面上。假设已知一个试探性的路径与最优的路径靠近，并且两条路径最多相距r 个单位，则应当将动态规划算法所要计算的节点限制在以试探性路径为中心、半径为r 的区域内，这样做无疑可以提高算法的效率。
我们不能盲目地搜索整个空间，应该利用人工智能空间搜索策略的剪枝技术，根据问题本身的特殊性将搜索空间限定在一个较小的区域范围内。若问题是搜索一条得分最高（或代价最小）的路径，那么，在搜索时，如果当前路径的得分低于某个下限（或累积代价已经超过某个上限），则对当前路径进行剪枝，即不再搜索当前路径的后续空间。
在序列两两比对中, Fickett和Ukkonen 设计了一种称为定界约束的方法来缩小搜索空间、减少计算量，其中得分矩阵的上界和下界可以预先确定或动态变化。为了在多维空间上使用动态规划算法，Carrillo 和Lipman 将这种思想引入到多重序列比对，即先进行初步的序列双重比对，以便限制进一步作多重序列全面比对所需要的多维空间，从而减少计算量。这样可以克服多序列的维数、空间以及运算量之间的矛盾。
这里介绍Carrillo-Lipman 的优化计算方法，该算法是基于下述关系的，即
多重序列比对与向两个序列投影之间的关系，也就是通过公式（3-36）将SP 得分与序列两两比对的得分联系起来。
设k 条序列s 1, s2, ... ,sk 的长度分别为n 1, n2, ... , nk ，按照SP 得分模型计算这些序列的最优比对。依然采用动态规划方法，但并不计算超晶格空间中所有的节点，而是仅处理“相关”的节点。但是，哪些节点是相关的呢？这需要观察节点在两条序列所在平面上的投影。
假设α是关于k 条序列s 1, s2, ... ,sk 的最优多重比对。注意：这并不意味着α向某两条序列的投影是这两条序列的最优比对。我们可以按下述方法测试超晶格空间中的一个节点是否为相关节点：从某个节点向任何两条序列所在的平面投影，如果该投影是这两条序列两两最优比对的一部分（前面一部分），则该节点是相关节点。
在说明具体的优化计算方法之前，先介绍一种计算两条序列经过特定断点的最优比对算法。设有两条序列s 、t ，已知它们的两个断点分别是i 、j （即s 和t 分别在i 和j 处一分为二），则s 、t 对于经过特定断点（i 、j ）的最优比对可分为两个部分，一是对应于两条序列前缀0:s:i 与0:t:j 的最优比对，另一个是对应于两条序列后缀i :s:m 与j :t:n 的最优比对, m、n 分别为两条序列的长度。实际上，经过特定断点的最优比对是两个比对。为了得到特定断点的最优比对，用两个矩阵A 和B ，其值为：
a i,j = sim(0:s:i , 0:t:j )
b i,j = sim(i :s:m , j :t:n )
对于矩阵A 的计算与标准算法一样，而对于矩阵B 的计算则是反方向的，即先对B 的最后一行和最后一列进行初始化，然后反向推进到（0，0）。这样，矩阵A 与B 的和（A+B=C）包含了经过任意特定断点（i 、j ）的最优比对得分。我们称C 矩阵为总得分矩阵，而A 、B 分别是前缀和后缀的得分矩阵。
根据矩阵C ，我们可以迅速求出两个序列的最优比对，而不需要像上一节那样用反向跟踪的方法求出最优比对所对应的路径。图3.19列出了根据某个打分模型计算得到的矩阵A 和C
。可
以看出，根据C 的最大值，可非常容易地找出最优比对所对应的路径。最优路径通过一系列断点，经过这一系列断点的所有最优比对得分值相同，实际上这个得分值就是两条序列比对的得分。换句话说，这些断点都处于最优路径上。
虽然最优多重比对的投影不一定是两两最优比对，但是我们可以为投影的得分设立一个下限。设α是关于s 1, s2, ... ,sk 的最优比对，如果SP-score(α) ≥ L，则可以证明
(3-39)
其中score(αi,j ) 是α在s i , 和s j 所在平面投影的得分，
(3-40)
这里，L 实际上是最优多重比对的一个下限，而L i,j 是序列s i 和序列s j 比对得分的一个下限。现在，需要判断超晶格中的一个节点i = (i1, i2, „, ik )是否在L 的限制下与最优比对相关。
简单地说， 如果一个节点满足 (3-39)式 的条件，则该节点是相关的；若条件不满足，即score(αi,j ) 小，则不可能是相关的，因此，i 肯定不会处于最优路径上。当然，上述条件只是一个必要条件，但不是充分条件。满足条件只是说明i 可能处于最优路径，但不一定处于最优路径。条件的作用是限制搜索空间，提高算法的实施效率。
将判断条件进一步具体化，则对于所有的1≤ x
C x ，y [ix , iy ] ≥ L x ，y (3-41)
则i 是相关的。这里，C x ，y 是序列s x 和s y 的总得分矩阵，C x ，y [ix , i y ]表示在点[ix , i y ]处的值。这个具体的条件是根据前面的叙述推演出来的。假设最优比对α所对应的路径通过节点（i 1, i 2,„,ix , .., iy ,„,ik ），则C x ，y [ix , iy ] ≥ score(αi,j ) ≥ L x ，y 。 因此，如果C x ，y [ix , iy ]
为了得到一个合理的下限L ，我们可以任选一个包含所有序列的多重比对，计算该多重比对的得分，以此作为L 。若选取的L 接近于最优值，算法速度将大大提高。注意：L 的值不能超过最优值，否则算法出错。
在实现上述优化计算方法时，必须非常仔细。不可能对超晶格中的每一个节点都进行上述判断，否则，算法的计算时间不会有多大的减少。因此，我们需要一种剪枝方法，它一次能够将许多不相关节点删除。
从超晶格的零点0 = (0, 0, „ , 0) 出发，此节点总是相关的。逐步扩展节点，直到终点(n1, „, nk ) 。在扩展过程中，仅分析相关节点。以数组a[]保存各节点的计算结果。称一个节点i 影响另一个节点j ，如果在计算a[j ]时用到i 。这又称为j 依赖于i ，每一个节点依赖于其它2 k -1个节点。为了便于处理，设置一个缓冲区，该缓冲区内仅存放相关节点的后续节点。首先将0放入其中。当一个节点i 进入缓冲区时，其对应的值a[i ]被初始化，然后a[i ]的值在随后的阶段中被更新。当该节点离开缓冲区时，其值即为该点真正的值，并用于其他节点（依赖于此节点）的计算。当然，其后续节点是否要计算，还取决于i 是否为相关节点，若不是，则不再计算其后续的其他节点。
具体的过程如下：设j 是一个依赖于i 的相关节点，如果j 不在缓冲区内，则将其放入缓冲区，并计算
a[j ] ← a[i ] + SP-score(Column(s , i , b ))
如果早已在缓冲区中，则按下式更新
a[j ] ← max(a[j ], a[i ] + SP-score(Column(s , i , b )))
Carrilo-Lipman 算法要求待比较的多个序列具有较大的相似性，并且序列数不能太多。幸运的是，在实际工作中，许多生物分子序列的比较满足这两个条件。
3.3.4 星形比对
如果采用标准的动态规划算法计算最优的多重序列比对，会需要很长的时间，即使采用前面所介绍的方法，也一样需要大量的时间。所以，必须考虑其他的方法，首选的就是所谓启发式方法。启发式方法不一定保证最终能得到最优解，但在大多数情况下，其计算结果接近于最优结果；并且，重要的一点是，这类方法能够大大减少所需的计算时间，加快计算速度。目前所用的算法大部分将序列多重比对转化为序列两两比对，逐渐将两两比对组合起来，最终形成完整的多序列比对。这种方式又称为渐进法。可以按照星形结构或者树形结构组合两两序列比对。
星形比对是一种启发式方法，首先由Gusfield 提出。星形比对的基本思想是：在给定的若干条序列中，选择一个核心序列，通过该序列与其他序列的两两比对，形成所有序列的多重比对α，从而使得α在核心序列和任何一个其他序列方向的投影是最优的两两比对。
设s 1, s2, ... ,sk 是k 条待比对的序列（暂时不考虑如何选择核心序列）。假设已知核心序列是s c ，c 介于1到k 之间，则可以利用标准的动态规划算法求出所有s i 和s c 的最优两两比对。这个工作需要时间为O （kn ）（假设所有序列的长度为n ）。将这些序列的两两比对聚集起来，并采用“只要是空位，则永远是空位”的原则。聚集过程从某一个两两比对开始，如s 1和s c ，然后逐步加上其他的两两比对。在这个过程中，逐步增加s c 中的空位字符，以适应其他的比对，但决不删除s c 中已存在的空位字符。假设在上述过程中的某一时刻，有一个由s c 指导的、已经建立好的部分序列的多重比对，接下来就是加入一个新的、与s c 两两比对的序列，如果需要，2
则插入新的空位字符。这个过程一直进行到所有的两两比对都加入以后结束，所需的计算量为O(kn) 。因此，算法总的时间复杂度为O(kn+kn) 。 222
那么，如何选择核心序列呢？一种方法就是尝试将每一个序列分别作为核心序列，按上述过程进行，取结果最好的一个。另一种方法是计算所有的两两比对，取下式值最大的一个：
(3-42) 例如，有5个序列，
s 1 = ATTGCCATT
s 2 = ATGGCCATT
s 3 = ATCCAATTTT
s 4 = ATCTTCTT
s 5 = ACTGACC
根据它们之间的相似性度量建立如图3.20所示的表格，得分函数见公式（3-2）。从表格中可以看出，选s c =s1时，式（3-42）取最大值。接下来，计算s 1与其他序列的最优两两比对。假设得到下述各两两比对：
ATTGCCATT ATTGCCATT-- ATTGCCATT ATTGCCATT
ATGGCCATT ATC-CAATTTT ATCTTC-TT ACTGACC--
图3.20 两两比对的得分
在这个例子中，核心序列s 1中的空位字符是由于s 3序列在末端引入的两个空位字符。最后多重序列比对的结果如下：
ATTGCCATT--
ATGGCCATT--
ATC-CAATTTT
ATCTTC-TT--
ACTGACC----
星形比对是一种多重序列比对近似的方法，然而是一种比较好的近似方法。如果用代价来评判来多重序列的比对结果，则可以证明，用该方法所得到的多重序列比对的代价不会大于最优多重序列比对代价的两倍。
3.3.5 树形比对
假设有k 条待比对的序列，同时有一棵具有k 个叶节点的树，其中每个叶节点对应一条序列。如果将序列也赋予树的内部节点，则可以计算树中每个分支的权值。权值代表对应分支连接的两条序列之间的相似度。所有权值的和就是这棵树（对应于将序列赋予内部节点的特定方式）的得分。树形比对的问题就是寻找一种树的内部节点序列赋予方式，使得树的得分最大。星形比对可以认为是树形比对的一种特例。
例如，图3.21中有4条简单序列，将CT 、CG 、CT 分别赋予节点x 、y 、z ，则树的得分为8。这里，假设a=b时，p(a,b)=1，否则p(a,b)=0，p(a,-)= p(-，b)= -1。
树形比对问题实际上是组合问题，是一个NP 完全问题，算法的复杂度是指数函数。在进行树形比对时一般采用启发式方法，以距离代替权值，而不是用相似度。
在将多重序列比对分解为两两序列比对的过程中，需要有一个算法合并两个比对, 即比对的比对算法（Alignment of alignments, AA）。
假设有两个多重序列比对α1和α2，α1代表序列s 1,s 2,„,si 的多重比对，α2代表序列
t 1,t 2,„,tj 的多重比对，并且（s 1，s 2，„，s i ）⋂（t 1，t 2，„，t j ）=φ，α代表s 1和t 1的两两比对，则计算与α相一致的α1和α2比对的算法如下：
① 标定α1的各列，如果s 1在比对中对应位置的编辑操作不是插入或删除，则这些列分别标记为
s 1对应位置上的字符a 1,a 2,„,als1（ls 1为序列s 1的长度）；
② 标定α2的各列，如果t 1在比对中对应的位置编辑操作不是插入或删除，则这些列分别标记为
t 1对应位置上的字符b 1,b 2,„,bl2（lt 1为序列t 1的长度）；
③ 对a 1,a 2,„,al1和b 1,b 2,„,bl2进行比对；
④ 在所得到的比对中，对于α1、α2和α中原来有插入或删除操作的位置，恢复其原有的实际字
符或空位字符“-”。
例：
α1：s 1 -H-LVV α2：t 1 L-HCLV α：s 1 -H-LVV
s 2 G-VLVC t2 VLHCL- t1 LHCLV-
s 3 GN-LVV
AA 算法的输出为
--H--LVV
-G--VLVG
-GN--LVV
L-HC-LV-
V-HC-L--
分别对第1、2列和4、5列进行压缩，则最后结果为
-H —LVV
G —VLVG
GN —LVV
LHCLV-
VHCL--
对于n 条序列的进行树形比对的一种最直接算法的计算过程如下：
① 初始化，对于每条序列，生成一个叶节点；
② 利用AA 算法合并两个节点，形成一个新节点，合并的结果放在新节点中，原来的两个节点作
为新节点的子节点；
③ 反复执行②，直到形成n 个叶节点的树根为止，根节点中的序列即为最终的多重比对结果。 与两个序列的比较一样，对于多重序列比较，也需要进行显著性的检验，以判断一个多重序列的比对是否具有真正的意义。实际上，任何一个序列比较过程都应该包含两个步骤，即序列比对和显著性检查。
3.3.6 其他多重序列比对算法
在多重序列近似比对算法方面，除渐进方法之外，还有许多其他方法，如用遗传算法、模拟退火算法、隐马尔柯夫模型等解决序列的多重比对问题。多重序列比对程序有ClustalW 、MAP 、MSA 、PILEUP 等，这些程序用于进行全局比对。可进行局部比对的程序有PIMA 、Block Maker、MEME 、MACAW 、SAM 等，另一个程序是Dialign ，该程序注重于寻找多条序列中相似的区域。
目前，使用最广泛的多重序列比对程序是ClustalW 。ClustalW 是一种渐进的比对方法，先将多个序列进行两两比对，基于这些比较，计算得到一个距离矩阵，该矩阵反映每对序列之间的
关系。根据距离矩阵计算产生系统发生树，从最紧密的两条序列开始，逐步引入邻近的序列，并不断重新构建比对，直到所有序列都被加入为止。
ClustalW 的程序可以自由使用，在任何主要的计算机平台上都可以运行。在美国国家生物技术信息中心NCBI 的FTP 服务器上可以找到下载的软件包，在欧洲生物信息学研究所EBI 的主页还提供了基于Web 的ClustalW 服务，用户可以把序列和各种要求通过表单提交到服务器上，服务器把计算的结果用Email 返回给用户。EBI 的ClustalW 网址是：http://www.ebi.ac.uk/clustalw/。ClustalW 对用户输入序列的格式和输出格式的选择比较灵活，可以是FASTA 格式，也可以是其他格式。
例：
3.3.7 统计特征分析
对于所得到的多重序列比对，我们往往需要进行归纳分析，总结这些序列的特征，或者给出这些序列共性的表示。一种是保守序列表示方式，表示序列每个位置上最可能出现的字符（或者所有可能出现的字符）。另一种表示方式是所谓的特征统计图谱（profile ）表示方式，或特征统计矩阵（matrix ）表示方式。
给定一个长度为L 的多重序列比对α，令P=(P1,P 2,„,PL ) ，P 表示在α的每一列上各种字符出现的概率分布。令P j =(pj,0,p j,1,„,pj,|A|) ，A 代表字母表，P j,k 代表字母表A 中第k 个字符在第j 列出现的概率。需要说明的是，第0个字符是特殊的空位符号“-”。称P 为多重序列比对α的特征统计图谱，或特征统计矩阵，它反映α各个列的特征。
利用特征统计矩阵可以判断一个序列是否满足一定的特征。给定一个序列s=a1a 2„am ，定义字符a 在第j 位的代价为
（3-43） 其中，w 是代价函数，|A|代表字母表A 的长度，A k 代表A 的第k 个字符，特别地A 0代表空缺字符“-”。如果用特征统计矩阵衡量整个序列s ，其代价为
（3-44）
将一条序列与特征统计矩阵相对照，如果代价值小，说明该序列具有相应的特征，否则，该序列不具备相应的特征。在序列数据库中，往往将各个序列按照同源关系进行分类，形成一系列的家族。对于每一个家族，构建对应的特征统计矩阵。这样，在进行数据库查询时，可以利用特征统计矩阵判断一条输入的序列是否属于特定的家族。反过来，也可以利用已知的特征统计矩阵搜索序列数据库，在研究家族关系时，这比用单个序列搜索数据库更有意义。当用单个序列搜索数据库时，我们希望找到与查询序列相似的序列；而用特征统计矩阵搜索数据库时，我们希望考察家族的成员关系。
返回总目录
第三章 序列比较
3.3 序列多重比对
与序列两两比对不一样，序列多重比对（Multiple Alignment ）的目标是发现多条序列的共性。如果说序列两两比对主要用于建立两条序列的同源关系和推测它们的结构、功能，那么，同时比对一组序列对于研究分子结构、功能及进化关系更为有用。例如，某些在生物学上有重要意义的相似性只能通过将多个序列对比排列起来才能识别。同样，只有在多序列比对之后，才能发现与结构域或功能相关的保守序列片段。对于一系列同源蛋白质，人们希望研究隐含在蛋白质序列中的系统发育的关系，以便更好地理解这些蛋白质的进化。在实际研究中，生物学家并不是仅仅分析单个蛋白质，而是更着重于研究蛋白质之间的关系，研究一个家族中的相关蛋白质，研究相关蛋白质序列中的保守区域，进而分析蛋白质的结构和功能。序列两两比对往往不能满足这样的需要，难以发现多个序列的共性，必须同时比对多条同源序列。
图3.14是从多条免疫球蛋白序列中提取的8个片段的多重比对。这8个片段的多重比对揭示了保守的残基（一个是来自于二硫桥的半胱氨酸，另一个是色氨酸）、保守区域（特别是前4个片段末端的Q-PG ）和其他更复杂的模式，如1位和3位的疏水残基。实际上，多重序列比对在蛋白质结构的预测中非常有用。
多重比对也能用来推测各个序列的进化历史。从图3.14可以看出，前4条序列与后4条序列可能是从两个不同祖先演化而来，而这两个祖先又是由一个最原始的祖先演化得到。实际上，其中的4个片段是从免疫球蛋白的可变区域取出的，而另4个片段则从免疫球蛋白的恒定区域取出。当然，如果要详细研究进化关系，还必须取更长的序列进行比对分析。
对于多重序列比对的定义，实际上是两个序列的推广。设有k 个序列s 1, s 2, ... ,sk ，每个序列由同一个字母表中的字符组成，k 大于2；通过插入操作，使得各序列s 1, s2, ... ,sk 的长度一样，从而形成这些序列的多重比对。如果将各序列在垂直方向排列起来，则可以根据每一列观察各序列中字符的对应关系，如图3.14。
通过序列的多重比对，可以得到一个序列家族的序列特征。当给定一个新序列时，根据序列特征，可以判断这个序列是否属于该家族。对于多序列比对，现有的大多数算法都基于渐进比对的思想，在序列两两比对的基础上逐步优化多序列比对的结果。进行多序列比对后，可以对比对结果进行进一步处理，例如构建序列的特征模式，将序列聚类，构建分子进化树等。
3.3.1 SP模型
SP 模型（Sum-of-Pairs ，逐对加和）是一种多重序列比对的评价模型。在多重比对中，首先要对所得到的比对进行评价，以确定其优劣。例如，对图3.14中的8条序列进行比对，可以得到另外两种结果，如图3.15所示。那么，这样的三个多重比对，哪一个更好呢？这就需要有一种方法来评价一个多重比对。
评价一个多重序列比对比评价序列两两比对结果更复杂。这里，我们假设得分（代价）函数具有加和性，即多重比对的得分是各列得分总和。那么，我们首先考虑如何给比对的每一列打分，然后将各列的和加起来，成为一个总得分。
在处理每一列时，自然的处理方式是寻找一个具有k 个变量的打分函数（k 是参与多重比对的序列的个数），而每一个变量或者是一个来自特定字母表中的字符，或者是一个空位。我们很难得到这样一种具有k 个变量的表达式函数。另一方面，这种隐式函数不具有统一的形式，随着k 的变化，函数的表现形式也发生变化，不利于计算机处理。可以考虑使用显式函数，在具体实现显式函数时，用一个k 维数组来表示该显式函数（类似于打分矩阵），指定对应于k 个变量各种组合的函数值。这带来一个新问题，即所需的数组空间很大，而且随着k 的变化，数据结构也要随之动态变化。
我们所期望的函数在形式上应该简单，具有统一的形式，不随序列的个数而发生形式变化。根据得分函数的意义，函数值应独立于各参数的顺序，即与待比较的序列先后次序无关。另外，对相同的或相似字符的比对，奖励的得分值高，而对于不相关的字符比对或空白，则进行惩罚（得分为负值）。满足上述条件的一个函数就是常用的逐对加和SP 函数。SP 函数定义为一列中所有字符对得分之和：
(3-34)
其中，c1,c2,„,ck是一列中的k 个字符，p 是关于一对字符相似性的打分函数。对于p 可采用不同的定义。下面一个是SP 函数计算例子。
总得分根据字符两两比较得分计算演化而来，即逐对计算p(1,2)，p (1,3)，... ，p (1,8)，p (2,3)，p(2,4)，... ，p (2,8)，... ，p (7,8) 的所有得分，再加和得到结果。在上述计算中，我们假定：如果两个对比的字符相同，则得分为0，否则，得分为-1。所以，上述一列的得分为：（-7-6-5-4-3-2-1）+2= -26。
将一个多重比对所有列的得分全部加起来，其和即为该多重序列比对的得分。
注意：在进行多重序列比对时，可能会出现两个空位字符的比对，因此需要扩充函数p 的定义域，即增加p(-，-) 的定义。通常的定义是：
p(- , -) =
0 （3-35）
这似乎有点意外，因为一般只要是遇到空位字符，其得分就应该是负的，所以两个空位字符的比对应得到更多的负分。但是，这样处理是有道理的。在序列多重比对中，我们往往在得到整体比对的基础上进一步分析两条序列的对应关系。例如，根据图3.14的比对结果，取出最后两条比对的序列，见图3.16（a ）。这里存在空位字符对比的列，相当于这两条序列都进行了插入操作。但是由于插入位置相同[如图3.16（a ）箭头所指位置]，这两条序列本身在此位置上是完全相同的，所以，此位置上的编辑代价为零，或者得分为0。因而在分析这两条序列时，可以去掉这些空位字符对比的列，得到由多重序列比对结果推演的序列两两比对，如图3.16（b ）所示。其结果又称为多重序列比对在两条特定序列上的投影（projection ）。
若先处理每一个序列对，而在处理序列对时，逐个计算字符对，最后加和，则SP 得分模型的计算公式如下：
(3-36)
其中，α是一个多重比对，αi,j 是由α推演出来的序列s i 和s j的两两比对。
具体计算SP 时，我们可以先对多重序列比对a 的每一列进行计算，然后将每一列的得分值相加；也可以先计算每一对推演出来的两两序列比对的得分，然后再加和。但是注意，上述两种计算过程仅在p(- , -)=0这一条件成立下才等价。
3.3.2 多重比对的动态规划算法
多重序列比对的最终目标是通过处理得到一个得分最高（或代价最小）的序列对比排列，从而分析各序列之间的相似性和差异。如同处理序列两两比对一样，依然可以用动态规划算法。 回想前一节中穿越得分矩阵的路径（见图3.9）。得分矩阵相当于二维平面，而对于三条序列，每一种可能的比对可类似地用三维晶格中的一条路径表示，而每一维对应于一条序列，如图
3.17所示。路径的起点为晶格的左下角，路径的终点为晶格的右上后角。图3.17中的路径记为 (0,0,0)，(1,0,0)，(2,1,0)，(3,2,0)，(3,3,1)，(4,3,2)。如果存在多条序列，则形成的空间是超晶格（Hyperlattice ）。
假设在三维晶格的顶端、前面、后面各有一平行光源，这些光源将代表多重序列比对的路径投影到相对的侧面，即投影到一对序列所在的平面（图3.17仅画出了右边的光源）。将多重比对投影到两个序列所在的平面在加速优化计算中将起着重要的作用。
一个多重序列比对投影的结果可能比原来要短，例如，下面的多重序列比对
G---SNS
GN----S
GNAVSNS
如果在前两条序列方向进行投影，则投影结果为：
G-SNS
GN--S
如果路径的某一段沿投影方向进行，那么该段路径就不可能产生投影。图3.17中的第一段沿着第一条序列向前进，推进方向垂直于其他两条序列，则平行光源对这一段不产生直线投影，而仅仅产生一个投影点。
图3.17三条序列所对应的三维晶格
序列两两比对的动态规划算法经改进后可直接用于序列的多重比对。就二重序列而言，将动态规划算法的计算过程看成是在二维平面上按一定顺序访问每一个节点，访问节点的先后顺序取决于节点之间的关系，如图3.7所示。二维平面（或二维晶格）与前面所介绍的超晶格相似，上一节得分矩阵中的位置与每一个晶格节点相对应。在计算过程中，对每个节点逐一计算其得分（或代价）。每个节点的得分代表两个序列前缀的最优比对的得分，而晶格最后一个节点（右下角节点）的得分即为两条完整序列的比对得分。
在超晶格中，序列的放置与上一节有所不同，计算是从左下角进行的，而不是得分矩阵中的左上角。如图3.18所示，计算过程从超晶格的坐标点（0,0,„,0）开始，按节点之间的依赖关系向右上后方推进，直到计算完最后一个节点。实际计算时，可以采用类似二维的计算方式，即先低维后高维（对应于先行后列）。在图3.18中，当前点的得分计算取决于与它相邻的7条边，分别对应于匹配、替换或引入空位等三种编辑操作。计算各操作的得分（包括前趋节点的得分），选择一个得分最大的操作，并将得分和存放于该节点。在三维或超晶格中，计算公式与上一节相似，但计算一个节点依赖于更多的前趋节点。在三维情况下要考虑7个前趋节点，在k 维情况下要考虑2-1个前趋节点。
k
假设以k 维数组A 存放超晶格，则计算过程如下：
① a[ 0, 0, „ ,0 ] = 0
② a[ i ] = max {a[ i - b ] + SP-score[Column(s , i , b )]}
这里i 是一个向量，代表当前点，b 是具有k 个元素的非零二进向量，代表i 与前一个点的相对位置差（例如，在二维的情况下，b =（1，1）、（1，0）或（0，1）），s 代表待比对的序列集合，而
其中，s j [ij ]表示第j 条序列在第i j 位的字符，SP-score(Column(s , i , b )) 代表SP 模型的得分值。计算过程是一个递推的过程，在计算每个晶格节点得分的时候，将其各前趋节点的值分别加上从前趋节点到当前点的SP 得分，然后，取最大值作为当前节点的值。
随着待比对的序列数目增加，计算量和所要求的计算空间猛增。 对于k 条序列的比对，动态规划算法需要处理k 维空间里的每一个节点，计算量自然与晶格中的节点数成正比，而节点数等于各序列长度的乘积。另外，计算每个节点依赖于其前趋节点的个数。在二维情况下，前趋节点个数等于3，三维等于7，四维等于15。对于k 维超晶格，前趋节点的个数等于2- 1。因此，用动态规划方法计算多重比对的时间复杂度为O(2∏i=1,...,k ∣s i ∣) 。这个计算量是巨大的，而动态规划算法对计算空间的要求也是很大的。 k k
称一个问题在多项式时间内可解，当且仅当存在一个时间复杂度为O(n) 的求解算法，其中c 是常数，n 是问题的输入规模（如序列的长度、序列的个数）。例如前面介绍的序列两两比较动态规划算法的时间复杂度为O(n) 。假设有k 个长度n 的序列，则利用SP 模型寻找最优多重序列比对算法的时间复杂度为O(2n ) ，这里k 和n 都是问题的规模。 k k 2c
可以证明，利用SP 模型寻找最优多重序列比对是一个NP-完全问题。1971年，Cook 给出NP-完全问题的定义。P 类问题为多项式界的问题；而NP-完全问题是这样一类问题，如果其中的某个问题存在多项式界的算法，则NP 类中的每一个问题都存在一个多项式界算法。NP-完全问题通常被认为是一些人们难以在有限的时间、空间内对问题求出最佳解的问题，几乎所有专家都认为不可能在多项式时间内准确求解NP-完全问题。到目前为止，已经证明属于NP-完全问题的达两千多个，这些问题同旅行售货员问题、有向图哈密顿回路等问题是一样难的，换句话说，要么大家都有多项式时间算法，要么都没有。目前，人们对NP-完全问题通常采用以下二种方法进行
求解：舍去寻找最优解的要求，寻找对一般问题比较接近最优解的近似解；利用非常规的求解技术求解，例如，利用神经网络、遗传算法等方法进行问题求解。
对于生物信息学的问题，有许多是NP-完全问题。在实际应用中，有一些变通的方法可以用来求解NP-完全问题。概括如下：（1）只求解规模比较小的问题；（2）利用动态规划、分支约束等技术减小搜索空间，提高求解问题的效率；（3）针对具体问题的特点，根据实际输入情况，设计实用求解算法，这样的算法虽然在最坏的情况下其时间复杂度是非多项式的，但是算法执行的平均效率和复杂度与多项式的算法相当；（4）采用近似算法或者启发式方法，如局部搜索、模拟退火、遗传算法等。
对基于SP 模型寻找最优多重序列比对这样一个问题，可以用近似的方法求解，其算法的时间复杂度可用多项式表示。下面分别介绍几种实用的多重序列比对算法。
3.3.3 优化计算方法
如果待比对的序列很多，多重序列比对所形成的超晶格空间将会非常大，算法不可能访问所有的节点，而且也没有必要这样做。正如序列两两比对一样，可以想象代表最优的路径处于主对角线附近，至少不会在超晶格中的某个平面上。假设已知一个试探性的路径与最优的路径靠近，并且两条路径最多相距r 个单位，则应当将动态规划算法所要计算的节点限制在以试探性路径为中心、半径为r 的区域内，这样做无疑可以提高算法的效率。
我们不能盲目地搜索整个空间，应该利用人工智能空间搜索策略的剪枝技术，根据问题本身的特殊性将搜索空间限定在一个较小的区域范围内。若问题是搜索一条得分最高（或代价最小）的路径，那么，在搜索时，如果当前路径的得分低于某个下限（或累积代价已经超过某个上限），则对当前路径进行剪枝，即不再搜索当前路径的后续空间。
在序列两两比对中, Fickett和Ukkonen 设计了一种称为定界约束的方法来缩小搜索空间、减少计算量，其中得分矩阵的上界和下界可以预先确定或动态变化。为了在多维空间上使用动态规划算法，Carrillo 和Lipman 将这种思想引入到多重序列比对，即先进行初步的序列双重比对，以便限制进一步作多重序列全面比对所需要的多维空间，从而减少计算量。这样可以克服多序列的维数、空间以及运算量之间的矛盾。
这里介绍Carrillo-Lipman 的优化计算方法，该算法是基于下述关系的，即
多重序列比对与向两个序列投影之间的关系，也就是通过公式（3-36）将SP 得分与序列两两比对的得分联系起来。
设k 条序列s 1, s2, ... ,sk 的长度分别为n 1, n2, ... , nk ，按照SP 得分模型计算这些序列的最优比对。依然采用动态规划方法，但并不计算超晶格空间中所有的节点，而是仅处理“相关”的节点。但是，哪些节点是相关的呢？这需要观察节点在两条序列所在平面上的投影。
假设α是关于k 条序列s 1, s2, ... ,sk 的最优多重比对。注意：这并不意味着α向某两条序列的投影是这两条序列的最优比对。我们可以按下述方法测试超晶格空间中的一个节点是否为相关节点：从某个节点向任何两条序列所在的平面投影，如果该投影是这两条序列两两最优比对的一部分（前面一部分），则该节点是相关节点。
在说明具体的优化计算方法之前，先介绍一种计算两条序列经过特定断点的最优比对算法。设有两条序列s 、t ，已知它们的两个断点分别是i 、j （即s 和t 分别在i 和j 处一分为二），则s 、t 对于经过特定断点（i 、j ）的最优比对可分为两个部分，一是对应于两条序列前缀0:s:i 与0:t:j 的最优比对，另一个是对应于两条序列后缀i :s:m 与j :t:n 的最优比对, m、n 分别为两条序列的长度。实际上，经过特定断点的最优比对是两个比对。为了得到特定断点的最优比对，用两个矩阵A 和B ，其值为：
a i,j = sim(0:s:i , 0:t:j )
b i,j = sim(i :s:m , j :t:n )
对于矩阵A 的计算与标准算法一样，而对于矩阵B 的计算则是反方向的，即先对B 的最后一行和最后一列进行初始化，然后反向推进到（0，0）。这样，矩阵A 与B 的和（A+B=C）包含了经过任意特定断点（i 、j ）的最优比对得分。我们称C 矩阵为总得分矩阵，而A 、B 分别是前缀和后缀的得分矩阵。
根据矩阵C ，我们可以迅速求出两个序列的最优比对，而不需要像上一节那样用反向跟踪的方法求出最优比对所对应的路径。图3.19列出了根据某个打分模型计算得到的矩阵A 和C
。可
以看出，根据C 的最大值，可非常容易地找出最优比对所对应的路径。最优路径通过一系列断点，经过这一系列断点的所有最优比对得分值相同，实际上这个得分值就是两条序列比对的得分。换句话说，这些断点都处于最优路径上。
虽然最优多重比对的投影不一定是两两最优比对，但是我们可以为投影的得分设立一个下限。设α是关于s 1, s2, ... ,sk 的最优比对，如果SP-score(α) ≥ L，则可以证明
(3-39)
其中score(αi,j ) 是α在s i , 和s j 所在平面投影的得分，
(3-40)
这里，L 实际上是最优多重比对的一个下限，而L i,j 是序列s i 和序列s j 比对得分的一个下限。现在，需要判断超晶格中的一个节点i = (i1, i2, „, ik )是否在L 的限制下与最优比对相关。
简单地说， 如果一个节点满足 (3-39)式 的条件，则该节点是相关的；若条件不满足，即score(αi,j ) 小，则不可能是相关的，因此，i 肯定不会处于最优路径上。当然，上述条件只是一个必要条件，但不是充分条件。满足条件只是说明i 可能处于最优路径，但不一定处于最优路径。条件的作用是限制搜索空间，提高算法的实施效率。
将判断条件进一步具体化，则对于所有的1≤ x
C x ，y [ix , iy ] ≥ L x ，y (3-41)
则i 是相关的。这里，C x ，y 是序列s x 和s y 的总得分矩阵，C x ，y [ix , i y ]表示在点[ix , i y ]处的值。这个具体的条件是根据前面的叙述推演出来的。假设最优比对α所对应的路径通过节点（i 1, i 2,„,ix , .., iy ,„,ik ），则C x ，y [ix , iy ] ≥ score(αi,j ) ≥ L x ，y 。 因此，如果C x ，y [ix , iy ]
为了得到一个合理的下限L ，我们可以任选一个包含所有序列的多重比对，计算该多重比对的得分，以此作为L 。若选取的L 接近于最优值，算法速度将大大提高。注意：L 的值不能超过最优值，否则算法出错。
在实现上述优化计算方法时，必须非常仔细。不可能对超晶格中的每一个节点都进行上述判断，否则，算法的计算时间不会有多大的减少。因此，我们需要一种剪枝方法，它一次能够将许多不相关节点删除。
从超晶格的零点0 = (0, 0, „ , 0) 出发，此节点总是相关的。逐步扩展节点，直到终点(n1, „, nk ) 。在扩展过程中，仅分析相关节点。以数组a[]保存各节点的计算结果。称一个节点i 影响另一个节点j ，如果在计算a[j ]时用到i 。这又称为j 依赖于i ，每一个节点依赖于其它2 k -1个节点。为了便于处理，设置一个缓冲区，该缓冲区内仅存放相关节点的后续节点。首先将0放入其中。当一个节点i 进入缓冲区时，其对应的值a[i ]被初始化，然后a[i ]的值在随后的阶段中被更新。当该节点离开缓冲区时，其值即为该点真正的值，并用于其他节点（依赖于此节点）的计算。当然，其后续节点是否要计算，还取决于i 是否为相关节点，若不是，则不再计算其后续的其他节点。
具体的过程如下：设j 是一个依赖于i 的相关节点，如果j 不在缓冲区内，则将其放入缓冲区，并计算
a[j ] ← a[i ] + SP-score(Column(s , i , b ))
如果早已在缓冲区中，则按下式更新
a[j ] ← max(a[j ], a[i ] + SP-score(Column(s , i , b )))
Carrilo-Lipman 算法要求待比较的多个序列具有较大的相似性，并且序列数不能太多。幸运的是，在实际工作中，许多生物分子序列的比较满足这两个条件。
3.3.4 星形比对
如果采用标准的动态规划算法计算最优的多重序列比对，会需要很长的时间，即使采用前面所介绍的方法，也一样需要大量的时间。所以，必须考虑其他的方法，首选的就是所谓启发式方法。启发式方法不一定保证最终能得到最优解，但在大多数情况下，其计算结果接近于最优结果；并且，重要的一点是，这类方法能够大大减少所需的计算时间，加快计算速度。目前所用的算法大部分将序列多重比对转化为序列两两比对，逐渐将两两比对组合起来，最终形成完整的多序列比对。这种方式又称为渐进法。可以按照星形结构或者树形结构组合两两序列比对。
星形比对是一种启发式方法，首先由Gusfield 提出。星形比对的基本思想是：在给定的若干条序列中，选择一个核心序列，通过该序列与其他序列的两两比对，形成所有序列的多重比对α，从而使得α在核心序列和任何一个其他序列方向的投影是最优的两两比对。
设s 1, s2, ... ,sk 是k 条待比对的序列（暂时不考虑如何选择核心序列）。假设已知核心序列是s c ，c 介于1到k 之间，则可以利用标准的动态规划算法求出所有s i 和s c 的最优两两比对。这个工作需要时间为O （kn ）（假设所有序列的长度为n ）。将这些序列的两两比对聚集起来，并采用“只要是空位，则永远是空位”的原则。聚集过程从某一个两两比对开始，如s 1和s c ，然后逐步加上其他的两两比对。在这个过程中，逐步增加s c 中的空位字符，以适应其他的比对，但决不删除s c 中已存在的空位字符。假设在上述过程中的某一时刻，有一个由s c 指导的、已经建立好的部分序列的多重比对，接下来就是加入一个新的、与s c 两两比对的序列，如果需要，2
则插入新的空位字符。这个过程一直进行到所有的两两比对都加入以后结束，所需的计算量为O(kn) 。因此，算法总的时间复杂度为O(kn+kn) 。 222
那么，如何选择核心序列呢？一种方法就是尝试将每一个序列分别作为核心序列，按上述过程进行，取结果最好的一个。另一种方法是计算所有的两两比对，取下式值最大的一个：
(3-42) 例如，有5个序列，
s 1 = ATTGCCATT
s 2 = ATGGCCATT
s 3 = ATCCAATTTT
s 4 = ATCTTCTT
s 5 = ACTGACC
根据它们之间的相似性度量建立如图3.20所示的表格，得分函数见公式（3-2）。从表格中可以看出，选s c =s1时，式（3-42）取最大值。接下来，计算s 1与其他序列的最优两两比对。假设得到下述各两两比对：
ATTGCCATT ATTGCCATT-- ATTGCCATT ATTGCCATT
ATGGCCATT ATC-CAATTTT ATCTTC-TT ACTGACC--
图3.20 两两比对的得分
在这个例子中，核心序列s 1中的空位字符是由于s 3序列在末端引入的两个空位字符。最后多重序列比对的结果如下：
ATTGCCATT--
ATGGCCATT--
ATC-CAATTTT
ATCTTC-TT--
ACTGACC----
星形比对是一种多重序列比对近似的方法，然而是一种比较好的近似方法。如果用代价来评判来多重序列的比对结果，则可以证明，用该方法所得到的多重序列比对的代价不会大于最优多重序列比对代价的两倍。
3.3.5 树形比对
假设有k 条待比对的序列，同时有一棵具有k 个叶节点的树，其中每个叶节点对应一条序列。如果将序列也赋予树的内部节点，则可以计算树中每个分支的权值。权值代表对应分支连接的两条序列之间的相似度。所有权值的和就是这棵树（对应于将序列赋予内部节点的特定方式）的得分。树形比对的问题就是寻找一种树的内部节点序列赋予方式，使得树的得分最大。星形比对可以认为是树形比对的一种特例。
例如，图3.21中有4条简单序列，将CT 、CG 、CT 分别赋予节点x 、y 、z ，则树的得分为8。这里，假设a=b时，p(a,b)=1，否则p(a,b)=0，p(a,-)= p(-，b)= -1。
树形比对问题实际上是组合问题，是一个NP 完全问题，算法的复杂度是指数函数。在进行树形比对时一般采用启发式方法，以距离代替权值，而不是用相似度。
在将多重序列比对分解为两两序列比对的过程中，需要有一个算法合并两个比对, 即比对的比对算法（Alignment of alignments, AA）。
假设有两个多重序列比对α1和α2，α1代表序列s 1,s 2,„,si 的多重比对，α2代表序列
t 1,t 2,„,tj 的多重比对，并且（s 1，s 2，„，s i ）⋂（t 1，t 2，„，t j ）=φ，α代表s 1和t 1的两两比对，则计算与α相一致的α1和α2比对的算法如下：
① 标定α1的各列，如果s 1在比对中对应位置的编辑操作不是插入或删除，则这些列分别标记为
s 1对应位置上的字符a 1,a 2,„,als1（ls 1为序列s 1的长度）；
② 标定α2的各列，如果t 1在比对中对应的位置编辑操作不是插入或删除，则这些列分别标记为
t 1对应位置上的字符b 1,b 2,„,bl2（lt 1为序列t 1的长度）；
③ 对a 1,a 2,„,al1和b 1,b 2,„,bl2进行比对；
④ 在所得到的比对中，对于α1、α2和α中原来有插入或删除操作的位置，恢复其原有的实际字
符或空位字符“-”。
例：
α1：s 1 -H-LVV α2：t 1 L-HCLV α：s 1 -H-LVV
s 2 G-VLVC t2 VLHCL- t1 LHCLV-
s 3 GN-LVV
AA 算法的输出为
--H--LVV
-G--VLVG
-GN--LVV
L-HC-LV-
V-HC-L--
分别对第1、2列和4、5列进行压缩，则最后结果为
-H —LVV
G —VLVG
GN —LVV
LHCLV-
VHCL--
对于n 条序列的进行树形比对的一种最直接算法的计算过程如下：
① 初始化，对于每条序列，生成一个叶节点；
② 利用AA 算法合并两个节点，形成一个新节点，合并的结果放在新节点中，原来的两个节点作
为新节点的子节点；
③ 反复执行②，直到形成n 个叶节点的树根为止，根节点中的序列即为最终的多重比对结果。 与两个序列的比较一样，对于多重序列比较，也需要进行显著性的检验，以判断一个多重序列的比对是否具有真正的意义。实际上，任何一个序列比较过程都应该包含两个步骤，即序列比对和显著性检查。
3.3.6 其他多重序列比对算法
在多重序列近似比对算法方面，除渐进方法之外，还有许多其他方法，如用遗传算法、模拟退火算法、隐马尔柯夫模型等解决序列的多重比对问题。多重序列比对程序有ClustalW 、MAP 、MSA 、PILEUP 等，这些程序用于进行全局比对。可进行局部比对的程序有PIMA 、Block Maker、MEME 、MACAW 、SAM 等，另一个程序是Dialign ，该程序注重于寻找多条序列中相似的区域。
目前，使用最广泛的多重序列比对程序是ClustalW 。ClustalW 是一种渐进的比对方法，先将多个序列进行两两比对，基于这些比较，计算得到一个距离矩阵，该矩阵反映每对序列之间的
关系。根据距离矩阵计算产生系统发生树，从最紧密的两条序列开始，逐步引入邻近的序列，并不断重新构建比对，直到所有序列都被加入为止。
ClustalW 的程序可以自由使用，在任何主要的计算机平台上都可以运行。在美国国家生物技术信息中心NCBI 的FTP 服务器上可以找到下载的软件包，在欧洲生物信息学研究所EBI 的主页还提供了基于Web 的ClustalW 服务，用户可以把序列和各种要求通过表单提交到服务器上，服务器把计算的结果用Email 返回给用户。EBI 的ClustalW 网址是：http://www.ebi.ac.uk/clustalw/。ClustalW 对用户输入序列的格式和输出格式的选择比较灵活，可以是FASTA 格式，也可以是其他格式。
例：
3.3.7 统计特征分析
对于所得到的多重序列比对，我们往往需要进行归纳分析，总结这些序列的特征，或者给出这些序列共性的表示。一种是保守序列表示方式，表示序列每个位置上最可能出现的字符（或者所有可能出现的字符）。另一种表示方式是所谓的特征统计图谱（profile ）表示方式，或特征统计矩阵（matrix ）表示方式。
给定一个长度为L 的多重序列比对α，令P=(P1,P 2,„,PL ) ，P 表示在α的每一列上各种字符出现的概率分布。令P j =(pj,0,p j,1,„,pj,|A|) ，A 代表字母表，P j,k 代表字母表A 中第k 个字符在第j 列出现的概率。需要说明的是，第0个字符是特殊的空位符号“-”。称P 为多重序列比对α的特征统计图谱，或特征统计矩阵，它反映α各个列的特征。
利用特征统计矩阵可以判断一个序列是否满足一定的特征。给定一个序列s=a1a 2„am ，定义字符a 在第j 位的代价为
（3-43） 其中，w 是代价函数，|A|代表字母表A 的长度，A k 代表A 的第k 个字符，特别地A 0代表空缺字符“-”。如果用特征统计矩阵衡量整个序列s ，其代价为
（3-44）
将一条序列与特征统计矩阵相对照，如果代价值小，说明该序列具有相应的特征，否则，该序列不具备相应的特征。在序列数据库中，往往将各个序列按照同源关系进行分类，形成一系列的家族。对于每一个家族，构建对应的特征统计矩阵。这样，在进行数据库查询时，可以利用特征统计矩阵判断一条输入的序列是否属于特定的家族。反过来，也可以利用已知的特征统计矩阵搜索序列数据库，在研究家族关系时，这比用单个序列搜索数据库更有意义。当用单个序列搜索数据库时，我们希望找到与查询序列相似的序列；而用特征统计矩阵搜索数据库时，我们希望考察家族的成员关系。
返回总目录