RDKit 2D/3D 分子处理全流程：从 SMILES 到构象优化

做分子性质预测、虚拟筛选或者化合物结构分析的时候我最常用也最离不开的一个工具就是 RDKit。不管是手里拿到的还只是一行 SMILES还是要处理已经带有三维坐标的文件绝大部分工作最后都能落到这个工具包上。简单说RDKit 就是化学信息学领域处理 2D 结构式和 3D 构象的一把瑞士军刀它能帮你完成从分子读入、清洗、画图、子结构匹配、指纹相似度计算到三维构象生成、力场优化、分子叠加对齐这一整条链路。这篇文章主要是想聊聊我平时用 RDKit 处理 2D 和 3D 分子的完整工作流包括每一步为什么要这么做、参数到底怎么选、以及那些文档里不会明说但实际一定会踩的坑。不管你是刚接触化学信息学的学生还是已经在做药物设计、材料筛选的开发人员这套流程应该都能直接拿去用。1. 环境准备与分子读入从一行 SMILES 开始的完整流程1.1 装环境为什么优先选 conda-forgeRDKit 的安装方式其实有不少pip 也能装但我个人强烈建议用 conda 从 conda-forge 通道装。原因很简单RDKit 的底层依赖了 Boost、NumPy 这些比较敏感的编译环境用 pip 装预编译包偶尔会遇到 OpenMP 版本不匹配或者某些绘图功能缺失的问题。conda-forge 上的版本更新及时而且会把 cairo、boost、numpy 这些配套依赖一起解决掉尤其是后面要用到 2D 描绘功能时没有 cairo 支持会特别难受。conda create -n rdkit-env -c conda-forge rdkit python3.11 conda activate rdkit-env装好之后我通常建议先跑一个自检脚本确认关键模块都能正常导入同时看一眼版本号方便后续排查问题from rdkit import Chem from rdkit.Chem import AllChem, Draw, rdMolDraw2D print(Chem.RDKFingerprint)如果你的环境已经有其他深度学习依赖也可以直接把 rdkit 加进去不用刻意隔离。只是要注意一点不要把 conda-forge 和 defaults 通道混得太乱否则解算依赖时容易卡半天。1.2 读入分子时的常见坑位RDKit 读入 2D 分子最基础的方式就是Chem.MolFromSmiles一行 SMILES 就能生成一个Mol对象。比如乙醇是CCO苯是c1ccccc1布洛芬大概是CC(C)Cc1ccc(C(C)C(O)O)cc1。但你如果直接拿这个函数去处理一批数据很快就会发现它有个让人头疼的地方——解析失败时不会抛异常而是默默返回None。这就意味着你写循环批量读 SMILES 的时候不能想当然地认为每个分子都成功解析了一定要做空值判断并且把失败记录打出来不然程序跑到一半才发现某个分子丢了排查起来很麻烦。from rdkit import Chem smiles_list [CCO, c1ccccc1, invalid_smiles] for smi in smiles_list: mol Chem.MolFromSmiles(smi) if mol is None: print(f解析失败: {smi}) continue print(f成功: {smi}, 原子数 {mol.GetNumAtoms()})还有一个新手容易忽略的点RDKit 内部的Mol对象本质上是一张图原子是节点键是边。它默认不显式保存氢原子所以你去数原子数时乙醇CCO只有 3 个重原子而不是 9 个原子。这一点在后续做 3D 构象生成时尤其重要因为三维坐标计算时通常需要先通过Chem.AddHs把氢原子显式加回来。2. 2D 层面到底能做什么子结构匹配、指纹与相似度2.1 子结构搜索不是简单的字符串匹配处理 2D 分子最频繁的需求就是子结构搜索。比如你想在一批化合物里筛出所有含有羧基的分子天然想法可能是直接在 SMILES 里找C(O)O这种字符串但这种方法很快会出问题因为同一个官能团在不同的分子里写出来可能长得完全不一样。RDKit 的方式是让你定义 SMARTS 模式它是一个描述“原子类型 连接关系 环境”的模式语言。举个例子比较稳健的羧基 SMARTS 可以写成[CX3](O)[OX2H1]。这里CX3表示碳原子连接了 3 个其他原子O是羰基氧OX2H1是连了 2 个原子的羟基氧。比单纯找字符串靠谱得多。carboxyl Chem.MolFromSmarts([CX3](O)[OX2H1]) mol Chem.MolFromSmiles(CC(C)Cc1ccc(C(C)C(O)O)cc1) match mol.GetSubstructMatch(carboxyl) print(match) # 返回匹配到的重原子索引元组实际操作中还有一个非常实用的细节如果同一个模式在一个分子里出现多次你希望全部找出来而不是只取第一组就要用GetSubstructMatches它返回的是所有匹配结果。做构象叠加或骨架对齐时这个功能特别有用因为一个分子里可能有两个甚至更多的候选匹配位点需要枚举后取最优的那一个。2.2 指纹计算怎么选 ECFP 还是 MACCS相似度计算是 2D 层面另一个核心操作。做虚拟筛选、聚类、去重都离不开分子指纹。RDKit 里最常用的两类指纹一类是 Morgan 指纹对应大家常说的 ECFP另一类是 MACCS keys。Morgan 指纹的本质是把每个原子周围的环境做“圆形展开”半径越大囊括的原子环境和键合信息越远。常用的 ECFP4 其实就是 Morgan 指纹半径等于 2、再按哈希方式映射到固定位长的位向量。我的经验是做药化分子相似性搜索和活性预测时Morgan 指纹通常表现最好因为它对分子的局部环境描述细腻而 MACCS keys 只有 166 位更多用于快速粗筛、骨架多样性分析。from rdkit import DataStructs from rdkit.Chem import AllChem mol1 Chem.MolFromSmiles(CCO) mol2 Chem.MolFromSmiles(CCOC) fp1 AllChem.GetMorganFingerprintAsBitVect(mol1, radius2, nBits2048) fp2 AllChem.GetMorganFingerprintAsBitVect(mol2, radius2, nBits2048) sim DataStructs.TanimotoSimilarity(fp1, fp2) print(sim)批量计算相似度时千万别写两层 for 循环然后逐个调TanimotoSimilarity数据量一大就慢得离谱。RDKit 提供了DataStructs.BulkTanimotoSimilarity一次算出一个指纹和一堆指纹的相似度列表速度提升非常明显。2.3 2D 描绘让别人一眼看懂你的分子很多时候你不需要把分子讲得多复杂直接画一张 2D 结构图就够了。RDKit 的rdMolDraw2D模块支持生成 PNG 或 SVG。我最常用的是MolDraw2DCairo因为它能直接输出位图放在论文和 PPT 里很方便。from rdkit.Chem import rdMolDraw2D mol Chem.MolFromSmiles(CC(C)Cc1ccc(C(C)C(O)O)cc1) drawer rdMolDraw2D.MolDraw2DCairo(400, 300) rdMolDraw2D.PrepareAndDrawMolecule(drawer, mol) drawer.FinishDrawing() drawer.WriteDrawingText(ibuprofen.png)如果目的不是示意图而是讲解机制我会把子结构匹配到的原子高亮出来或者把原子索引号画出来这样在一些远程协作场景里描述“第 7 号原子发生了什么反应”就会清楚得多。3. 3D 构象生成的完整链路从加氢到坐标输出3.1 为什么不能直接把 2D 结构拉成 3D很多刚接触的人以为 3D 结构就是给 2D 坐标加一个 Z 值就行这是个很大的误解。分子的三维构象不是平面结构简单抬高而是要在满足键长、键角、二面角优先倾向以及非键原子不重叠的前提下找到一套合理的空间坐标。RDKit 生成 3D 构象的核心逻辑是距离几何法。它先根据 2D 分子图的连接关系估算出原子两两之间的距离上下界然后在一个高维距离矩阵里随机采样再把采样结果投影到三维空间。这一步做完后分子的原子坐标已经有了但键长键角通常不够理想所以还要用力场优化一下。一个容易踩的坑做 3D 之前一定先调用Chem.AddHs(mol)把氢原子补上。因为 2D 表示里的氢通常是隐式的如果不加氢就生成构象你会得到一个没有氢、看起来好像也合理但实际物理上很怪的结构。而且后面很多力场参数是依赖氢原子的少了氢就没法正常算。3.2 EmbedMolecule 的关键参数与随机种子RDKit 里实际生成单构象的函数是AllChem.EmbedMolecule。推荐的做法是使用带经验知识约束的 ETKDG 方法也就是设置useExpTorsionAnglePrefsTrue和useBasicKnowledgeTrue它会参考已知的扭转角偏好和常见环构象库生成的构象合理度明显更高。mol Chem.MolFromSmiles(CCO) mol Chem.AddHs(mol) AllChem.EmbedMolecule(mol, randomSeed42, useExpTorsionAnglePrefsTrue, useBasicKnowledgeTrue) AllChem.MMFFOptimizeMolecule(mol) conformer mol.GetConformer() pos conformer.GetAtomPosition(0) print(pos.x, pos.y, pos.z) # 单位是埃这里有个很重要的经验一定要设置randomSeed。如果不设置程序内部会取系统随机数结果就是每次运行生成的构象可能都不一样。做实验、写论文、跑流程的时候如果你不固定随机种子别人根本没法复现你的结果。我一般习惯固定为 42 或者某个跟数据集编号有关的值。MMFFOptimizeMolecule这一行是力场优化默认选用 MMFF94 力场。MMFF94 对有机分子尤其是类药分子的几何优化效果通常不错。如果 MMFF 因为元素类型或参数缺失报错可以退一步用UFFOptimizeMolecule通用性更强但几何精度略差。3.3 多构象生成并不是越多越好有时候你需要的不只是一个构象而是尽量覆盖分子可能采取的多种空间形状比如做构象搜索或者 3D 相似性分析。这时可以用EmbedMultipleConfs一次性生成多个构象然后逐一优化。confs AllChem.EmbedMultipleConfs(mol, numConfs50, randomSeed1, useExpTorsionAnglePrefsTrue, useBasicKnowledgeTrue) res AllChem.MMFFOptimizeMoleculeConfs(mol)但我要提醒一句生成多少个构象取决于你的下游任务。虚拟筛选的 3D 对接每个配体生成 50 到 200 个构象通常就够用了如果你只是想要一个合理的展示结构生成一个就完了。我见过有人一上来就生成 500 个构象最后发现大部分都聚成一类白白浪费时间。生成之后还有一个重要问题你怎么判断这些构象是不是真的“不一样”这时候不能只看坐标而要算两两之间的 RMSD。RMSD 低于某个阈值比如 0.5 埃的两个构象就可以认为是重复的后续用聚类去冗余才有意义。3.4 3D 结构如何写出到文件得到 3D 构象后最常用的输出格式是 SDF。SDF 的好处是不仅能保存分子连接关系还能保存每个原子的三维坐标以及分子的属性字段。from rdkit.Chem import SDWriter mol Chem.MolFromSmiles(CCO) mol Chem.AddHs(mol) AllChem.EmbedMolecule(mol, randomSeed42) AllChem.MMFFOptimizeMolecule(mol) writer SDWriter(ethanol_3d.sdf) writer.write(mol) writer.close()写文件之前一定要确保分子确实带有坐标。判断方法很简单mol.GetNumConformers() 0就说明有构象。如果分子没有任何构象写出来也会报错或丢失 3D 信息。4. 让 2D 和 3D 协同工作子结构对齐与构象比较4.1 子结构匹配如何转成 3D 操作前面说 2D 子结构匹配是找原子索引而一旦加了氢、生成了构象原子的编号顺序就会变原来 2D 里第 0 号原子是碳加氢之后第 1 号原子可能是氢。很多人在这里栽跟头原因就是没搞清楚索引映射。正确做法是在加氢之前先做子结构匹配记录重原子索引再对Mol的副本做 3D 操作。如果你要在生成 3D 结构后继续用 2D 的匹配结果最好先把匹配到的原子索引存下来因为那组索引在AddHs之前是稳定的之后就不一定了。mol Chem.MolFromSmiles(CC(C)Cc1ccc(C(C)C(O)O)cc1) substructure Chem.MolFromSmarts([OH1]) match mol.GetSubstructMatch(substructure) mol3d Chem.AddHs(mol) AllChem.EmbedMolecule(mol3d, randomSeed42)这里你会发现match里保存的是重原子的索引在mol和mol3d上如果原子顺序没变前几个重原子的索引还是可以对应的。但如果你之后要移除氢、再做凯库勒化之类的操作索引可能就不可靠了所以最稳妥的策略是优先处理完所有 2D 逻辑之后再统一转向 3D。4.2 分子对齐与 RMSD 计算3D 层面最常用的操作之一是把两个分子的构象放到同一坐标系下比较。比如你想知道分子 A 的活性构象和分子 B 的某种构象在空间上差多少就可以找一个公共骨架做对齐然后算 RMSD。RDKit 的AlignMol可以完成这件事你提供参考分子和目标分子以及两组原子索引的映射关系它会通过最小二乘法旋转平移坐标让对应原子的空间距离尽量小然后返回 RMSD 值。from rdkit.Chem import AllChem core Chem.MolFromSmarts(c1ccccc1) mol1 Chem.MolFromSmiles(c1ccc(C)cc1) mol2 Chem.MolFromSmiles(c1ccc(C(O)O)cc1) for m in [mol1, mol2]: m Chem.AddHs(m) AllChem.EmbedMolecule(m, randomSeed42) AllChem.MMFFOptimizeMolecule(m) match1 mol1.GetSubstructMatch(core) match2 mol2.GetSubstructMatch(core) rmsd AllChem.AlignMol(mol1, mol2, atomMaplist(zip(match1, match2))) print(rmsd)这里我要特别强调一个细节如果公共骨架在分子里出现不止一次比如分子里有两个苯环GetSubstructMatch只返回第一组匹配不一定是最优的。更稳妥的做法是用GetSubstructMatches枚举所有匹配分别算对齐后的 RMSD取最小值。4.3 构象聚类的去冗余思路做构象搜索后经常需要对几十甚至几百个构象做去冗余。RMSD 是天然的距离度量有了构象间的 RMSD 矩阵之后可以用 RDKit 自带的 Butina 聚类或者直接用层次聚类。Butina 聚类的思路是选择周围邻居最多的构象作为类代表然后不断移除已聚类的点。这个算法速度很快而且不要求你预先指定类别数适合构象去重。实际操作上你只需要从每个构象里取坐标计算所有构象对之间的 RMSD。批量计算 RMSD 可以用AllChem.GetConformerRMSD或者GetBestRMSD。后者会先做对齐再算 RMSD更适合比较柔性分子的不同构象因为你去掉整体旋转和平移的影响后才能真正看出构象本身的差别。5. 高频报错与排查实录5.1 SMILES 解析失败这是最容易遇到也最好解决的问题。MolFromSmiles返回None时第一反应不应该是怀疑 RDKit 坏了而是检查 SMILES 本身的合法性。常见原因包括元素符号写错、价态不满足、方括号里的原子写法有问题。如果你想看到具体错误信息可以把sanitizeFalse传进去然后手动调用Chem.SanitizeMol它会告诉你到底是哪一步出了问题是价态不对、还是芳香性定义有冲突、还是某个原子没有正确闭合。mol Chem.MolFromSmiles(C1CCCC1, sanitizeFalse) try: Chem.SanitizeMol(mol) except Chem.AtomValenceException as e: print(e)还有一种情况是高通量筛选时遇到大量 SMILES 失败这时最好的策略是直接跳过并记录不要试图当场修复每个字符串。等流程跑完后再单独看失败列表很多问题其实是同一个原因造成的批量修起来更快。5.2 3D 嵌入失败EmbedMolecule偶尔会返回 -1意思是在距离几何采样阶段没能找到合理构象。常见场景是分子含有苛刻的刚性环系、金属配位环境或者某些原子的半径约束相互冲突。解决思路不是死磕一个函数而是多用EmbedMultipleConfs加randomSeed换种子多试几次。距离几何本质是随机采样换个种子往往就成功了。如果多试几次还是失败再考虑把问题分子用 UFF 力场重试优化或者检查 SMILES 是否包含不常见元素导致距离约束计算失败。5.3 加氢去氢后索引错乱这是 2D/3D 协同处理里最让人头疼的问题。加了氢之后原子数量变多索引变化是必然的。要避免混淆最关键的原则是想清楚你到底需要哪个阶段的信息。如果你需要的是“某个官能团的位置”就在 2D 阶段先匹配并保存索引如果你需要的是“氢原子在 3D 空间里的位置”就等加氢和构象生成之后再去遍历原子。不要在加氢之后还试图用原来的重原子索引去直接访问也不要反过来。5.4 力场优化报错MMFFOptimizeMolecule报错通常跟元素类型有关系特别是含有硼、硅、某些金属原子时MMFF94 的参数不齐全。这时候要么改UFFOptimizeMolecule要么检查分子是否需要先做电荷计算再优化而不是硬扛。6. 实用心得流程设计比单点操作更重要最后聊一点经验。处理 2D 和 3D 分子单点 API 用熟只是第一步真正考验人的是把整个流程串起来。比如一个典型任务给定一批 SMILES需要生成 3D 构象、计算指纹相似性、找出与参考分子最接近的候选。你需要明确每一步的数据格式转换SMILES 转MolMol加氢加氢后算指纹还是生成构象生成构象后怎么评估。只要一个环节输出格式不匹配后面再强也会被卡住。我自己在项目里常用的顺序是先用 2D 操作完成清洗和筛选比如剔除无效 SMILES、按子结构淘选感兴趣分子再对筛选后的分子统一加氢、生成构象、力场优化如果需要比较构象先对齐再算 RMSD避免直接比较绝对坐标所有过程固定随机种子并记录处理日志方便追溯。这个流程看起来很朴素但真的能救大命。尤其是数据集很大的时候提前在 2D 层面把不想要的分子筛掉能省掉后面大量 3D 构象生成和优化的时间。再提一个小技巧很多人画 2D 结构图时喜欢直接用默认布局但 RDKit 的 2D 坐标生成结果跟分子初始表示有关。如果你发现某个分子画出来特别奇怪、键长乱得没法看可以试试先把分子做一次Chem.MolFromSmiles再转一次 canonical SMILES很多时候重新读入一次就能拿到更规整的 2D 坐标。RDKit 这套工具真正强大之处不在于某个函数有多神奇而在于它把 2D 的结构逻辑和 3D 的空间逻辑无缝地接在了一起。理解这一点之后你再遇到“这个分子有没有羧基”“它的构象长什么样”“两个分子空间上像不像”这样的问题基本都能在几分钟之内写出可复用的代码解决。