RDKit 快速上手：分子读写、分子可视化、子结构匹配、分子描述符计算
RDKit 是化学信息学领域的瑞士军刀。我最早接触它是在做药物虚拟筛选的项目里，当时手动解析 SDF 文件差点疯掉，后来发现 RDKit 一行代码就能搞定——嗯，从那天起我就再也没离开过它。这一章我们直奔最常用的四个技能：读写分子、画出来看看、搜某段结构、算各种描述符。学完你就能在 Python 里流畅地操作分子数据了。
1. 分子读写 —— 从文件到 Python 对象
先说怎么把分子捞进来。RDKit 统一用一种叫 Mol 的对象表示分子，不管你从 SMILES、SDF、还是 Mol2 读入，最后都是 Mol。我个人最常用的是 SMILES 字符串，简单轻量。
from rdkit import Chem
# 从 SMILES 读入
mol = Chem.MolFromSmiles('CCO')
print(mol.GetNumAtoms()) # 输出 3MolFromSmiles 如果失败会返回 None。我在项目中遇到过用户上传了一堆乱写的 SMILES，程序直接崩了。后来我养成习惯：
mol = Chem.MolFromSmiles(smiles_str)
if mol is None:
print(f"⚠ 无效 SMILES: {smiles_str}")
continue读 SDF 文件 更常见，尤其是从 ZINC、ChEMBL 下载的数据。RDKit 用 SDMolSupplier 按顺序读取所有分子：
suppl = Chem.SDMolSupplier('molecules.sdf')
for mol in suppl:
if mol is None:
continue
# 处理 mol写入也很直接。写 SMILES 用 MolToSmiles，写 SDF 用 SDWriter：
smiles_out = Chem.MolToSmiles(mol)
writer = Chem.SDWriter('output.sdf')
writer.write(mol)
writer.close()SDWriter，导致文件只写了一半。建议用 with 上下文管理器，或者手动 close()。另外，SDF 里可能包含多个分子的属性数据，mol.GetPropsAsDict() 可以一次性取出来。2. 分子可视化 —— 看一眼你的分子长什么样
光有数据不够，你得能“看见”分子长什么样。RDKit 的 Draw 模块可以快速生成 PNG 图片，在 Jupyter 笔记本里直接显示：
from rdkit.Chem import Draw
from IPython.display import display
mol = Chem.MolFromSmiles('CC(C)NC1=C(C(=O)NC1=O)C2=CC=CC=C2')
img = Draw.MolToImage(mol, size=(400, 300))
display(img)画多个分子并排对比也很实用：
mols = [Chem.MolFromSmiles(s) for s in ['CCO', 'CCN', 'CCC']]
img_grid = Draw.MolsToGridImage(mols, molsPerRow=3, subImgSize=(250,200))
display(img_grid)我自己的经验：当你要汇报分子设计结果时，把分子画成带原子编号、带立体构型、带 highlight 的图会非常有说服力。RDKit 支持给特定原子着色，比如标记出匹配的子结构位置——这部分在下一小节会用到。
3. 子结构匹配 —— 查查分子里有没有某个片段
你有一个候选分子，想查它是否包含某个药效团。或者你在筛库，想找出所有含苯环的分子。子结构匹配就是干这个的。
RDKit 用 MolFromSmarts 定义查询模式（SMARTS 语言），然后用 GetSubstructMatches 返回所有匹配位置的原子索引：
mol = Chem.MolFromSmiles('CC(=O)OC1=CC=CC=C1C(=O)O')
# 查羧基 -C(=O)OH
pattern = Chem.MolFromSmarts('C(=O)[OH]')
matches = mol.GetSubstructMatches(pattern)
print(matches) # 输出 ((7, 8, 9),)如果只要判断有没有匹配（不管几次），用 HasSubstructMatch：
if mol.HasSubstructMatch(pattern):
print("分子中含有羧基")我曾在做代谢物预测时犯过一个错：SMARTS 里忘记加原子类型限制，匹配出了一堆假阳性。后来我习惯先拿几个已知阳性和阴性的分子测试一下模式再大规模跑。
Draw.MolToImage 的 highlightAtoms 参数画出命中区域，效果非常直观。4. 分子描述符计算 —— 把分子变成数字
机器学习需要数值特征，分子描述符就是那把“翻译”钥匙。RDKit 的 Descriptors 模块内置了几百种描述符，从最简单的分子量到复杂的拓扑指数都有。我通常先算这几个最常用的：
| 描述符 | 含义 | 代码 |
|---|---|---|
| MolWt | 分子量 | Descriptors.MolWt(mol) |
| MolLogP | 脂水分配系数 | Descriptors.MolLogP(mol) |
| TPSA | 拓扑极性表面积 | Descriptors.TPSA(mol) |
| NumHDonors | 氢键供体数 | Descriptors.NumHDonors(mol) |
| NumHAcceptors | 氢键受体数 | Descriptors.NumHAcceptors(mol) |
| NumRotatableBonds | 可旋转键数 | Descriptors.NumRotatableBonds(mol) |
一次性算全部描述符可以用 CalcMolDescriptors，它返回一个字典：
from rdkit.ML.Descriptors import MoleculeDescriptors
calc = MoleculeDescriptors.MolecularDescriptorCalculator(['MolWt','MolLogP','TPSA'])
results = calc.CalcDescriptors(mol)
print(dict(zip(calc.GetDescriptorNames(), results)))为什么我推荐先熟悉这六个？它们对应着经典的 类药五规则（Lipinski's Rule of Five）。我在做早期筛选时，常常先用这五个描述符做快速过滤，跑一轮能干掉 60%-70% 的候选分子，后面再用复杂模型精筛。
5. 实战小案例 —— 批量计算并筛选
把上面的技能串起来。假设你有一个 library.sdf，里面存了 1000 个候选分子，你想找出那些同时满足 “分子量 < 500” 且 “LogP ≤ 5” 的分子，并保存为新的 SDF。
from rdkit import Chem
from rdkit.Chem import Descriptors, SDWriter
suppl = Chem.SDMolSupplier('library.sdf')
writer = SDWriter('filtered.sdf')
for mol in suppl:
if mol is None:
continue
mw = Descriptors.MolWt(mol)
logp = Descriptors.MolLogP(mol)
if mw < 500 and logp <= 5:
writer.write(mol)
writer.close()
print("筛选完成")这个脚本在我的日常工作中反复出现。我一般在筛选之前还会加一句 Chem.SanitizeMol(mol) 来确保分子结构没有异常，避免后续计算描述符时出错。