第2章 分子文件格式与读取：PDB、mmCIF、SDF
做结构比对，第一步就是把手头的分子文件读进内存。这步看似简单，但坑特别多。我在项目中遇到过，不同格式的数据混在一起，读取逻辑不统一，最后比对结果全乱套。
本章我们重点啃三个主流格式：PDB、mmCIF 和 SDF。我们还要用 Python 的 Bio.PDB 库，把原子坐标和链信息老老实实抽出来。
2.1 三大格式的来龙去脉
先搞清楚它们是什么，适合干什么。
| 格式 | 全称 | 典型用途 | 特点 |
|---|---|---|---|
| PDB | Protein Data Bank | 蛋白质/核酸 X-ray 结构 | 固定列宽、行号区分数据 |
| mmCIF | Macromolecular Crystallographic Information File | PDB 官方新标准 | 键-值结构、支持大结构 |
| SDF | Structure Data File | 小分子/配体/化合物库 | 多记录、包含属性字段 |
我入行时读的第一批数据就是 PDB。这东西历史悠久，但格式确实有点“古板”。每一行都按列对齐，从 1 数到 80，少一个字符都不行。后来 PDB 官方自己都推不动了，mmCIF 就接班了。
2.2 PDB 格式——老兵不死
PDB 文件是纯文本，每行开头有个记录名，比如 ATOM、HETATM、TER、END。最重要的就是 ATOM 记录，里面藏着每个原子的三维坐标。
我们来看一条典型的 ATOM 行：
ATOM 1 N ALA A 1 7.354 0.509 18.683 1.00 0.00 N
这里面每列代表什么？
- 第 1-6 列：记录名（ATOM）
- 第 7-11 列：原子序号（1）
- 第 13-16 列：原子名（N）
- 第 18-20 列：氨基酸名（ALA）
- 第 22 列：链标识（A）
- 第 23-26 列：残基序号（1）
- 第 31-38 列：X 坐标（7.354）
- 第 39-46 列：Y 坐标（0.509）
- 第 47-54 列：Z 坐标（18.683）
2.3 mmCIF 格式——新王登基
mmCIF 的读法和 PDB 完全不同。它是键-值对结构，很像一个带层级关系的配置文件。
loop_
_atom_site.group_PDB
_atom_site.id
_atom_site.type_symbol
_atom_site.Cartn_x
_atom_site.Cartn_y
_atom_site.Cartn_z
ATOM 1 N 23.456 45.678 12.345
ATOM 2 CA 24.567 46.789 13.456
它用 loop_ 开头定义表格，每个字段名都带 _atom_site. 前缀。说到底就是更规范了，也塞得下更大的蛋白结构。
2.4 SDF 格式——小分子专属
SDF 是 MDL 公司定的格式，一个文件可以装多个分子，每个分子之间用 $$$$ 分隔。
一个记录的例子：
LIGAND_NAME
SDfile
18 17 0 0 0 0 0 0 0 0999 V2000
1.2345 2.3456 3.4567 C 0 0 0 0 0 0 0 0 0 0 0 0
2.3456 3.4567 4.5678 N 0 0 0 0 0 0 0 0 0 0 0 0
...
M END
> <ID>
LIG-001
> <SMILES>
CC(=O)OC1=CC=CC=C1C(=O)O
$$$$
第 4 行是原子计数和坐标块；M END 表示坐标结束；后面 > 开头的部分是属性字段，可以放 ID、SMILES、活性值等。
我在做虚拟筛选项目时，从 PubChem 下了几百万个 SDF，就是靠属性字段过滤活性的。这招非常实用。
2.5 用 Python 读取 PDB 文件
我推荐用 Bio.PDB 库。它由生信大神 Thomas Hamelryck 开发，能优雅地处理 PDB、mmCIF 和 mmTF 格式。
安装 Bio.PDB
pip install biopython
读取一个 PDB 文件
from Bio.PDB import PDBParser
# 创建解析器
parser = PDBParser(QUIET=True) # QUIET=True 屏蔽警告
# 解析文件，返回 Structure 对象
structure = parser.get_structure("1ubq", "1ubq.pdb")
QUIET=True 这个参数很关键。一开始我没加，结果跑脚本时满屏的格式警告信息，看得头皮发麻。
遍历模型、链、残基、原子
Bio.PDB 的对象层级是：Structure → Model → Chain → Residue → Atom
for model in structure:
for chain in model:
print(f"Chain ID: {chain.get_id()}")
for residue in chain:
# 跳过水分子
if residue.get_id()[0] == " ":
for atom in residue:
print(f" Atom: {atom.get_name()} | Coord: {atom.get_coord()}")
这样就能拿到每个原子的名字和三维坐标，后面做构象比对就直接用这些数据。
提取链级别的序列和结构
for model in structure:
for chain in model:
chain_id = chain.get_id()
chain_seq = ""
chain_coords = []
for residue in chain:
if residue.get_id()[0] == " ":
chain_seq += residue.get_resname()
for atom in residue:
if atom.get_name() == "CA": # 只取 alpha 碳
chain_coords.append(atom.get_coord())
print(f"Chain {chain_id}: {len(chain_coords)} CA atoms")
这里只取了 CA 原子，做构象比对时，CA 是最常用的主干原子。
2.6 读取 mmCIF 和 SDF
mmCIF 读取
from Bio.PDB.MMCIFParser import MMCIFParser
parser = MMCIFParser(QUIET=True)
structure = parser.get_structure("1ubq", "1ubq.cif")
后面的遍历逻辑和 PDB 完全一样，因为 MMCIFParser 也会返回 Structure 对象。这对我们用户来说最友好，接口统一。
SDF 读取（用 RDKit）
from rdkit import Chem
supplier = Chem.SDMolSupplier("ligands.sdf")
for mol in supplier:
if mol is None:
continue
conf = mol.GetConformer()
for i in range(mol.GetNumAtoms()):
pos = conf.GetAtomPosition(i)
print(f"Atom {i}: {pos.x} {pos.y} {pos.z}")
RDKit 是化学信息学的瑞士军刀。不过要注意，SDF 中多记录之间用 $$$$ 分隔，Chem.SDMolSupplier 已经帮我们处理好了。
list() 直接全部加载，内存会爆。我习惯用生成器逐条处理。
2.7 本章知识框架
下面这张图是整个流程的简化版本，帮你梳理 PDB→Bio.PDB→原子坐标的路径：
流程很清楚：三种格式的文件经过各自的解析器，最终统一输出为原子坐标和链信息。这个结构在后面做构象比对时，会反复用到。