4. PDB文件解析（一）：原子记录、坐标、元素信息提取
好，咱们今天正式开始碰PDB文件的骨头。我个人做了这么多年蛋白质结构分析，发现一个规律：只要是跟结构数据打交道的活儿，八成以上的bug都出在最开始的解析阶段。说白了，你后面算得再花哨，如果原子坐标读进来就是错的，那整个分析就是垃圾进垃圾出。
这一节，我们专门把PDB文件里最核心的记录——原子记录，给剥开了看清楚。
4.1 原子记录长什么样？
你随便打开一个PDB文件，拉到中间，会看到一大片以“ATOM”开头的行。这就是原子记录。如果你看到的是“HETATM”，别慌，它俩格式几乎一样，只是ATOM代表标准残基的原子，HETATM代表配体、水分子或非标准残基。
我刚开始学的时候，觉得这俩没区别，结果有一次跑分子对接，死活找不到配体对接口袋。排查了半天，发现是我的解析脚本把HETATM行全过滤掉了。嗯，那滋味可不好受。
| 记录类型 | 含义 | 常见场景 |
|---|---|---|
| ATOM | 标准残基原子 | 蛋白质主链、侧链 |
| HETATM | 非标准原子 | 配体、水分子、金属离子 |
一条原子记录长这样：
ATOM 1 N ALA A 1 -0.525 1.362 0.000 1.00 0.00 N这一行里信息量巨大。我用一个表格把这列位置和含义给你标清楚，你自己写解析脚本的时候，直接按这个偏移量切字符串就行：
| 列 | 起始 | 结束 | 字段 | 示例 |
|---|---|---|---|---|
| 1-6 | 1 | 6 | 记录名 | ATOM |
| 7-11 | 7 | 11 | 原子序号 | 1 |
| 13-16 | 13 | 16 | 原子名称 | N |
| 17 | 17 | 17 | 备用位置指示符 | |
| 18-20 | 18 | 20 | 残基名称 | ALA |
| 22 | 22 | 22 | 链标识符 | A |
| 23-26 | 23 | 26 | 残基序号 | 1 |
| 27 | 27 | 27 | 插入码 | |
| 31-38 | 31 | 38 | X坐标 | -0.525 |
| 39-46 | 39 | 46 | Y坐标 | 1.362 |
| 47-54 | 47 | 54 | Z坐标 | 0.000 |
| 55-60 | 55 | 60 | 占据率 | 1.00 |
| 61-66 | 61 | 66 | 温度因子(B-factor) | 0.00 |
| 77-78 | 77 | 78 | 元素符号 | N |
4.2 坐标信息才是核心
X、Y、Z这三列是PDB文件最值钱的数据。单位是Å（埃，10⁻¹⁰米）。
这里有个小坑：坐标值可能是正数也可能是负数，但格式必须是右对齐，预留8个字符宽度。假如你把坐标读成整数，那就完蛋了。我建议你理解一下PDB文件的列宽限制，这是为了兼容早期的Fortran程序而设计的。
坐标读取的Python代码，我个人习惯用切片加strip：
def parse_atom_line(line):
rec_name = line[0:6].strip()
if rec_name not in ('ATOM', 'HETATM'):
return None
atom_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_id = line[21:22].strip()
res_seq = int(line[22:26].strip())
i_code = line[26:27].strip()
x = float(line[30:38].strip())
y = float(line[38:46].strip())
z = float(line[46:54].strip())
occupancy = float(line[54:60].strip())
b_factor = float(line[60:66].strip())
element = line[76:78].strip()
if not element:
element = atom_name[0]
return {
'serial': atom_serial,
'name': atom_name,
'residue': res_name,
'chain': chain_id,
'resid': res_seq,
'x': x, 'y': y, 'z': z,
'occupancy': occupancy,
'b_factor': b_factor,
'element': element
}4.3 元素信息的推断与验证
有些PDB文件很规范，元素符号写得清清楚楚。有些文件，尤其是从旧数据库导出的，元素列可能是空的，或者写成了“CA”（钙）跟“CA”（Cα）混在一起。你想想看，如果不做区分，把钙离子当成了碳原子，那计算出来的静电势就全错了。
我曾经帮一个学生调试过类似问题。他跑MD模拟预处理，老是报原子类型不一致。我让他把解析出来元素符号打印出来一看，好家伙，所有金属离子都变成了第一个字母。后来我教他用原子名称加上残基名称来双校验，这才解决问题。
这是我的改进版元素推断逻辑：
def infer_element(atom_name, residue_name):
# 如果是常见配体残基中的金属，直接返回
metal_map = {
'ZN': 'ZN', 'MG': 'MG', 'CA': 'CA',
'FE': 'FE', 'MN': 'MN', 'CU': 'CU'
}
if residue_name in metal_map:
return metal_map[residue_name]
# 常规蛋白原子：取第一个字母
element_map = {
'N': 'N', 'C': 'C', 'O': 'O', 'S': 'S',
'P': 'P', 'H': 'H'
}
first_char = atom_name[0]
return element_map.get(first_char, 'X')4.4 原子记录的整体数据结构
你解析完一整条原子记录后，最好设计一个结构体或字典来装这些信息。我习惯把所有原子存成一个列表，每个元素都是一个字典。这样后续按残基分组、按链切片都很方便。
下面这张SVG图帮你搭建一下整个原子记录解析的思维框架：
解析完一条原子记录后，最终的数据结构推荐这样组织：
[
{
'serial': 1,
'name': 'N',
'residue': 'ALA',
'chain': 'A',
'resid': 1,
'x': -0.525, 'y': 1.362, 'z': 0.000,
'occupancy': 1.00,
'b_factor': 0.00,
'element': 'N'
},
...
]4.5 实战小建议
你如果要自己写PDB解析器，记住这几点：
- 不要用split()切行：PDB格式是固定列宽的，空格位置有含义。split会把空白合并，导致你拿不到正确的列值。
- 备用位置指示符要留着：同一个残基可能有多个构象，指示符是A、B还是空，会影响后续建模。
- 插入码别丢掉：有的PDB文件里残基序号会重复（比如122A、122B），如果不处理插入码，你的残基索引就会乱套。
好了，原子记录的基本解析逻辑已经清楚了。从下一段开始，我们会把解析出来的数据组织成更高级的结构，方便后续做距离计算、序列比对、或者机器学习特征提取。目前先把地基打牢，后面才不容易翻车。