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 下的质子化态才更合理。
标准化流程可以概括为一条可复用的处理链:读取 → 去盐去溶剂 → 去重 → 键级与芳香性修正 → 质子化态调整 → 加氢 → 坐标检查 → 输出。下面这张图展示了该流程中各类输入文件经过的关键处理节点与输出判断,便于在实际脚本中对照实现。
四、可执行的清洗检查清单
为了让清洗过程可复现、可审计,建议在脚本中固定输出一份检查记录。以下清单可直接作为函数返回值或日志字段使用:
| 检查项 | 判断条件 | 处理动作 |
|---|---|---|
| 原子数一致性 | 计数行原子数 = 原子块行数 | 不一致则丢弃并记录 |
| 空记录 | 原子数为 0 或全部坐标为 0 | 丢弃 |
| 盐离子/溶剂 | 残基名在剔除列表且无共价连接 | 删除该片段 |
| 重复分子 | InChIKey 前 14 位相同 | 按策略保留一条 |
| 芳香性 | 环系原子是否被标记为芳香 | 重建键级后重新感知 |
| 质子化态 | 羧酸/胺基/咪唑是否符合目标 pH | 调整电荷与氢数 |
| 氢原子 | 重原子价键是否饱和 | 补全氢原子 |
| 三维坐标 | 坐标是否退化或缺失 | 重新生成构象 |
clean_molecule(mol, ph=7.4) 函数,返回清洗后的分子对象与一份字典形式的检查报告。这样在批量处理数千条记录时,可以快速定位是哪一类问题导致失败率偏高,而不是逐条人工排查。
完成上述清洗后,分子才具备进入构象工程环节的基本条件:原子类型明确、电荷合理、芳香性正确、氢原子完整、坐标可用于后续优化。下一阶段我们将在此基础上讨论构象搜索与能量评估,但前提是这一步的数据质量已经过关。