3. 结构质量检查：残基缺失、原子缺失、异常键长键角、Ramachandran图

拿到一个PDB文件，我最习惯做的事不是直接看原子坐标，而是先跑一遍质量检查。你想想看，如果模型里残基都缺了一串，或者键长快到断裂了，后面不管是做分子动力学还是对接，都是白费功夫。我在项目里遇到过好多次，就是因为开头没做这步检查，后面跑了几天的模拟结果全是废的。所以这一章，我们就把这几个核心检查项挨个捋一遍。

说白了，结构质量检查就四件事：残基有没有缺、原子有没有丢、键长键角是不是在合理范围内、主链二面角落在不落在允许区。下面我从工程实战的角度一个个讲，每个检查我都会告诉你常见的坑，以及我自己的处理方式。

3.1 残基缺失检查

残基缺失是最容易被忽略的问题。很多PDB文件里的loop区、或是末端残基直接被删掉了，但PDB header里根本没注明。我最早做结构分析的时候，拿过一个X-ray结构，链A从1到250，链B却只有10到200——中间缺了十几号残基，而我用了好几天才发现。

检查残基缺失，核心方法是用PDB文件的SEQRES记录和实际原子记录的ATOM/HETATM序列做比对。我给一个常用的Python片段：

from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure("1abc", "1abc.pdb")
for model in structure:
    for chain in model:
        residues = list(chain.get_residues())
        res_ids = [res.get_id()[1] for res in residues]
        gaps = []
        for i in range(1, len(res_ids)):
            if res_ids[i] - res_ids[i-1] != 1:
                gaps.append((res_ids[i-1]+1, res_ids[i]-1))
        if gaps:
            print(f"Chain {chain.id} gaps: {gaps}

这里要注意：PDB文件里残基编号可能跳过某些数，比如从50直接跳到55。这不一定是缺失，可能是故意留的间隔用来放插入片段。我建议先读SEQRES里的完整序列，再跟实际残基编号比对，这样能排除假阳性。

3.2 原子缺失检查

原子缺失比残基缺失更隐蔽。有时候残基编号还在，但少了几个侧链原子。比如一个赖氨酸只有Cα、Cβ，后面四个原子全没了。这就是典型的电子密度不够导致的“部分残基”。

我检查原子缺失的惯用思路是：遍历每个残基，跟标准氨基酸原子列表比对。Biopython里的standard_aa_names字典可以帮忙。下面是我常用的检查函数：

from Bio.PDB.Polypeptide import standard_aa_names

def check_missing_atoms(structure):
    for model in structure:
        for chain in model:
            for residue in chain:
                if residue.get_resname() in standard_aa_names:
                    expected_atoms = set(standard_aa_names[residue.get_resname()])
                    actual_atoms = set(a.name for a in residue.get_atoms())
                    missing = expected_atoms - actual_atoms
                    if missing:
                        print(f"Missing atoms in {residue.get_resname()}{residue.get_id()[1]}: {missing}

不过这里有个坑：有些标准库把氧原子命名为O，而有些PDB里主链羰基氧可能叫O或O1。我建议统一用atom.name.strip()清理空格再比较。另外，非标准残基（比如配体、修饰氨基酸）不要用这个办法，否则会报一大堆假阳性。

3.3 异常键长键角检查

键长键角是结构质量的直接体现。正常的Cα-C键长约1.53Å，C-N肽键约1.33Å。如果偏差超过0.05Å，说明这个区域可能有几何冲突。我记得刚入行时用过Engh & Huber的键长键角参数表，那是做X-ray修结构的老标准，现在很多工具还用这套参数。

用Biopython算键长很简单：

from Bio.PDB.vectors import calc_bond_length

def bond_length_check(structure, max_dev=0.1):
    for model in structure:
        for chain in model:
            for residue in chain:
                atoms = list(residue.get_atoms())
                for i in range(len(atoms)-1):
                    for j in range(i+1, len(atoms)):
                        if atoms[i].is_disordered() or atoms[j].is_disordered():
                            continue
                        d = calc_bond_length(atoms[i].get_vector(), atoms[j].get_vector())
                        # 这里需要查标准键长表简单判断
                        if abs(d - 1.5) > max_dev and abs(d - 1.33) > max_dev and abs(d - 1.54) > max_dev:
                            print(f"Abnormal bond: {atoms[i].name}-{atoms[j].name} {d:.2f}Å")

这个代码很粗糙，实战中建议用ProCheck或MolProbity的几何参数。我自己的习惯是先将PDB转成mmCIF格式，然后直接用mmCIF自带的_struct_conn表里的键长数据做统计。如果无键长信息可用，就用CCDC的CSD平均键长做参照。

键角检查也是类似思路。但键角计算需要三个原子成组，而且涉及二面角。我更推荐用现成的tool——比如pdb-tools里的pdb_validate_angles。下面是我用的一个简化版：

# 用ProCheck的命令行检查
procheck -pdb my_structure.pdb -outdir procheck_output
# 然后从procheck_output/bands.ls提取z-score异常的残基

3.4 Ramachandran图分析

Ramachandran图是结构质量的“照妖镜”。它展示的是主链二面角φ和ψ的分布。一个高质量的晶体结构，90%以上的残基应该落在“最允许区”，而完全不允许区的残基占比应小于2%。

实战中最简单的方式是用Rampage（验证工具，可网页运行）或者用Biopython画一个。我自己写过一个脚本，调用rampage命令行：

rampage -f my_structure.pdb -o rampage_output
# 读取输出文件中的favoured、allowed、outlier统计

如果你非要自己算二面角，Biopython有现成接口：

from Bio.PDB.vectors import calc_dihedral
from Bio.PDB.PDBParser import PDBParser

structure = PDBParser(QUIET=True).get_structure("test", "test.pdb")
for model in structure:
    for chain in model:
        residues = list(chain.get_residues())
        for i in range(1, len(residues)-1):
            try:
                n = residues[i-1]["N"].get_vector()
                ca = residues[i]["CA"].get_vector()
                c = residues[i]["C"].get_vector()
                n_next = residues[i+1]["N"].get_vector()
                phi = calc_dihedral(c, n, ca, c)
                psi = calc_dihedral(n, ca, c, n_next)
                print(f"{residues[i].get_resname()}{residues[i].get_id()[1]} phi={phi:.2f} psi={psi:.2f})
            except KeyError:
                continue

注意一点：Phi和Psi的合理范围会受氨基酸类型影响。比如甘氨酸的允许区更大，而脯氨酸的φ角几乎锁死在-60°左右。我见过有人把甘氨酸超出的点也标记为outlier，这就不对了。所以一定要按氨基酸类型做分类统计。

下面这张SVG图可以帮你快速理解质量检查的完整流程——从PDB文件读入，到输出检查报告：

3.5 实战避坑与自动化

一次完整的质量检查，我一般在预处理流水线里写成脚本，每次出新结构都自动跑一遍。这里我给出一个简易的Python脚本流程：

# quality_check.py  – 配合Biopython和外部工具
import sys, subprocess, tempfile, os
from Bio.PDB import PDBParser

pdb_file = sys.argv[1]
parser = PDBParser(QUIET=True)
structure = parser.get_structure("target", pdb_file)

# 1. 残基缺失检查
gaps = []
for chain in structure[0]:
    res_ids = [res.get_id()[1] for res in chain if res.get_id()[0] == ' ']
    for i in range(1, len(res_ids)):
        if res_ids[i] - res_ids[i-1] > 1:
            gaps.append((chain.id, res_ids[i-1]+1, res_ids[i]-1))

# 2. 原子缺失检查（略，见前文）

# 3. 键长检查：调用ProCheck
subprocess.run(["procheck", "-pdb", pdb_file, "-outdir", "procheck_out"])

# 4. Ramachandran检查：调用Rampage
subprocess.run(["rampage", "-f", pdb_file, "-o", "rampage_out"])

# 5. 输出汇总
with open("qa_report.txt", "w") as f:
    f.write(f"Gaps found: {len(gaps)}\n")
    for g in gaps:
        f.write(f"Chain {g[0]}: residues {g[1]}-{g[2]} missing\n")
    # 读取procheck/rampage的报告追加
print("Done. See qa_report.txt")

最后提醒一句：所有检查结果都要写入审计跟踪日志。我习惯用JSON格式记录，包含检查时间、工具版本、阈值、每个异常的详细信息。这样几周后反查时，能快速知道当初是怎么处理的。这点我们在后面第8章会详细展开。