2. 分子文件与数据清洗实战

在真实项目中,我们拿到的分子数据往往来自不同渠道:商业数据库导出的 SDF、对接软件生成的 MOL2、晶体结构解析得到的 PDB,甚至还有从文献补充材料里手工整理的 SMILES 列表。这些文件在进入构象搜索、药效团建模或机器学习训练之前,必须先经过一轮系统化的清洗与标准化。本节围绕 SDF、MOL2、PDB 三类常见格式,讲清楚“读进来之后到底要处理什么、怎么处理、处理到什么程度算合格”。

一、三类格式的读取要点与常见陷阱

SDF 是最常用的二维/三维分子交换格式,一个文件可包含多条记录,每条记录由标题行、程序行、注释行、计数行、原子块、键块以及可选的属性块(如 > <PUBCHEM_...>)组成,记录之间以 $$$$ 分隔。读取时最容易出问题的地方有三处:一是计数行中的原子数与实际原子块行数不一致,导致解析错位;二是属性块中字段名大小写不统一,后续按字段取值时匹配失败;三是部分记录缺少三维坐标,只有二维坐标甚至全零坐标,若直接用于构象生成会得到退化结构。用 RDKit 读取时,推荐显式关闭 sanitize 再逐步处理:

from rdkit import Chem

suppl = Chem.SDMolSupplier("input.sdf", sanitize=False, removeHs=False)
for mol in suppl:
    if mol is None:
        continue
    # 先做基础检查,再决定是否 sanitize
    if mol.GetNumAtoms() == 0:
        continue
    try:
        Chem.SanitizeMol(mol)
    except Exception as e:
        print("sanitize failed:", e)
        continue

MOL2 格式由 Tripos 提出,常用于对接与力场计算,其原子类型(如 C.3、N.am、O.co2)携带了杂化与化学环境信息,但不同软件写出的 MOL2 在原子类型命名、子结构块(@<TRIPOS>MOLECULE、@<TRIPOS>ATOM、@<TRIPOS>BOND)的完整性上差异很大。读取 MOL2 时要特别关注:键级是否缺失(很多对接输出只给单键)、电荷字段是否为 0、以及是否存在重复原子坐标。若键级缺失,后续芳香性感知会失败,需要根据原子类型和几何距离重建键级。

PDB 格式面向生物大分子,但小分子配体常以 HETATM 记录出现。读取 PDB 时,CONECT 记录不一定完整,氢原子通常缺失,且同一配体可能因晶体学无序被拆成多个 altLoc。处理策略是:先按 HETATM 提取目标配体,检查 CONECT 是否覆盖全部重原子,若缺失则用距离阈值(如 1.9 Å 以内视为成键)补全连接表,再统一加氢。

二、盐离子、溶剂与重复结构的处理

盐离子和溶剂是数据清洗中最常见的“噪声”。Na⁺、K⁺、Cl⁻、Ca²⁺、Mg²⁺、SO₄²⁻ 等通常不是目标分子的一部分,除非研究的是金属酶配位或离子通道。溶剂分子如 HOH、EDO、GOL、DMS、PEG 片段也应剔除。处理时不能简单按元素过滤,因为有些配体本身含氯或含硫。更稳妥的做法是维护一个“保留/剔除”白名单与黑名单,结合残基名、原子数、是否与目标分子有共价连接来判断。

重复结构有两种含义:一是同一文件内多条记录表示同一分子(如不同质子化态或不同构象),二是同一分子被重复导入。判断重复不能只看标题或 ID,而应基于标准化后的 InChIKey 或 canonical SMILES。对于同一分子的多构象记录,若目的是构象集成,应保留;若目的是去重建库,则只保留能量最低或坐标最完整的一条。

三、芳香性感知、质子化状态与标准化流程

芳香性感知是标准化中最容易出错的环节。RDKit 默认使用 Chem.SanitizeMol 中的 SANITIZE_SETAROMATICITY 步骤,依据 Hückel 规则和环系电子数判断。但若输入分子来自 MOL2 且键级全为单键,芳香环会被误判为脂肪环。此时应先根据原子类型和键长重建双键,再执行芳香性感知。对于含氮杂环,还要注意吡啶型氮与吡咯型氮的区别,前者不贡献孤对电子,后者贡献。

质子化状态直接影响分子电荷、氢键供受体数量和后续对接打分。生理 pH 7.4 下,羧酸通常去质子化为 COO⁻,脂肪胺通常质子化为 NH₃⁺,咪唑环以中性形式为主。处理时可用 rdkit.Chem.MolStandardize 中的 Uncharger 与 Normalizer,或调用 OpenBabel 的 --partialcharge 与 -p 7.4 选项。需要强调的是,质子化状态不是“越标准越好”,而应服务于具体任务:若研究的是胃酸环境下的药物吸收,pH 1.5 下的质子化态才更合理。

标准化流程可以概括为一条可复用的处理链:读取 → 去盐去溶剂 → 去重 → 键级与芳香性修正 → 质子化态调整 → 加氢 → 坐标检查 → 输出。下面这张图展示了该流程中各类输入文件经过的关键处理节点与输出判断,便于在实际脚本中对照实现。

分子文件读取与标准化处理链 SDF 输入 多记录 / 属性块 / $$$$ MOL2 输入 原子类型 / 键级可能缺失 PDB 输入 HETATM / CONECT 不全 读取与解析 原子块 / 键块 坐标 / 电荷 / 属性 去盐 / 去溶剂 Na⁺ K⁺ Cl⁻ HOH EDO 保留共价连接片段 重复结构判定 InChIKey / canonical SMILES 多构象按需保留 键级与芳香性 重建双键 / 环系电子数 吡啶氮 vs 吡咯氮 质子化状态 pH 7.4 生理条件 COO⁻ / NH₃⁺ / 中性 加氢与坐标检查 补全氢原子 检测全零 / 二维坐标 生成三维构象 标准化输出 SDF / SMILES / InChIKey 电荷与氢原子完整 异常记录 解析失败 / 空原子 sanitize 报错 解析异常 关键判断依据 输入格式差异:SDF 看计数行与属性块,MOL2 看原子类型与键级,PDB 看 CONECT 完整性 去盐去溶剂:按残基名与共价连接判断,不盲目按元素过滤 芳香性与质子化:先修键级再感知芳香性,质子化态按目标 pH 选择

四、可执行的清洗检查清单

为了让清洗过程可复现、可审计,建议在脚本中固定输出一份检查记录。以下清单可直接作为函数返回值或日志字段使用:

检查项判断条件处理动作
原子数一致性计数行原子数 = 原子块行数不一致则丢弃并记录
空记录原子数为 0 或全部坐标为 0丢弃
盐离子/溶剂残基名在剔除列表且无共价连接删除该片段
重复分子InChIKey 前 14 位相同按策略保留一条
芳香性环系原子是否被标记为芳香重建键级后重新感知
质子化态羧酸/胺基/咪唑是否符合目标 pH调整电荷与氢数
氢原子重原子价键是否饱和补全氢原子
三维坐标坐标是否退化或缺失重新生成构象
实际项目中,建议把上述检查项封装为一个 clean_molecule(mol, ph=7.4) 函数,返回清洗后的分子对象与一份字典形式的检查报告。这样在批量处理数千条记录时,可以快速定位是哪一类问题导致失败率偏高,而不是逐条人工排查。
芳香性感知与质子化状态之间存在耦合:先调整质子化态可能改变环系电子数,从而影响芳香性判断。推荐顺序是先去盐去溶剂、去重,再修键级与芳香性,最后调整质子化态并加氢。若顺序颠倒,可能出现羧酸去质子化后咪唑环芳香性丢失的异常情况。

完成上述清洗后,分子才具备进入构象工程环节的基本条件:原子类型明确、电荷合理、芳香性正确、氢原子完整、坐标可用于后续优化。下一阶段我们将在此基础上讨论构象搜索与能量评估,但前提是这一步的数据质量已经过关。