第4章 解析PDB文件内容——把原子坐标读进内存

说到解析PDB文件，我觉得这是整个质控体系的基础。你要是连文件内容都读不准，后面什么验证、修复、统计都是空中楼阁。

我刚开始做结构解析那阵子，曾试着用一行Python脚本来读PDB。结果发现，这些看起来规规矩矩的文本文件，实际解析起来充满陷阱。

今天咱们就深入看看，怎么正确地把PDB文件里那些数据扒出来。

4.1 PDB格式的骨架

PDB文件本质上是个固定宽度的文本格式。每一行80个字符，每个字段的位置都是定死的。说白了，就是个纯文本的数据库表。

我个人习惯把PDB的记录类型分成三组：

记录类别典型记录作用
头信息HEADER, TITLE, KEYWDS, EXPDTA描述实验信息、来源、作者等元数据
坐标核心ATOM, HETATM, TER, END原子坐标、残基链、B因子等结构数据
附属记录CRYST1, REMARK, SSBOND, SITE晶胞、注释、化学键、活性位点等额外信息

你在解析的时候，必须清楚每种记录的字段偏移。少一个空格，读出来的数据就全乱了。

经验之谈：
处理PDB文件之前，先确认它的版本。PDB v3.0开始对字段格式做了调整，特别是REMARK记录。我习惯在解析前用grep "^HEADER"看一眼版本号，再决定用哪个解析模板。

4.2 坐标核心部分——占80%的工作量

ATOM记录是整个解析的重头戏。每行都代表一个原子，字段排列是这样的：

记录类型:  ATOM
序列号:    1-6    (右对齐)
原子名称:  13-16  (左对齐)
残基名称:  18-20  (左对齐)
链ID:      22     (单字符)
残基号:    23-26  (右对齐)
插入码:    27     (可选项)
X坐标:     31-38  (右对齐, 三位小数)
Y坐标:     39-46
Z坐标:     47-54
占有率:    55-60  (一位小数)
B因子:     61-66  (两位小数)
元素:      77-78  (右对齐)
电荷:      79-80  (可选项)

你想想看，这就意味着你不能用简单的空格分割来解析。如果你的代码是line.split()，那十有八九会读错那些带空格的原子名。我在项目中遇到过好几次，" CA "和"CA "被当成同一个东西，导致后续的残基匹配全崩了。

4.3 关键附属记录

除了坐标，有几个记录对质控特别重要。我挑几个重点说说。

4.3.1 CRYST1——晶胞参数

这条记录告诉你晶体的对称性和晶胞尺寸。很多验证和重计算都得靠它。

CRYST1   41.810   41.810   38.000  90.00  90.00 120.00 P 31 2 1    6

字段解析：

  • a, b, c：晶胞边长，6-15列，右对齐
  • α, β, γ：晶胞角度，16-25列
  • 空间群：26-40列，例：P 31 2 1
注意坑：
不同版本的PDB文件在CRYST1记录里，空间群的格式并不统一。有些写P 31 2 1，有些写P3121。我建议在加载完成后，统一用标准空间群符号映射表做一次标准化。

4.3.2 REMARK——自由文本的沼泽

REMARK记录是给人类看的，不是给机器看的。但有些质控指标会藏在里面，比如分辨率和R因子。

我见过最麻烦的解析任务是提取REMARK 200里分辨率的信息。它经常写成这样的格式：

REMARK   2 RESOLUTION. 2.50 ANGSTROMS.

你能直接取数字吗？不行。有些文件写成2.5，有些写2.50，甚至前面还会有杂乱的注释。我常用的办法是用正则表达式去抓(\d+\.?\d*)\s*ANGSTROMS，然后做容错转换。

4.4 解析流程——别自己造轮子

核心观点：
我强烈建议不要从零写PDB解析器。直接复用BioPython的PDB解析模块，再封装自己的质控逻辑。这样稳定，还省事。

不过，自己写一套轻量级解析方案，对理解格式本身帮助很大。下面是一段我自己项目里用过的简易解析代码：

def parse_atom_record(line):
    """从固定宽度的ATOM记录中提取字段"""
    if not line.startswith(('ATOM', 'HETATM')):
        return None
    
    record = {}
    record['type'] = line[0:6].strip()
    record['serial'] = int(line[6:11].strip())
    record['name'] = line[12:16].strip()
    record['resname'] = line[17:20].strip()
    record['chain'] = line[21].strip()
    record['resnum'] = int(line[22:26].strip())
    record['x'] = float(line[30:38].strip())
    record['y'] = float(line[38:46].strip())
    record['z'] = float(line[46:54].strip())
    record['occupancy'] = float(line[54:60].strip())
    record['bfactor'] = float(line[60:66].strip())
    return record

这段代码能处理大部分标准的PDB文件。但我曾经遇到过一批老旧文件，字段不到80列，或者有异常的填充字符。这时候你就要加点容错逻辑，比如用.strip()配合if len(line) > 30的检查。

4.5 常见坑和实战建议

  • 空行和注释行：解析时记得过滤掉空行以及以ANISOU开头的各向异性温度因子记录，除非你关注这些数据。
  • 链ID缺失：有些旧PDB文件链ID列是空格。我建议在加载时自动赋'A'作为默认链。
  • 残基编号重复：同一链里出现两个1号残基的情况，在NMR模型里很常见。这时候得用insertion_code区分，或者启用model_index。
  • B因子异常：我见过B因子是负数或者超过100的情况。解析时标记这类异常，但不要直接抛弃，留到后续质控再判断。
避坑指南：
我曾经负责一个批次的PDB文件清洗任务，300个文件里有23个的ATOM记录不是80列。原因是一些旧软件在写入时用了混合长度。解法很简单：在解析前先对每行做line.ljust(80)补齐。

好，你知道了怎么解析PDB文件，还了解了常见坑。下一步就是把这些解析出来的数据，套进质控规则里做自动化检测。不过那是后面章节的事了，今天先到这里。

← 上一章 📖 返回目录 下一章 →