3、分子坐标文件格式：PDB、XYZ、MOL2等常见格式的结构与读取方式

做RMSD计算，第一步就是搞定分子坐标文件。这些文件格式看着多，其实核心就三样东西：原子符号（或元素）、三维坐标、以及可能的连接信息（键）。我入行那会儿被各种格式搞晕过，后来发现只要抓住每类格式的“骨架”，读起来就顺手了。

今天咱们就把PDB、XYZ、MOL2这三个最常见的老朋友捋一遍。你想想看，市面上90%的分子模拟工具都在用它们。我会把每种格式的结构拆开，再给一段Python读取代码，这样你拿到任意一个文件都能快速上手。

3.1 PDB格式——生物大分子界的“普通话”

PDB（Protein Data Bank）格式是蛋白质、核酸结构的标配。它有一定的历史包袱，但胜在标准化。每个PDB文件由若干“记录类型”组成，比如ATOM、HETATM、TER、END。我最常用的是ATOM记录，它里面藏着每个原子的坐标、残基名称、链标识等。

ATOM记录是固定宽度的文本，每行80个字符。字段位置是硬编码的：

起始列结束列字段名描述
16记录类型ATOM 或 HETATM
711原子序号整数，从1开始
1316原子名称如 CA、CB、N
1820残基名称如 ALA、GLY
2222链标识单字符，如 A
2326残基序号整数
3138x坐标浮点数，单位埃
3946y坐标浮点数
4754z坐标浮点数
5560占有率默认1.0
6166温度因子B值
7778元素符号右对齐
我的习惯：写读取函数时，不要用字符串切片去硬解析每列，而是用行长度判断。有些PDB文件的行末尾有多余空格，但核心列位置不变。我一般先检查行是否以 "ATOM" 或 "HETATM" 开头，再取子串。

下面是一个从头解析ATOM坐标的Python函数。我故意没用第三方库，方便你理解原理：

def parse_pdb_atoms(pdb_path):
    atoms = []
    with open(pdb_path, 'r') as f:
        for line in f:
            if line.startswith(('ATOM', 'HETATM')):
                # 固定列解析
                x = float(line[30:38].strip())
                y = float(line[38:46].strip())
                z = float(line[46:54].strip())
                atom_name = line[12:16].strip()
                res_name  = line[17:20].strip()
                chain     = line[21].strip()
                res_num   = int(line[22:26].strip())
                element   = line[76:78].strip()
                atoms.append({
                    'x': x, 'y': y, 'z': z,
                    'atom_name': atom_name,
                    'res_name': res_name,
                    'chain': chain,
                    'res_num': res_num,
                    'element': element
                })
    return atoms
曾经踩过的坑：有的PDB文件里，元素符号在76-78列，但原子名称在13-16列。如果你直接用元素符号去查原子半径或质量，可能会捡到“CA”（碳）或者“CH”（碳氢？不存在）。我建议优先用元素符号字段，如果为空再fallback到原子名称的第一个字母。

3.2 XYZ格式——极简主义，一把梭

XYZ格式简单到令人感动。第一行是原子数，第二行是注释（可以随便写，一般写步数或能量），后面每行一个原子：元素符号 + x y z。没了，真的没了。它不存储成键信息，也不区分残基、链之类的。适合小分子、团簇、未成键体系。

举个例子：

5
water cluster snapshot
O     0.0000   0.0000   0.1179
H     0.7572   0.0000  -0.4716
H    -0.7572   0.0000  -0.4716
O     2.5000   0.0000   0.1179
H     3.2572   0.0000  -0.4716

读取XYZ简直就像在念小学生作文：

def parse_xyz(xyz_path):
    with open(xyz_path, 'r') as f:
        n_atoms = int(f.readline().strip())
        comment = f.readline().strip()
        atoms = []
        for i in range(n_atoms):
            parts = f.readline().split()
            if len(parts) < 4:
                continue
            element = parts[0]
            x, y, z = float(parts[1]), float(parts[2]), float(parts[3])
            atoms.append({
                'element': element,
                'x': x, 'y': y, 'z': z
            })
    return atoms, comment

你看，核心逻辑不超过10行。但正因为简单，它缺少很多元信息。如果你想做蛋白质受体的RMSD，XYZ格式会丢掉残基和链的归属，这时候就得用PDB。所以，我把XYZ称为“小而美的裸坐标格式”。

3.3 MOL2格式——键、电荷、原子类型都给你

MOL2是Tripos公司推出的格式，在药物设计领域用得很多。它的特点是分块（section）：比如@<TRIPOS>MOLECULE、@<TRIPOS>ATOM、@<TRIPOS>BOND、@<TRIPOS>SUBSTRUCTURE等等。每个块处理不同的信息。

ATOM块里每行格式：原子ID + 原子名称 + x + y + z + 原子类型 + 子结构编号 + 子结构名称 + 电荷。注意原子类型不是元素符号，而是专门定义的SYBYL类型（比如C.3、C.2、N.3、O.3等）。我当初因为这个吃过亏——直接用原子类型去匹配元素周期表，发现C.3和C.2都是碳，但类型不同。

BOND块则记录成键关系：键ID + 起始原子ID + 终点原子ID + 键型。键型有1（单键）、2（双键）、3（三键）、ar（芳香键）、am（酰胺键）等。

要读取MOL2，我推荐按"@<TRIPOS>"标记分割段落：

def parse_mol2(mol2_path):
    atoms = []
    bonds = []
    current_section = None
    with open(mol2_path, 'r') as f:
        for line in f:
            line = line.strip()
            if line.startswith('@<TRIPOS>'):
                current_section = line.split('>')[1].strip()  # 截取块名
                continue
            if not line:
                continue
            if current_section == 'ATOM':
                parts = line.split()
                if len(parts) >= 9:
                    atom_id = int(parts[0])
                    atom_name = parts[1]
                    x, y, z = float(parts[2]), float(parts[3]), float(parts[4])
                    atom_type = parts[5]
                    subst_id = int(parts[6])
                    subst_name = parts[7]
                    charge = float(parts[8]) if len(parts) >= 9 else 0.0
                    atoms.append({
                        'id': atom_id,
                        'name': atom_name,
                        'x': x, 'y': y, 'z': z,
                        'type': atom_type,
                        'subst_id': subst_id,
                        'subst_name': subst_name,
                        'charge': charge
                    })
            elif current_section == 'BOND':
                parts = line.split()
                if len(parts) >= 4:
                    bond_id = int(parts[0])
                    atom1 = int(parts[1])
                    atom2 = int(parts[2])
                    bond_type = parts[3]
                    bonds.append({
                        'id': bond_id,
                        'atom1': atom1,
                        'atom2': atom2,
                        'type': bond_type
                    })
    return atoms, bonds

MOL2的一大优点是自带键连接信息，这在RMSD计算中很少用到（因为RMSD只看坐标），但在做结构比对、药效团匹配时，键信息能让拓扑对齐更准确。

3.4 三种格式核心差异对比

特性PDBXYZMOL2
原子数量ATOM记录行数第一行显式给出MOLECULE块定义
坐标精度固定宽度，三位小数自由格式，可高精度自由格式，通常六位小数
成键信息无直接存储（需CONECT）无有BOND块
残基/链信息包含无通过子结构间接表示
电荷无全局电荷（CRYST1可选）无原子级电荷
读取难度中（固定列，易错位）低中（分段解析）
常见应用蛋白质、核酸小分子、轨迹帧配体、药物分子
一句话选择建议：做蛋白质RMSD用PDB；算小分子或团簇的RMSD用XYZ；既要坐标又要键和电荷信息用MOL2。我个人会在项目里同时准备PDB和XYZ两套格式，PDB给下游可视化，XYZ给快速数值计算。

3.5 知识体系一览

下面这张图概括了三种格式在RMSD计算中的定位和读取要点：

分子坐标格式核心结构 PDB格式 XYZ格式 MOL2格式 记录类型 ATOM / HETATM 固定宽度 80列 残基链信息 极简结构 行1: 原子数 行2: 注释 坐标行: 元素 x y z 分段结构 @ATOM @BOND SYBYL原子类型 读取关键：行解析 + 字段提取 蛋白质受体 RMSD → PDB 小分子/团簇 RMSD → XYZ 药物分子+键拓扑 → MOL2

3.6 避坑总结

  • PDB列对齐问题：有的软件生成PDB时会用制表符或全角空格，导致解析错位。我习惯先用len(line)检查行长度，低于78的自动补空格。
  • XYZ负数坐标：元素符号可能在一个字符串内，比如“C -0.5 1.2 -3.0”，split后还是能正常工作。但小心元素为“Cu”时跟在列首不会被误读。
  • MOL2的@大小写：标准的WRITER生成的是大写的TRIPOS，但有些软件可能写小写。建议用line.upper().startswith('@')来统一。
  • 缺失电荷：MOL2的电荷字段在早期版本中可选。如果解析时发现这一列缺失（行长度不够），我会默认电荷为0.0。

嗯，到这里三种格式的骨骼你已经摸清了。在接下来的课程中，我们会直接用这些读取函数把坐标喂给RMSD计算核心。放心，代码我会保持和今天一样——朴实、无冗余、直接可用。

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