第二季 · 第 51 章

寻找生命的相似之处

《算法的故事》第二季 · 第五篇：破解生命的密码——从桑格的胰岛素、戴霍夫的蛋白质图谱，到Needleman–Wunsch与Smith–Waterman序列比对算法

1950年代，英国剑桥。一位生物化学家正面对一项困难得令人难以想象的工作。他想弄清楚一种蛋白质的内部结构。这种蛋白质并不陌生。它叫作胰岛素。人们已经知道，胰岛素与人体血糖调节密切相关。但是，要真正理解这种蛋白质，科学家还需要回答一个更加基础的问题。组成胰岛素的氨基酸，究竟按照什么顺序连接在一起？一种蛋白质包含许多氨基酸。它们不像一箱可以随意调换位置的零件。连接顺序的变化，可能改变蛋白质的结构、性质与功能。因此，科学家不能只知道胰岛素里面有哪些氨基酸。还必须知道：第一个是什么，第二个是什么，第三个又是什么？这位研究者叫作：弗雷德里克·桑格（Frederick Sanger）。经过多年研究，他与同事终于确定了胰岛素完整的氨基酸序列及其链间连接关系。1955年前后，这项工作已经完整建立。1958年，桑格因蛋白质结构，特别是胰岛素结构方面的工作，获得诺贝尔化学奖。 ^{[PMC][MRC Laboratory of Molecular Biology]}对于分子生物学来说，这项成就不仅意味着一种激素的结构终于得到揭示。它还让人们更加清楚地认识到：生命中的某些重要分子，可以被表示成一串具有严格顺序的符号。就像一句话。就像一段程序。也像一种使用特殊字母书写的语言。随着越来越多蛋白质的序列被测定，一项新的研究问题逐渐出现。假设我们已经知道两种蛋白质的氨基酸序列。它们看起来有些相似。有些位置完全相同。有些位置换成了另一种氨基酸。还有一些位置，似乎多出或者缺少了某个氨基酸。那么：我们怎样判断，两条序列究竟有多相似？如果它们长度不同，又该怎样比较？如果某些位置发生过插入和删除，能不能仍然找出对应关系？更进一步：我们能不能利用这些对应关系，帮助研究不同蛋白质之间可能存在的共同进化来源？这道问题，终于从生物化学走向计算机科学。1970年，两位在美国芝加哥工作的研究者发表了一篇论文。A General Method Applicable to the Search for Similarities in the Amino Acid Sequence of Two Proteins《一种适用于寻找两种蛋白质氨基酸序列相似性的一般方法》。作者是：索尔·尼德曼（Saul B. Needleman）。以及：克里斯蒂安·文施（Christian D. Wunsch）。

他们提出：将两条蛋白质序列放入一个二维矩阵。让矩阵中的每个位置，表示两条序列里一对氨基酸的比较。然后，在这个矩阵中寻找得分最高的合法路径。这项方法后来被称为：Needleman–Wunsch算法。又过了十一年。1981年，两位研究者进一步提出：如果两条很长的序列整体并不相似，但其中各有一小段高度相似，该怎么办？是不是应该允许算法放弃那些完全不相关的部分，只寻找最值得比较的局部片段？他们是：坦普尔·史密斯（Temple F. Smith）。迈克尔·沃特曼（Michael S. Waterman）。他们的研究被称为：Smith–Waterman算法。它与Needleman–Wunsch算法一起，构成现代序列比对理论中最经典的两条道路。一条追求整体对应。另一条寻找最有意义的局部相似片段。两者共同使用了一种我们已经多次见过的数学工具：动态规划（Dynamic Programming）。不过，这一次，算法不再寻找道路、通信路径或资源分配方案。它要在生命序列的字符之间，寻找可能被漫长进化过程掩盖的对应关系。

第一幕：1955年，生命第一次能够被完整地写成一串氨基酸字母

今天，我们已经很习惯用字母表示生命信息。例如，一段DNA序列可以写成：

ACGTTGCA

一条蛋白质序列可以写成：

MALWMRLLPLLALLALWGPDPAAA

但在20世纪中叶，这些看起来非常普通的字符序列，并不是科学家一开始就已经完整掌握的资料。

蛋白质由什么构成？

蛋白质由氨基酸残基连接而成。生物体的常见蛋白质使用20种标准氨基酸作为主要构建单元。为了便于记录，研究者可以用单个英文字母代表不同氨基酸。例如：

字母	氨基酸
A	丙氨酸（Alanine）
C	半胱氨酸（Cysteine）
G	甘氨酸（Glycine）
L	亮氨酸（Leucine）
M	甲硫氨酸（Methionine）
S	丝氨酸（Serine）
T	苏氨酸（Threonine）
V	缬氨酸（Valine）

这样，一条蛋白质便可以表示为：

S=s_1s_2\cdots s_n

其中每个：

s_i

是一种氨基酸残基的符号。不过，必须特别强调：蛋白质序列是研究蛋白质结构与功能的重要基础，却不是蛋白质三维结构本身。蛋白质怎样折叠、与哪些分子相互作用，以及最终具有怎样的功能，还涉及其他层次的信息。

桑格为什么选择胰岛素？

胰岛素是一种已知具有重要生理作用的蛋白质。它也相对较小。成熟胰岛素由两条多肽链构成。A链有21个氨基酸残基。B链有30个氨基酸残基。合计：

21+30=51

个残基。两条链之间由二硫键连接，A链内部还存在一个二硫键。桑格团队需要逐步确定：每条链的氨基酸顺序。不同片段之间的重叠关系。以及链与链之间的连接方式。这项工作不是利用今天的测序仪，一次性输出51个字符。他们必须使用化学试剂、蛋白质水解、纸色谱、电泳以及其他实验方法，逐步识别片段，再将线索拼接起来。这是一项极其艰苦的实验研究。^[PMC]

但是，当两条蛋白质序列都被测出来以后，新的问题才刚刚开始

假设两种动物体内的某种蛋白质，已经分别测得氨基酸序列。科学家希望知道：它们是否具有较高的序列相似性？某些对应位置上的氨基酸变化，是否与进化过程有关？是否有一些位置，长期保持不变？这些都是重要的研究问题。但仅靠肉眼观察，难以稳定处理越来越长、越来越多的序列。于是，分子生物学开始需要一种新能力。不仅要能测定序列。还要能比较序列。这便是算法登场的时刻。

第二幕：1965年，玛格丽特·戴霍夫开始建立生命的“字典”

1950年代末到1960年代，科学家能够测定的蛋白质序列逐渐增加。但是，这些序列分散在不同研究论文和实验报告中。每一篇论文可能只研究一种或少量蛋白质。如果研究者想比较不同动物、不同组织或不同功能的蛋白质，便需要查阅许多资料。而当数据越来越多时，人工整理会变得非常困难。这时候，出现了一位后来被称为生物信息学先驱的研究者。玛格丽特·奥克利·戴霍夫（Margaret O. Dayhoff）。

她做了一件当时并不显眼，却十分重要的工作

戴霍夫与Richard Eck、Robert Ledley等研究者合作，将散落在科学文献中的蛋白质序列资料，逐渐整理成可供计算分析的集合。1965年，他们出版了：Atlas of Protein Sequence and Structure《蛋白质序列与结构图谱》。它是早期具有里程碑意义的蛋白质序列资料集，也是计算机辅助序列整理发展的重要成果。^[PubMed]

1965年《蛋白质序列与结构图谱》封面。戴霍夫团队的工作帮助蛋白质序列从分散的实验报告，逐步转变为可以集中比较和计算的数据。

为什么建立序列资料库如此重要？

因为当研究者只有两条序列时，人工比较也许还勉强可行。但如果有：100条。1000条。甚至数百万条序列。问题就完全不同。首先，我们需要可靠地保存这些序列。然后，需要使用一致的符号表示方式。还要处理来源、物种、功能和测定资料。最后，需要算法在其中寻找可能的对应关系。这实际上形成了现代生物信息学的一条完整工作链。

图51-1：从实验测序到算法比较

实验测定获得蛋白质或核酸的分子序列

资料整理使用统一的符号和来源信息保存序列

序列比较寻找相同、替换、插入和缺失的可能对应关系

生物学解释结合统计证据与其他资料研究进化、结构和功能

戴霍夫还留下了一项重要思想：不同氨基酸替换，不一定具有同样的意义

假设两个蛋白质对应位置上，一个是：

A

另一个是：

V

这两种氨基酸虽然不同，但某些理化性质相似。又假设一个位置从：

W

变成：

G

这种替换可能在某些情境下具有不同的生化影响。因此，序列比较并不一定适合永远使用：相同就加1。不同就扣1。更精细的蛋白质序列分析，可以使用氨基酸替换评分矩阵。戴霍夫及其合作者后来对氨基酸替换概率和PAM矩阵的发展作出了重要贡献。我们将在本章后面进一步讨论。但到1960年代末，最迫切的问题仍然是：如果允许两条序列在某些位置插入空位，怎样在大量可能的排列中寻找最佳对应？这个问题，把我们带到1969年的芝加哥。

第三幕：1969年7月21日，一篇来自芝加哥的论文递交了期刊

1969年7月21日。一份关于蛋白质序列比较的研究论文，送达《Journal of Molecular Biology》。论文的作者是：Saul B. Needleman。Christian D. Wunsch。两人的工作单位包括：美国西北大学生物化学系。以及芝加哥的一所退伍军人事务研究医院。这项细节可以直接从原始论文首页核查。^{[doc.aporc.org]}

他们面对的困难，究竟是什么？

假设有两条蛋白质序列：

A B C D E F

以及：

A B D E F

两条序列的长度只相差一个字符。如果我们坚持逐位直接比较：

A B C D E F
A B D E F

从第三位开始，对应关系就会发生偏移。但如果我们允许第二条序列插入一个空位：

A B C D E F
A B - D E F

我们就能看到：除了一处插入或删除所对应的缺口，其余字符能够很好地对齐。这项操作叫作：序列比对（Sequence Alignment）。它不是改变真实序列的排列顺序。而是在适当位置加入表示对应缺失的：Gap（空位）。使两条序列可以在同一组列中进行比较。

为什么不能随意重新排序？

因为序列的顺序具有生物学意义。如果一条序列是：

A C G T

我们不能仅仅因为另一条序列拥有同样的字符，就把它们随意交换位置。例如：

T G C A

虽然使用了同样四种字符，但顺序完全不同。序列比对要求：保留每条原始序列内部的先后顺序。它允许在某些位置加入空位。但不允许为了追求高分而任意重新排列字符。这便是算法所必须遵守的第一项结构约束。

暴力尝试所有空位位置，会有多困难？

假设两条序列各有100个字符。为了进行完整比对，我们可以在两条序列中加入各种不同数量和位置的空位。这些可能性组合会迅速增长。如果直接列举每一种比对方式，再计算得分，计算量可能非常巨大。因此，Needleman与Wunsch的核心问题就是：能否不枚举全部候选，而是通过某种矩阵结构，高效排除那些不可能形成最佳结果的路径？

原始论文的真正贡献

1970年3月28日，论文正式发表。第48卷第3期。第443—453页。原文在摘要中明确说明：方法将两条蛋白质序列中可能的氨基酸配对，组织到二维数组中。不同序列对应关系，可以转化成数组中的合法路径。再通过累积最优匹配值，寻找最佳对应关系。^{[ScienceDirect][doc.aporc.org]}原始论文首先使用一种简单的评分：相同氨基酸记1分。不同氨基酸记0分。并且讨论了可以加入Gap惩罚与更加复杂的评分函数。^{[doc.aporc.org]}这点很重要。因为今天我们常在教材里看到的三方向递推公式，虽然准确体现了Needleman–Wunsch型全局比对的动态规划思想，但并不是对1970年论文中每一项原始计算细节的逐字复刻。原始论文使用更一般的矩阵路径和跳跃关系，讨论允许不同长度的缺口。后来，Sellers、Waterman、Gotoh等研究者又发展了不同的形式化距离模型和更高效的计算方法。因此，我们将清楚区分：1970年真实提出的矩阵最大匹配思想。以及：今天最常用于教学的线性Gap罚分三方向动态规划实现。只有这样，才能同时尊重历史与数学。

第四幕：鲸鱼的肌红蛋白，为什么会与人类的血红蛋白放在一起比较？

这一次，我们不使用虚构的两条蛋白质作为历史主线。直接翻开1970年论文的实验部分。Needleman与Wunsch选择了两组具有不同研究目的的蛋白质。第一组：鲸的肌红蛋白（Myoglobin）。以及：人的血红蛋白β链（β-Hemoglobin）。第二组：牛胰腺核糖核酸酶（Ribonuclease）。以及：鸡蛋清溶菌酶（Lysozyme）。他们希望通过这些比较，研究序列评分方法能否区分预期具有共同进化关系的蛋白质，以及没有相同预期的蛋白质组合。^{[doc.aporc.org]}

为什么选择肌红蛋白与血红蛋白？

肌红蛋白主要参与肌肉中的氧结合和储存。血红蛋白则是血液中的重要氧运输蛋白。两者功能并不完全相同。但它们都与氧结合相关，也具有重要的结构联系。因此，研究者有理由预期它们的序列包含反映共同演化来源的相似性。如果算法能够识别这种相似性，就有机会为分子进化研究提供证据。但是，仅仅找到一些相同字符还不够。因为即使随机生成两条蛋白质序列，也可能碰巧包含相同的氨基酸。尤其在序列很长时，随机匹配并不罕见。

作者怎样判断匹配是否具有意义？

他们不仅计算真实蛋白质的匹配分数。还把其中一条蛋白质的氨基酸顺序随机打乱。然后重新进行序列比较。并将随机序列产生的匹配得分，与真实序列的匹配得分进行比较。这是一个非常重要的实验设计。因为研究者希望区分：真正反映序列组织关系的高匹配分数。以及仅仅由氨基酸组成和偶然匹配造成的高分。原始论文记录，他们在实验中对相关序列进行了十次随机重排比较，并估计了随机匹配分数的均值与标准差。作者还明确承认：随机样本只有十份，是受限于当时可用计算时间。因此，相关统计估计存在不小的不确定性。^{[doc.aporc.org]}这项细节非常值得保留。因为它让我们看见1970年的真实计算条件。今天，在一台普通计算机上进行十次短序列比对，可能只需要极短时间。但在当时，反复处理蛋白质序列、生成随机对照并比较得分，已经涉及不容忽视的计算资源。

这与算法有什么关系？

它提醒我们：序列比对有两个不同层次的问题。第一个是计算问题：在给定评分规则下，怎样找到得分最高的比对？第二个是科学解释问题：这个最高得分，究竟是否具有生物学或统计学意义？动态规划主要解决第一个问题。第二个问题还需要随机背景模型、统计分析、蛋白质结构、进化关系及其他生物学证据。尤其要避免一个常见误解：序列相似性（similarity）不等于已经证明共同祖先关系（homology）。在现代生物学的严格用法中，同源性表示共同进化来源，是一种关系判断。它不是简单的“相似百分比”。我们可以说：两条序列有70%的序列一致性。但不应直接说：“它们具有70%的同源性。”因为共同祖先关系不能像长度比例那样定义成一个连续的百分比。这一点，对于后面学习BLAST和系统发育树都非常重要。

第五幕：四个字母与三个字母，为什么能构成一张数学地图？

现在，我们开始真正构造序列比对算法。先暂时放下复杂的蛋白质。使用两条非常短的DNA教学序列。第一条：

\boxed{A=\text{ACGT}}

第二条：

\boxed{B=\text{AGT}}

可以发现，第二条序列比第一条少一个C。如果希望完整比较两条序列，一种自然的比对方式是：

A C G T
A - G T

这里：-表示空位。整个比对共有四列。但为什么这就是最佳比对？如果我们改变评分规则，结果还会一样吗？如果有许多可能的空位位置，怎样确保找到得分最高的一种？这正是动态规划要解决的问题。

第一步：建立评分规则

为了方便计算，我们约定：

比对情况	得分
字符相同	+2
字符不同	−1
一个字符与Gap对齐	−2

例如：

s(A,A)=2

s(A,G)=-1

而一个字符与空位对齐，扣2分。于是：

A C G T
A - G T

的得分为：

2+(-2)+2+2

即：

\boxed{4}

注意，这是一套教学评分规则。它不是DNA变化概率的完整生物学模型。

第二步：定义动态规划状态

令：

F(i,j)

表示：第一条序列前

i

个字符，与第二条序列前

j

个字符进行全局比对时，能够获得的最高分数。例如：

F(1,1)

表示：第一条序列的前一个字符A，与第二条序列的前一个字符A，比对所得的最优分数。所以：

F(1,1)=2

而：

F(4,3)

表示完整的：ACGT和：AGT之间，能够达到的最佳全局比对分数。我们最终希望得到的就是：

\boxed{F(4,3)}

第三步：为什么每个格子只需要考虑三个方向？

想象我们正在计算：

F(i,j)

也就是两条序列各自前一部分的最佳对应关系。最后一列比对，只可能属于以下三种情况。

图51-2：每个比对状态的三种前驱选择

左上方两个字符对齐A / A

上方A的字符对GapC / –

左方Gap对B的字符– / G

如果最后一列对齐两个真实字符，那么此前对应：

F(i-1,j-1)

如果最后一列是A中的字符，对齐B里的空位，那么此前对应：

F(i-1,j)

如果最后一列是A里的空位，对齐B中的字符，那么此前对应：

F(i,j-1)

因此，所有可能的完整比对，按照最后一列的类型，都能被划分到这三种情况。这正是动态规划的关键。我们不必列举所有完整比对。只需要知道三个较短前缀的最佳结果。

第四步：递推公式

于是，得到经典的线性Gap罚分全局比对递推：

\boxed{ F(i,j)=\max \begin{cases} F(i-1,j-1)+s(a_i,b_j)\\ F(i-1,j)-d\\ F(i,j-1)-d \end{cases} }

其中：

d>0

是每个Gap字符的罚分。本例：

d=2

而：

s(a_i,b_j)= \begin{cases} 2,&a_i=b_j\\ -1,&a_i\ne b_j \end{cases}

注意，我们用的是最大化得分。所以Gap罚分写成：

-d

这与上一章某些路径算法采用的“最小化代价”表示形式不同，但数学结构是一样的。

第五步：边界应该怎样初始化？

如果第一条序列为空，而第二条序列包含三个字符：AGT那么，要进行完整的全局比对，只能是：

- - -
A G T

需要三个Gap。每个Gap扣2分。所以：

F(0,3)=-6

同样：

F(4,0)=-8

因此，初始化为：

\boxed{F(i,0)=-2i}

\boxed{F(0,j)=-2j}

以及：

\boxed{F(0,0)=0}

这意味着：全局比对不允许随意跳过序列开头或结尾，所有字符都必须在最终比对中获得一个位置。

第六幕：一格一格填完矩阵，最佳答案竟然出现在右下角

现在，我们将刚才的公式用于：

A=\text{ACGT}

B=\text{AGT}

第一个格子：比较A与A

F(1,1)

有三个候选。左上方：

F(0,0)+2=2

上方：

F(0,1)-2=-4

左方：

F(1,0)-2=-4

因此：

\boxed{F(1,1)=2}

第二个字符C，遇到第二条序列中的G

现在考虑：

F(2,2)

对应：AC与：AG左上方的两个字符：C与G。不相同。因此：

F(1,1)-1=1

上方：

F(1,2)-2=0-2=-2

左方：

F(2,1)-2=0-2=-2

所以：

\boxed{F(2,2)=1}

关键位置：G与G

当我们计算：

F(3,2)

时，对应：ACG与：AG最后一位G与G相同。左上方：

F(2,1)+2=0+2=2

上方：

F(2,2)-2=1-2=-1

左方：

F(3,1)-2=-2-2=-4

因此：

\boxed{F(3,2)=2}

最后一个位置：T与T

最终：

F(4,3)

对应：ACGT与：AGTT与T相同。左上方：

F(3,2)+2=2+2=4

上方：

F(3,3)-2=0-2=-2

左方：

F(4,2)-2=0-2=-2

所以：

\boxed{F(4,3)=4}

完整矩阵如下。

图51-3：Needleman–Wunsch型全局比对矩阵

	–	A	G	T
–	0	-2	-4	-6
A	-2	2	0	-2
C	-4	0	1	-1
G	-6	-2	2	0
T	-8	-4	0	4

行对应序列ACGT的前缀，列对应AGT的前缀。采用匹配+2、不匹配−1、每个Gap−2。右下角4是最佳全局比对得分。

但得分4，还不是完整的比对结果

我们已经知道：最佳得分是4。但怎样知道空位应该放在哪里？这就需要另一项动态规划技巧：回溯（Traceback）。从矩阵右下角开始。查看当前格子的最佳分数来自哪个方向。如果来自左上：说明两个真实字符被对齐。如果来自上方：说明A的一个字符对齐了Gap。如果来自左方：说明Gap对齐了B的一个字符。

从右下角逆向回溯

首先：

F(4,3)=4

来自：

F(3,2)+2

因此：T与T对齐。接着：

F(3,2)=2

来自：

F(2,1)+2

因此：G与G对齐。再看：

F(2,1)=0

来自：

F(1,1)-2

因此：C与Gap对齐。最后：

F(1,1)=2

来自：

F(0,0)+2

因此：A与A对齐。我们从后往前得到：

T / T
G / G
C / -
A / A

反转顺序：

A C G T
A - G T

于是，算法不仅算出了最优分数：

\boxed{4}

也恢复出了对应的最佳全局比对。

这与第49章的维特比算法有什么联系？

上一章，维特比算法在有限状态的网格图中，保存抵达每个状态的最佳前驱。本章，序列比对算法在二维矩阵中，保存抵达每个前缀状态的最佳前驱。一个处理：随时间变化的隐藏状态。另一个处理：两条序列的前缀位置。但它们都使用：最优子结构、状态复用与前驱回溯。这就是动态规划为什么能够跨越不同领域的原因。

第七幕：1981年，为什么有人认为“整条序列都对齐”反而会掩盖真正的发现？

十一年过去了。时间来到1981年。那时，分子生物学研究者已经逐渐积累了更多的序列资料。但一个新的问题越来越明显。并不是所有值得比较的蛋白质或核酸序列，都应该从第一个字符一直比对到最后一个字符。

假设两条长序列中间拥有一段共同片段

第一条：

T T A C G T A A

第二条：

G G A C G T T

如果从第一位开始比较：T与G不同。第二位仍然不同。如果要求整条序列都进入最终比对，那么前后不相关的区域会引入额外的失配与Gap罚分。但如果我们只看中间：

A C G T
A C G T

就会发现一段完全相同的局部区域。它们在各自序列中的位置并不相同。而且两端还存在不相关的字符。可是，这段共同片段，可能具有值得进一步研究的生物学意义。

1970年的方法是否完全不能发现局部相似性？

不是。这一点需要澄清。Needleman与Wunsch的1970年原始研究，并不只是今天教学里那种固定的全局三方向递推。它讨论过矩阵路径、最大匹配和不同Gap安排。而序列比对方法在1970年代也经历了持续发展。因此，不能简单地将1970年算法描述成“完全不知道局部相似区域”。1981年Smith与Waterman的突破，是在已有序列距离与动态规划研究基础上，给出了能够明确优化最佳局部片段得分的递推方法。它允许最优路径在序列内部开始。也允许在内部结束。而不必被迫覆盖两条完整序列。

史密斯与沃特曼的真实论文

1981年3月25日，《Journal of Molecular Biology》发表：Identification of Common Molecular Subsequences《共同分子子序列的识别》。作者：Temple F. Smith。Michael S. Waterman。第147卷第1期。第195—197页。整篇正文只有三页。^[PubMed]

论文开头直接讨论：在长分子序列中，寻找最大程度相似的局部片段，是一个重要的问题。作者回顾了：Needleman与Wunsch在1970年的矩阵方法。Sellers在1974年发展的序列距离度量。以及Waterman与合作者在1976年对不同长度插入、删除的研究。他们并不是从零开始建立全部序列比对数学。而是在这些研究基础上，提出了更加明确的局部最优片段识别方法。^{[布朗大学计算机科学系]}

1974年的Sellers为什么值得记住？

Peter H. Sellers研究了使用最少编辑操作将一条序列转变成另一条序列的数学方法。例如：一个字符替换。插入一个字符。删除一个字符。这些操作可以拥有不同的代价。于是，我们可以问：从序列A变成序列B，最低需要多少代价？这与最大化相似得分，看起来是两个不同问题。但在适当的评分参数关系下，它们可以互相转换。Smith、Waterman与Walter Fitch在1981年的另一篇论文中，明确研究了Needleman–Wunsch与Sellers方法之间的等价条件。^[PubMed]因此，序列比对的历史还包含一项数学上的演进：从寻找最大匹配。到定义序列距离。再到局部片段的最优相似性。

Smith–Waterman真正改变了哪里？

它没有推翻动态规划。反而保留了动态规划。也没有改变序列必须保持顺序这一基本约束。最关键的变化是：允许某个动态规划状态从零重新开始。这个看似微小的修改，会让整个优化目标发生根本变化。下一幕，我们就从这一项变化出发，亲手推导局部比对算法。

第八幕：只在递推公式里增加一个0，算法为什么就从全局走向局部？

回忆上一幕中的经典全局比对公式：

F(i,j)= \max \begin{cases} F(i-1,j-1)+s(a_i,b_j)\\ F(i-1,j)-d\\ F(i,j-1)-d \end{cases}

我们必须从矩阵左上角开始。最后在右下角结束。因为算法负责覆盖两条完整序列。但是，如果我们希望只找到两条序列内部最相似的一个片段，就需要允许路径从中间开始。Smith–Waterman的核心做法可以表示为：

\boxed{ H(i,j)=\max \begin{cases} 0\\ H(i-1,j-1)+s(a_i,b_j)\\ H(i-1,j)-d\\ H(i,j-1)-d \end{cases} }

请注意第一项：

\boxed{0}

它就是关键。

为什么要允许得分重新变成0？

假设我们正在寻找一段相似的局部片段。但前面已经连续遇到许多不匹配字符。累计得分变成了：

-5

现在，我们来到一段全新的、可能高度相似的区域。如果必须保留前面的负5分，那么后面的好匹配就必须先补偿过去的损失。但如果我们真正的目标是：找出最佳局部相似区域。那么，前面那些完全不相关的字符，不应该被强制加入这段局部比对。因此，只要某条延伸路径的得分低于零，就可以选择：从这里重新开始。将得分设为0。这是一个非常简洁却深刻的思想。对于局部匹配来说，如果过去已经成为负担，最佳选择可能是放弃过去。当然，这只是动态规划问题中的数学操作，并非对生物进化过程的真实模拟。

边界条件也改变了

对于Smith–Waterman局部比对：

\boxed{H(i,0)=0}

\boxed{H(0,j)=0}

所有第一行与第一列都设为0。因为我们不要求从两条序列的第一个字符开始。

最终答案又在哪里？

全局比对从右下角：

F(n,m)

读取最佳分数。但局部比对不一样。最佳片段可能结束在矩阵中的任何位置。因此，我们需要寻找：

\boxed{ \max_{0\le i\le n,\ 0\le j\le m}H(i,j) }

也就是整张矩阵里的最大值。然后，从最大值所在位置开始回溯。当回溯遇到：

H(i,j)=0

就停止。因为那意味着：局部最优片段已经到达它应该开始的位置。这正是1981年原始论文描述的方法。论文在递推式中包含零值，并明确规定：先找到矩阵最大元素，再沿最优前驱回溯，直到遇到得分为零的元素。^{[布朗大学计算机科学系]}

两种算法的关键区别

图51-4：全局比对与局部比对

比较项目	Needleman–Wunsch型	Smith–Waterman型
目标	完整序列的最优对应	最佳局部片段
第一行、列	累计Gap罚分	全部为0
递推候选	三个方向	三个方向加0
回溯起点	右下角	矩阵最大分值
回溯终点	左上角	第一个零分状态

再看之前的局部比对例子

第一条：

T T A C G T A A

第二条：

G G A C G T T

采用：匹配：

+2

不匹配：

-1

Gap：

-2

那么，共同片段：

A C G T
A C G T

得分为：

2+2+2+2

即：

\boxed{8}

这就是当前例子里的最佳局部比对分数。而且，这段局部相似区域并没有位于两条序列的相同绝对位置之后再完整延伸到末尾。它只是两条较长序列内部的一段高度一致片段。

为什么一个0就能发挥这样的作用？

因为它改变了动态规划状态的语义。在全局比对中：

F(i,j)

表示两条前缀完整覆盖的最优得分。而在局部比对中：

H(i,j)

表示所有恰好在位置 $i,j$ 结束的局部对齐片段，能够获得的最佳非负得分。它不是简单地在整个矩阵上追求“匹配越多越好”。而是针对每一个可能的局部终点，计算最佳片段。最终再从这些终点中选出最大者。这才是Smith–Waterman的真正数学含义。

第九幕：互动实验室——亲手改变序列，看最优路径怎样出现

现在，我们将两套算法放进同一张可以操作的动态规划矩阵。你可以选择：短序列中的一个插入或缺失。两条较长序列之间的共同局部片段。或者两条包含更多替换的DNA序列。然后切换：全局比对。局部比对。还可以修改Gap罚分。每改变一个条件，矩阵与最优路径都会重新计算。最重要的是，我们并不直接跳到最终答案。你可以逐格填表。观察每个状态依赖哪些已经算好的位置。最后再显示回溯路径。

实验一：全局比对与局部比对共用一张矩阵

图51-5：序列比对动态规划实验台使用匹配 +2、不匹配 −1、可调线性Gap罚分。点击“下一格”，按从上到下、从左到右的顺序填表。

选择序列

选择比对模式

每个Gap字符的罚分-2

序列AACGT

序列BAGT

已计算内部格子0 / 12

橙色：刚计算的格子。绿色：全部完成后的一条最优回溯路径。负分允许出现在全局模式；局部模式不会保留负分。

建议亲手完成四次实验

**第一次：**选择短序列、全局比对，直接计算全部。你会得到：

A C G T
A - G T

最佳得分：

4

这是我们刚才手工推导的结果。**第二次：**保持相同序列，改成局部比对。比较全局与局部的结果。因为局部算法不必覆盖全部字符，它可能找到得分相同但片段更短的最佳结果，也可能出现多个同分最优片段。**第三次：**选择“局部片段”，使用Smith–Waterman模式。你会看到：

A C G T
A C G T

最佳得分：

8

它能忽略两端不相关的部分。**第四次：**把Gap罚分从：

-2

改成：

-4

观察最优路径是否发生变化。这项实验会让你认识到：最佳序列比对不是单纯由两条字符串决定的，还取决于我们如何定义不同编辑事件的分数和代价。这也是生物信息学中，选择合适评分模型非常重要的原因。

第十幕：1982年，Gotoh重新思考一个问题——为什么连续缺失十个字符，不应等同于十次独立缺失？

到目前为止，我们一直使用非常简单的Gap罚分。每插入一个空位：

-2

如果连续插入三个空位：

-6

这种方式叫作：线性Gap罚分（Linear Gap Penalty）。它的数学形式是：

\boxed{g(L)=dL}

其中：

L

表示连续空位的长度。

d

表示每个空位字符的代价。但研究者逐渐意识到：在生物序列中，连续的插入或删除，有时可以来自一次长度较长的进化事件。这与多个彼此独立、相距较远的插入或删除事件，并不一定具有相同的生物学解释。

一个简单的例子

假设第一条序列：

A C G T A C G T

第二条序列：

A C - - - - G T

中间有一段连续的四字符空位。如果使用线性罚分：

d=2

那么：

g(4)=8

现在考虑另一种情况。如果四个空位分散在四个不同区域，按照同样的线性罚分，总代价还是：

4\times2=8

但从某些进化模型的角度看：一次连续删除四个残基。与四次独立删除一个残基。未必应该支付完全相同的代价。这就引出了：仿射Gap罚分（Affine Gap Penalty）。

一次打开Gap，后面延长Gap

我们定义：

\boxed{g(L)=g_o+(L-1)g_e}

其中：

g_o

是开启一个Gap的代价。

g_e

是每增加一个额外空位的延长代价。通常：

g_o>g_e>0

例如：

g_o=5

g_e=1

那么连续四个空位的总代价是：

g(4)=5+3\times1=8

但如果是四个彼此独立的单字符Gap，则需要开启四次Gap。总代价：

4\times5=20

于是，评分模型能够区分：一次较长的连续缺失。与多次独立缺失。

图51-6：线性罚分与仿射罚分

查看图表数据

length	线性罚分：2L	仿射罚分：5+(L−1)
1	2	5
2	4	6
3	6	7
4	8	8
5	10	9
6	12	10
7	14	11
8	16	12

这里显示的是扣分的正数大小，不是实际加到比对得分里的负值。参数仅为教学示例，不代表某类生物序列必须采用这一组数值。

历史并没有在1981年停止

1982年12月15日，日本研究者：后藤修（Osamu Gotoh）发表：An Improved Algorithm for Matching Biological Sequences《一种改进的生物序列匹配算法》。刊登于《Journal of Molecular Biology》第162卷第3期，第705—708页。他提出了针对相关Gap代价形式的高效算法，使原本可能需要更复杂跳跃搜索的比对问题，可以在：

\boxed{O(nm)}

时间内求解。^{[ScienceDirect]}这里需要保留一项非常重要的历史区别。我们在第五幕使用的：

\max(\text{左上},\text{上},\text{左})

三方向递推，本身已经能够以：

O(nm)

时间处理线性Gap罚分。但如果Gap代价取决于一整段连续空位的长度，就不能简单地每次仅扣同一个固定数值。对于一般长度相关Gap代价，若在每个状态里朴素枚举各种长度，计算量就会增加。Gotoh的关键贡献，是利用更精细的状态表示，在仿射Gap模型下仍然保持高效递推。

Gotoh怎样做到？

直观地说，我们不再只用一个矩阵记录全部情况。而是区分：当前以两个真实字符对齐结束。当前以A中的字符对Gap结束。当前以Gap对B中的字符结束。常见实现使用三个动态规划状态矩阵。这样，程序就能知道：这一格是在开启一个新的Gap？还是在延长已经打开的Gap？两种操作采用不同的代价。这再次体现出动态规划的一条重要原则：如果未来的代价依赖过去的某种状态，就需要将这项历史信息包含进动态规划状态。这与第49章维特比算法为什么要记住编码器的历史状态，具有非常相似的数学思想。

第十一幕：1992年，为什么蛋白质序列比对还需要更聪明的评分表？

现在，我们已经学会了：怎样调整Gap罚分。但还有一个问题。对于蛋白质序列，我们是否应该将全部不同氨基酸的替换，都看成相同程度的不匹配？例如：A变成V。与：W变成G。它们是否必须一律扣1分？不一定。因为不同氨基酸拥有不同的：大小。电荷。极性。疏水性。以及其他理化性质。在蛋白质进化过程中，不同种类的氨基酸替换，也可能具有不同的统计倾向。因此，研究者开始建立更加精细的：替换评分矩阵（Substitution Matrix）。

从PAM矩阵到BLOSUM矩阵

戴霍夫及其合作者发展了PAM替换模型。PAM与根据进化替换规律估计的氨基酸变化概率有关。后来的研究者还发展了其他统计评分方法。1992年11月15日：史蒂文·亨尼科夫（Steven Henikoff）。以及：乔嘉·亨尼科夫（Jorja G. Henikoff）。发表：Amino Acid Substitution Matrices from Protein Blocks《从蛋白质保守片段构造氨基酸替换矩阵》。这项工作发展了：BLOSUM矩阵。其中，BLOSUM62成为蛋白质序列比对中非常常用的一种评分矩阵。^[PubMed]

它怎样改变动态规划？

注意，动态规划公式其实不需要重新发明。原来我们使用：

s(a_i,b_j)= \begin{cases} 2,&a_i=b_j\\ -1,&a_i\ne b_j \end{cases}

现在，我们可以让：

s(a_i,b_j)

直接从BLOSUM62等替换矩阵中查询。于是，两个不同的氨基酸组合，可以拥有不同得分。但全局或局部比对的递推结构仍然保持不变。这说明：动态规划负责寻找给定模型下的最优比对；评分矩阵负责表达什么样的字符对应更值得奖励或惩罚。这是两项不同职责。

为什么不直接说“相似分数越高，就越能证明共同祖先”？

因为不同序列的长度、氨基酸组成和随机背景不同。一条非常长的序列，即使与另一条无关，也可能偶然出现得分较高的局部片段。因此，研究者还需要判断：当前分数与随机背景相比，到底有多不寻常？这涉及：统计显著性。背景分布。序列组成。以及具体数据库搜索条件。在下一章BLAST中，我们会进一步研究：为什么只有一个最高比对得分，还不够支持高质量的生物序列搜索？

第十二幕：算法复杂度——为什么两条短序列可以精确比较，整个基因组却不能总用同样的方法硬算？

我们再回到计算机科学。假设两条序列长度分别为：

n

和：

m

现代教材中的线性Gap全局比对，以及相应线性Gap局部比对，都建立一个：

(n+1)\times(m+1)

的二维动态规划矩阵。因此，需要计算的内部状态数量是：

\boxed{nm}

每个状态只比较固定数量的候选。所以，时间复杂度为：

\boxed{O(nm)}

如果还需要保存整张分数表和回溯指针：

\boxed{O(nm)}

空间也很常见。

一百个字符对一百个字符

n=100,\quad m=100

内部状态数量：

100\times100=10\,000

这对于现代计算机非常容易处理。

一万对一万

10\,000\times10\,000=100\,000\,000

需要一亿个内部状态。虽然某些优化实现仍然可以处理这样的规模，但内存和计算量已经明显增加。

一百万对一百万

10^6\times10^6=10^{12}

即：一万亿个状态。如果每个状态仅保存一个4字节整数：

4\times10^{12}\text{ bytes}

就已经约为4 TB。而实际完整比对程序还可能需要保存方向信息、额外矩阵和其他状态。因此，直接保存完整矩阵并不现实。

只计算最高分，能不能减少内存？

可以。因为计算：

F(i,j)

时，只需要：上一行相关位置。以及当前行左侧位置。因此，如果我们只关心最终得分，而不要求输出完整回溯路径，可以通过滚动数组，只保存必要的相邻行。空间从：

O(nm)

减少为：

\boxed{O(\min(n,m))}

当然，时间复杂度仍然通常是：

O(nm)

这是一种空间优化。并没有减少需要计算的状态总数。

如果仍然需要恢复完整比对呢？

1975年，计算机科学家Dan Hirschberg提出了一项重要的线性空间算法思想。通过分治与动态规划结合，可以在某些序列问题中恢复最优结果，同时将额外空间需求保持在线性量级。1988年，Eugene Myers和Webb Miller进一步发展了适合生物序列比对、包含仿射Gap罚分的线性空间算法。他们的论文：Optimal Alignments in Linear Space于1988年3月发表。^[PubMed]但必须区分：减少内存，并不必然意味着计算时间也会大幅下降。分治方法可能需要重新计算部分动态规划状态，以换取较低的存储消耗。

这又引出下一章的问题

即使我们已经将精确比对算法优化得非常好，如果需要：拿一条序列。与数据库中的一千万条序列逐个比较。那么，大规模精确动态规划仍然可能十分昂贵。因此，生物信息学研究者开始寻找新的算法设计。不一定对数据库中每个候选都执行完整精确比对。而是先寻找可能具有显著相似性的短片段。再重点扩展和验证这些候选。这就是后来BLAST研究的核心方向之一。但在进入下一章之前，我们还要把本章的两种经典算法写成真正可运行的程序。

第十三幕：Python实战——从零实现两种经典序列比对算法

现在，我们将本章的全部核心算法整合为一份Python程序。它可以在两种模式下运行：Global： Needleman–Wunsch型全局比对。Local： Smith–Waterman型局部比对。程序只使用Python标准库。不依赖Biopython。也不需要任何网络服务。它会返回：最佳比对分数。对齐后的两条序列。比对中相同的位置。以及局部比对在原始序列中的位置范围。

第一部分：通用的动态规划比对函数


def align(
    a: str,
    b: str,
    mode: str = "global",
    match: int = 2,
    mismatch: int = -1,
    gap: int = -2,
):
    """
    双序列动态规划比对。

    mode:
        "global" -> Needleman-Wunsch型全局比对
        "local"  -> Smith-Waterman型局部比对

    gap为单个空位字符的得分，通常为负数。

    使用线性Gap罚分，而非仿射Gap罚分。
    """

    if mode not in ("global", "local"):
        raise ValueError("mode必须为global或local")

    n = len(a)
    m = len(b)

    # 分数矩阵
    dp = [
        [0] * (m + 1)
        for _ in range(n + 1)
    ]

    # 保存最佳前驱方向
    parent = [
        [None] * (m + 1)
        for _ in range(n + 1)
    ]

    if mode == "global":
        for i in range(1, n + 1):
            dp[i][0] = dp[i - 1][0] + gap
            parent[i][0] = "up"

        for j in range(1, m + 1):
            dp[0][j] = dp[0][j - 1] + gap
            parent[0][j] = "left"

    # 局部比对的最佳得分与终点
    best_score = 0
    best_position = (0, 0)

    # 逐格计算动态规划矩阵
    for i in range(1, n + 1):
        for j in range(1, m + 1):

            pair_score = (
                match
                if a[i - 1] == b[j - 1]
                else mismatch
            )

            diagonal = (
                dp[i - 1][j - 1]
                + pair_score
            )

            up = dp[i - 1][j] + gap
            left = dp[i][j - 1] + gap

            # 按对角、上、左的固定顺序处理平局
            if mode == "local":
                value = 0
                direction = None
            else:
                value = float("-inf")
                direction = None

            for candidate, candidate_direction in (
                (diagonal, "diag"),
                (up, "up"),
                (left, "left"),
            ):
                if candidate > value:
                    value = candidate
                    direction = candidate_direction

            dp[i][j] = value
            parent[i][j] = direction

            if mode == "local" and value > best_score:
                best_score = value
                best_position = (i, j)

    if mode == "global":
        best_score = dp[n][m]
        best_position = (n, m)

    # 从最佳终点开始回溯
    i, j = best_position

    end_a, end_b = i, j

    aligned_a = []
    aligned_b = []

    while i > 0 or j > 0:

        # 局部模式：到达零分就停止
        if mode == "local" and dp[i][j] == 0:
            break

        direction = parent[i][j]

        if direction == "diag":
            aligned_a.append(a[i - 1])
            aligned_b.append(b[j - 1])

            i -= 1
            j -= 1

        elif direction == "up":
            aligned_a.append(a[i - 1])
            aligned_b.append("-")

            i -= 1

        elif direction == "left":
            aligned_a.append("-")
            aligned_b.append(b[j - 1])

            j -= 1

        else:
            break

    aligned_a.reverse()
    aligned_b.reverse()

    result_a = "".join(aligned_a)
    result_b = "".join(aligned_b)

    # 用于显示字符相同和不同的位置
    middle = "".join(
        "|" if x == y and x != "-" else " "
        for x, y in zip(result_a, result_b)
    )

    return {
        "mode": mode,
        "score": best_score,
        "aligned_a": result_a,
        "aligned_b": result_b,
        "middle": middle,
        "range_a": (i, end_a),
        "range_b": (j, end_b),
        "matrix": dp,
    }


def print_alignment(result):
    print("模式:", result["mode"])
    print("得分:", result["score"])
    print()
    print(result["aligned_a"])
    print(result["middle"])
    print(result["aligned_b"])
    print()
    print("A区间:", result["range_a"])
    print("B区间:", result["range_b"])

这里有一项值得关注的设计。对于全局比对，我们从：

(n,m)

开始回溯。而对于局部比对，我们需要在计算矩阵的过程中，持续记录：

\max H(i,j)

对应的位置。因此，虽然两者可以共用绝大部分程序逻辑，但它们的起点、终点与边界条件并不相同。

第二部分：全局比对测试


result = align(
    "ACGT",
    "AGT",
    mode="global",
    match=2,
    mismatch=-1,
    gap=-2,
)

print_alignment(result)

assert result["score"] == 4

assert result["aligned_a"] == "ACGT"
assert result["aligned_b"] == "A-GT"

print("全局比对测试通过")

预期的主要输出：

模式: global
得分: 4

ACGT
| ||
A-GT

A区间: (0, 4)
B区间: (0, 3)

这里区间采用Python切片习惯：起点包含。终点不包含。所以：

(0, 4)

对应原始序列的第1至第4个字符。

第三部分：局部比对测试


result = align(
    "TTACGTAA",
    "GGACGTT",
    mode="local",
    match=2,
    mismatch=-1,
    gap=-2,
)

print_alignment(result)

assert result["score"] == 8

assert result["aligned_a"] == "ACGT"
assert result["aligned_b"] == "ACGT"

print("局部比对测试通过")

预期的主要输出：

模式: local
得分: 8

ACGT
||||
ACGT

A区间: (2, 6)
B区间: (2, 6)

这证明：即使两条序列开头不同、结尾不同，Smith–Waterman仍然能够找出内部高度相似的局部区域。

第四部分：一条序列与自己比较

如果：

ACGT

与：

ACGT

比较。每个位置都完全相同。没有Gap。没有不匹配。那么：

4\times2=8

所以：

\boxed{\text{score}=8}

我们可以同时验证全局与局部比对：


for mode in ("global", "local"):
    result = align(
        "ACGT",
        "ACGT",
        mode=mode,
    )

    assert result["score"] == 8
    assert result["aligned_a"] == "ACGT"
    assert result["aligned_b"] == "ACGT"

print("完全相同序列测试通过")

第五部分：不相似序列的局部比对

假设：

AAAA

与：

CCCC

比较。每个对应字符都不同。而我们的评分规则规定：不匹配扣1分。Gap也扣分。因此，如果要求找出最高的非负局部比对分数，最优选择可能是不选取任何正得分局部片段。所以：

\boxed{H_{\max}=0}

这不是程序没有找到路径。而是局部比对模型允许：没有任何值得保留的正得分片段。


result = align(
    "AAAA",
    "CCCC",
    mode="local",
)

assert result["score"] == 0
assert result["aligned_a"] == ""
assert result["aligned_b"] == ""

print("无正分局部片段测试通过")

第六部分：检查算法的时间复杂度

如果序列A长度为：

n

序列B长度为：

m

那么，程序内部的两层循环是：

for i in range(1, n + 1):
    for j in range(1, m + 1):

总共执行：

nm

次状态更新。每次更新只执行固定数量的比较。因此：

\boxed{T(n,m)=O(nm)}

而：

dp

与：

parent

都保存：

(n+1)(m+1)

个状态。所以空间复杂度同样为：

\boxed{S(n,m)=O(nm)}

这份实现的优点不是已经达到工业级速度。而是：代码与数学公式几乎一一对应。我们能够非常清楚地看到：初始化。递推。记录前驱。寻找终点。以及回溯。这些步骤，构成了最经典的动态规划求解过程。

第十四幕：一项容易被误解的科学结论——算法找到的“相似”，究竟意味着什么？

现在，我们已经可以让Python程序输出：

ACGT
||||
ACGT

也能够得到一个明确的最高分数。但是，这还不是生物学研究的终点。甚至不能自动证明：两条序列具有共同的进化祖先。

一个随机序列，也可能碰巧出现很好的局部匹配

想象我们生成两条很长的DNA序列。每个字符都随机选择：A。C。G。T。即使两条序列完全没有生物学上的共同来源，也有可能碰巧包含：

ACGT

这样的相同片段。如果序列非常长，相同短片段出现的机会还会增加。因此：局部匹配存在，不等于一定存在共同的进化关系。这也是为什么，Needleman与Wunsch在1970年的研究中，已经使用随机打乱序列的方法进行对照。他们不是只输出一个最高匹配得分，就立即宣布两个蛋白质之间存在确定的生物学关系。

序列一致性与序列相似性也不完全相同

例如：

A C G T
A C G T

四个位置完全一致。序列一致性：

100\%

但蛋白质分析中，还有一种更广义的“相似”。假设两个对应的氨基酸不同。但它们具有相似的理化性质。或者在某个经过统计训练的替换矩阵中，获得较高的替换得分。这时，我们可能认为它们属于值得考虑的保守替换。因此：**Identity（序列一致性）**关注对应位置字符是否完全相同。**Similarity（序列相似性）**则可能依赖更复杂的评分定义。不同软件给出的相似性统计，可能采用不同的矩阵和阈值。所以，比较结果时还需要核查：评分系统。Gap罚分。所比较区域的长度。以及统计显著性。

更不能把比对直接等同于蛋白质功能预测

如果两条蛋白质序列高度相似，确实可以成为研究共同结构或功能的有用线索。但功能推断还需要考虑：序列覆盖范围。关键活性位点。蛋白质结构。物种与进化背景。以及实验验证。尤其当两条蛋白质只在一个很短的局部片段中相似时，不能直接宣称它们具有完全相同的功能。NCBI Bookshelf的进化基因组学教材明确指出，同源性指共同祖先关系；仅凭一项可能偶然产生的序列比对，不能自动推断生物学功能。^{[NCBI Bookshelf]}

今天的算法，还需要面对更大的研究对象

我们已经知道：两条100字符序列。两条1000字符序列。都可以利用动态规划比较。但现代生物信息学的序列规模，可以远远超过1970年研究者最初面对的情况。例如：一个完整的基因组。大量测序读段。千万级的蛋白质数据库。以及不断增长的公共序列资料集。因此，研究者还需要：高效的索引。启发式搜索。数据库过滤。以及并行计算。并不是所有任务都适合对全部候选直接执行完整的：

O(nm)

精确动态规划。这便为下一章BLAST的出现，准备了真正的历史与算法背景。

尾声：当生命成为一串字母，数学开始寻找它们之间的历史

1950年代，剑桥。桑格和他的同事们，在实验室里逐步解开胰岛素的氨基酸顺序。他们首先证明了一件重要的事。蛋白质并不是由一组任意混合的氨基酸构成。氨基酸的连接顺序，是分子结构的重要组成部分。这种顺序可以被测定。也可以被记录。随着实验技术发展，越来越多蛋白质的序列进入科学文献。1965年，戴霍夫与合作者出版《蛋白质序列与结构图谱》。她试图让散落在不同研究报告中的序列，成为能够集中保存和比较的科学资料。于是，一个新的问题出现了。如果我们已经能够读出生命分子的序列，下一步应该怎样理解不同序列之间的关系？1969年7月21日，芝加哥。Needleman与Wunsch向期刊提交研究论文。他们没有选择用肉眼逐一比较所有可能的氨基酸对应。而是将序列比较放进二维矩阵。让每个位置代表一对氨基酸。让每条合法路径代表一种保持原始顺序的对应方案。通过累积最佳匹配分数，排除大量不必要的候选。1970年3月，这项研究正式发表。科学家第一次拥有了一种能够系统搜索两条蛋白质序列最佳对应关系的重要计算方法。它后来成为经典的Needleman–Wunsch算法。时间又过去十一年。1981年，Smith与Waterman提出：如果两条长序列整体并不相似，但内部存在重要的共同片段，就不应该强迫算法覆盖所有字符。于是，他们让动态规划得分可以重新从零开始。再从整张矩阵中的最高分位置，回溯到零。这便形成了经典的局部序列比对方法。后来，Gotoh继续改进Gap罚分的计算。戴霍夫的替换模型和亨尼科夫夫妇的BLOSUM矩阵，为蛋白质序列比较提供了更细致的评分依据。Hirschberg、Myers、Miller等研究者，则进一步探索如何在有限内存中恢复最优比对。这是一条跨越数十年的研究道路。并不是一位研究者独自完成了全部工作。

本章真正学会的两项算法

第一项：Needleman–Wunsch型全局比对。它适合研究需要完整覆盖的两条序列之间，怎样获得最优对应关系。经典线性Gap版本的递推为：

\boxed{ F(i,j)=\max \begin{cases} F(i-1,j-1)+s(a_i,b_j)\\ F(i-1,j)-d\\ F(i,j-1)-d \end{cases} }

最终答案来自：

\boxed{F(n,m)}

第二项：Smith–Waterman型局部比对。它允许算法在任意位置重新开始，从两条序列中寻找最优局部对应片段。递推为：

\boxed{ H(i,j)=\max \begin{cases} 0\\ H(i-1,j-1)+s(a_i,b_j)\\ H(i-1,j)-d\\ H(i,j-1)-d \end{cases} }

最终从：

\boxed{\max_{i,j}H(i,j)}

开始回溯，直到遇到零分状态。两个算法具有几乎相同的矩阵结构。却解决了不同的优化问题。而这种区别，正是算法设计最重要的思想之一。有时候，决定答案的不是计算步骤多么复杂，而是我们究竟怎样定义一个“好答案”。

为什么这一章要放在维特比算法之后？

在第49章，我们研究：怎样从大量可能的隐藏历史中，找出累计分数最高的一条状态路径。现在，我们研究：怎样从大量可能的字符对应关系中，找出评分最高的一条比对路径。它们都使用了：最优子结构。动态规划状态。局部转移。以及前驱回溯。但是，两项算法的状态定义不同。维特比算法的状态通常是：某个时间点上的编码器状态或隐藏状态。序列比对算法的状态则是：两条序列分别已经处理了多少个字符。这再次说明：动态规划并不是只能处理背包、最短路径或某一种固定题型的技巧。它是一种组织具有重叠子问题和最优子结构的计算问题的方法。

算法找到了相似性，生命的故事却没有因此结束

假设两条蛋白质序列，经过比对后显示出明显的共同片段。研究者可能据此继续探索：它们是否具有共同的进化来源？某些保守位置是否与蛋白质功能有关？某次插入或缺失是否影响了分子结构？它们是否属于某个蛋白质家族？这些问题，都可以从序列比对结果开始。但不能只凭一项评分就直接得出全部结论。因为生命的历史，比一张得分矩阵更加复杂。同样的短序列，可能偶然出现。不同位置上的替换，可能具有不同的统计意义。而真正的进化关系，还需要结合更广泛的证据。因此，算法所做的，不是代替科学家宣布最终的生物学事实。而是将原本难以系统搜索的海量对应关系，转化成一套可以计算、比较和检验的候选结构。它让科学家能够在更加明确的数学条件下，提出下一步问题。生命将历史写在了分子的顺序里。而序列比对算法，使我们第一次能够更加系统地比较这些顺序。一条序列，告诉我们一种分子由什么组成。两条序列之间的对应关系，则可能成为理解生命如何变化的一条线索。

历史档案：从蛋白质测序到现代序列比对

年份	科学家与事件	历史意义
1955年前后	Frederick Sanger等人完成胰岛素序列与连接关系研究	建立蛋白质完整氨基酸顺序的实验基础
1958年	Sanger获得诺贝尔化学奖	表彰其蛋白质结构研究
1965年	Margaret Dayhoff等人出版蛋白质图谱	推动序列资料集中整理与计算分析
1969年7月21日	Needleman与Wunsch论文收到	提出矩阵最大匹配方法
1970年3月28日	Needleman–Wunsch论文发表	序列比对动态规划历史的重要里程碑
1974年	Peter H. Sellers发展序列距离方法	将序列差异形式化为可优化的编辑代价
1976年	Waterman等人推广序列距离模型	支持更一般的插入、删除长度代价
1981年3月25日	Smith–Waterman论文发表	明确提出最佳局部片段的动态规划识别
1982年12月15日	Osamu Gotoh发表改进算法	高效处理仿射Gap罚分
1988年3月	Myers与Miller发表线性空间比对研究	推动节省内存的最优比对实现
1992年11月15日	Henikoff夫妇发表BLOSUM研究	发展蛋白质氨基酸替换评分体系

主要参考文献与一手资料

桑格与胰岛素序列研究MRC Laboratory of Molecular Biology — 1958 Fred Sanger英国MRC分子生物学实验室的历史资料，核查桑格测定胰岛素51个氨基酸残基及其1958年诺贝尔奖的研究背景。

1965年 · 戴霍夫的蛋白质图谱Collecting, Comparing, and Computing Sequences — Historical Study研究Margaret Dayhoff、Richard Eck和Robert Ledley如何建立早期蛋白质序列资料集，以及计算机在序列整理中的作用。

1970年 · Needleman–Wunsch原始论文A General Method Applicable to the Search for Similarities in the Amino Acid Sequence of Two Proteins — PDF Journal of Molecular Biology — 期刊书目核心一手史料。建议重点阅读第443—445页的矩阵路径思想，以及第447—449页关于鲸肌红蛋白、人血红蛋白β链和随机序列对照的实验。

1981年 · Smith–Waterman原始论文Identification of Common Molecular Subsequences — PDF Journal of Molecular Biology — 期刊书目建议重点阅读第195页含零值的递推式，以及第196页关于从最高分回溯到零的解释。这是本章局部比对算法最直接的一手依据。

1982年 · Gotoh仿射Gap算法Osamu Gotoh — An Improved Algorithm for Matching Biological Sequences用于核查Gap长度相关代价的高效计算及1982年的重要算法改进。

1988年 · 线性空间最优比对Myers & Miller — Optimal Alignments in Linear Space解释Hirschberg分治思想如何进一步用于生物序列比对，以及内存优化为何对长序列非常重要。

1992年 · BLOSUM替换矩阵Henikoff & Henikoff — Amino Acid Substitution Matrices from Protein Blocks蛋白质序列比对评分模型的重要研究，说明为什么不同氨基酸替换需要不同得分。

生物学解释：序列一致性与共同祖先关系NCBI Bookshelf — Evolutionary Concept in Genetics and Genomics用于区分序列相似性、一致性与同源性，避免将算法分数直接解释为共同祖先关系的百分比。

**历史与教学重构说明：**桑格的胰岛素研究、戴霍夫的序列图谱、Needleman–Wunsch的投稿与发表日期、1970年实际比较的蛋白质，以及Smith–Waterman在1981年的递推方法，均有相应史料依据。本章的实验室氛围和科学问题演进采用历史叙事写法，不代表存在未经证实的特定人物对话。本章使用的ACGT、AGT与TTACGTAA等序列，是为数学教学构造的例子。经典三方向、线性Gap罚分的Needleman–Wunsch型实现，是现代教学中广泛使用的简化形式，不应误写成1970年原始论文全部计算过程的精确复刻。互动实验与Python代码没有包含仿射Gap、复杂氨基酸替换矩阵或统计显著性估计，因此不应直接用作正式的生物学证据判定系统。

下一章 · 第52章

《在浩瀚基因数据库里寻找一个线索：BLAST的诞生》

1990年，美国马里兰州，贝塞斯达。分子生物学正在进入一个新的时代。科学家能够测定的DNA和蛋白质序列越来越多。这些序列被保存在不断扩大的数据库里。一名研究者刚刚获得一段新的蛋白质序列。他想知道：数据库里是否已经存在某种相似的蛋白质？它可能有什么功能？是否拥有值得研究的共同进化来源？从理论上说，我们已经拥有Smith–Waterman算法。它能够找到两条序列之间的最佳局部对应。但如果数据库拥有几十万甚至更多条序列，逐一进行完整动态规划，将需要大量计算。于是，一群研究者开始寻找另一条道路。不必先精确比较全部候选。能否首先寻找短而得分较高的共同片段？然后只扩展最有希望的候选区域？1990年10月5日，五位研究者在《Journal of Molecular Biology》发表论文：Basic Local Alignment Search Tool作者包括：Stephen F. Altschul、Warren Gish、Webb Miller、Eugene W. Myers与David J. Lipman。这项工具后来广为人知，简称：BLAST。它的原始论文第215卷第3期，第403—410页，既讨论快速近似局部比对，也讨论如何分析匹配结果的统计显著性。 ^[PubMed]下一章，我们将回到1990年的NCBI研究环境，详细追踪BLAST怎样通过短词匹配、候选扩展与统计模型，大幅提高大规模序列数据库的搜索效率。还将亲手构造一个微型DNA数据库，使用Python实现简化的种子搜索与延伸算法，并与本章的Smith–Waterman精确局部比对进行比较。我们将看到：精确算法告诉我们最佳答案在哪里；优秀的启发式搜索，则让我们有机会在海量数据中，以可承受的计算代价尽快找到最值得研究的候选。

——《算法的故事》第二季 · 第五篇《破解生命的密码》 · 第51章完