3、PDB文件解析：使用Biopython读取PDB文件的基本方法

PDB文件，说白了就是蛋白质的「身份证」。里面记录了每个原子的坐标、残基信息、链标识等等。我刚开始接触生物信息学那会儿，面对动辄几万行的PDB文件，说实话，头皮发麻。但你用Biopython去读，一切就简单多了。

我个人习惯，在开始任何结构分析之前，第一步永远是「读进来再说」。Biopython的Bio.PDB模块就是干这个的。

3.1 核心模块导入与基础读取

Biopython提供了两种读取方式：PDBParser和MMCIFParser。PDB格式大家最熟悉，mmCIF格式则是PDB现在的官方标准格式。大部分公开数据已经迁移到mmCIF了，但老数据还是PDB居多。我自己的做法是：优先检测扩展名，然后选择对应的解析器。

先看看最基础的导入和读取：

from Bio.PDB import PDBParser
from Bio.PDB import PDBList
# 创建一个解析器
parser = PDBParser(QUIET=True)
# 读取本地PDB文件
structure = parser.get_structure("1ubq", "1ubq.pdb")

这里有个小细节：QUIET=True参数是我强烈建议加上的。不加的话，遇到文件里有些格式不太规范的记录，解析器会疯狂输出警告信息。我在项目中遇到过，一次处理几百个结构的时候，那个警告日志刷得我直接怀疑人生。加了这个参数，只有真正的错误才会被打印。

3.2 结构层次模型：从Structure到Atom

Biopython把PDB文件组织成了一个四层级的对象模型。你想想看，这就像文件系统的目录结构一样：

• Structure —— 最顶层，就像根目录
• Model —— 第二层，通常只有一个（NMR可能有多个）
• Chain —— 第三层，蛋白质链（A、B、C...）
• Residue —— 第四层，氨基酸残基
• Atom —— 最底层，单个原子

为什么会这样设计？因为实际工作中，你可能需要精确控制到任意一个层次。比如我要校验批次一致性，常常需要拿到两个不同批次解析结果里的Chain对象，然后逐残基去比较。这种层次结构让我能精准定位。

# 从structure开始遍历
structure = parser.get_structure("1ubq", "1ubq.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:
            hetflag, resseq, icode = residue.get_id()
            print(f"    Residue: {residue.get_resname()} {resseq}")
            for atom in residue:
                print(f"      Atom: {atom.get_name()} -> {atom.get_vector()}")

这段代码看着简单，但实际跑起来输出量非常大。对于1ubq这种小蛋白（76个残基），完整遍历一次的Atom数量也在600个左右。

3.3 常用属性与坐标获取

拿到Atom之后，最常做的事就是获取它的三维坐标。Biopython用Atom.get_vector()返回一个Vector对象。我自己比较喜欢用它来进行空间距离计算：

from Bio.PDB import Vector

# 获取两个原子的坐标
atom1 = residue["CA"]  # CA是Cα原子
atom2 = residue["CB"]  # CB是Cβ原子

coord1 = atom1.get_vector()   # Vector(x, y, z)
coord2 = atom2.get_vector()

# 计算两个原子之间的距离
dist = (coord1 - coord2).norm()
print(f"CA-CB distance: {dist:.3f} Å")

除了坐标，每个Atom还有很多关键属性——温度因子、占据率、元素类型等。在批次一致性校验中，温度因子和占据率经常用来判断结构质量是否稳定。

3.4 二级结构与残基类型判断

PDB文件本身通常不包含二级结构信息。但Biopython提供了一些好用的工具来判断残基类型：

from Bio.PDB.vectors import calc_dihedral
from Bio.PDB.Polypeptide import PPBuilder

# 提取一条链中的所有标准氨基酸
ppb = PPBuilder()
for pp in ppb.build_peptides(structure[0]["A"]):
    print(f"Peptide length: {len(pp)}")
    for residue in pp:
        print(f"  {residue.get_resname()} {residue.get_id()[1]}")

这个PPBuilder帮了大忙。它会把非标准残基（比如HETATM记录的配体）自动过滤掉，只留下标准氨基酸。我记得在做一个多结构比对的项目时，如果不先用这个把水分子和配体过滤掉，后面的比对算法会直接爆炸。

# 判断残基是否属于给定的氨基酸列表
common_aa = {"ALA", "ARG", "ASN", "ASP", "CYS", ...}
if residue.get_resname() in common_aa:
    print("Standard residue")
else:
    print("Non-standard (maybe ligand or water)")

3.5 完整的读写流程

到这里，我们已经能把一个PDB文件读进来，提取出所有原子坐标、名字、残基信息。但如果有多个批次的结构文件，要怎么做一致性校验？

我的套路是这样：

1. 用PDBParser分别读取每个批次的PDB文件
2. 取第一个Model（通常就是唯一的一个）
3. 遍历Chain，对每个残基的Cα坐标做相似度比较
4. 记录下偏差超过阈值的残基位置

说白了，这跟比对两个字符串的差异没什么本质不同。但关键在于，你得先把PDB文件里藏着的「元数据」都扒出来。

# 一个简单的批次比较骨架
def compare_structures(pdb_a, pdb_b, tolerance=1.0):
    parser = PDBParser(QUIET=True)
    struct_a = parser.get_structure("batch1", pdb_a)
    struct_b = parser.get_structure("batch2", pdb_b)
    
    chain_a = struct_a[0]["A"]
    chain_b = struct_b[0]["A"]
    
    for res_a, res_b in zip(chain_a, chain_b):
        if res_a.get_resname() != res_b.get_resname():
            print(f"Residue mismatch at {res_a.get_id()}")
            continue
        # 比较Cα位置
        ca_a = res_a["CA"].get_vector()
        ca_b = res_b["CA"].get_vector()
        diff = (ca_a - ca_b).norm()
        if diff > tolerance:
            print(f"Deviated: {res_a.get_id()} diff={diff:.3f}Å")

这个函数虽然简陋，但已经能抓到很多问题。我自己就曾经用它发现过一个批次中某个残基的侧链方向整个反了180度——这种情况如果不校验，直接拿去跑分子动力学模拟，结果基本就是废的。

这张图把这五个层级的关系画得很清楚了。从Structure一路向下，你可以通过.get_list()或直接下标索引快速地访问任意一层。这种层次结构，本质上就是为逐级精细校验而生的。