1. 从一张画满荧光笔的分子图说起做化学信息学的人迟早会碰到一个绕不开的需求手里有一堆分子想快速知道某个片段到底出现在哪些结构里而且最好能一眼看出来它长在哪个位置。我第一次遇到这个场景是要从一批候选化合物里筛出所有含苯并咪唑骨架的分子当时傻乎乎地一个个用肉眼比对几十个分子下来眼睛都快瞎了还漏掉了好几个。后来才意识到这种活儿本来就该交给 RDKit 去做——它不仅能帮你把匹配的子结构找出来还能直接在分子图上用高亮标出来省时省力还不容易出错。这篇内容就是围绕“搜索并高亮展示分子中匹配的子结构”这件事展开的。核心工具是 RDKit一个在化学信息学领域被广泛使用的开源工具包。它提供了子结构搜索Substructure Search和高亮绘制Highlight Rendering两大能力前者负责判断“这个片段在不在分子里”后者负责把“在哪儿”可视化出来。适合的读者包括做虚拟筛选的药物化学研究者、处理化合物库的数据科学家、以及任何需要批量分析分子结构的人。哪怕你之前只用 RDKit 画过简单的分子式这篇内容也能让你把子结构搜索和高亮这套组合拳打顺。我下面会从最基础的匹配原理讲起然后一步步拆到代码层面把 SMARTS 的写法、匹配函数的区别、高亮参数的调优、以及批量处理时的性能坑都过一遍。中间会穿插我自己踩过的坑和实测有效的技巧尽量让你少走弯路。2. 子结构搜索的底层逻辑RDKit 到底在比什么2.1 图匹配而非字符串匹配很多人第一次听到“子结构搜索”下意识会以为是在做字符串包含判断就像 Python 里abc in abcdef那样。但分子结构的本质是图Graph原子是节点化学键是边。RDKit 做子结构搜索实际上是在做**子图同构Subgraph Isomorphism**的判断——问的是“查询片段这个图能不能嵌入到目标分子这个图里”。这个区别很关键。字符串匹配是线性的、位置敏感的而图匹配要考虑原子的连接关系、环的闭合方式、以及化学环境。举个例子查询片段是“苯环”目标分子里有一个苯环但苯环上还连着其他基团。字符串匹配可能会因为“多了字符”而失败但图匹配只关心苯环这个子图是否存在不关心它外面还挂了什么。这就是为什么子结构搜索必须用图算法而不能用简单的文本比对。RDKit 底层用的是 VF2 算法的变体来做子图同构判断。VF2 是一种经典的回溯搜索算法通过逐步扩展匹配节点、剪枝不满足条件的候选最终找到所有可能的嵌入方式。你不需要自己实现这个算法但理解它的存在有助于你明白为什么有些查询会返回多个匹配结果因为同一个片段可能在分子里有多种嵌入方式以及为什么某些复杂查询会变慢因为回溯搜索的状态空间可能很大。2.2 SMARTS描述子结构的专用语言要让 RDKit 知道你想找什么片段你得用一种它能读懂的语言来描述。这个语言就是SMARTSSMiles ARbitrary Target Specification。SMARTS 是 SMILES 的扩展SMILES 用来描述一个具体分子而 SMARTS 用来描述一类结构模式。举个最直观的例子。苯的 SMILES 是c1ccccc1这是一个具体的分子。如果你想表达“任意一个苯环”SMARTS 可以写成c1ccccc1但更通用的写法是[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1其中[#6]表示任意碳原子:表示芳香键。当然实际使用中很少有人写这么啰嗦通常直接用c1ccccc1就够了因为小写c在 SMARTS 里就代表芳香碳。SMARTS 的强大之处在于它支持各种通配符和逻辑操作。比如[#6]匹配任意碳[C,N]匹配碳或氮[!O]匹配非氧原子*匹配任意原子。还有针对化学环境的修饰符比如[C;H3]表示恰好连了三个氢的碳[C;R]表示处于环上的碳。这些语法让你能精确地描述“我要找的是什么”。我个人的经验是刚开始用 SMARTS 时最容易犯的错是把 SMILES 当 SMARTS 用。比如你想找羧基写了C(O)O这在很多情况下确实能匹配到但它也会匹配到酯基里的C(O)O部分因为酯基也包含这个子图。如果你只想找游离的羧酸得写成[CX3](O)[OX2H1]明确指定羟基上的氧连了一个氢。这种细节在批量筛选时特别重要写错一个修饰符结果可能差出几百个分子。2.3 匹配函数的家族HasSubstructMatch 与 GetSubstructMatchesRDKit 提供了两个层次的子结构搜索函数用途不同返回值也不同。HasSubstructMatch返回一个布尔值只告诉你“有没有匹配”不告诉你“匹配在哪里”。这个函数适合做快速筛选比如你有一个十万级的化合物库想先粗筛出所有含某个片段的分子用这个函数就够了速度快内存占用小。GetSubstructMatches返回一个元组的元组每个内层元组是一组原子索引表示一次匹配中查询片段的各个原子分别对应目标分子的哪些原子。这个函数适合需要精确定位的场景比如你要在高亮图上标出匹配位置就必须知道具体是哪些原子参与了匹配。还有一个GetSubstructMatch注意没有复数 s只返回第一个匹配的原子索引元组。如果你只关心“找到一个就行”用这个比GetSubstructMatches更省资源因为它找到第一个匹配就停了不会继续搜索所有可能的嵌入方式。选择哪个函数取决于你的下游需求。如果只是统计有多少分子含有某片段HasSubstructMatch足够如果要画高亮图必须用GetSubstructMatches或GetSubstructMatch。我见过有人用HasSubstructMatch筛完之后又对每个阳性分子调GetSubstructMatches来画图这其实多了一次搜索不如直接调后者然后判断返回是否为空。3. 高亮绘制的参数调优让匹配位置一眼可见3.1 高亮原理原子索引到颜色的映射RDKit 的绘制模块rdkit.Chem.Draw提供了Draw.MolToImage和Draw.MolToFile等函数它们都接受一个highlightAtoms参数。这个参数是一个原子索引的序列RDKit 会把这些原子以及它们之间的键用高亮颜色画出来。关键点在于高亮的是原子键的高亮是附带的。当你传入一组原子索引时RDKit 会检查这些原子之间的键如果一条键的两个端点都在高亮集合里这条键也会被高亮。所以如果你只高亮了一个原子它周围的键不会变色如果你高亮了整个苯环的六个碳环上的六条键也会跟着高亮。这个机制意味着GetSubstructMatches返回的原子索引元组可以直接喂给highlightAtoms不需要额外转换。比如匹配结果返回(0, 1, 2, 3, 4, 5)你直接把这个元组传给highlightAtoms苯环就会被完整高亮。3.2 颜色与样式的自定义默认的高亮颜色是黄色说实话在白色背景上对比度一般尤其是当分子本身有颜色或者打印出来的时候黄色高亮经常看不清。RDKit 允许你通过highlightColor参数自定义颜色传一个 RGB 元组就行。我通常用(0.8, 0.2, 0.2)这种偏红的颜色在白底上非常醒目。如果你需要同时高亮多个不同的子结构比如一个分子里既有苯环又有羧基你想用不同颜色区分那就得用highlightAtomColors参数。这个参数接受一个字典键是原子索引值是 RGB 元组。你可以把第一个匹配的原子染成红色第二个匹配的染成蓝色以此类推。还有一个容易被忽略的参数是highlightBonds它允许你单独指定要高亮的键。通常情况下你不需要用它因为原子高亮会自动带动键高亮。但在某些特殊场景下比如你想高亮一条键但不想高亮它的端点原子就得手动指定。不过这种需求很少见大部分时候用原子高亮就够了。3.3 图例与标注的添加画高亮图时经常需要在图旁边加个图例说明“红色高亮的是苯环蓝色高亮的是羧基”。RDKit 的Draw.MolToImage支持legend参数可以传一个字符串作为图例。如果你用Draw.MolsToGridImage画多个分子的网格图legends参数可以传一个字符串列表每个分子对应一个图例。图例的字体大小和位置 RDKit 没有直接暴露参数来调但你可以通过Draw.DrawingOptions来全局调整。比如Draw.DrawingOptions.legendFontSize可以改图例字号Draw.DrawingOptions.legendFraction可以改图例占整张图的高度比例。这些参数在批量出图时很有用尤其是当分子图很小、图例文字挤在一起的时候。我自己的习惯是如果只是自己看不加图例也行高亮颜色本身就说明了问题。但如果是给别人看或者放进报告里图例一定要加而且要写清楚 SMARTS 表达式这样别人能复现你的筛选条件。4. 从单分子到批量完整实操链路4.1 环境准备与依赖安装RDKit 的安装方式取决于你的操作系统和 Python 环境。最省事的方式是用 condaconda install -c conda-forge rdkit如果你用 pip也可以装pip install rdkit但要注意pip 版本的 RDKit 在某些平台上可能缺少一些可选依赖比如 Cairo 绘图后端。如果你发现画出来的图是空白的或者报错说找不到绘图后端大概率是这个问题。conda 版本通常更完整推荐优先用 conda。安装完之后验证一下from rdkit import Chem from rdkit.Chem import Draw print(Chem.__version__ if hasattr(Chem, __version__) else ok)如果没报错说明环境没问题。我建议在虚拟环境里装 RDKit因为它依赖的库比较多跟其他科学计算库混在一起容易出冲突。4.2 单分子搜索与高亮的完整代码下面这段代码是我平时最常用的模板从 SMILES 读入分子搜索子结构然后画高亮图from rdkit import Chem from rdkit.Chem import Draw # 目标分子 mol Chem.MolFromSmiles(c1ccc2c(c1)nc(n2)c3ccccc3) # 查询片段苯并咪唑骨架 query Chem.MolFromSmarts(c1ccc2c(c1)nc(n2)) # 搜索所有匹配 matches mol.GetSubstructMatches(query) # 画高亮图 if matches: img Draw.MolToImage(mol, highlightAtomsmatches[0], highlightColor(0.8, 0.2, 0.2), size(400, 300)) img.save(highlighted.png) else: print(未找到匹配)这段代码里matches是一个元组的元组。如果分子里只有一个苯并咪唑matches就只有一个元素如果有多个就会有多个。我通常只高亮第一个匹配因为多个匹配同时高亮会让图变得很乱。如果你确实需要展示所有匹配可以用highlightAtoms传一个扁平化的列表把所有匹配的原子都放进去但颜色就没法区分了。4.3 批量处理与网格图输出当你有一批分子要处理时逐个画图再拼起来太麻烦。RDKit 的Draw.MolsToGridImage可以直接画网格图from rdkit import Chem from rdkit.Chem import Draw smiles_list [ c1ccc2c(c1)nc(n2)c3ccccc3, c1ccc2c(c1)nc(n2)CCCC, c1ccccc1, c1ccc2c(c1)nc(n2)C(O)O ] query Chem.MolFromSmarts(c1ccc2c(c1)nc(n2)) mols [] highlights [] legends [] for smi in smiles_list: mol Chem.MolFromSmiles(smi) if mol is None: continue match mol.GetSubstructMatch(query) if match: mols.append(mol) highlights.append(match) legends.append(Match) else: mols.append(mol) highlights.append(()) legends.append(No match) img Draw.MolsToGridImage(mols, molsPerRow2, highlightAtomListshighlights, legendslegends, subImgSize(300, 250)) img.save(grid.png)这里有个细节highlightAtomLists参数接受一个列表每个元素对应一个分子的高亮原子索引。对于没有匹配的分子传空元组()就行RDKit 会正常画出分子但不加高亮。legends参数同理每个分子对应一个图例字符串。实测下来MolsToGridImage在分子数量超过 50 个时生成的图片会非常大加载和渲染都会变慢。如果你要处理几百个分子建议分批出图每批 20 到 30 个或者直接用MolsToFile逐个保存成单独的图片文件。5. 那些文档里不会写的踩坑经验5.1 芳香性感知导致的匹配失败这是我最常遇到的坑没有之一。RDKit 在读入分子时会做芳香性感知Aromaticity Perception决定哪些环是芳香环。如果你的 SMARTS 里用了小写字母如c表示芳香碳但目标分子在 RDKit 眼里那个环不是芳香环匹配就会失败。举个例子吡啶的 SMILES 是c1ccncc1RDKit 会正确识别为芳香环。但如果你从某些文件格式读入分子或者分子经过了某些化学修饰RDKit 可能没有正确感知芳香性导致c1ccncc1这个 SMARTS 匹配不上。解决办法是用Chem.SanitizeMol强制做一次完整的化学感知或者在 SMARTS 里用大写字母加芳香性修饰符比如[C;R]表示环上的碳不依赖芳香性判断。我个人的习惯是在批量处理前先对每个分子调一次Chem.SanitizeMol确保 RDKit 对分子的理解是一致的。虽然这会稍微增加处理时间但能避免很多莫名其妙的匹配失败。5.2 匹配数量爆炸与性能陷阱有些 SMARTS 查询会返回大量的匹配结果尤其是当查询片段很小或者很通用的时候。比如你用[#6]去匹配一个含 50 个碳的分子会返回 50 个匹配每个匹配就是一个碳原子。如果你用GetSubstructMatches并且没有设置uniquifyFalseRDKit 会返回所有可能的嵌入方式数量可能是指数级的。GetSubstructMatches有一个uniquify参数默认是True表示对结果去重。去重的规则是如果两组匹配的原子集合相同只保留一组。但即使去重之后匹配数量仍然可能很大。还有一个maxMatches参数可以限制返回的最大匹配数默认是 1000。如果你只关心“有没有匹配”或者“前几个匹配”把maxMatches设小一点能显著提速。另外GetSubstructMatches还有一个useChirality参数默认是False。如果你需要区分手性比如只想匹配特定构型的片段得把它设成True。但开启手性匹配会变慢因为算法需要额外比较原子的立体化学信息。5.3 高亮图的原子索引错位这个问题比较隐蔽通常出现在你从文件读入分子之后。RDKit 在读入分子时原子索引的顺序可能跟你预期的不同。比如你从 SDF 文件读入一个分子SDF 里原子的顺序是 A、B、C但 RDKit 内部可能会重新排列原子索引。如果你用外部工具生成的原子索引去高亮就会标错位置。解决办法是始终用 RDKit 自己返回的索引。GetSubstructMatches返回的索引就是 RDKit 内部的原子索引直接传给highlightAtoms不会错位。如果你需要跟外部数据对照先用mol.GetAtomWithIdx(i).GetSymbol()确认一下索引对应的原子是不是你想要的。还有一个相关的问题是当你对分子做了Chem.AddHs或者Chem.RemoveHs之后原子索引会变。加氢会插入新的氢原子索引全部后移去氢会删除氢原子索引前移。如果你在高亮之前做了这些操作一定要重新获取匹配结果不能用之前的索引。5.4 绘图后端的兼容性问题RDKit 的绘图功能依赖后端默认在 Jupyter 里用IPython.display显示在脚本里用PIL或Cairo保存文件。如果你在无图形界面的服务器上跑代码可能会遇到ImportError或者OSError提示找不到显示设备。解决办法是明确指定用PIL后端from rdkit.Chem.Draw import rdMolDraw2D drawer rdMolDraw2D.MolDraw2DCairo(400, 300) drawer.DrawMolecule(mol, highlightAtomsmatches[0]) drawer.FinishDrawing() with open(output.png, wb) as f: f.write(drawer.GetDrawingText())MolDraw2DCairo是纯文件输出的不依赖图形界面在服务器上跑很稳。MolDraw2DSVG则是输出 SVG 矢量图适合放进网页或者需要无损缩放的场景。我一般在服务器上批量出图时都用MolDraw2DCairo从来没出过问题。6. 进阶玩法让子结构搜索更贴合实际需求6.1 多片段联合搜索有时候你要找的不是单个片段而是“同时含有 A 和 B 两个片段”的分子。最直接的做法是分别搜索然后取交集query_a Chem.MolFromSmarts(c1ccccc1) query_b Chem.MolFromSmarts(C(O)O) match_a mol.GetSubstructMatches(query_a) match_b mol.GetSubstructMatches(query_b) if match_a and match_b: # 两个片段都存在 all_atoms list(match_a[0]) list(match_b[0]) img Draw.MolToImage(mol, highlightAtomsall_atoms)但这样高亮出来两个片段是同一种颜色分不清哪个是哪个。更好的做法是用highlightAtomColors给不同片段染不同颜色colors {} for idx in match_a[0]: colors[idx] (0.8, 0.2, 0.2) # 红色 for idx in match_b[0]: colors[idx] (0.2, 0.2, 0.8) # 蓝色 img Draw.MolToImage(mol, highlightAtomslist(colors.keys()), highlightAtomColorscolors)这样苯环是红色羧基是蓝色一目了然。如果你有更多片段可以继续加颜色但建议不要超过三种否则图会变得很花。6.2 用 SMARTS 做官能团筛选子结构搜索最常见的应用场景之一就是官能团筛选。RDKit 官方文档里有一份常用的 SMARTS 列表我摘几个我经常用的官能团SMARTS说明羧酸[CX3](O)[OX2H1]注意要指定羟基氢酯基[CX3](O)[OX2][#6]氧另一端连碳酰胺[CX3](O)[NX3]氮可以是伯胺、仲胺、叔胺磺酸基[SX4](O)(O)[OX2H1]硫连两个双键氧和一个羟基硝基[NX3](O)O氮连两个双键氧芳香环c1ccccc1苯环更通用的是c1ccccc1或[#6;R]1:[#6;R]:[#6;R]:[#6;R]:[#6;R]:[#6;R]:1这些 SMARTS 我实测过大部分但要注意不同的 RDKit 版本对芳香性的处理可能有细微差异同一个 SMARTS 在不同版本上匹配结果可能不完全一致。如果你要做严格的筛选建议在目标版本上先做一轮验证。6.3 匹配结果的可视化统计当你处理完一批分子后通常需要统计一下匹配情况。比如“这批分子里有多少含羧酸多少含酰胺”。这个用 Python 的collections.Counter就能做from collections import Counter results Counter() for mol in mols: if mol.HasSubstructMatch(carboxylic_acid): results[carboxylic_acid] 1 if mol.HasSubstructMatch(amide): results[amide] 1 print(results)如果你想把统计结果可视化可以用 matplotlib 画个柱状图。但要注意RDKit 本身不提供统计图表功能得靠外部库。我一般用 matplotlib 或者 seaborn简单画个 bar chart 就够了。还有一个技巧是把匹配结果导出成 CSV方便后续在 Excel 里做透视表。RDKit 的PandasTools模块可以直接把分子和属性转成 DataFrame配合 pandas 做统计分析非常方便。7. 我个人的几条实操建议第一SMARTS 写完一定要验证。我习惯拿几个已知含目标片段的分子和几个已知不含的分子做正负对照确认 SMARTS 的匹配行为符合预期。这一步花不了几分钟但能避免后面批量筛选时出现系统性错误。第二高亮图不要追求一次画完所有匹配。如果一个分子里有多个匹配全部高亮会让图变得难以辨认。我通常只高亮第一个匹配或者用不同颜色高亮前两三个剩下的在文字里说明。第三批量处理时注意内存。RDKit 的分子对象虽然不算大但如果你一次加载几万个分子并全部保留在内存里还是会吃不少资源。我一般用生成器逐行读 SDF 或 SMILES 文件处理完一个就释放一个避免内存堆积。第四绘图尺寸要跟用途匹配。如果是放进论文或报告用size(600, 450)或更大保证打印清晰如果只是在终端里快速看一眼size(300, 200)就够了。SVG 格式适合需要缩放的场景PNG 适合直接嵌入文档。第五保留中间结果。子结构搜索的结果匹配的原子索引建议保存下来不要每次重新算。尤其是当你的 SMARTS 比较复杂、分子比较大的时候重新搜索可能很耗时。把匹配结果存成 JSON 或 pickle下次直接读能省不少时间。这套流程我在多个项目里反复用过从几十个分子的手工验证到几万个分子的批量筛选基本都能覆盖。RDKit 的子结构搜索和高亮绘制虽然只是它众多功能中的一小块但把这一块吃透能解决化学信息学里相当一部分日常需求。
