如何给 AlphaFold 预测打分RMSD 和 lDDT 的正确打开方式【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold你刚跑完一批 AlphaFold 预测输出目录里躺着model_1.npy到model_5.npy五个坐标文件形状都是[残基数, 37, 3]——每个残基 37 种原子类型37 这个数字来自 alphafold/common/residue_constants.py 里的atom_types残基没有某种原子时对应位置补零。你的任务对实验结构给每个模型打分排序找出最好的那个。打分靠两个数RMSD均方根偏差管整体构象lDDT局部距离差异测试管局部质量。下面我们就用这两个数把排序任务完整跑一遍边算边讲每个数的来历。想自己动手验证的话先把仓库拉下来git clone https://gitcode.com/GitHub_Trending/al/alphafold。先把对齐做对Kabsch 旋转的 30 行核心逻辑RMSD 这一侧的活儿不复杂但它的第一步是 90% 的人第一次都会踩错的地方算之前必须先把两个结构转到同一个坐标系里去。为什么必须旋转之后再算说白了就是把两个结构放进同一个坐标系后逐原子算差、平方、平均、开根号。如果预测结构只是在空间里转了几十度逐原子距离会炸到很大RMSD 直接飙到 10 Å 以上——可结构本身一点没变。所以它的严格定义其实是找一组最优旋转 $R$ 和平移 $t$让原子偏差最小$$RMSD(R, t) \sqrt{\frac{1}{N} \sum_{i1}^{N} | \mathbf{x}_i - (R\mathbf{x}_i t) |^2}$$这个公式问的是最佳叠加之后每个原子平均离它的对应原子有多远。注意平方项——大偏差比如甩飞的末端会被放大这是 RMSD 怕局部翻车的根源后面会专门说。用 SVD 求最优旋转好消息是最优解有闭式解对就是那个 Kabsch 算法。技巧是先把两个结构都平移到各自质心平移 $t$ 就此消失再对交叉协方差矩阵做 SVD 求 $R$。写成代码核心二十行以内def kabsch_rmsd(pred, target): # pred / target: [N, 3]一一对应的原子坐标取 Cα 即可 c_p, c_t pred.mean(axis0), target.mean(axis0) p, t pred - c_p, target - c_t # SVD 给出最优旋转 u, _, vh np.linalg.svd(p.T t) d np.linalg.det(u vh) r u np.diag([1.0, 1.0, d]) vh # d0 时翻转最小奇异方向防镜像 aligned p r return float(np.sqrt(np.mean(np.sum((aligned - t)**2, axis1))))这段代码干的其实就四件事中心化、SVD、防镜像、算 RMSD。注意p r的乘法方向行向量在右侧乘旋转矩阵等价于列向量左乘 $R^T$换成列向量习惯写法代码长相会不同但数学完全一样。手动验证旋转矩阵是否右手系⚠️ 整段 Kabsch 里最容易出事的是一行d np.linalg.det(u vh)。当两组点云互为近似镜像时实际里更常见的是预测缺了某个结构域或原子对错了序列SVD 会直接解出一个 det -1 的反射矩阵。照用不误的话等于把一个结构翻进了肚子里面算出来的 RMSD 是个假值。diag([1, 1, d])这个技巧就是翻转最小奇异值方向、把反射强扭成正规旋转det 1。你完全可以手动验证这一步检查np.linalg.det(r)是否约等于 1、r.T r是否约等于单位阵。这个坑第一次跑几乎必踩症状就是数字怎么都对不上但就是查不出原因。Cα 还是全原子实际跑的时候建议先只取 Cα稳、快、且足够代表主链构象。Cα 的索引直接查atom_orderatom_order[CA]是 1所以coords[:, 1, :]就是整条 Cα 轨迹。要算全原子 RMSD 的话用原子掩码把一边缺原子的位置剔除别让零值坐标参与平均。这个仓库里 relax 阶段用的就是同款思路alphafold/relax/relax.py 第 70 行用np.sqrt(np.sum((start_pos - min_pos)**2) / N)比较能量最小化前后的坐标把结果作为收敛指标存进debug_data[rmsd]——它比较的是同一坐标系的两次采样所以省掉了 SVD但逐原子差、平方、平均、开根号的内核一模一样。到这里整体分就齐了。但五个模型跑完后你大概率会撞见怪事两个模型 RMSD 分别是 2.1 Å 和 2.3 Å挨得很近打开结构一看一个核心域打包得漂漂亮亮另一个末端却是一团乱麻。这时候就得请出第二个指标了。为什么 RMSD 说 2.1ÅlDDT 却说只有 0.55lDDT 问的是另一个问题 两个指标最关键的区别RMSD 问每个原子在哪必须对齐后才能算lDDT 问成对距离的分布对不对完全不需要对齐。lDDT 的实现在仓库里就是 alphafold/model/lddt.py 的一个函数。训练主流程里它作为损失项lDDT_CA的一部分参与监督但你对预测结果做事后评估时直接调用它就行。这个函数的核心逻辑非常直白# 摘自 model/lddt.py 的核心片段输入形状 (batch, length, 3) 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)) # 只统计在真实结构里足够近的点对cutoff 默认 15 Å并排除自配对 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]))) 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)) norm 1. / (1e-10 jnp.sum(dists_to_score, axis(-2, -1))) return norm * (1e-10 jnp.sum(dists_to_score * score, axis(-2, -1)))距离分箱 0.5/1/2/4Å 的来历看打分公式。对任意一个原子对设两个结构里它俩的距离差为 $\Delta d |d_{\text{true}} - d_{\text{pred}}|$得分是$$s 0.25\left( \mathbb{1}[\Delta d 0.5] \mathbb{1}[\Delta d 1.0] \mathbb{1}[\Delta d 2.0] \mathbb{1}[\Delta d 4.0] \right)$$翻译成人话距离差在 4 Å 以内至少拿 0.25差得越小叠加的档越多0.5 Å 以内直接满分 1.0。为什么是这四档、各占 0.25这是 lDDT 原始论文Mariani et al., 2013, Bioinformatics的设计——把连续误差变成阶梯式达标0.5 Å 是主链热涨落的量级1 Å 是侧链呼吸2 Å 算构象变化4 Å 基本就是这个接触已经断了。分箱比裸误差对结构噪声的容忍度好得多。有个细节要心里有数源码注释明确说了这是近似 lDDT原论文里对物理可行性比如键长违背的修正项没算进去所以它和 CASP 官方报的 lDDT 数值不会严格相等。末端甩飞了为什么它不爆表用一个迷你场景就能理解两个数为什么会打架你的预测里 N 端前 20 个残基整体甩飞了 5 Å但 120 个残基的核心域和实验结构贴合得严丝合缝。RMSD 把 20 个残基 5 Å 的偏差平方后摊进 140 个残基数字掉到 1.8 Å 附近——看着还不错可末端其实是坏的lDDT 看成对接触核心域内部绝大部分点对距离差都小于 0.5 Å拿满分少数跨末端的接触断了记 0 分。最终 0.87 告诉你核心可信末端存疑。反过来也成立如果某个结构域相对实验结构整体转了 2 Å域内接触完好、只有域间距离变了lDDT 只在跨域接触上扣分RMSD 却会被整体拖低。所以看到两个数打架时先问自己一句两个结构到底差在哪——RMSD 告诉你整体偏了lDDT 告诉你哪类接触偏了。再看lddt()里的掩码设计只有两端原子在真实结构中都存在true_points_mask、且真实距离小于 cutoff默认 15 Å的点对才进分母。这让它对预测缺原子实验结构缺侧链天然免疫——这正是 RMSD 的软肋它得先手动做掩码两个结构长度不同时还得先做序列比对找对应。批量跑 500 个模型时该信哪个数如果规模从 5 个模型变成 500 个 target、每个 5 个模型就得把前面两节的流程固化成一条流水线。从原始坐标到最终评分整条链路长这样什么任务信 RMSD什么任务信 lDDT做全局构象比较时信 RMSD比较变构蛋白的开放/闭合态、聚类 MD 轨迹、监控 relax 能量最小化是否收敛alphafold/relax/relax.py 干的就是这件事。它简单、快对齐后就是 O(N)单位是 Å物理意义直观。评估局部质量时信 lDDT结构有缺失区域、想知道哪些残基可信、或者要比较长度不同的模型。它免对齐、对缺原子鲁棒加per_residueTrue还能拿到逐残基质量曲线和预测自带的 pLDDT 置信曲线正好配套看两个模型 pLDDT 打平时用 lDDT 决胜负预测末端 pLDDT 本来就低、lDDT 末端也低那末端就别信。 一条实操经验RMSD 排大形状lDDT 查细节。RMSD 高但 lDDT 高的模型多半是全局取向或结构域排布的问题打开结构看一眼就知道RMSD 低但 lDDT 低的模型是局部细节成片出问题别选。另外别忘了成本差异lDDT 要建完整距离矩阵长蛋白质和复合物上开销明显更大批量场景先 RMSD 快筛、再 lDDT 复核也更省机器。还有个容易忽略的坑逐残基 lDDT 的平均值并不等于全局 lDDT【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 如何给 AlphaFold 预测打分？RMSD 和 lDDT 的正确打开方式
如何给 AlphaFold 预测打分RMSD 和 lDDT 的正确打开方式【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold你刚跑完一批 AlphaFold 预测输出目录里躺着model_1.npy到model_5.npy五个坐标文件形状都是[残基数, 37, 3]——每个残基 37 种原子类型37 这个数字来自 alphafold/common/residue_constants.py 里的atom_types残基没有某种原子时对应位置补零。你的任务对实验结构给每个模型打分排序找出最好的那个。打分靠两个数RMSD均方根偏差管整体构象lDDT局部距离差异测试管局部质量。下面我们就用这两个数把排序任务完整跑一遍边算边讲每个数的来历。想自己动手验证的话先把仓库拉下来git clone https://gitcode.com/GitHub_Trending/al/alphafold。先把对齐做对Kabsch 旋转的 30 行核心逻辑RMSD 这一侧的活儿不复杂但它的第一步是 90% 的人第一次都会踩错的地方算之前必须先把两个结构转到同一个坐标系里去。为什么必须旋转之后再算说白了就是把两个结构放进同一个坐标系后逐原子算差、平方、平均、开根号。如果预测结构只是在空间里转了几十度逐原子距离会炸到很大RMSD 直接飙到 10 Å 以上——可结构本身一点没变。所以它的严格定义其实是找一组最优旋转 $R$ 和平移 $t$让原子偏差最小$$RMSD(R, t) \sqrt{\frac{1}{N} \sum_{i1}^{N} | \mathbf{x}_i - (R\mathbf{x}_i t) |^2}$$这个公式问的是最佳叠加之后每个原子平均离它的对应原子有多远。注意平方项——大偏差比如甩飞的末端会被放大这是 RMSD 怕局部翻车的根源后面会专门说。用 SVD 求最优旋转好消息是最优解有闭式解对就是那个 Kabsch 算法。技巧是先把两个结构都平移到各自质心平移 $t$ 就此消失再对交叉协方差矩阵做 SVD 求 $R$。写成代码核心二十行以内def kabsch_rmsd(pred, target): # pred / target: [N, 3]一一对应的原子坐标取 Cα 即可 c_p, c_t pred.mean(axis0), target.mean(axis0) p, t pred - c_p, target - c_t # SVD 给出最优旋转 u, _, vh np.linalg.svd(p.T t) d np.linalg.det(u vh) r u np.diag([1.0, 1.0, d]) vh # d0 时翻转最小奇异方向防镜像 aligned p r return float(np.sqrt(np.mean(np.sum((aligned - t)**2, axis1))))这段代码干的其实就四件事中心化、SVD、防镜像、算 RMSD。注意p r的乘法方向行向量在右侧乘旋转矩阵等价于列向量左乘 $R^T$换成列向量习惯写法代码长相会不同但数学完全一样。手动验证旋转矩阵是否右手系⚠️ 整段 Kabsch 里最容易出事的是一行d np.linalg.det(u vh)。当两组点云互为近似镜像时实际里更常见的是预测缺了某个结构域或原子对错了序列SVD 会直接解出一个 det -1 的反射矩阵。照用不误的话等于把一个结构翻进了肚子里面算出来的 RMSD 是个假值。diag([1, 1, d])这个技巧就是翻转最小奇异值方向、把反射强扭成正规旋转det 1。你完全可以手动验证这一步检查np.linalg.det(r)是否约等于 1、r.T r是否约等于单位阵。这个坑第一次跑几乎必踩症状就是数字怎么都对不上但就是查不出原因。Cα 还是全原子实际跑的时候建议先只取 Cα稳、快、且足够代表主链构象。Cα 的索引直接查atom_orderatom_order[CA]是 1所以coords[:, 1, :]就是整条 Cα 轨迹。要算全原子 RMSD 的话用原子掩码把一边缺原子的位置剔除别让零值坐标参与平均。这个仓库里 relax 阶段用的就是同款思路alphafold/relax/relax.py 第 70 行用np.sqrt(np.sum((start_pos - min_pos)**2) / N)比较能量最小化前后的坐标把结果作为收敛指标存进debug_data[rmsd]——它比较的是同一坐标系的两次采样所以省掉了 SVD但逐原子差、平方、平均、开根号的内核一模一样。到这里整体分就齐了。但五个模型跑完后你大概率会撞见怪事两个模型 RMSD 分别是 2.1 Å 和 2.3 Å挨得很近打开结构一看一个核心域打包得漂漂亮亮另一个末端却是一团乱麻。这时候就得请出第二个指标了。为什么 RMSD 说 2.1ÅlDDT 却说只有 0.55lDDT 问的是另一个问题 两个指标最关键的区别RMSD 问每个原子在哪必须对齐后才能算lDDT 问成对距离的分布对不对完全不需要对齐。lDDT 的实现在仓库里就是 alphafold/model/lddt.py 的一个函数。训练主流程里它作为损失项lDDT_CA的一部分参与监督但你对预测结果做事后评估时直接调用它就行。这个函数的核心逻辑非常直白# 摘自 model/lddt.py 的核心片段输入形状 (batch, length, 3) 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)) # 只统计在真实结构里足够近的点对cutoff 默认 15 Å并排除自配对 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]))) 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)) norm 1. / (1e-10 jnp.sum(dists_to_score, axis(-2, -1))) return norm * (1e-10 jnp.sum(dists_to_score * score, axis(-2, -1)))距离分箱 0.5/1/2/4Å 的来历看打分公式。对任意一个原子对设两个结构里它俩的距离差为 $\Delta d |d_{\text{true}} - d_{\text{pred}}|$得分是$$s 0.25\left( \mathbb{1}[\Delta d 0.5] \mathbb{1}[\Delta d 1.0] \mathbb{1}[\Delta d 2.0] \mathbb{1}[\Delta d 4.0] \right)$$翻译成人话距离差在 4 Å 以内至少拿 0.25差得越小叠加的档越多0.5 Å 以内直接满分 1.0。为什么是这四档、各占 0.25这是 lDDT 原始论文Mariani et al., 2013, Bioinformatics的设计——把连续误差变成阶梯式达标0.5 Å 是主链热涨落的量级1 Å 是侧链呼吸2 Å 算构象变化4 Å 基本就是这个接触已经断了。分箱比裸误差对结构噪声的容忍度好得多。有个细节要心里有数源码注释明确说了这是近似 lDDT原论文里对物理可行性比如键长违背的修正项没算进去所以它和 CASP 官方报的 lDDT 数值不会严格相等。末端甩飞了为什么它不爆表用一个迷你场景就能理解两个数为什么会打架你的预测里 N 端前 20 个残基整体甩飞了 5 Å但 120 个残基的核心域和实验结构贴合得严丝合缝。RMSD 把 20 个残基 5 Å 的偏差平方后摊进 140 个残基数字掉到 1.8 Å 附近——看着还不错可末端其实是坏的lDDT 看成对接触核心域内部绝大部分点对距离差都小于 0.5 Å拿满分少数跨末端的接触断了记 0 分。最终 0.87 告诉你核心可信末端存疑。反过来也成立如果某个结构域相对实验结构整体转了 2 Å域内接触完好、只有域间距离变了lDDT 只在跨域接触上扣分RMSD 却会被整体拖低。所以看到两个数打架时先问自己一句两个结构到底差在哪——RMSD 告诉你整体偏了lDDT 告诉你哪类接触偏了。再看lddt()里的掩码设计只有两端原子在真实结构中都存在true_points_mask、且真实距离小于 cutoff默认 15 Å的点对才进分母。这让它对预测缺原子实验结构缺侧链天然免疫——这正是 RMSD 的软肋它得先手动做掩码两个结构长度不同时还得先做序列比对找对应。批量跑 500 个模型时该信哪个数如果规模从 5 个模型变成 500 个 target、每个 5 个模型就得把前面两节的流程固化成一条流水线。从原始坐标到最终评分整条链路长这样什么任务信 RMSD什么任务信 lDDT做全局构象比较时信 RMSD比较变构蛋白的开放/闭合态、聚类 MD 轨迹、监控 relax 能量最小化是否收敛alphafold/relax/relax.py 干的就是这件事。它简单、快对齐后就是 O(N)单位是 Å物理意义直观。评估局部质量时信 lDDT结构有缺失区域、想知道哪些残基可信、或者要比较长度不同的模型。它免对齐、对缺原子鲁棒加per_residueTrue还能拿到逐残基质量曲线和预测自带的 pLDDT 置信曲线正好配套看两个模型 pLDDT 打平时用 lDDT 决胜负预测末端 pLDDT 本来就低、lDDT 末端也低那末端就别信。 一条实操经验RMSD 排大形状lDDT 查细节。RMSD 高但 lDDT 高的模型多半是全局取向或结构域排布的问题打开结构看一眼就知道RMSD 低但 lDDT 低的模型是局部细节成片出问题别选。另外别忘了成本差异lDDT 要建完整距离矩阵长蛋白质和复合物上开销明显更大批量场景先 RMSD 快筛、再 lDDT 复核也更省机器。还有个容易忽略的坑逐残基 lDDT 的平均值并不等于全局 lDDT【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考