3. PDB文件读取与解析：用Biopython把结构“拆开”
拿到一个PDB文件后，第一件事就是把它读进内存，然后一层层拆开：模型→链→残基→原子。说白了，就像剥洋葱。我个人习惯用Biopython的PDBParser，它稳定、文档全，社区里踩过的坑基本都能搜到答案。
3.1 安装与导入
Biopython是必须装的。你在终端跑一行就行：
pip install biopython然后在你脚本里这样导入核心模块：
from Bio.PDB import PDBParser
from Bio.PDB import PDBList # 如果需要从PDB在线下载
from Bio.PDB import Polypeptide # 处理肽链的工具3.2 创建解析器 & 读取文件
Biopython的PDBParser默认会做一些结构验证。你可以选择是否让解析器“容忍”PDB文件的错误。我一般这样写：
parser = PDBParser(QUIET=True) # QUIET=True 忽略非致命警告
structure = parser.get_structure("1ubq", "1ubq.pdb")第一个参数是结构ID（随便起），第二个是文件路径。返回的structure对象就是整个结构的根。为什么要设QUIET=True？因为原始PDB里经常有奇怪的注释行、格式小偏差，你不让它安静，它会在终端刷一大堆警告，看着心烦。
QUIET=False时会抛出异常。我曾经在一个老旧晶体结构上中过招——解析到一半报错，查了半天才发现是残基编号缺了几个。所以工业级脚本里我一般先QUIET=True读进来，再手动校验。
3.3 遍历层级结构：Model → Chain → Residue → Atom
Biopython把PDB结构映射成树形对象。走一遍你就清楚了：
for model in structure:
print(f"Model ID: {model.get_id()}")
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:
print(f" Atom: {atom.get_name()} 坐标: {atom.get_coord()}")这段代码会输出类似：
Model ID: 0
Chain: A
Residue: MET 1
Atom: N 坐标: [ 1.234 2.345 3.456]
Atom: CA 坐标: [ 1.567 2.678 3.789]
Residue: GLN 2
...你想想看，这种显式的遍历特别适合Debug。每次拿到新PDB，我第一件事就是跑一遍这个循环，确认文件结构是否正常。
3.4 常用对象属性与方法
我整理了一张速查表，平时写脚本直接翻：
| 对象 | 常用方法 / 属性 | 说明 |
|---|---|---|
| Structure | get_id(), get_models() | 结构ID、模型迭代器 |
| Model | get_id(), get_chains() | 模型序号（通常0或1） |
| Chain | get_id(), get_residues() | 链字母（A/B/C…） |
| Residue | get_resname(), get_id(), get_atoms() | 残基名称（如ALA）、元组(hetflag, seq, insertion) |
| Atom | get_name(), get_coord(), get_bfactor() | 原子名（CA/CB/N…）、坐标、温度因子 |
这里要注意：residue.get_id()返回的是(hetflag, residue_number, insertion_code)，比如(' ', 123, ' ')。其中hetflag如果是W表示水分子，H_开头的是配体。我们做清洗时经常需要跳过水或配体，用这个字段过滤就可以。
Polypeptide.is_aa(residue)。它比手动比对残基名列表靠谱得多。我有个项目里要提取所有氨基酸的主链原子，直接用了这段代码：
from Bio.PDB.Polypeptide import is_aa
for residue in chain:
if is_aa(residue):
# 处理主链原子
pass3.5 快速获取原子坐标的捷径
有时候我们不需要遍历所有层级，只需要某个残基的某个原子。Biopython提供了Selection.unfold_entities或直接索引的方法。不过我更喜欢用“向量化”的方式——先把所有原子聚在一起：
from Bio.PDB import Selection
atoms = Selection.unfold_entities(structure, 'A') # A 代表原子层级
coords = [atom.get_coord() for atom in atoms]
print(f"总原子数: {len(coords)}")这样能快速拿到所有原子的坐标数组，方便后面做坐标变换或RMSD计算。
3.6 避坑指南：PDB文件里常见的“坑”
我曾经在解析一个来自PDB的晶体结构时，发现链编号竟然是小写字母"a"而不是大写"A"。Biopython默认按大写处理，结果chain.get_id()返回了"a"，我写了个比较if chain_id == 'A'死活匹配不上，浪费了20分钟。后来养成习惯：先打印所有链ID看看。
另外，有些PDB文件里包含多个模型（比如NMR结构）。这时候structure.get_models()能返回多个模型对象。遍历时别忘了指定你要哪个模型：
first_model = structure[0] # 像字典一样按索引获取
# 或者
for model in structure:
if model.get_id() == 0:
# 处理第一个模型
breakresidue.get_id()[1]返回的整数可能会引起误解。我一般用str(residue.get_id()[1]) + (residue.get_id()[2] or '')来构造一个稳定标识。
3.7 本章知识结构图
下面这张SVG图总结了从PDB文件到原子坐标的层级关系，画得挺直观：
这张图从文件到原子，四层嵌套，你写遍历代码时脑子里的结构就跟它一模一样。记住：Structure -> Model -> Chain -> Residue -> Atom，每一层都可以用get_xxx()拿下级对象。