AlphaFold结构比较实战:RMSD和lDDT到底该怎么选
【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold
两个结构叠在一起,肉眼看着几乎一样,评估分数却差了0.3。问题出在哪?
蛋白质结构相似度评估不是"数字越小越好"这么简单。在AlphaFold结构比较中,RMSD(对应原子的均方根偏差,衡量两套坐标差多远)是每个人先算的数,也是最容易误读的数。它由你选了哪些原子、用什么方式对齐共同决定。
为什么不能只看一个数
一次模型评审上,两个数吵了起来:
甲:预测和实验的骨架RMSD才1.8Å,这模型质量不错。 乙:同一个模型,活性位点里60%的关键残基偏了2Å以上。
两边都对。这个数字是全局平均:局部几个大错,会被几百个正确残基稀释掉;反过来,处处偏1Å的模型,分数比"核心区完美、末端拧了"的模型还差。说白了,平均数表达不了"错在哪"。
所以AlphaFold里是双指标配合:RMSD管整体形状,lDDT(local Distance Difference Test,局部距离差异测试,按局部距离表差多少打分)管局部贴合。下面这张CASP14预测图很直观:绿色是实验结构,蓝色是计算预测,底部标的是GDT分数。
RMSD:坑藏在对齐这一步
仓库里这个数字只在一处直接出现:alphafold/relax/relax.py中Amber弛豫结束后,比较弛豫前后的坐标算偏差,写进debug信息,看优化把结构挪动多少。注意语境:两个结构来自同一条弛豫轨迹,同一个坐标系,不需要对齐。
但拿预测结构比实验结构时,情况完全不同。两个结构的原点和朝向毫无关系,第一步必须先做叠加(superposition,把一点云扣到另一点云上)。叠加不是可选项,它决定数字大小。
完整流程就四步:
pred_ca = pred_pos[:, 1, :] # 取Cα(37原子表里的第1列) true_ca = true_pos[:, 1, :] mask = atom_mask[:, 1] # 排除缺失残基 pred_c = pred_ca[mask] - pred_ca[mask].mean(0) # 质心移到原点 true_c = true_ca[mask] - true_ca[mask].mean(0) R = kabsch(pred_c, true_c) # 解最优旋转矩阵 rmsd = rms(pred_c @ R - true_c) # 对应原子差的均方根这段伪代码把从原始坐标到最终数字的管线压缩了。仓库里最简的版本在alphafold/relax/relax.py:同坐标系,只剩最后一步开根号。
对齐分平移和旋转。平移把质心挪到原点,旋转用Kabsch算法(本质上就是求一个最优旋转矩阵,让两个点云重合得最好)解出来。对齐完,算数字就很简单:
$$\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left|x_i - x_i'\right|^2}$$
白话翻译:每对对应原子的差先平方、求和、除以N、开根号,就是"平均每个原子挪动了多远"。x_i和x_i'是对齐后两个结构里同一个残基的坐标,N是参与计算的原子数。
为什么单独拎出Cα?完整原子表有37列(见alphafold/common/residue_constants.py,这里定义了标准原子顺序),大部分是侧链原子,很多残基根本没有。Cα是每个残基都有的原子,串起来就是蛋白主链:
![AlphaFold
【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold
创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考