📚 课程目录第2章全局与局部比对

全局与局部比对：Needleman-Wunsch与Smith-Waterman算法原理

序列比对这件事，说白了就是找两段序列的相似之处。但你要问“怎么找”，就有两个思路了：一个是把整条序列从头比到尾，另一个是只盯着最像的那一小段。前者叫全局比对，后者叫局部比对。我给课题组搭过几次比对流程，这两个选择直接影响结果走向——选错了，你可能会错过关键的同源结构域。

全局比对：Needleman-Wunsch 算法

Needleman-Wunsch算法（NW算法）是1970年提出的，算得上序列比对的“老祖宗”。它的核心思想很朴素：用一张二维矩阵记录每次匹配、插入、删除的得分，最后回溯找到最优路径。我习惯把整个过程拆成三步：初始化、填充矩阵、回溯路径。

第一步：初始化矩阵

假设我们有两条序列：SeqA = GATTACA和SeqB = GCATGCU。矩阵的行数 = len(SeqA) + 1，列数 = len(SeqB) + 1。第一行和第一列用来表示空位（gap），得分按空位罚分累加。比如空位罚分设为-2，那么第一行第j列就填 -2*j，第一列同理。

         G   C   A   T   G   C   U
     0  -2  -4  -6  -8 -10 -12 -14
G   -2
A   -4
T   -6
T   -8
A   -10
C   -12
A   -14

嗯，这里要注意：初始化不只是为了填数字，它还规定了比对的对齐起点。如果你忘了处理第一行第一列，后面的动态规划就没法正确传播。

第二步：填充矩阵

每个格子 F(i,j) 的得分来自三个方向：

左上角（对角线）：匹配或替换得分 F(i-1,j-1) + S(SeqA[i], SeqB[j])
左方（插入空位）：F(i-1,j) + gap_penalty
上方（删除空位）：F(i,j-1) + gap_penalty

取这三个值中的最大值。匹配得分我一般设+2，错配-1，空位罚分-2。举个例子，计算 F(1,1)：SeqA[1]=G，SeqB[1]=G，匹配得+2。左上角是0，所以对角得分=0+2=2；左边是-2，加空位罚分=-4；上边是-2，加空位罚分=-4。最大值是2。这样填下去，直到右下角。

         G   C   A   T   G   C   U
     0  -2  -4  -6  -8 -10 -12 -14
G   -2   2   0  -2  -4  -6  -8 -10
A   -4   0   1   2   0  -2  -4  -6
T   -6  -2  -1   0   4   2   0  -2
T   -8  -4  -3  -2   2   3   1  -1
A   -10 -6  -5  -1   0   1   2   0
C   -12 -8  -4  -3  -2  -1   3   1
A   -14 -10 -6  -2  -4  -3   1   2

右下角的得分就是全局最优比对的得分，这里是2。

第三步：回溯路径

从右下角开始，每次往得分最高的前驱方向移动。如果某个格子来自对角线，就表示匹配/替换；来自左边或上边，就表示插入空位。回溯完一条完整路径，我们就得到了比对结果。我曾经给一个全长基因做比对，手动走回溯走了半个多小时——后来赶紧写了个脚本。嗯，别学我，这种机械活儿交给代码就好。

局部比对：Smith-Waterman 算法

NW算法适合两条长度相近的完整序列。但如果你只想找一段保守结构域，比如跨膜区域，全局比对会把不相关的区域也硬对齐上，结果得分被拉低。这时候就需要Smith-Waterman算法（SW算法）。它只有一个关键改动，但这个改动让整个方法不一样。

改动在于：填充矩阵时，如果某个方向来的得分加完是负数，就直接写0。举个例子，F(i-1,j) + gap 如果是负的，那这个格子就填0，而不是把它传下去。这意味着局部比对可以“重新开始”，就像在做滑动窗口一样，只找高分区域。

SW算法三个改动点：

矩阵中不允许出现负值，负值直接归零
回溯时从得分最高的格子开始，而不是从右下角
遇到0就停止回溯，这意味着我们只取一段连续的高分片段

SW填充示例

还是刚才两条序列，用同样的评分矩阵，但加一条规则：每个格子计算完先和0比，取max。结果矩阵就不一样了：

         G   C   A   T   G   C   U
     0   0   0   0   0   0   0   0
G    0   2   0   0   0   2   0   0
A    0   0   1   2   0   0   1   0
T    0   0   0   0   4   2   0   0
T    0   0   0   0   2   3   1   0
A    0   0   0   2   0   1   2   0
C    0   0   2   0   0   0   3   1
A    0   0   0   1   0   0   1   2

最高得分是4（在(3,4)位置）。从那里回溯，得到局部比对片段：AT vs AT？不对，仔细看路径：SeqA[3:4]=AT，SeqB[4:5]=TG？等一下，这里得分路线是(3,4)来自(2,3)（匹配A-A），再来自(1,2)（匹配C？）。嗯，实际上这个例子不完美——我换一条更经典的序列来演示吧。

避坑指南：我曾经在比对短序列时直接用NW算法，结果把两段不相干的区域强行对齐，得到一个伪高分。后来改成SW，一下子就找到了真正的保守motif。所以，什么时候用全局、什么时候用局部，看你的生物学问题：全长同源比对用NW，motif查找/结构域比对用SW。

两种算法的对比

维度	Needleman-Wunsch	Smith-Waterman
比对范围	全序列	局部高分段
负分处理	允许负值传递	负值归零
回溯起点	右下角	矩阵中的最大值
回溯终点	左上角	遇见0为止
时间复杂度	O(n*m)	O(n*m)
典型应用	全基因比对、多序列比对对齐	数据库搜索（BLAST核心）

从计算复杂度看，两者都是O(n*m)，但SW在回溯时只访问局部高分区域，实际速度会快一些。不过，NW和SW都是精确算法——如果你有足够的时间和内存，它们一定能找到最优解。现实中我们用BLAST，就是把SW的精确性做了近似加速。

SVG图：算法流程对比

实际应用中的选择

如果你是做全长cDNA比对，或者构建多序列比对需要的对齐骨架，毫无疑问用NW。但如果你在做数据库搜索，或者想找某个短模序（比如锌指结构域），SW的局部比对更合适。BLAST本质上是加速的SW——它先用种子词索引定位可能匹配的区域，再用SW精确比对。

参数陷阱：我曾经在SW中把gap罚分设得太高（-10），导致所有比对结果都避开了空位，把本来需要插入的片段压缩成了错配。另外，match/mismatch的比值也要合理：一般match取正，mismatch取负，gap罚分大概是对应错分配分的2~3倍。你可以用BLOSUM62或PAM250这类现成的评分矩阵，比自己设靠谱得多。

还有一个细节：实际代码中，NW和SW的矩阵填充循环是完全一样的，就差在负值处理和回溯规则上。所以如果你能写一个通用动态规划函数，再暴露两个参数（allow_negative和backtrack_mode），就可以同时复用两种算法。我一般这样写：

def fill_matrix(seq_a, seq_b, score_matrix, gap_penalty, allow_negative=False):
    n, m = len(seq_a), len(seq_b)
    F = [[0] * (m+1) for _ in range(n+1)]
    for i in range(1, n+1):
        F[i][0] = F[i-1][0] + gap_penalty
    for j in range(1, m+1):
        F[0][j] = F[0][j-1] + gap_penalty
    for i in range(1, n+1):
        for j in range(1, m+1):
            match = F[i-1][j-1] + score_matrix[(seq_a[i-1], seq_b[j-1])]
            delete = F[i-1][j] + gap_penalty
            insert = F[i][j-1] + gap_penalty
            F[i][j] = max(match, delete, insert)
            if not allow_negative and F[i][j] < 0:
                F[i][j] = 0
    return F