RDKit入门：分子读取、可视化、基本操作与子结构搜索
RDKit是什么？说白了一个开源的化学信息学工具包。我当年刚接触计算化学的时候，最烦的就是处理分子文件格式。一个SMILES写法能有好几种变体，PDB里原子命名乱得一批。RDKit把这些破事全给你兜底了。这一章我就带你快速上手，把分子读进来、画出来、做点基本操作，再搞一下子结构搜索——这几个能力几乎是你后续做分子属性预测的基础。
安装与基本环境
先搞定环境。我个人习惯用conda，干净利落。
conda create -n rdkit-env python=3.9
conda activate rdkit-env
conda install -c conda-forge rdkit装完后跑个测试，看看能不能创建分子对象。
from rdkit import Chem
from rdkit.Chem import Draw
mol = Chem.MolFromSmiles('CCO')
print(mol.GetNumAtoms())输出应该是3。如果报错，八成是你环境里混了不同版本的numpy。嗯，这里要注意：RDKit对numpy版本有隐性要求，避坑指南——千万别用pip硬装rdkit-pypi，那玩意有坑，我踩过。
分子读取：多种格式一站搞定
RDKit支持的输入格式多得离谱。SMILES、SDF、Mol、PDB、MOL2，甚至还能读mmCIF。你在实际项目里遇到什么格式，它基本都能吞。
从SMILES读
SMILES是最常见的方式，一行字符串表达一个分子。但有个坑：同一个分子可以有不同SMILES写法。
mol1 = Chem.MolFromSmiles('CCO')
mol2 = Chem.MolFromSmiles('OCC')
# 上面两个都是乙醇，但原子顺序不一样我记得去年有个项目，两个团队对同一个数据集跑特征工程，结果特征数量对不上。查了一天，最后发现是SMILES写法不一致导致的。所以后来我统一用Chem.MolToSmiles(mol)做规范化，把分子转成标准SMILES再往下传。
从文件读
SDF是最常用的分子文件格式之一。它能存多个分子，是虚拟筛选数据集的标准容器。
supplier = Chem.SDMolSupplier('ligands.sdf')
for mol in supplier:
if mol is not None:
print(mol.GetNumHeavyAtoms())注意SDMolSupplier返回的是迭代器。有个细节：如果某个分子解析失败，它返回None而不是抛异常。所以一定要判空，不然你后面调用GetNumAtoms会当场崩。
从PDB读
PDB是蛋白质结构的主流格式，但RDKit读PDB有个默认行为：会帮你加氢。
mol = Chem.MolFromPDBFile('protein.pdb', removeHs=False)
# 用removeHs=False来保留原始状态我建议：除非你要做力学计算，否则别让RDKit自动加氢。蛋白PDB里的氢原子有时候是结晶水推断出来的，加错了反而影响后续子结构匹配。
可视化：一张图顶一千个SMILES
搞分子属性预测，不看结构瞎调参数，那叫盲人摸象。RDKit的可视化接口非常直接。
from rdkit.Chem import Draw
mol = Chem.MolFromSmiles('c1ccccc1C(=O)O')
img = Draw.MolToImage(mol, size=(400, 400))
img.save('benzoic_acid.png')如果你想批量画分子，使用MolToGridImage。
smiles_list = ['CCO', 'CC=O', 'CCN']
mols = [Chem.MolFromSmiles(smi) for smi in smiles_list]
img = Draw.MolsToGridImage(mols, molsPerRow=3, subImgSize=(250, 250))
img.save('grid.png')可视化在调试时极其有用。我经常这样用：读入一批分子，先画出前10个看看结构是否合理——有时候数据源造假的分子，主链都断开了，特征算出来全是NaN。
基本操作：增删改查四板斧
分子不是死的，你可以对它做各种变换。下面这几个操作是我在项目中用得最多的。
原子与键的遍历
mol = Chem.MolFromSmiles('CCO')
for atom in mol.GetAtoms():
print(atom.GetSymbol(), atom.GetAtomicNum(), atom.GetDegree())
for bond in mol.GetBonds():
print(bond.GetBeginAtomIdx(), bond.GetEndAtomIdx(), bond.GetBondType())原子序号从0开始。我遇到过一个坑：顺序不保证与文件输入一致。尤其在处理SDF中多分子时，每个分子的原子索引是独立的，你不能硬编码原子号。
加氢与去氢
mol = Chem.MolFromSmiles('CCO')
mol_H = Chem.AddHs(mol) # 显式加氢
mol_noH = Chem.RemoveHs(mol_H) # 去掉氢做力场优化或3D构象生成前必须加氢。但做子结构搜索时，显式氢反而容易干扰匹配。我通常的做法是：先加氢做计算，去氢后再做搜索。
分子规范化
mol = Chem.MolFromSmiles('[Na]OC(=O)c1ccccc1')
mol_norm = Chem.MolFromSmiles(Chem.MolToSmiles(mol))
# 这会把钠盐形式统一成规范写法规范化在数据集整合阶段特别重要。不同来源的分子写法五花八门，不规范化会导致重复计算或特征维度爆炸。
分子属性设置
mol.SetProp('_Name', 'Benzoic Acid')
mol.SetIntProp('total_atoms', mol.GetNumAtoms())
print(mol.GetProp('_Name'))我习惯在分子对象上挂载自定义属性，比如实验值、计算后的特征向量。RDKit的分子对象本质上是个dict，随你怎么存。但注意：序列化到SDF时只会保留字符串属性，数值属性会丢失。
子结构搜索：从大海里捞化学针
子结构搜索是化学信息学的核心操作。你想问一个分子里有没有某个官能团？搜一下就知道。你想在百万级化合物库中筛选含羧基的所有分子？也是搜。
SMARTS模式匹配
SMARTS是SMILES的升级版，支持通配符和逻辑算子。比如羧基的SMARTS模式是[CX3](=O)[OX1H1]——它表示一个碳原子（C），配位数为3（X3），双连氧，单连羟基氧。
mol = Chem.MolFromSmiles('CC(=O)O')
pattern = Chem.MolFromSmarts('[CX3](=O)[OX1H1]')
matches = mol.GetSubstructMatches(pattern)
print(matches)返回的是匹配到的原子索引元组。对于包含多个官能团的分子，GetSubstructMatches会返回所有出现位置。
是否包含子结构
mol = Chem.MolFromSmiles('c1ccccc1C(=O)O')
pattern = Chem.MolFromSmarts('[OH]')
print(mol.HasSubstructMatch(pattern)) # True这个返回值是布尔型，适合做条件筛选。我在项目中常常这样用：把SMARTS写进一个列表，循环判断每个分子是否符合某个官能团定义，然后生成特征向量——0/1，简单暴力但有效。
移除匹配的子结构
mol = Chem.MolFromSmiles('CC(=O)O')
pattern = Chem.MolFromSmarts('[CX3](=O)[OX1H1]')
rms = Chem.DeleteSubstructs(mol, pattern)
# 剩下乙基部分这个操作在分解分子骨架时很有用。比如你先移除所有甲基，再移除所有羟基，剩下的就是母核。
替换子结构
mol = Chem.MolFromSmiles('CC(=O)O')
pattern = Chem.MolFromSmarts('[CX3](=O)[OX1H1]')
replacement = Chem.MolFromSmiles('C')
new_mol = Chem.ReplaceSubstructs(mol, pattern, replacement)[0]这是做虚拟衍生物的利器。把骨架上的某个基团替换成其他基团，批量生成候选分子。注意ReplaceSubstructs返回的是列表，因为同一个分子可能有多个匹配位置。
本章知识图谱
常用操作速查表
| 功能 | 函数 | 典型返回 | 我踩过的坑 |
|---|---|---|---|
| 从SMILES创建分子 | Chem.MolFromSmiles | RDKit Mol对象 | 输入非法SMILES返回None，一定要判空 |
| 从文件读多分子 | Chem.SDMolSupplier | 迭代器 | 解析失败不抛异常，返回None |
| 可视化单分子 | Draw.MolToImage | PIL Image | 不给size参数默认200×200太小 |
| SMARTS子结构搜索 | Mol.GetSubstructMatches | 元组列表 | 有显式氢时容易漏匹配 |
| 替换子结构 | Chem.ReplaceSubstructs | Mol列表 | 返回列表不是单对象 |
| 分子规范化 | Chem.MolToSmiles | 标准SMILES字符串 | 不同输入可能得到同一SMILES |
避坑指南
我在实战中遇到最多的几个坑，一次性给你列出来：
- 原子序号不固定：同一分子从不同方式读入，原子索引顺序可能不同。做子结构匹配时不要硬编码原子索引，用GetSubstructMatches返回值代替。
- SDF属性丢失：往分子上设整数、浮点属性，导出SDF时只有字符串属性保留。如果数据要持久化，统一用
SetProp把数值转成字符串。 - 加氢导致子结构搜索偏慢：加氢后的分子节点数暴增，SMARTS匹配复杂度上升。我的习惯：搜索前用
Chem.RemoveHs去氢。 - PDB的多模型问题：
Chem.MolFromPDBFile默认只读第一个模型。如果你的PDB含多个构象，需要用PDBMultilineMolSupplier或自己轮询。
一个小项目串联本章知识
把今天讲的东西串起来。假设你收到一个SDF文件，里面有一堆候选分子，你需要在其中找到所有含羧基的分子，并把羧基替换成甲基，最后输出新分子的SMILES列表。
from rdkit import Chem
supplier = Chem.SDMolSupplier('candidates.sdf')
carboxyl = Chem.MolFromSmarts('[CX3](=O)[OX1H1]')
methyl = Chem.MolFromSmiles('C')
output = []
for mol in supplier:
if mol is None:
continue
if mol.HasSubstructMatch(carboxyl):
new_mols = Chem.ReplaceSubstructs(mol, carboxyl, methyl)
for new_mol in new_mols:
smi = Chem.MolToSmiles(new_mol)
output.append(smi)
print(f"共生成 {len(output)} 个衍生物分子")这个例子把分子读取、空值检查、子结构匹配、替换、SMILES输出整个串起来了。你在真实项目里，可能还要加个可视化确认一下——嗯，别嫌麻烦，看一眼结构能省你很多debug时间。