第4章 搞定PDB文件：用Biopython把结构数据抓出来

做 RMSD 计算，第一步就是读文件。PDB 格式看着密密麻麻，其实规律很简单。我刚开始搞分子模拟那会儿，对着 PDB 文档啃了半天，后来发现用 Biopython 几行代码就能把原子坐标全扒出来。这一章咱就把 PDB 格式拆开看，再手把手用 Biopython 把它变成 Python 对象。

4.1 PDB 格式长什么样？—— 一行一行拆给你看

PDB 是纯文本，每行 80 个字符，用关键字开头标识记录类型。最常见的就几种：

ATOM 行的格式是固定的，拿一个典型行举例：

ATOM      1  N   ALA A   1      -1.524  21.642  11.244  1.00  9.17           N  

从左往右分别是：记录名、序号、原子名、残基名、链ID、残基序号、x、y、z 坐标、占有率、温度因子、元素符号。嗯，坐标就在 7-9 列，其他够用就行。我记得有次调试一个脚本，死活读不对坐标，最后发现是 PDB 文件里用了非标准空格，Biopython 一把就解决了。

4.2 用 Biopython 解析 PDB —— 一句话搞定

你没必要自己写解析器。Biopython 的 Bio.PDB 模块就是干这个的。安装很简单：

pip install biopython

用法：

from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)  
structure = parser.get_structure("my_protein", "1ubq.pdb")

这里 QUIET=True 表示不打印警告信息。我个人习惯都加上，避免输出被日志淹没。 get_structure 返回一个 Structure 对象，它里面嵌套了 Model、Chain、Residue、Atom 四级结构。说白了就像一个俄罗斯套娃：

Structure → Model(0) → Chain(A) → Residue(1) → Atom(N)

为什么要分 Model？因为 NMR 结构往往存了多个构象，每个构象是一个 Model。X-ray 晶体结构通常只有一个 Model。我建议你始终指定 model = 0，除非你知道自己在处理多模型文件。

4.3 从结构里提取原子坐标 —— 实战环节

要算 RMSD，先得拿到参考和待比较的坐标。比如提取所有 Cα 原子的坐标：

coords = []
for model in structure:
    for chain in model:
        for residue in chain:
            if residue.has_id("CA"):
                atom = residue["CA"]
                coords.append(atom.get_vector())

这里 atom.get_vector() 返回一个 Vector 对象，可以转成 numpy 数组。为什么强调 Cα？因为 RMSD 通常只比骨架原子，Cα 是主链的代表。你想想看，如果比所有重原子，计算量大不说，侧链构象差异也会让 RMSD 变很大。

提取出来的坐标列表，直接用 numpy.array 做成矩阵：

import numpy as np
coords_np = np.array([v.get_array() for v in coords])

4.4 踩坑与技巧

• 多模型陷阱：NMR 结构有多个 Model，默认遍历会得到所有模型。记得指定 structure[0] 取第一个模型。
• 残基编号可能是负的：某些 PDB 文件会用负数表示 N 端插入的残基，Biopython 直接用整数存，注意不要当成过滤条件。
• HETATM 和 ATOM 要一起处理：如果只遍历 ATOM，会漏掉配体或修饰氨基酸。建议用 chain.get_residues()，它默认包含所有残基。
• 坐标精度：PDB 坐标是三位小数，Biopython 读取后存为浮点数，足够 RMSD 计算用。

我当年第一次用 Biopython 时，忘了加 QUIET=True，终端刷满警告信息，差点以为程序崩了。后来养成习惯，解析前先读几行原始文件确认格式，再用 Parser 一把过。

下面这张图总结了 PDB 解析的完整流程，从原始文本到坐标数组，一目了然。