4、使用Python解析PDB文件：读取原子坐标、链、残基信息
我个人经验告诉我，解析PDB文件这事，确实是个基本功。搞RMSD计算，第一步就是把结构文件里的坐标、链信息、残基编号这些数据给提出来。别看PDB文件就是个文本文件，里面的门道还是不少的。
做RMSD评估，说白了就是对比两组原子坐标。但你想想看，坐标从哪来？结构里是不是只有A链？水分子要不要保留？这些都得在解析这一步处理好。我在项目里遇到过，同一个晶体结构，不同实验室修出来的版本，链标识和残基编号可能都不一样——嗯，这里要注意。
PDB文件长什么样
一个标准PDB文件是一行一行的文本记录。每行以特定的记录类型开头，比如ATOM、HETATM、TER、END。我们最关心的就是ATOM记录。
典型的ATOM记录格式（固定列宽，注意是列列列，不是逗号分隔）：
ATOM 1 N ALA A 1 8.342 16.044 22.408 1.00 0.00 N
ATOM 2 CA ALA A 1 9.062 15.309 21.326 1.00 0.00 C
ATOM 3 C ALA A 1 8.162 15.028 20.129 1.00 0.00 C
ATOM 4 O ALA A 1 7.314 14.175 20.228 1.00 0.00 O每一列对应什么？我习惯按起始位置来记（从1开始数）：
| 起始列 | 结束列 | 字段 | 说明 |
|---|---|---|---|
| 1 | 6 | 记录类型 | ATOM 或 HETATM |
| 7 | 11 | 原子序号 | 整型，从1递增 |
| 13 | 16 | 原子名称 | 如 CA、N、O |
| 17 | 17 | 备选位置标识 | A/B 或空格 |
| 18 | 20 | 残基名称 | 如 ALA、ILE |
| 22 | 22 | 链标识 | 如 A、B |
| 23 | 26 | 残基序号 | 整型 |
| 31 | 38 | x 坐标 | 单位：埃 |
| 39 | 46 | y 坐标 | 单位：埃 |
| 47 | 54 | z 坐标 | 单位：埃 |
用Python自写解析器
解析PDB文件这事儿，我推荐两个路子。一个是直接用Biopython，省事。另一个是自己写个简单的解析器，轻量可控。
先看自己写的版本。我个人偏爱这个方式，因为不用多装一个包，而且逻辑完全在掌控之下。
def parse_pdb(filepath):
atoms = []
with open(filepath, 'r') as f:
for line in f:
if line.startswith(('ATOM', 'HETATM')):
record_type = line[0:6].strip()
serial = int(line[6:11].strip())
atom_name = line[12:16].strip()
alt_loc = line[16:17].strip()
res_name = line[17:20].strip()
chain = line[21:22].strip()
res_seq = int(line[22:26].strip())
x = float(line[30:38].strip())
y = float(line[38:46].strip())
z = float(line[46:54].strip())
atoms.append({
'serial': serial,
'atom_name': atom_name,
'res_name': res_name,
'chain': chain,
'res_seq': res_seq,
'x': x,
'y': y,
'z': z
})
return atoms用Biopython解析
如果你项目里已经装了Biopython，那解析起来就轻松多了。我从Biopython开始接触结构生物信息，确实方便。
from Bio.PDB import PDBParser
parser = PDBParser(QUIET=True) # QUIET=True 可以少吐一堆警告
structure = parser.get_structure('my_protein', 'example.pdb')
for model in structure:
for chain in model:
print(f"Chain: {chain.get_id()}")
for residue in chain:
print(f" Residue: {residue.get_resname()} {residue.get_id()[1]}")
for atom in residue:
coord = atom.get_vector()
print(f" Atom: {atom.get_name()} ({coord[0]:.3f}, {coord[1]:.3f}, {coord[2]:.3f})")同样的解析，用Biopython做出来的是一个嵌套结构：Structure → Model → Chain → Residue → Atom。我一般通过这个层级来提取坐标。
| 方法 | 优点 | 缺点 |
|---|---|---|
| 自写解析器 | 无依赖、代码可控、轻量 | 需要自己处理边界情况 |
| Biopython | 功能完整、处理紊乱更稳健 | 额外依赖、封装较厚 |
提取链和残基信息
很多时候我们不只要原子坐标，还需要按链和残基来组织数据。比如做RMSD计算时，通常是对齐两个结构的Cα原子，而且只对齐主链，还要排除水分子。
下面这个函数，提取指定链、指定残基类型（只取骨架原子）的坐标：
def extract_ca_coords(atoms, chain_id='A', exclude_water=True):
coords = []
for atom in atoms:
if atom['chain'] != chain_id:
continue
if exclude_water and atom['res_name'] == 'HOH':
continue
if atom['atom_name'] == 'CA':
coords.append((atom['x'], atom['y'], atom['z']))
return coords这里有个坑：PDB里的水分子残基名是 HOH，有时也写成 WAT。我在项目里碰到过一个PDB混合了两种写法。所以做RMSD之前，一定先去确认数据结构。你的对比对象里是不是混进了溶剂？
解析流程全景图
下面这张SVG图，是我做解析时脑中过的一遍整体流程。从PDB文件开始，按行读入，筛选ATOM记录，再分别提取链、残基、原子坐标，最后按RMSD计算需求做二次筛选。
避坑指南 & 实战经验
- 多模型PDB： 有些NMR结构包含多个模型。Biopython里会返回多个Model对象，你得指定取第一套坐标。自写解析器就要小心了，MODEL记录要正确处理。
- 备选位置： 同一个原子可能有A和B两种构象。RMSD计算时通常只取A或主占位位置。我一般会过滤 alt_loc == ' ' 或依次取A、B。
- 残基插入码： 某些残基序号后面带了插入码（比如 10A、10B），残基序号字段不是纯整数了。解析时建议读成
形式。 - 文件编码： 大部分PDB是ASCII，但偶尔遇到UTF-8 BOM头，open 的时候用 encoding='utf-8' 比较稳。