4. 化学信息学工具：RDKit安装与基本使用、分子读取与操作

聊到分子指纹，RDKit 是绕不过去的工具。说实话，我当初刚接触化学信息学的时候，第一个装的库就是它。RDKit 在药物筛选圈子里几乎是标配——处理分子、算指纹、做相似性比较，它都能搞定。我们这一章就先把 RDKit 装好，然后上手读写分子文件，再熟悉一些基础的分子操作方法。

4.1 安装 RDKit

安装 RDKit 我推荐直接用 conda。它的底层依赖（比如 Boost）比较麻烦，用 conda 能省掉一堆编译的麻烦。我个人习惯是专门为 RDKit 建一个环境，避免和其他库打架。

conda create -n rdkit_env python=3.9
conda activate rdkit_env
conda install -c conda-forge rdkit

等一下，如果你是 Windows 用户，装完后最好跑一下这个命令验证是否成功：

python -c "from rdkit import Chem; print(Chem.__version__)"

如果能打出版本号，恭喜你，环境搭好了。我后来碰到过一个坑：在 macOS 上如果用了 M1 芯片，conda 的默认 channel 可能会装到 x86 版本，效率大打折扣。建议指定 osx-arm64 的 channel，或者直接用 miniforge。

4.2 读取与写入分子文件

RDKit 支持的分子格式挺多的，日常最常用的就是 SDF 和 SMILES。我刚开始做项目时，手头有一批供应商的 SDF 文件，里面某些分子的原子坐标是乱的。嗯，之后被「坑」多了，就慢慢练出手感了。

4.2.1 从 SMILES 构建分子

from rdkit import Chem

# 从 SMILES 读取单个分子
smiles = "c1ccccc1C(=O)O"
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
    print("分子读取成功，原子数：", mol.GetNumAtoms())
else:
    print("SMILES 格式有误，请检查")

这里要注意，Chem.MolFromSmiles 如果解析失败会返回 None。我见过不少新手直接拿 None 往下算指纹，结果程序直接崩溃。所以你想想看，加个判空处理能省很多事。

4.2.2 从 SDF 文件读取

SDF 常用于存储 3D 构象或带属性的分子。RDKit 提供了批量读取的接口：

supplier = Chem.SDMolSupplier("ligands.sdf")
for i, mol in enumerate(supplier):
    if mol is None:
        continue
    print(f"分子 {i}: {mol.GetNumAtoms()} 个原子")

SDMolSupplier 是个生成器，懒加载，不会一下子把所有分子读进内存。这一点对于处理上百万个分子的大文件特别友好。不过，SDF 里有些分子格式不规范，比如氢原子数写错，会导致 mol 为 None。我建议加上跳过逻辑，避免只因为一个坏分子整个作业崩掉。

4.2.3 写入分子到文件

把处理好的分子写出去也是家常便饭：

# 写 SDF
writer = Chem.SDWriter("output.sdf")
for mol in mol_list:
    writer.write(mol)
writer.close()

# 写 SMILES
with open("output.smi", "w") as f:
    for mol in mol_list:
        smi = Chem.MolToSmiles(mol)
        f.write(smi + "\n")

写 SDF 时记得调用 close()，不然文件可能不完整。我在一次批量虚拟筛选后就吃过这个亏，凌晨跑完一看文件是空的，气得我直接加了 with 语句管理上下文。

4.3 分子基础操作

拿到分子对象后，我们常需要做一些操作，比如加氢、删氢、计算属性、子结构搜索。下面是我在项目中用得最多的几个。

4.3.1 显式氢与隐式氢

RDKit 默认的分子对象里氢原子是隐式的。如果需要做 3D 优化或对接，一般先要把氢显式化：

mol = Chem.MolFromSmiles("CCO")
mol_with_H = Chem.AddHs(mol)
print("隐式氢模型原子数：", mol.GetNumAtoms())      # 3
print("显式氢模型原子数：", mol_with_H.GetNumAtoms()) # 9

反过来，做完 3D 工作后，一般要删掉显式氢再存盘或算指纹：

mol_no_H = Chem.RemoveHs(mol_with_H)

4.3.2 计算分子属性

RDKit 自带一系列快速计算分子描述符的方法。这里列几个我在筛选 pipeline 里常用的：

from rdkit.Chem import Descriptors

mol = Chem.MolFromSmiles("c1ccccc1C(=O)O")
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
hbd = Descriptors.NumHDonors(mol)
hba = Descriptors.NumHAcceptors(mol)
print(f"MW: {mw:.1f}  LogP: {logp:.2f}  HBD: {hbd}  HBA: {hba}")

4.3.3 子结构搜索与匹配

药物筛选里常要判断分子是否包含某个药效团或警报结构（比如 PAINS）。RDKit 提供了子结构匹配：

pattern = Chem.MolFromSmarts("[$(C(=O)O)]")  # 羧基
if mol.HasSubstructMatch(pattern):
    print("分子中包含羧基")
    match_atoms = mol.GetSubstructMatch(pattern)
    print("匹配的原子索引：", match_atoms)

这里用的是 SMARTS 模式，比 SMILES 灵活很多。我常常用它来批量标记分子里的「问题片段」。比如 PAINS 过滤，把几十个警报 SMARTS 写成一个列表，一行代码扫过去，能做到秒级过滤。

4.4 核心流程与知识结构

下面这张 SVG 图展示了本章的 RDKit 操作流程，从读取分子到输出结果的完整路径：

4.5 实践中容易踩的坑

最后聊几个我实际经历过的教训，希望对你有帮助：

• 芳香性感知不一致：RDKit 默认的芳香性模型可能和第三方工具（比如 Open Babel）不同。同一根键可能在 RDKit 里是芳香键，在其他工具里是单双键交替。如果你需要跨工具交换数据，最好统一用 Kekulize 后的形式。
• 立体化学丢失：从 SDF 读取时，如果原子坐标不够精确或缺少顺反标记，RDKit 可能无法正确识别手性。我曾经遇到过一组手性分子，跑完筛选发现有一半的 R/S 反了，检查才发现源文件里已经丢掉了立体信息。我建议在读入后打印一下 Mol.GetStereo 做快速验证。
• 分子无效返回 None 但不报错：SDMolSupplier 对格式有问题的分子只是跳过，不会抛异常。我曾经有 5 万分子跑完筛选后，猛然发现实际只处理了 4.8 万。后来我加了一段计数代码：

total = 0
good = 0
for mol in supplier:
    total += 1
    if mol is not None:
        good += 1
print(f"总分子数：{total}，成功读取：{good}，丢失：{total - good}")

这么一查，就发现了那个带坏头的 SDF 文件。嗯，从那以后，我的 pipeline 里一直保留着这段日志。