Needleman-Wunsch算法

Needleman–Wunsch 算法（Needleman–Wunsch algorithm，NW）用动态规划（dynamic programming，DP）求两条字符串的最优全局比对（global pairwise alignment）：从第一个字符比到最后一个字符，允许插入/删除造成的空位，使整条比对在给定评分下得分最高（或编辑距离最小）。它是生物序列比对的起点，也是后续 Smith–Waterman、Gotoh、BLAST 等算法的共同祖先之一。

段末注释：全局指两端都必须吃进比对，不能只截最相似的一段；最优只相对你给定的匹配/错配/空位分数，换一套分数就换一条「最优」。

---

1. 要比的到底是什么

给两条序列 (A = a_1 a_2 \ldots a_n)、(B = b_1 b_2 \ldots b_m)（DNA / RNA / 蛋白质 / 任意字母表）。一次比对是把它们拉成等长，空位记为 (-)，且去掉空位后分别还原为 (A)、(B)。每一列只能是三种事件之一：

列类型	写法	生物学读法
匹配（match）	(a_i = b_j)	同一残基保留
错配（mismatch）	(a_i \neq b_j)	替换
空位（gap / indel）	(a_i) 对 (-)，或 (-) 对 (b_j)	插入或删除

比对的得分是各列分数之和。NW 在全部合法全局比对里取总分最优的那一条（并列最优可有多条）。

和「编辑距离」是同一张桌子的两面：最大化相似度 (S)，与最小化距离 (D)，在分数符号对偶时给出同一条路径（Sellers, 1974）。下文一律写最大化得分。

段末注释：indel = insertion or deletion，一条空位从序列 (A) 看是删除、从 (B) 看是插入，算法不区分方向。

---

2. 1970 原文和今天课堂上的 NW 不是同一份伪代码

Needleman & Wunsch（J. Mol. Biol. 1970, 48:443–453）要解决的是：两条蛋白质能否配出「最大匹配」，从而判断同源。原文允许空位不罚款（(d=0)），递推对更早的整行/整列取 (\max)，朴素实现是三次时间。Sankoff（1972）以及后来的教材，把「线性空位 + 三向邻居」收成今天人人会写的二次 DP。Wagner–Fischer（1974）在字符串校正问题上得到同一结构。

本文讲的是现代标准形：线性空位惩罚、(O(nm)) 填一张表。仿射空位（开口贵、延伸便宜）是 Gotoh（1982）的三矩阵推广，见第 8 节。

---

3. 评分：算法保证最优，不保证生物学正确

记 (s(x,y)) 为字符对分数，(d) 为线性空位分（通常 (d<0)）。一条比对的总分：

[
S = \sum_{\text{非空位列}} s(a,b) + d \cdot (#\text{空位字符})
]

教学常用：(s(x,x)=+1)，(s(x,y)=-1)（(x\neq y)），(d=-1)。蛋白质应换成替换矩阵（substitution matrix），如 BLOSUM62、PAM250；核酸常用 match/mismatch，或 NUC.4.4。

空位比错配「更贵」时，算法倾向替换而不是撕开序列；空位太便宜则会出现碎片化缺口。仿射模型 (g(\ell)=\alpha+(\ell-1)\beta)（(\alpha) 开口、(\beta) 延伸，(\lvert\alpha\rvert>\lvert\beta\rvert)）更贴近「一次断裂、一段连续 indel」的诱变图像，但那已经不是单矩阵 NW。

段末注释：BLOSUM 编号越大表示序列越近（BLOSUM80 近、BLOSUM45 远）；PAM 正好相反。默认蛋白质矩阵是 BLOSUM62。

---

4. 状态、递推、初始化

最优子结构：(A[1..i]) 与 (B[1..j]) 的最优全局比对，最后一列只可能是 ((a_i,b_j))、((a_i,-)) 或 ((-,b_j))；前缀必须已经最优。

定义 (F_{i,j}) 为 (A[1..i]) 与 (B[1..j]) 的最优得分。行对应 (A) 的前缀，列对应 (B) 的前缀。

边界（从空串一路空位走到当前前缀）：

[
F_{0,0}=0,\qquad F_{i,0}=i\cdot d,\qquad F_{0,j}=j\cdot d
]

递推（(i,j\ge 1)）：

[
F_{i,j}

\max
\begin{cases}
F_{i-1,j-1}+s(a_i,b_j) & \text{对角：对齐 }a_i\text{ 与 }b_j\
F_{i-1,j}+d & \text{向上：}a_i\text{ 对空位}\
F_{i,j-1}+d & \text{向左：}b_j\text{ 对空位}
\end{cases}
]

整表最优分在右下角 (F_{n,m})。比对本身不在格子里，而在回溯指针里：每个格子记下是谁给出了 (\max)（并列则都记，对应多条最优比对）。

段末注释：向上消耗 (A) 的一个字符、向左消耗 (B) 的一个字符；把行列对调时，「上/左」跟着对调，公式本身不变。

---

5. 数值例：GCATGCG vs GATTACA

取 (A=) GCATGCG，(B=) GATTACA，(s=+1/-1)，(d=-1)。这是教材里最常见的 7×7 例子。

5.1 填表

第一行、第一列按边界写成 (0,-1,\ldots,-7)。左上第一个实格子是 G vs G：

[
\max(0+1,;-1-1,;-1-1)=1
]

按行扫完得：

	(-)	G	A	T	T	A	C	A
(-)	0	−1	−2	−3	−4	−5	−6	−7
G	−1	1	0	−1	−2	−3	−4	−5
C	−2	0	0	−1	−2	−3	−2	−3
A	−3	−1	1	0	−1	−1	−2	−1
T	−4	−2	0	2	1	0	−1	−2
G	−5	−3	−1	1	1	0	−1	−2
C	−6	−4	−2	0	0	0	1	0
G	−7	−5	−3	−1	−1	−1	0	0

(F_{7,7}=0) 即最优全局得分。若干格子存在并列来源（例如 (F_{4,4}=1) 可来自对角或左邻），所以最优比对不唯一。

5.2 回溯

从 ((7,7)) 走回 ((0,0))。优先对角、再上、再左，得到一条最优路径：

[
(0,0)\xrightarrow{\text{对角 }G/G}
(1,1)\xrightarrow{\text{上 }C/-}
(2,1)\xrightarrow{\text{对角 }A/A}
(3,2)\xrightarrow{\text{左 }-/T}
(3,3)\xrightarrow{\text{对角 }T/T}
(4,4)\xrightarrow{\text{对角 }G/A}
(5,5)\xrightarrow{\text{对角 }C/C}
(6,6)\xrightarrow{\text{对角 }G/A}
(7,7)
]

写成比对（列分数：(+1,-1,+1,-1,+1,-1,+1,-1)）：

GCA-TGCG
G-ATTACA

总分 (0)，与 (F_{7,7}) 一致。另一条常见并列解是 GCATG-CG / G-ATTACA，分数同为 (0)。

回溯规则：

来源	(A) 写	(B) 写	指针
(F_{i-1,j-1}+s(a_i,b_j))	(a_i)	(b_j)	(i{-}1,;j{-}1)
(F_{i-1,j}+d)	(a_i)	(-)	(i{-}1,;j)
(F_{i,j-1}+d)	(-)	(b_j)	(i,;j{-}1)

段末注释：只算分、不回溯时，空间可压成两行 (O(\min(n,m)))；既要分又要比对且只要线性空间，用 Hirschberg（1975）分治，时间仍是 (O(nm))，常数大约翻倍。

---

6. 复杂度

每格 (O(1))，格数 ((n+1)(m+1))，时间、朴素空间都是 (O(nm))。两条长 (10^4) 的序列约 (10^8) 格，笔记本能算；两条染色体不行。Four Russians 可改进到 (O(nm/\log n))，实践中更常见的是 SIMD（parasail）或低分歧时的 Wavefront（WFA2），而不是改教科书递推。

---

7. 和邻居算法差在哪一句话

算法	改哪一句	适用
NW（本文）	边界强制空位罚到底；回溯从 ((n,m)) 到 ((0,0))	全长同源、长度接近
半全局（semi-global / overlap）	末端空位 (d=0)；从最后一行/列的最大格起回溯	读段对参考、引物、重叠群
Smith–Waterman（1981）	递推加 (\max(\ldots,0))；从全局最大格回溯到 (0)	局部同源、结构域、BLAST 的精确原型
Gotoh（1982）	三张表 (M,I,D) 区分开口/延伸	几乎所有生产级比对器的默认空位

Gotoh 的仿射递推（相似度形式，记号随教材略有出入）：

[
\begin{aligned}
M_{i,j} &= s(a_i,b_j)+\max(M_{i-1,j-1},,I_{i-1,j-1},,D_{i-1,j-1})\
I_{i,j} &= \max(M_{i-1,j}+\alpha,,I_{i-1,j}+\beta)\
D_{i,j} &= \max(M_{i,j-1}+\alpha,,D_{i,j-1}+\beta)
\end{aligned}
]

时间仍是 (O(nm))，只是每格常数变成约 3 倍。Biopython 的 PairwiseAligner(mode='global') 在开口≠延伸时走的就是 Gotoh，不是 1970 原文。

段末注释：长度差很大时硬跑 NW，短序列会被一长串末端空位「拉满」整条长序列，分数和生物学都无意义；改半全局或局部。

---

8. 实现：先手写线性空位，再生产用社区库

社区方案优先（应用场景 / 风险一并写）：

工具	场景	实现要点	风险
手写 NW（下面代码）	教学、改评分、单测对照	单矩阵线性空位	慢；不要当生产比对器
Biopython Bio.Align.PairwiseAligner	交互、<10 kb、脚本	C 后端 Gotoh；mode='global'	默认空位分为 0，配 BLOSUM 会疯狂插空；必须显式设 open_gap_score / extend_gap_score。Bio.pairwise2 已弃用
EMBOSS needle	要可审计的固定默认参数	经典全局比对 CLI	无 SIMD，批量慢
parasail（Daily, 2016）	高通量 NW/SW	SSE/AVX	优势在长序列批循环里才明显
edlib	只要编辑距离、DNA	Myers 位并行	不含替换矩阵
WFA2 / pywfa	长且分歧 <5%	波前，内存随差异长	远同源会退回 Gotoh 复杂度
BLAST / MMseqs2	库搜索、全对全	种子–延伸启发式	不保证全局最优；不能拿来当 NW 的替代证明

8.1 线性空位 NW（教学完整版）

from typing import List, Tuple


def needleman_wunsch(
    seq_a: str,
    seq_b: str,
    match: int = 1,
    mismatch: int = -1,
    gap: int = -1,
) -> Tuple[str, str, int]:
    """计算两条序列的一条最优全局比对（线性空位）。

    输入:
        seq_a, seq_b: 待比对字符串。
        match, mismatch, gap: 匹配分、错配分、线性空位分（gap 一般为负）。
    输出:
        (aligned_a, aligned_b, score): 一条最优比对及其总分。
    处理逻辑:
        填 (n+1)×(m+1) 得分表；回溯优先对角、再上、再左，并列时只取一条。
    """
    n, m = len(seq_a), len(seq_b)
    score_of = lambda x, y: match if x == y else mismatch

    f: List[List[int]] = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(n + 1):
        f[i][0] = i * gap
    for j in range(m + 1):
        f[0][j] = j * gap

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            f[i][j] = max(
                f[i - 1][j - 1] + score_of(seq_a[i - 1], seq_b[j - 1]),
                f[i - 1][j] + gap,
                f[i][j - 1] + gap,
            )

    aligned_a, aligned_b = [], []
    i, j = n, m
    while i > 0 or j > 0:
        if (
            i > 0
            and j > 0
            and f[i][j] == f[i - 1][j - 1] + score_of(seq_a[i - 1], seq_b[j - 1])
        ):
            aligned_a.append(seq_a[i - 1])
            aligned_b.append(seq_b[j - 1])
            i -= 1
            j -= 1
        elif i > 0 and f[i][j] == f[i - 1][j] + gap:
            aligned_a.append(seq_a[i - 1])
            aligned_b.append("-")
            i -= 1
        else:
            aligned_a.append("-")
            aligned_b.append(seq_b[j - 1])
            j -= 1

    return "".join(reversed(aligned_a)), "".join(reversed(aligned_b)), f[n][m]


if __name__ == "__main__":
    a, b, s = needleman_wunsch("GCATGCG", "GATTACA")
    print(a)  # GCA-TGCG
    print(b)  # G-ATTACA
    print(s)  # 0

蛋白质把 score_of 换成查表即可，递推不动。

8.2 生产路径：Biopython（显式空位）

from Bio.Align import PairwiseAligner, substitution_matrices

aligner = PairwiseAligner(mode="global", match_score=1, mismatch_score=-1)
aligner.open_gap_score = -1
aligner.extend_gap_score = -1  # 等于 open 时退化为线性空位，接近手写 NW
print(aligner.algorithm)        # 应看到 Needleman-Wunsch 或 Gotoh global ...

aln = next(aligner.align("GCATGCG", "GATTACA"))
print(aln.score)  # 0

prot = PairwiseAligner(
    mode="global",
    substitution_matrix=substitution_matrices.load("BLOSUM62"),
    open_gap_score=-11,
    extend_gap_score=-1,
)

aligner.score(a, b) 只算分、更快。最优比对可能极多，设 aligner.max_alignments 防止爆内存。

---

9. 什么时候不该用 NW

• 长度差大、只关心一段同源：SW 或半全局。
• 一条 query 对整个库：BLAST、DIAMOND、MMseqs2；NW 是 (O(nm)) 精确器，不是搜索引擎。
• 蛋白质同一度 <25%：序列 DP 进入暮光区，改 profile（HHsearch / HMMER）或结构（Foldseek）。
• 低复杂度 / 串联重复：未 mask 时分数虚高。
• 把「能出比对」当成「同源」：任何两条串 NW 都会吐出一条全局比对；同源要靠 bit score / E-value / 置换检验，不是靠路径存在。

同一条比对的百分同一度至少有四种分母（含不含内部空位、较短条、平均长），同一对序列可差十几个百分点。写报告时必须写清定义。

---

10. 要点

1. NW = 全局 + DP。状态是双前缀最优分，转移只有对角/上/左。
2. 课堂上的二次递推是 1970 之后定型的线性空位版；原文更一般、也更慢。
3. 分数是先验。最优保证的是算术，不是进化。
4. 要仿射空位用 Gotoh，要局部用 SW，要末端免费用半全局——改边界或加 (\max(0,\cdot))，不是另起一套哲学。
5. 手写只为看清指针；批量与长序列用 parasail / WFA2 / EMBOSS，库搜索用 BLAST 系。