打分矩阵与 Gap 罚分

快速概览

真实序列比对不能假设「所有错误代价相同」。打分矩阵通过量化氨基酸替换的可能性，Gap 罚分通过区分开启与延长代价，共同构建了模拟进化过程的评价体系。

• 理解 PAM (点突变) 与 BLOSUM (保守块) 矩阵的构建逻辑差异
• 掌握仿射 Gap 罚分（Affine Gap Penalty） 对连续 Indel 事件的偏好
• 学习如何根据任务目标（如远缘同源检测 vs. 精确定位）选择得分系统
• 建立从概率模型（Log-odds Ratio）到打分矩阵的数学联系

前置知识

• 编辑距离
• 生物学基础对象

所属板块 核心方法

索引、比对、组装与概率模型构成的核心算法主轴。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系，再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉，再进入本页。

1. 是什么

打分矩阵（Substitution Matrix / Scoring Matrix）和 Gap 罚分（Gap Penalty）共同构成了序列比对的评价体系。打分矩阵为每一对字符（碱基或氨基酸）的对齐指定一个得分，Gap 罚分则定义了插入或缺失（Insertion/Deletion, Indel）的代价。

在序列比对的动态规划算法中，每一步的最优决策都依赖于当前这一”列”对齐的得分。打分矩阵决定了”匹配好还是错配好，以及错配之间也有区别”，而 Gap 罚分决定了”开一个 Gap 还是延长已有 Gap 更划算”。一个好的打分体系应该反映真实的生物学进化过程。

2. 要解决什么生物信息学问题

编辑距离（Edit Distance）假设所有操作代价相同，但在真实的分子进化中：

• 替换不等价：亮氨酸（Leu, L）和异亮氨酸（Ile, I）都是疏水性氨基酸，结构相似，它们之间的替换在进化中更容易被接受。而亮氨酸和天冬氨酸（Asp, D）一个是疏水一个是带负电，这种替换通常受到强烈的负选择。
• Gap 不等价：一次长度为 5 的连续缺失（如一次滑动或重组事件）比 5 次独立的单碱基缺失更可能发生。进化中，大片段的插入或缺失往往是单次事件的结果。
• 蛋白质 vs 核酸：蛋白质有 20 种氨基酸，其生化性质的差异远大于 4 种碱基之间的差异，因此蛋白质比对需要更精细的打分矩阵。

打分矩阵和 Gap 罚分的目标，就是将这些生物学知识编码为数学量，使比对算法的优化目标与生物学真实一致。

3. 蛋白质替换矩阵

PAM (Percent Accepted Mutation)

由 Margaret Dayhoff 于 1978 年提出，基于进化模型外推。

构建过程

PAM 矩阵的构建分为以下几个步骤：

1. 收集数据：从高度相似（> 85% 一致性）的蛋白质序列对中，统计观察到的氨基酸替换事件。
2. 计算突变概率矩阵 $M$：$M_{a,b}$ 表示氨基酸 $a$ 在一个进化时间步内突变为 $b$ 的概率。
3. 外推：通过矩阵乘法 $PAM_n = M^n$ 模拟 $n$ 步进化后的替换概率。

PAM 数字的含义

• PAM 1：代表 1% 的氨基酸发生了”可接受”突变。这对应于约 1 亿年的进化时间。
• PAM 250：代表每个位置平均发生了 250 次可接受突变（注意：不是 250% 的氨基酸发生了变化，因为同一个位置可能经历多次突变，包括回复突变）。

如何选择 PAM 矩阵

PAM 编号	进化距离	适用场景
PAM 30	近缘	高度相似的序列（> 70% 一致性）
PAM 70	中等	中等相似度的同源序列
PAM 250	远缘	寻找远距离同源关系

直觉：PAM 数字越大，代表演化时间越长，矩阵中的得分差异越不明显（因为远缘序列中几乎任何替换都可能发生过），越适合找远缘同源序列。

BLOSUM (Blocks Substitution Matrix)

由 Henikoff 和 Henikoff 于 1992 年提出，基于实际观测到的保守区块统计。

构建过程

1. 收集保守区块：从 PROSITE 蛋白质家族数据库中，提取高度保守的蛋白质区域（称为 Blocks）。
2. 聚类：在每个 Block 中，将相似度超过阈值 $T\%$ 的序列聚类，用一条代表序列替代。
3. 统计：在聚类后的序列集合中，统计每对氨基酸在同一列中共现的频率。
4. 计算 Log-odds 得分：将频率转换为 Log-odds 分数。

BLOSUM 数字的含义

• BLOSUM 62：构建时将相似度 > 62% 的序列进行了聚类。这是最常用的替换矩阵。
• BLOSUM 80：构建时聚类阈值为 80%，保留了更多近缘序列的信息。

如何选择 BLOSUM 矩阵

BLOSUM 编号	聚类阈值	适用场景
BLOSUM 80	80%	近缘、高度保守的序列
BLOSUM 62	62%	通用场景（默认推荐）
BLOSUM 45	45%	远缘同源检测

直觉：BLOSUM 数字越大，代表用于构建矩阵的序列越相似，矩阵越”保守”，越适合找高度保守的区域。数字越小，矩阵中允许的替换越多，越适合检测远缘同源关系。

PAM vs BLOSUM 的核心区别

维度	PAM	BLOSUM
**构建方式**	进化模型外推	直接统计观测
**数据来源**	全长蛋白质序列对	保守区块（Blocks）
**数值含义**	进化时间步数（越大越远）	聚类阈值（越大越近）
**优点**	理论基础清晰，可外推	直接反映真实数据，更实用

在实际使用中，BLOSUM 62 是最广泛使用的默认选择，因为它在大多数场景下表现稳健。

4. Gap 罚分模型

线性 Gap 罚分（Linear Gap Penalty）

最简单的 Gap 模型：每个 Gap 字符的代价相同。

$\text{Cost}(k) = k \cdot g$

其中 $g$ 是常数罚分值，$k$ 是 Gap 长度。这种模型简单但不符合生物学直觉——它假设每个 Gap 字符是独立事件。

仿射 Gap 罚分（Affine Gap Penalty）

为了反映”连续 Indel 更易发生”的直觉，我们使用：

$\text{Cost}(k) = g_o + k \cdot g_e$

• $g_o$ (Gap Open)：开启一个新缺口的高额代价。
• $g_e$ (Gap Extend)：每增加一个碱基的低廉代价。

算法影响：在动态规划中，这要求我们维护三个平行的矩阵（$M, I, D$），分别记录匹配、插入和删除的状态，从而追踪 Gap 是否正在延长。详见 仿射 Gap 罚分。

常见参数设置

应用场景	Gap Open ($g_o$)	Gap Extend ($g_e$)	说明
蛋白质近缘比对	-10 ~ -11	-1 ~ -2	保守默认
蛋白质远缘比对	-12 ~ -15	-1 ~ -2	较大的开启罚分
核酸比对	-5 ~ -6	-2 ~ -3	较小的开启罚分

5. 从概率到得分：Log-odds Ratio

矩阵中的得分 $s(a, b)$ 并不是拍脑袋定的，它通常是：

$s(a, b) = \lambda \cdot \log \frac{P_{observed}(a, b)}{P_{expected}(a, b)}$

其中：

• $P_{observed}(a, b)$ 是在同源序列中观察到 $a$ 和 $b$ 配对的概率。
• $P_{expected}(a, b) = P(a) \cdot P(b)$ 是在随机序列中 $a$ 和 $b$ 配对的期望概率。
• $\lambda$ 是缩放因子，通常取 $\frac{1}{\ln 2}$ 或 $\frac{10}{\ln 10}$ 以将得分转换为整数。
• $\log$ 通常取自然对数或以 2 为底的对数。

得分的生物学解读

• 分数为正：观测到的替换频率高于随机概率，$s(a,b) > 0$，暗示演化上的正选择或中性漂变。例如 BLOSUM62 中 $s(W, W) = 11$（色氨酸-色氨酸），因为色氨酸是最大且最罕见的氨基酸，它的保守性极高。
• 分数为负：观测频率低于随机，$s(a,b) < 0$，暗示生化性质不相容或演化上的负选择。例如 $s(W, P) = -4$（色氨酸-脯氨酸），因为色氨酸是疏水的，而脯氨酸会破坏二级结构。
• 对角线最高：每个氨基酸与自己配对的得分通常是该行/列中的最大值。

6. Worked Example：BLOSUM62 矩阵片段解读

以下是 BLOSUM62 矩阵的一个子集（展示部分氨基酸）：

	A	R	N	D	C	Q	E	G	H	I
A	4	-1	-2	-2	0	-1	-1	0	-2	-1
R	-1	5	0	-2	-3	1	0	-2	0	-3
N	-2	0	6	1	-3	0	0	0	1	-3
D	-2	-2	1	6	-3	0	2	-1	-1	-3

解读关键条目

• $s(A, A) = 4$：丙氨酸（Ala, A）是小的疏水氨基酸，自身保守性较高。
• $s(D, E) = 2$：天冬氨酸（Asp, D）和谷氨酸（Glu, E）都是带负电的氨基酸，生化性质相近，替换可接受。
• $s(R, D) = -2$：精氨酸（Arg, R）带正电，天冬氨酸带负电，它们之间的替换会破坏电荷平衡，不利于蛋白质功能。
• $s(N, N) = 6$：天冬酰胺（Asn, N）的自身得分很高，因为它在许多功能位点（如糖基化位点）中是保守的。

在比对中的应用

假设我们比对一个蛋白质序列片段，某列中一个物种是 D（Asp），另一个物种是 E（Glu）。使用 BLOSUM62，这一对的得分为 +2，说明这是一个可以接受的保守替换。但如果一个是 R（Arg），另一个是 D（Asp），得分为 -2，说明这是一个不太可能的功能性替换。

7. 核酸打分矩阵

对于核酸序列（DNA/RNA），由于只有 4 种碱基，打分矩阵相对简单。常见的方案包括：

• 简单匹配/错配：匹配 +1，错配 -1。
• 转换/颠换模型：区分转换（Transition, A↔G, C↔T）和颠换（Transversion, 嘌呤↔嘧啶）。转换的频率约为颠换的 2 倍，因此可以设置转换错配 -1，颠换错配 -2。
• 矩阵名称：如 BLAST 的核酸矩阵 NUC.4.4（匹配 +1，错配 -3）。

转换-颠换的差异源于分子机制：嘌呤之间或嘧啶之间的替换只需要一次碱基的化学变化（脱氨），而嘌呤与嘧啶之间的替换需要两次变化。

8. 复杂度分析

打分矩阵本身是 $|\Sigma| \times |\Sigma|$ 的查找表（Lookup Table），其中 $\Sigma$ 是字母表大小：

• 蛋白质：$|\Sigma| = 20$，矩阵大小为 $20 \times 20 = 400$ 个条目。
• 核酸：$|\Sigma| = 4$，矩阵大小为 $4 \times 4 = 16$ 个条目。

在比对算法中，每次查表的复杂度为 $O(1)$。打分矩阵的选择不会改变比对算法本身的渐进复杂度（仍为 $O(mn)$），但会显著影响比对结果的生物学质量。

适用前提

• 替换矩阵的构建依赖于已有的进化数据。如果研究的是人工合成序列或高度变异的病毒序列，标准矩阵可能不适用。
• Gap 罚分的参数需要根据具体任务调整。没有一组参数在所有场景下都是最优的。

9. 与真实工具的连接

• BLAST 默认使用 BLOSUM62 作为蛋白质比对的替换矩阵，可通过 -matrix 参数切换。详见 BLAST。
• ClustalW / MUSCLE / MAFFT 等 MSA 工具均支持自定义替换矩阵和 Gap 罚分参数。
• BWA 等短序列比对工具通常使用简单的匹配/错配得分（如匹配 +1，错配 -3，Gap Open -6，Gap Extend -1），因为核酸的替换模式比蛋白质简单。
• 替换矩阵数据库：NCBI 提供了完整的 PAM 和 BLOSUM 系列矩阵，可在 NCBI FTP 站点下载。

常见误区

• BLOSUM62 永远是最好的选择：
• 不是。BLOSUM62 是一个通用默认值，但在特定场景下可能不是最优的。如果搜索远缘同源，BLOSUM45 或 PAM250 可能更合适；如果比较高度保守的序列，BLOSUM80 可能更好。
• PAM 和 BLOSUM 是完全独立的体系：
• 实际上它们在某种程度上是等价的。PAM 250 和 BLOSUM 45 大致对应相同的进化距离。研究表明，在相同进化距离下，两者产生的比对结果高度一致。选择哪一种更多是历史和使用习惯的问题。
• 得分为正就意味着生物学上有意义：
• 不一定。正得分只说明该替换比随机期望更频繁，但不一定意味着具有功能意义。得分只是一个统计指标，生物学解释还需要结合蛋白质结构、功能位点等信息。

相关页面

全局比对算法

应用得分矩阵的标准 DP 实现

延伸阅读

多序列比对（MSA）

基于 SP 得分的群体比较

延伸阅读

BLAST

启发式搜索中的得分阈值应用

延伸阅读