RDKit 分子操作：读取写入、子结构搜索与分子绘图

讲完了 SMILES 的基本概念和批次加载，咱们该动手碰碰真正的分子对象了。RDKit 在整个分子信息学里就像一把瑞士军刀——能读、能写、能搜、能画。这一章我们就挨个儿把这些功能走一遍。

我个人习惯把“读取”和“写入”看作一个闭环。你从文件读进来，处理完再写出去，中间才能去做各种分析。子结构搜索则是 RDKit 最核心的“找茬”能力——给定一个片段，看分子里有没有。绘图嘛，是为了让你肉眼确认结果，别光看数字。

3.1 分子读取：从 SMILES 到 RDKit Mol 对象

最基础的操作是把 SMILES 字符串转成 rdkit.Chem.rdchem.Mol 对象。一个 Mol 对象里不光有原子和键，还有立体化学、环信息、氢原子计数等等。RDKit 在背后做了很多清理工作。

核心函数：Chem.MolFromSmiles(smiles)

from rdkit import Chem

smiles = "CC(=O)Oc1ccccc1C(=O)O"
mol = Chem.MolFromSmiles(smiles)

if mol is not None:
    print("成功读取分子，原子数:", mol.GetNumAtoms())
else:
    print("SMILES 解析失败")

注意，如果 SMILES 语法有问题或者 RDKit 无法理解，会返回 None。我在项目里吃过这个亏：有批 SMILES 里混了空格和换行，没加校验直接传进去，跑着跑着就崩了。所以后来我习惯在读取后立刻判空。

除了 SMILES，RDKit 还能从 Mol 文件、SDF 文件、InChI 字符串甚至分子名字（通过 PubChem 查找）创建 Mol 对象。不过我们目前聚焦在 SMILES 上。

3.2 分子写入：把 Mol 对象转回字符串或文件

处理完分子后，我们经常要把它写回 SMILES 或 SDF。写回 SMILES 时，RDKit 会“规范化”——比如把环的表示收拢成标准形式。

smiles_out = Chem.MolToSmiles(mol)
print(smiles_out)  # 例如输出: O=C(O)c1ccccc1OC(C)=O

你可能会发现输入和输出的 SMILES 不一样。别紧张，这很正常。RDKit 的规范化会把分子拍平成最简形式，比如把连在环上的羧基调整顺序。如果你的任务要求保留原始 SMILES，记得在用 MolFromSmiles 前备份一份。

写入文件（比如 SDF）也很常见，尤其是做虚拟筛选或对接的时候：

from rdkit import Chem
writer = Chem.SDWriter("output.sdf")
writer.write(mol)
writer.close()

有个小坑：SDF 文件里可以存多个分子，如果你要批量写入，记得在循环里用 writer.write() 一连写进去。我早期写代码忘了关 writer，导致文件锁住读不出来——后来我都用 with 语句了。

我的经验：批量读写时用 Chem.SDMolSupplier 和 Chem.SDWriter 配合，效率很高。

3.3 子结构搜索：找分子里是否包含某个片段

子结构搜索是 RDKit 最常用的功能之一。你给一个查询片段（比如羟基、苯环、羧基），它能判断目标分子里有没有这个片段，还能定位在哪里。

核心步骤两步走：

  1. 创建查询分子：用 Chem.MolFromSmiles(query_smiles) 生成查询片段。
  2. 执行匹配：mol.HasSubstructMatch(query_mol) 返回 True/False，或者用 mol.GetSubstructMatches(query_mol) 得到所有匹配的原子索引。
from rdkit import Chem

# 目标分子：阿司匹林
aspirin = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")

# 查询：羧基 -C(=O)OH
carboxyl = Chem.MolFromSmiles("C(=O)O")

if aspirin.HasSubstructMatch(carboxyl):
    matches = aspirin.GetSubstructMatches(carboxyl)
    print(f"找到 {len(matches)} 个羧基，原子索引: {matches}")
else:
    print("无羧基")

你知道我当年做药物虚拟筛选时，最头疼的就是非精确匹配的问题。比如你想找苯环，但分子里有取代的苯环（比如对氯苯），普通子结构搜索不会自动考虑“任意取代”的模式。这时候就要用到 SMARTS 模式了。

SMARTS 是 SMILES 的扩展，支持通配符、原子属性限制等。比如任意芳香环上的碳，可以用 [c]。而“任意芳香杂环”可以用 [c,n,o,s] 枚举。

# 查询任意五元芳香杂环（含一个N）
query = Chem.MolFromSmarts("[c,n]1[c,n][c,n][c,n][c,n]1")
mol_demo = Chem.MolFromSmiles("c1cc[nH]c1")  # 吡咯
print(mol_demo.HasSubstructMatch(query))  # True

注意：SMARTS 写起来比 SMILES 灵活但容易出错。我建议先用 Chem.MolFromSmarts 测试一下，看看匹配结果对不对。

3.4 分子绘图：一图胜千言

RDKit 自带的绘图功能非常强大——不需要安装额外包就能把分子结构画出来。核心函数是 Draw.MolToImage(mol) 或 Draw.MolToFile(mol, "output.png")。

但注意：如果你在 Jupyter Notebook 里，直接用 Draw.MolToImage 就能显示。在脚本里运行时，需要指定保存路径。

from rdkit import Chem
from rdkit.Chem import Draw

aspirin = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")
img = Draw.MolToImage(aspirin, size=(300, 300))
img.save("aspirin.png")  # 保存到当前目录

也可以绘制多分子网格图，方便对比数据集：

mols = [Chem.MolFromSmiles(smi) for smi in ["CCO", "CCN", "CC(=O)O"]]
img_grid = Draw.MolsToGridImage(mols, molsPerRow=3, subImgSize=(200,200))
img_grid.save("grid.png")

如果你要给分子标记原子索引或者标注子结构匹配的结果，用 Draw.MolToImage 的 highlightAtoms 参数就好：

# 高亮羧基
mol = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O")
query = Chem.MolFromSmiles("C(=O)O")
matches = mol.GetSubstructMatches(query)
# 展平所有匹配原子索引
highlight_atoms = [atom for match in matches for atom in match]
img = Draw.MolToImage(mol, highlightAtoms=highlight_atoms, size=(400, 300))
img.save("highlighted.png")

技巧：在高亮时可以用 highlightBonds 参数同步高亮匹配的键，让可视化更清晰。

3.5 本章知识体系总览

下面这张图梳理了 RDKit 分子操作的核心脉络：从 SMILES 出发，经过读取/写入、子结构搜索和绘图，形成闭环。

RDKit 分子操作核心流程 SMILES 字符串 MolFromSmiles → Mol MolToSmiles / SDWriter 子结构搜索 / SMARTS Draw.MolToImage 结果文件 / 图像 应用场景：虚拟筛选、 分子验证、结构绘图

你能看出来，整个流程是单向但分叉的：先读取，再根据需求选择写入、搜索或绘图。当然也可以组合使用，比如先搜索出匹配的子结构，然后高亮绘制出来。

3.6 综合案例：从 SMILES 文件批量读取、搜取代谢物并绘图

咱们来个实战。假设你有一个 CSV 文件（molecules.csv），第一列是 SMILES，第二列是名字。你想找出所有含羧基的分子，并把它们画到一张网格图里。

import pandas as pd
from rdkit import Chem
from rdkit.Chem import Draw

# 读取数据
df = pd.read_csv("molecules.csv")
smiles_list = df["SMILES"].tolist()
names_list = df["Name"].tolist()

# 子结构搜索
carboxyl = Chem.MolFromSmiles("C(=O)O")
matched_mols = []
matched_names = []

for smi, name in zip(smiles_list, names_list):
    mol = Chem.MolFromSmiles(smi)
    if mol is None:
        continue
    if mol.HasSubstructMatch(carboxyl):
        matched_mols.append(mol)
        matched_names.append(name)

# 绘图
if matched_mols:
    img = Draw.MolsToGridImage(
        matched_mols,
        legends=matched_names,
        molsPerRow=4,
        subImgSize=(250, 250)
    )
    img.save("carboxyl_mols.png")
    print(f"找到 {len(matched_mols)} 个含羧基分子，已保存图片")
else:
    print("未发现含羧基分子")

这个小脚本几乎涵盖了本章的所有知识点：SMILES 读取、子结构搜索、批量绘图。我在实际项目中经常用这种模式快速筛查活性片段。

关键点回顾

  • 读取用 MolFromSmiles，记得判空
  • 写入用 MolToSmiles 或 SDWriter
  • 子结构搜索用 HasSubstructMatch 或 GetSubstructMatches，复杂条件用 SMARTS
  • 绘图单分子 MolToImage，多分子 MolsToGridImage，高亮用 highlightAtoms

好了，这一章基本就是这样。RDKit 的分子操作其实很直观——你只要记住每个步骤的对应的函数名，剩下的就是多写多试。我刚开始用的时候也经常忘传参数，画出来的图歪歪扭扭。不过别怕，多调试几次就会有手感。

下一次我们就要进入差异分析的实质了——如何比较两个批次之间的分子特征。到时候我们会用到指纹和相似度，但这些都是建立在今天学的分子对象基础之上的。

← 上一章 📖 返回目录 下一章 →