跑完 AlphaFold 之后:PDB/MMCIF 输出文件里到底该找什么

跑完 AlphaFold 之后:PDB/MMCIF 输出文件里到底该找什么 跑完 AlphaFold 之后PDB/MMCIF 输出文件里到底该找什么【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold你刚跑完一次 AlphaFold 单体预测输出目录里躺着ranked_0.pdb、unrelaxed_model_1.pdb、confidence_model_1.json、pae_model_1.json……打开终端一堆文件不知道该从哪下手到底该看哪个 PDBpLDDT 分数藏在哪一列用 Python 怎么把 MMCIF 打开这篇就是回答这三个问题的。先说结论日常分析就看ranked_0.pdb质量指标看 B 因子和 JSON 文件格式转换和模型比较用 Biopython 几行代码搞定。如果你只关心预测结果可不可信直接跳到质量指标怎么看那一节。拿到你要的结构文件输出目录里文件不少但真正常用的就下面这几个文件一句话说明ranked_0.pdb置信度最高的结构默认已过松弛最常用relaxed_model_*.pdb各模型经 Amber 能量松弛后的结构unrelaxed_model_*.pdb模型直接输出的原始结构ranked_*.cif每个ranked_*.pdb对应的 MMCIF 版本confidence_model_*.json逐残基 pLDDT 分数及置信度分类pae_model_*.jsonPAE 残基对误差矩阵unrelaxed/relaxed/ranked三者的关系其实很简单模型先吐出 5 个原始预测unrelaxed其中置信度最好的会被 Amber 做能量松弛relaxed默认只松弛最优的那个--models_to_relaxbest最后所有模型按置信度重新排序写出ranked_0.pdb~ranked_4.pdb——所以ranked_0.pdb默认就是松弛后的最优模型日常分析拿它就对了。想让输出更符合你的需求run_alphafold.py里有几个参数值得知道参数作用--output_dir结果保存目录--model_presetmonomer/multimer等模型预设--models_to_relaxall/best/none控制松弛范围--num_multimer_predictions_per_model多聚体模式下每个模型的预测次数PDB 和 MMCIF到底该用哪个先搞清楚一件事两种格式装的是同一套原子坐标差别在壳。维度PDBMMCIF结构固定列宽文本键值对表格扩展性弱难加新字段强可容纳任意元数据元数据有限丰富作者、软件、全局 pLDDT 等工具兼容可视化软件普遍支持蛋白质数据库当前标准格式手动编辑方便费劲在 AlphaFold 的输出里两者是严格一一对应的每个ranked_*.pdb旁边都有一个同名.cif比如ranked_0.pdb↔ranked_0.cif。选择建议本地分析、可视化、对接前处理用 PDB要写进论文、提交数据库、或者需要完整元数据比如 AlphaFold 写进 MMCIF 的全局 pLDDT 指标用 MMCIF。PDB 里每条 ATOM 记录长这样注意标出来的几处关键字段ATOM 1 N ALA A 1 23.12 18.45 12.78 1.00 98.76 N原子名是 N残基 ALA、链 A、残基号 1接着是 xyz 坐标。最后一个数 98.76 是 B 因子——在 AlphaFold 的输出里B 因子被复用来存该残基的 pLDDT 分数这是后面提取质量指标的关键。用 Python 打开结构文件Biopython 是这类工作的标配pip install biopython就能用。下面这段代码打开ranked_0.pdb看一下结构里有多少残基from Bio.PDB import PDBParser # 解析 PDB 文件并查看结构概况 parser PDBParser(QUIETTrue) structure parser.get_structure(target, ranked_0.pdb) model next(iter(structure)) # 取第一个模型 chain next(iter(model)) # 取第一条链 residues list(chain.get_residues()) print(f残基数: {len(residues)}) print(前 5 个残基:, [r.get_resname() for r in residues[:5]])再来看最有用的操作提取每个残基的 pLDDT 并画分布图。pLDDT 就在 CA 原子的 B 因子里一行列表推导就能拿到import matplotlib.pyplot as plt from Bio.PDB import PDBParser parser PDBParser(QUIETTrue) structure parser.get_structure(target, ranked_0.pdb) # 取每个残基 CA 原子的 B 因子作为 pLDDT plddt [r[CA].get_bfactor() for r in structure.get_residues() if r.get_id()[2] and CA in r] plt.plot(range(1, len(plddt) 1), plddt, b-, labelpLDDT) plt.axhline(90, colorr, ls--, label高置信阈值 90) plt.axhline(50, colororange, ls--, label低置信阈值 50) plt.ylabel(pLDDT); plt.xlabel(Residue) plt.title(AlphaFold 逐残基 pLDDT 分数); plt.legend(); plt.show()注意这里CA in r这个过滤不能省低置信度区域的残基可能压根没有 CA 原子。质量指标怎么看pLDDT 是逐残基的预测置信度0–100 分官方代码里的分档逻辑是90 高置信结构基本可信、50–90 中等骨架大致对侧链位置别太当真、50 低置信通常是无序区、柔性环或结合界面。PAEpae_model_*.json则是残基对级别的误差估计主要回答两个结构区域之间的相对位置有没有把握对多结构域和多聚体界面判断特别有用。想知道我的结果可不可信跑一遍这段统计就够了import numpy as np from Bio.PDB import PDBParser parser PDBParser(QUIETTrue) structure parser.get_structure(target, ranked_0.pdb) plddt np.array([r[CA].get_bfactor() for r in structure.get_residues() if CA in r]) # 均值、中位数和三个置信度档位占比 print(f均值 pLDDT: {plddt.mean():.1f}, 中位数: {np.median(plddt):.1f}) print(f高置信 (90): {(plddt 90).mean() * 100:.1f}%) print(f中置信 (50-90): {((plddt 50) (plddt 90)).mean() * 100:.1f}%) print(f低置信 (50): {(plddt 50).mean() * 100:.1f}%) 经验法则均值低于 70 或者低置信占比超过 20%先怀疑预测本身别急着拿结构去做下游分析。confidence_model_*.json里也存了同样的分数和 D/L/M/H 分类不想自己算时直接读它。格式转换PDB ↔ MMCIFBiopython 内置了双向转换PDB 转 MMCIF 就这几行from Bio.PDB import PDBParser from Bio.PDB.mmcifio import MMCIFIO def pdb_to_mmcif(pdb_file, cif_file): 把 PDB 坐标文件转成 MMCIF 格式。 structure PDBParser(QUIETTrue).get_structure(target, pdb_file) io MMCIFIO() io.set_structure(structure) io.save(cif_file) pdb_to_mmcif(ranked_0.pdb, ranked_0_converted.cif)反向的 MMCIF 转 PDB 同理换成MMCIFParser解析、PDBIO写出即可。但注意这是坐标层面的转换AlphaFold 原版.cif里那些元数据见下一节的坑换不了。多模型比较与结构叠加ranked_0到ranked_4到底差多少用 Biopython 的Superimposer算 CA 原子 RMSD 即可from Bio.PDB import PDBParser from Bio.PDB.Superimposer import Superimposer parser PDBParser(QUIETTrue) m1 parser.get_structure(m1, ranked_0.pdb) m2 parser.get_structure(m2, ranked_1.pdb) ca1 [a for a in m1.get_atoms() if a.get_name() CA] ca2 [a for a in m2.get_atoms() if a.get_name() CA] superimposer Superimposer() superimposer.set_atoms(ca1, ca2) print(fRMSD: {superimposer.rms:.2f} Å)RMSD 怎么解读小于 1 Å 说明两个模型几乎一致结构很稳1–2 Å 有局部差异常见于柔性区超过 2 Å 就值得认真看看差在哪了——可能某段在两个模型里整体移走了。想看细节叠加之后逐残基算位移import numpy as np superimposer.apply(ca2) # 把 m2 叠加到 m1 坐标系下 # 叠加后每个残基 CA 的位移 displ [np.linalg.norm(a1.get_coord() - a2.get_coord()) for a1, a2 in zip(ca1, ca2)] print(f平均位移: {np.mean(displ):.2f} Å, 最大: {np.max(displ):.2f} Å)位移大的残基如果正好对应 pLDDT 低谷说明模型差异来自不确定的区域放心若高置信区也在大幅移动那这个区域的结构要打个问号。常见坑与快速排错问题pLDDT 读出来全是 0 或位数不对→ 原因读错了原子或者没过滤缺 CA 的残基。 → 解决只取CA in r的残基用r[CA].get_bfactor()对照confidence_model_*.json里的数值验证。问题ranked_0.pdb 和 model_1.pdb 对不上号→ 原因ranked_*是按置信度重排的ranked_0未必对应model_1且ranked_0默认是松弛后的结构坐标与unrelaxed有小幅差异。 → 解决查ranking_debug.json里面有每个ranked_i对应的原始模型名和排序依据。问题用 Biopython 转出来的 MMCIF 少了东西→ 原因MMCIFIO只写坐标块不会写 AlphaFold 原版.cif里的_ma_qa_metric全局 pLDDT、_ma_software_group、引用信息等元数据。 → 解决对外发布或入库直接用 AlphaFold 原生的ranked_*.cifBiopython 转换只用于纯坐标场景。问题解析 MMCIF 时链标识变成A或丢了→ 原因旧版 Biopython 对label_asym_id的处理有坑。 → 解决升级 biopython仍不放心就用grep ^_atom_site.label_asym_id ranked_0.cif | sort -u直接核对原始文件。问题PDB 里的残基数比序列短→ 原因pLDDT 50 的无序区残基原子缺失属正常现象不是文件损坏。 → 解决序列长度以输入的 FASTA 为准别拿 PDB 残基数当全长。一句话总结记住一条主线日常分析看ranked_0.pdb置信度最高的松弛结构pLDDT 藏在 CA 原子的 B 因子里Biopython 负责解析、转换和模型比较发布级元数据直接用 AlphaFold 原生的.cif。延伸阅读README.md官方仓库文档输出文件结构完整说明docs/technical_note_v2.3.0.mdAlphaFold 2.3 技术说明notebooks/AlphaFold.ipynb官方演示 Notebook含结果可视化流程【免费下载链接】alphafoldOpen source code for AlphaFold 2.项目地址: https://gitcode.com/GitHub_Trending/al/alphafold创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考