本文详解 连锁不平衡（Linkage Disequilibrium，LD） 矩阵的计算原理与实现方法，重点对应 AiCE（AI-informed Constraints for protein Engineering，人工智能约束蛋白工程）框架 AiCEmulti 模块中的 02.caculated_ld.py 流水线。与 酶改造-modelpaper-AiCE 中 §4.4 互补：该文侧重框架总览，本篇专述 LD 的数学定义、编码策略与组合打分。
段末注释：AiCE 为 Cell 2025 发表的通用蛋白突变设计框架，详见 酶改造-modelpaper-AiCE；AiCEmulti 为其组合突变提名模块。
1. 什么是 LD 矩阵
连锁不平衡（Linkage Disequilibrium，LD） 描述两个位点上等位基因（allele）在样本间非随机共现的程度。在群体遗传学中，若位点 A 的等位基因与位点 B 的等位基因在个体间倾向于一起出现（或一起不出现），则称二者存在 LD。
AiCE 将这一概念借用于蛋白工程：把逆折叠采样得到的多序列比对（Multiple Sequence Alignment，MSA）中每条序列视为一个「个体」，把各位点的氨基酸（经伪 DNA 编码后）视为「等位基因型」，用 LD 衡量哪些位点的突变在采样序列中倾向于协同出现。高 LD 的组合突变，负向上位效应（negative epistasis） 风险相对更低。
段末注释：等位基因指同一基因座上的不同变异形式；AiCE 中「基因座」被映射为序列位点。MSA 为多条同源序列按位点对齐的表示。
2. 经典 LD 的数学定义
对两个双等位基因位点，记：
- 位点 A：等位基因 (A/a)，(A) 的频率为 (p_A)
- 位点 B：等位基因 (B/b)，(B) 的频率为 (p_B)
- 单倍型 (AB) 的联合频率为 (p_{AB})
2.1 连锁不平衡系数 (D)
[
D = p_{AB} - p_A \cdot p_B
]
- (D = 0)：两地位点统计独立（无 LD）
- (D \neq 0)：存在非随机关联
2.2 标准化 (D’)（Lewontin 标准化）
[
D’ = \frac{D}{D_{\max}}, \quad D_{\max} = \min\bigl(p_A p_b,; p_a p_B\bigr)
]
(D’) 取值范围为 ([-1, 1])，便于在不同等位基因频率下比较 LD 强度。
2.3 平方相关系数 (r^2)（AiCE 实际使用）
[
r^2 = \frac{D^2}{p_A(1-p_A), p_B(1-p_B)}
]
(r^2 \in [0, 1])：
| (r^2) | 含义 |
|---|---|
| (= 1) | 两地位点完全共变 |
| (\approx 0) | 两地位点近似独立 |
AiCE 官方实现通过 PLINK 的 --r2 square 输出 (r^2) 方阵，即 .ld 文件中的数值。
段末注释：PLINK 为群体遗传学常用基因型分析工具；Hardy–Weinberg 平衡（HWE） 等群体遗传假设在 AiCE 场景下不适用，故流水线中关闭 HWE 检验。
3. AiCE 中 LD 矩阵的计算流程
AiCE 的 LD 计算在官方仓库 scripts/02.caculated_ld.py 中实现，整体为四步流水线：
1 | flowchart LR |
命令：
1 | python scripts/02.caculated_ld.py <seq_dir> <output_ld_dir> |
对每个 *.fa 输入，输出同前缀的 .vcf 与 .ld 文件。
3.1 Step 1：蛋白 MSA → 伪 DNA
输入：逆折叠采样得到的 .fa 文件。首条序列为野生型（wild type，WT）参考，其余 (N) 条为采样序列（默认 (N = 1000)）。
方法：每条氨基酸序列按固定最优密码子表逐残基翻译为 DNA。这不是真实基因组编码，仅为 PLINK 提供等位基因型表示：
| 氨基酸 | 伪 DNA 密码子 | 氨基酸 | 伪 DNA 密码子 |
|---|---|---|---|
| A | GCG | M | ATG |
| R | CGT | F | TTT |
| N | AAC | P | CCG |
| D | GAT | S | AGC |
| C | TGC | T | ACC |
| Q | CAG | W | TGG |
| E | GAA | Y | TAT |
| G | GGC | V | GTG |
| H | CAT | L | CTG |
| I | ATT | K | AAA |
| - / X | — |
设计意图：PLINK 处理的是碱基型 SNP 数据，不能直接读蛋白 MSA。最优密码子映射是格式适配——同一氨基酸始终对应同一密码子，因此氨基酸层面的共现结构完整传递到核苷酸层面。
蛋白长度为 (L) 残基时，伪 DNA 长度为 (3L) bp。
3.2 Step 2：DNA FASTA → VCF
参考序列：首条 DNA 序列（对应 WT）。
逐碱基扫描（convert_fasta_to_vcf）：
- 对每个碱基位点 (pos)，统计所有样本（除参考外）的碱基分布
- 若存在与 REF 不同的碱基 → 写入一条 VCF 记录
- 每个样本基因型编码为二倍体纯合：
0|0（REF）或1|1（ALT）等
VCF 示例：
1 | #CHROM POS REF ALT ... FORMAT sample1 sample2 ... |
变异过滤：无变异的位点跳过；gap（-）与插入/缺失简化为 <INS> / <DEL> 标记。
重要细节：由于每个氨基酸对应 3 个碱基，一个氨基酸位点的变异通常会在伪 DNA 上产生最多 3 个高度相关的 SNP（同一密码子内共变）。最终 VCF 中的 SNP 数 (M) 一般小于 (3L)，且不一定等于氨基酸位点数 (L)。
3.3 Step 3：PLINK 计算 (r^2) 矩阵
1 | plink --vcf prefix.vcf \ |
| 参数 | 含义 |
|---|---|
--r2 square |
输出 (M \times M) 的 (r^2) 对称方阵 |
--maf 0.00001 |
极低最小等位基因频率阈值，保留稀有变异 |
--hwe 0 |
关闭 Hardy–Weinberg 检验（合成序列不适用） |
--geno 1 / --mind 1 |
允许 100% 样本/位点保留 |
输出：prefix.ld，(M \times M) 矩阵，(L_{ij} = r^2_{ij})。
3.4 Step 4：格式转换与清理
- 用
sed将.ld中 tab 分隔转为逗号分隔，供com_mut_prediction.py读取 - 删除中间文件（
.log、.nosex、_dna.fasta等）
4. 矩阵元素的生物学（工程）含义
在 AiCE 语境下，(L_{ij} = r^2_{ij}) 表示：
在 (N) 条逆折叠采样序列中，SNP (i) 与 SNP (j) 的等位基因型共变强度。
| (r^2) 范围 | 解读（AiCEmulti 视角） |
|---|---|
| 接近 1 | 两位点突变在采样序列中几乎总是一起出现 → 组合协同性好 |
| 0.5–0.8 | 中等共现；AiCE 默认阈值 ≥ 0.5 即推荐 |
| 接近 0 | 两位点独立变化 → 组合可能遭遇负向上位 |
与 SCA 的分工：
| 方法 | 刻画对象 | 数据来源 |
|---|---|---|
| LD | 采样序列中的统计共现 | 伪 DNA → VCF → PLINK |
| SCA（Statistical Coupling Analysis，统计耦合分析） | 进化耦合结构 | pySCA 对蛋白 MSA 的协方差分解 |
AiCEmulti 两者联合使用，LD 作为数据驱动的共现过滤器，SCA 作为进化约束过滤器。
段末注释：SCA 原用于从天然 MSA 识别功能 sector；AiCE 将其与 LD 一并应用于逆折叠合成 MSA。
5. 组合突变的下游打分
对候选位点组合 (\mathcal{S} = {p_1, p_2, \ldots, p_k})（1-based 氨基酸位点编号），从 LD 矩阵提取子矩阵 (\mathbf{L}_{\mathcal{S}})，由 com_mut_prediction.py 计算以下导出指标：
5.1 Mean Pairwise LD（默认筛选用）
[
\overline{\mathrm{LD}} = \frac{2}{k(k-1)} \sum_{i < j} L_{p_i p_j}
]
取子矩阵上三角非对角元素的算术均值。
5.2 Log Mean Pairwise LD
[
\overline{\log\mathrm{LD}} = \mathrm{mean}{i<j}\bigl[\log(L{p_i p_j} + 10^{-10})\bigr]
]
对数变换压缩极值，使分布更平滑。
5.3 Multilocus LD（连乘形式）
[
\mathrm{MLD} = \prod_{i < j} L_{p_i p_j}
]
强调所有位点对同时高关联的组合。
5.4 推荐规则
输出格式（.ld.result）：
1 | Mutation Type Mean Pairwise score Log Mean Pairwise score Logical Flag (0/1) |
- 默认：Mean Pairwise LD ≥ 0.5 → Logical Flag = 1（推荐该组合）
- 阈值可通过
04.com_mut_prediction.sh的-t参数自定义
完整 AiCEmulti 命令示例：
1 | # 1. 构建 LD 矩阵 |
6. 实现细节与使用注意
6.1 伪 DNA 编码的必要性
PLINK 设计用于真实基因型数据（VCF 格式），不能直接读取蛋白 MSA。最优密码子翻译是接口适配层，不改变「哪些序列在哪些位点共变」的统计结构。
6.2 矩阵维度与位点索引
- LD 矩阵行/列对应 VCF 中有变异的核苷酸 SNP（(M) 个），而非直接的氨基酸位点数 (L)
- SCA 矩阵为 (L \times L)（氨基酸位点级）
- 组合打分时，
.comb文件中的位点编号为氨基酸 1-based 索引；使用时需确保与 LD 矩阵维度对齐。若 (M \neq L)，需检查 VCF 构建逻辑或做显式位点映射
6.3 样本量对估计稳定性的影响
(r^2) 估计依赖采样序列数 (N)（ProteinMPNN 默认 num_seq_per_target=1000）：
- (N) 过小 → LD 估计方差大、假阳性/假阴性增多
- (N) 增大 → 共现模式更可靠，但计算成本线性增加
6.4 与群体遗传学用法的差异
| 维度 | 群体遗传学 | AiCE |
|---|---|---|
| 「个体」 | 真实生物样本 | 逆折叠采样序列 |
| 「等位基因」 | 天然 SNP | 伪 DNA 碱基型 |
| HWE 假设 | 通常检验 | 关闭（--hwe 0） |
| 生物学解释 | 连锁、选择、漂变 | 结构兼容序列的共变模式 |
7. 输入输出文件一览
| 文件 | 阶段 | 内容 |
|---|---|---|
prefix.fa |
输入 | 蛋白 MSA（参考 + 采样序列） |
prefix_dna.fasta |
中间 | 伪 DNA 序列（运行后删除） |
prefix.vcf |
中间/输出 | 变异位点与基因型 |
prefix.ld |
输出 | (M \times M) 的 (r^2) 矩阵 |
prefix.ld.result |
下游 | 组合突变的 LD 打分与推荐标志 |
8. 小结
| 环节 | 方法 | 输出 |
|---|---|---|
| 输入 | 逆折叠 MSA（.fa） |
参考 + (N) 条采样序列 |
| 编码 | 最优密码子表 → 伪 DNA | 每条序列 (3L) bp |
| 变异检测 | 逐碱基与 REF 比较 → VCF | 有变异位点 + 基因型 |
| LD 计算 | PLINK --r2 square |
(M \times M) 的 (r^2) 矩阵 |
| 组合筛选 | 子矩阵 Mean Pairwise LD ≥ 0.5 | .ld.result |
核心思想：AiCE 不直接用 LD 预测实验活性，而是用 LD 矩阵识别「在结构兼容的采样序列中倾向于一起出现的突变位点组合」，作为规避负向上位的统计共现过滤器。
9. 延伸阅读
- 酶改造-modelpaper-AiCE：AiCE 框架总览与全部指标链
- AiCE 官方仓库：https://github.com/ScorpioLea/AiCE
- 检索关键词：
linkage disequilibrium protein engineering,AiCE LD matrix,PLINK r2,pseudo-DNA MSA,combinatorial mutation epistasis