如何度量两个蛋白质结构的差距：AlphaFold 源码中 RMSD、lDDT 与 FAPE 三把尺子的实战指南

如何度量两个蛋白质结构的差距：AlphaFold 源码中 RMSD、lDDT 与 FAPE 三把尺子的实战指南

【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold

拿到 AlphaFold 的预测结构后，紧接着的问题几乎都是同一个：这个模型和参考结构到底差多远？也就是蛋白质结构比较（protein structure comparison）里的量化问题。这个仓库给出的答案不是单一指标，而是三把"尺子"：lDDT（local Distance Difference Test，局部距离差异测试）是源码里的一等公民，RMSD（均方根偏差，Root Mean Square Deviation）只在角落里露了一手，还有一把常被忽略的尺子叫 FAPE（Frame Aligned Point Error，框架对齐点误差）。下面直接按源码把它们的实现讲透，并给出动手可用的写法。

先想清楚：三种思路，三种用途

在翻代码之前，先把三者的本质差异摆出来，选错尺子比算错更糟。

• RMSD回答的是"把两个结构叠得最齐之后，每个对应原子平均挪了多远"。它必须先找一个最优刚体变换（平移+旋转），所以结果天然依赖对齐方式。
• lDDT完全绕开坐标对齐：分别算两个结构内部所有原子对的距离，再看这两张距离矩阵差多少。它只比较"结构内部的关系"，天然对整体旋转、平移免疫。
• FAPE介于两者之间：用一组局部坐标系（frame）把两个点云各自投影到局部参照系里再比误差，比全局对齐更抗局部畸变。

一句话选型：评局部质量用 lDDT，报整体差异用 RMSD（必须说清对齐方法），做训练或鲁棒比较时用 FAPE 的思路。

从 lddt.py 提取评分逻辑

lDDT 的实现全部集中在alphafold/model/lddt.py，整个文件只有 90 行左右，核心是一个函数。去掉断言后的主干逻辑如下：

# alphafold/model/lddt.py 核心逻辑（节选并加注释） # 1) 分别算两个结构的"距离矩阵"：结构内部所有点对的欧氏距离 dmat_true = jnp.sqrt(1e-10 + jnp.sum( (true_points[:, :, None] - true_points[:, None, :])**2, axis=-1)) dmat_predicted = jnp.sqrt(1e-10 + jnp.sum( (predicted_points[:, :, None] - predicted_points[:, None, :])**2, axis=-1)) # 2) 圈定要打分的距离对：只统计真实结构里距离小于 cutoff（默认 15 Å）的对， # 且两端原子都必须真实存在（mask），并排除自身对 dists_to_score = ( (dmat_true < cutoff).astype(jnp.float32) * true_points_mask * jnp.transpose(true_points_mask, [0, 2, 1]) * (1. - jnp.eye(dmat_true.shape[1]))) # 3) 每个距离对按"四档阶梯"打分：距离差 < 0.5 / 1 / 2 / 4 Å 各记 0.25 分 dist_l1 = jnp.abs(dmat_true - dmat_predicted) score = 0.25 * ((dist_l1 < 0.5).astype(jnp.float32) + (dist_l1 < 1.0).astype(jnp.float32) + (dist_l1 < 2.0).astype(jnp.float32) + (dist_l1 < 4.0).astype(jnp.float32))) # 4) 在"应评分的距离对"上做加权平均，得到 0~1 的分数 reduce_axes = (-1,) if per_residue else (-2, -1) norm = 1. / (1e-10 + jnp.sum(dists_to_score, axis=reduce_axes)) score = norm * (1e-10 + jnp.sum(dists_to_score * score, axis=reduce_axes))

有三个细节值得圈出来：

1. 筛选基准是真实结构：哪些距离对参与评分，看的是dmat_true < cutoff，而不是预测结构。这个不对称设计保证了评分集合是确定的。
2. 打分是离散的：不是用距离差的连续值，而是 0.5 / 1 / 2 / 4 Å 四档阶梯，每档 0.25 分。所以完全正确的结构拿 1.0 分，误差 0.4 Å 和 4.9 Å 拿 0 分——这是个"及格线"式的度量，对微小误差宽容、对大偏差惩罚干脆。
3. 它是"近似"的 lDDT：docstring 明确写了这里没有包含物理合理性修正项（如键长违反映在得分里的惩罚），拿它和评测体系里严格定义的 lDDT 直接对表会差一截。

另外注意per_residue=True时返回每个残基自己的分数。docstring 特意提醒：全局分数不等于各残基分数的算术平均，因为不同残基参与评分的距离对数量不一样。

lddt 在项目里的真实角色：监督信号，不只是评测工具

很多人以为 lddt.py 是给"事后评估"用的，其实它在管线里另有身份。看alphafold/model/modules.py的置信度模块（约 1093 行）：

• 取结构输出的 Cα 坐标列（37 原子表示中的第 1 列），对预测与真实结构逐残基算 lDDT；
• 用stop_gradient掐断梯度，把它当成标签；
• 把分数离散化到num_bins个桶里做 one-hot，和网络的 logits 算交叉熵。

也就是说，模型被训练成"预测每个残基的 lDDT 该落在哪个桶"，网络输出再经alphafold/common/confidence.py换算成最终报告里的逐残基置信度（pLDDT）。你在结果里看到的那个"某残基置信度 87"，源头就是这里这段距离差打分。

顺带一提，alphafold/model/all_atom.py里的find_optimal_renaming处理了一个更隐蔽的问题：部分氨基酸的原子命名有歧义（典型如缬氨酸的两个甲基 CG1/CG2），它会把两种命名方案各自算一遍 lDDT，取得分更高的那种。如果你的复现脚本里 lDDT 数值莫名偏低，先检查是不是撞上了这种命名歧义。

RMSD：源码里只有一处，Kabsch 对齐要自己补

在整个仓库里搜 RMSD，正经实现只有一处——alphafold/relax/relax.py第 70 行：

rmsd = np.sqrt(np.sum((start_pos - min_pos)**2) / start_pos.shape[0])

注意它的用途：这是 relax（Amber 能量最小化）流程的调试指标，衡量"弛豫把坐标挪动了多少"，属于位移量而非对齐后的结构差异，直接拿来做评估是不对的。

真正做评估用的"对齐后 RMSD"需要 Kabsch 算法求最优旋转，仓库本身没有现成函数（alphafold/model/geometry/提供的是旋转矩阵、四元数等基础件）。下面这个 20 行左右的实现是标准写法，可以直接放进你自己的评测脚本：

import numpy as np def kabsch_rmsd(pred: np.ndarray, true: np.ndarray) -> float: """pred / true 形状均为 (N, 3)，先做 Kabsch 最优刚体对齐，再返回 RMSD""" assert pred.shape == true.shape # 两个点云平移到质心为原点 p = pred - pred.mean(axis=0) t = true - true.mean(axis=0) # 对协方差矩阵做 SVD，奇异向量外积即最优旋转 u, _, vh = np.linalg.svd(p.T @ t) rot = vh.T @ u.T # 行列式为负说明 SVD 给的是"反射"，翻掉最后一列变回纯旋转 if np.linalg.det(rot) < 0: vh[-1, :] *= -1 rot = vh.T @ u.T aligned = p @ rot return float(np.sqrt(np.mean(np.sum((aligned - t) ** 2, axis=1))))

两个使用上的硬性约定：只比较一一对应的原子（先按序列比对锁定残基，再挑 Cα 或同一套原子表），以及缺失原子要用掩码剔除后再算——漏掉任何一条，数字都没有意义。

三把尺子对照表

动手前的几个常见坑

📌别把源码的"近似 lDDT"当标准分。如前所述，lddt.py的 docstring 声明其不含物理可行性修正，跨论文、跨评测系统比较时先确认对方用的是哪个口径。

📌RMSD 不报对齐方法就等于没报。只平移对齐、全原子对齐、Cα 对齐出来的数完全不同，写报告时注明"Kabsch 对齐、Cα 原子、N 个残基"是最低要求。

📌逐残基分数的平均 ≠ 全局分数。接触数不同的残基权重不同，per_residue=True的输出别直接mean()当整体结论。

📌O(N²) 的账要提前算。lDDT 的距离矩阵是length × length，链长到几千残基时显存会先于 CPU 时间出问题，批量评估时注意 batch 维度。

📌歧义原子命名会悄悄拉低分数。复现 lDDT 之前，确认处理了residue_constants里列出的原子重命名情形，或参考find_optimal_renaming的双命名取优策略。

收个尾

在 AlphaFold 这套代码里，"结构有多像"被拆成了三个层次：lDDT 管逐残基的局部几何质量，也是置信度的训练信号；RMSD 管对齐后的整体差异，仓库只留了一个弛豫监控的简化版，评估用的 Kabsch 版本需要按第二节的方式自行补齐；FAPE 则展示了"用局部参照系替代全局对齐"的折中路线。实际项目里三者配合使用：先跑 lDDT 定位低置信区域，再用 Kabsch-RMSD 给出一个可对外引用的整体数字，遇到大环区或柔性尾链时借助 FAPE 式的局部对齐避免被单点拖累。入口都在这份仓库里：run_alphafold.py负责完整预测流程，notebooks/AlphaFold.ipynb适合快速上手复现，指标相关的源码则集中在alphafold/model/目录下的lddt.py、all_atom.py和modules.py三个文件。

【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold