📚 课程目录第3章链与残基的提取与筛选

链与残基的提取与筛选：实现亚基分离

PDB 文件里装的可不止一个蛋白。很多时候，一个文件里塞了好几套链、配体、水分子，甚至还有核酸。你想想看，要是直接把整个结构扔进模型，出来的结果能对吗？

我个人习惯是，拿到 PDB 的第一件事——拆。把目标亚基单独拎出来，把水分子和杂原子清理掉，再按残基序号排好序。这一步做干净了，后面的特征对齐才能省心。

为什么必须做亚基分离？

我遇到过不少新手，直接从 PDB 里读所有原子，然后算 RMSD。结果呢？配体的坐标把整体偏差拉得老高，一个本该对齐得很好的结构，硬生生变成了废数据。

说白了，链与残基的筛选就是这个目的：

去掉干扰项——水、配体、离子、非标准残基
锁定目标亚基——只保留你要分析的那条链
规范残基命名——不同来源的 PDB 可能用不同命名，得统一
按序号对齐——残基编号偏移了，特征就全对不上

一句话总结：链和残基的筛选，决定了你后续特征对齐的基准是否正确。这一步省事，后面全是坑。

链 ID 的提取逻辑

PDB 里每条链都有一个字母或数字标识，写在 ATOM 记录的第五列（标准 PDB 格式）。比如链 A、链 B、链 C……但在 mmCIF 格式里，它是 auth_asym_id 字段。

我一般先用程序把链 ID 全扫一遍，看看文件里到底有哪些链。举个 Python 例子：

from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure("1abc", "1abc.pdb")

chains = {}
for model in structure:
    for chain in model:
        chain_id = chain.get_id()
        residues = [r.get_resname() for r in chain.get_residues()]
        chains[chain_id] = set(residues)

for cid, res_set in chains.items():
    print(f"链 {cid}: {len(res_set)} 种残基")

跑完这段，你就能看到每条链上都有什么残基。嗯，这时候你会发现，有些链其实只有几个水分子，或者是一小段配体。这些都可以直接去掉。

按残基名过滤

残基名是 PDB 里最常用的过滤条件。标准氨基酸是 20 种，核酸是 DA、DC、DG、DT、A、C、G、U，剩下的基本就是配体或水。

我在项目中常用的过滤规则是这样：

类别	残基名特征	处理方式
标准氨基酸	ALA, ARG, ASN, ASP, CYS, GLN, GLU, GLY, HIS, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP, TYR, VAL	保留
核酸	DA, DC, DG, DT, A, C, G, U	按需保留
水分子	HOH, WAT	丢弃
配体/离子	ATP, NAD, MG, ZN, CL …	丢弃或单独处理
非标准残基	MSE, SEP, TPO …	映射为标准残基或丢弃

小技巧：非标准残基 MSE（硒代甲硫氨酸）其实是 MET 的修饰版，在结构比对时可以直接映射为 MET。我曾经因为这个没处理，导致特征矩阵里多了一列奇怪的化学属性，排查了半天。

standard_aa = {
    'ALA','ARG','ASN','ASP','CYS',
    'GLN','GLU','GLY','HIS','ILE',
    'LEU','LYS','MET','PHE','PRO',
    'SER','THR','TRP','TYR','VAL'
}

def filter_standard_residues(chain):
    keep = []
    for residue in chain.get_residues():
        resname = residue.get_resname().strip()
        if resname in standard_aa:
            keep.append(residue)
    return keep

按残基序号范围筛选

同一个 PDB 里，残基序号可能不是从 1 开始的，也可能中间有跳跃。比如晶体结构里，有些 loop 区域电子密度不够，残基被跳过了。你直接按序号做特征对齐，肯定对不上。

我做过一个项目，两个 PDB 结构看起来是同一个蛋白，但一个的残基序号是从 2 开始的，另一个是从 1 开始。嗯，序号偏移了 1。这个问题如果不处理，后面所有特征位置全错位。

def extract_by_range(chain, start, end):
    selected = []
    for residue in chain.get_residues():
        seq_id = residue.get_id()[1]  # 残基序号
        if start <= seq_id <= end:
            selected.append(residue)
    return selected

序号范围筛选在亚基分离时特别有用。比如你要提取某个蛋白的催化域（残基 1-250），而调控域（残基 251-500）不需要，直接切掉就行。

注意：有些 PDB 文件里同一条链的残基序号不是严格递增的，或者有插入残基（如 1A, 1B）。处理时要把插入码一并考虑，不然会漏掉或错位。

亚基分离的完整流程

好了，把上面几个方法串起来，就是一条完整的亚基分离管线。我通常按这个步骤走：

读取 PDB——用 Biopython 或 Gemmi 加载结构
列出所有链——看有哪些链 ID 和对应的残基组成
筛选目标链——按链 ID 保留需要的亚基
过滤残基类型——只留标准氨基酸或核酸
按序号排序——确保残基顺序正确
序号归一化——如果有多结构比对，把序号统一偏移
输出清洗后的结构——保存为 PDB 或 mmCIF

下面是这个流程的示意图：

实战示例：从 6LU7 中提取主蛋白酶

拿 COVID-19 的主蛋白酶（PDB ID: 6LU7）来说吧。这个结构里其实有三条链：A（蛋白）、B（蛋白）、C（配体抑制剂）。我们要提取链 A 的蛋白部分。

from Bio.PDB import PDBParser, PDBIO

parser = PDBParser(QUIET=True)
structure = parser.get_structure("6lu7", "6lu7.pdb")

# 只保留链 A
chain_a = structure[0]['A']

# 过滤残基：只保留标准氨基酸
standard_aa = {
    'ALA','ARG','ASN','ASP','CYS',
    'GLN','GLU','GLY','HIS','ILE',
    'LEU','LYS','MET','PHE','PRO',
    'SER','THR','TRP','TYR','VAL'
}

class ResidueSelect:
    def accept_residue(self, residue):
        return residue.get_resname().strip() in standard_aa

# 写入新 PDB
io = PDBIO()
io.set_structure(chain_a)
io.save("6lu7_chainA_clean.pdb", ResidueSelect())

这段代码跑完，你得到的 PDB 就只有链 A 的标准氨基酸了。配体链 C 被自动丢弃，水分子也被过滤掉。

个人习惯：我一般在清洗后还会加一步——检查残基序号的连续性。如果发现跳跃超过 5，我会单独标记出来，免得后面做序列对齐时莫名其妙出错。

避坑指南

嗯，这里我说几个我踩过的坑：

链 ID 大小写敏感——'A' 和 'a' 是两条不同的链，千万别弄混
残基名有空格——PDB 格式里残基名右对齐，比如 ' ALA'，记得 strip()
插入码（insertion code）——有些残基序号带字母后缀，比如 '101A'，处理序号时要把插入码单独提取
ALT 定位（ alternate location）——同一个残基可能有多个坐标，选 occupancy 最高的那个
非标准残基的映射表——MSE→MET, SEP→SER, TPO→THR，这个映射表最好提前准备好

曾经有一次，我在处理一个膜蛋白结构时，没注意到链 ID 'A' 和 '1' 同时存在。结果把链 '1' 当成了链 'A' 的一部分，特征矩阵直接串位。排查了整整一天才发现是链 ID 混淆了。所以——先扫一遍所有链 ID，确认了再动手。

小结一下

链与残基的提取与筛选，说白了就是三个动作：选链、留残基、定序号。这三个动作做扎实了，亚基分离就水到渠成。后面你要做特征对齐、计算 RMSD、构建特征矩阵，才有可靠的基础。

记住一句话：清洗数据花的时间，会在后续分析里十倍地省回来。我个人深有体会。