做化学信息学的人几乎每天都在和分子结构打交道。化合物库去重、苗头化合物筛选、骨架聚类、ADMET性质过滤、药效团定义……这些任务看着方向完全不同最后却常常落到同一个问题这个分子里到底有没有我关心的片段在RDKit里回答这个问题最直接的工具就是SMARTS。SMARTS是Daylight在SMILES基础上发展出来的一套子结构查询语言用方括号描述原子属性、用单字符描述键类型、用逻辑运算符组合条件RDKit负责把这段模式文本解析成查询对象再通过子图同构匹配算法在目标分子里找出所有符合条件的原子集合。这篇文章我打算结合自己做化合物过滤和药效团建模的实际经历把RDKit对SMARTS的支持范围、语法细节、扩展方式和踩过的坑一次性讲清楚。内容包括RDKit到底支持哪些SMARTS写法、和Daylight规范有哪些出入、怎么从零写实用的子结构模式以及如何借助递归SMARTS、反应SMARTS和自定义逻辑把它的表达能力再往前推一步。刚接触RDKit的同学可以把它当入门地图已经写了不少SMARTS却总被匹配结果吓一跳的老手也能在这找到一些值得对照的细节。1. SMARTS在RDKit里的定位子结构查询的通用语言1.1 一张照片和一份拼图SMARTS到底在描述什么先说一个容易混淆的点SMILES和SMARTS根本是两种东西。SMILES描述的是一个确定分子的完整连接表比如CC(O)O就是醋酸原子数、键级、电荷全都写死了。SMARTS描述的是一类分子片段是一种模式。打个比方SMILES像一张通缉犯照片SMARTS像一份拼图描述这里有个戴眼镜的、穿深色外套的、身高大约一米七的人。照片只能和照片比是不是同一个人拼图却能在人群里把符合条件的嫌疑人全部捞出来。这个差异直接决定了RDKit里的函数设计。Chem.MolFromSmiles()构建的是一个具体分子对象后面跟着的是完整的分子图Chem.MolFromSmarts()构建的是一个查询对象它不在乎目标分子里有多少个匹配原子只关心有没有和在哪。这也是为什么[C]在SMARTS里能匹配所有脂肪碳而C在SMILES里只能代表带一个特定连接环境的碳原子——语法符号相同语义完全不同。RDKit的定位很明确它是化学信息学的瑞士军刀而SMARTS就是刀上最常用的那个刀头。凡是涉及按结构特征筛分子的逻辑最后几乎都能翻译成一个或一组SMARTS查询。1.2 四个最常用的API入口RDKit里和SMARTS相关的核心API其实不多我日常用到的主要就四个from rdkit import Chem # 1. 把SMARTS模式解析成查询对象 pattern Chem.MolFromSmarts(c1ccccc1) # 2. 判断分子是否包含该模式 mol Chem.MolFromSmiles(CC(O)Oc1ccccc1C(O)O) print(mol.HasSubstructMatch(pattern)) # True # 3. 拿到具体的匹配原子索引 matches mol.GetSubstructMatches(pattern) print(matches) # 比如 [(7, 8, 9, 12, 10, 11)]表示苯环上的6个原子 # 4. 批量场景下用子结构库 from rdkit.Chem import SubstructLibrary lib SubstructLibrary.SubstructLibrary() lib.AddMol(mol) print(lib.GetMatches(pattern))HasSubstructMatch适合做布尔过滤写条件判断时最顺手GetSubstructMatches适合需要知道具体原子编号的场景比如后面要画高亮图、要对匹配片段做替换或计算。SubstructLibrary则是处理大批量分子的利器它内部会构建索引匹配速度比逐分子调用HasSubstructMatch快得多后面第五节我会详细说。再举个真实业务里常见的例子。很多人做化合物库清洗时会用PAINS规则筛掉疑似假阳性片段RDKit的FilterCatalog直接内置了PAINS的SMARTS集合from rdkit.Chem import FilterCatalog params FilterCatalog.FilterCatalogParams() params.AddCatalog(FilterCatalog.FilterCatalogParams.FilterCatalogs.PAINS) catalog FilterCatalog.FilterCatalog(params) bad_mol Chem.MolFromSmiles(OC1NC(O)NC1O) print(catalog.HasMatch(bad_mol)) # True命中PAINS规则这条链路的底层就是成百上千条SMARTS模式在跑子结构匹配。所以别把SMARTS当成毕业论文里凑数的小技巧它是整个分子筛选体系的地基。2. RDKit支持的SMARTS语法全景与Daylight规范差异2.1 原子查询属性速查SMARTS的核心语法是方括号里的原子表达式。RDKit对Daylight规范里的原子属性支持得比较完整我把常用写法整理成了一张速查表写模式时对照着看基本够用写法含义示例场景[C]脂肪族碳匹配饱和碳原子[c]芳香族碳匹配苯环、吡啶环上的碳[#6]原子序数6任意碳不管芳香还是脂肪都匹配[N,n,O,S,P,F,Cl,Br,I]对应元素常规元素匹配[D2]显式连接度等于2限定取代程度[H1]总氢数等于1限定CH、NH等[v3]价态等于3限定价电子状态[X4]连接度含氢等于4限定杂化环境[R2]处于2个环中匹配桥头原子[r5]最小环尺寸为5匹配五元环[1]/[-1]形式电荷匹配带电原子[13C]同位素标记匹配特定同位素[C]/[C]手性方向匹配光学异构环境这里面最常踩的坑是[C]和[#6]的区别。[C]只匹配脂肪碳[#6]匹配所有碳。如果你用[C]去查一个含吡啶环的分子吡啶环上的5个碳是不会被匹配到的因为这些碳在RDKit的芳香性模型里属于芳香碳必须写成[c]才对。反过来如果你确实想不管芳香脂肪只要是碳就行那就应该写[#6]。2.2 键查询与逻辑组合原子条件搞定之后键的条件决定了片段的空间连接方式。RDKit支持的键查询写法如下写法含义-单键双键#三键:芳香键~任意键/\方向性键顺反异构环键注意和手性原子的区分键也可以用逻辑符号组合比如[C],:[O]表示碳和氧之间是双键或芳香键。但真正让SMARTS有表达力的是原子内部的逻辑运算符。逻辑运算有四类!表示否定,表示或表示与;也表示与但优先级比高。注意这个优先级陷阱[C;R;!r5]的意思是处于环中且不在五元环中的脂肪碳而[CR!r5]这种写法其实是非法或者容易产生歧义的。我的经验是条件多的时候统一用分号隔开别混用和;可读性和准确性都更好。配合递归SMARTS还能写更复杂的条件比如不是羰基相连接的伯胺primary_amine Chem.MolFromSmarts([NX3;H2;!$(NCO)])这里的$()就是递归SMARTS后面专门讲。2.3 RDKit与Daylight规范的差异RDKit对Daylight SMARTS的兼容度很高但高不等于一模一样。我实测下来有三处需要特别注意。第一是芳香性模型不同。Daylight最早用的是1988年定义的芳香性规则RDKit实现的是增强模型对稠环、大环、含杂原子的芳香体系判断会有出入。同一个芳香SMARTS在两个体系里可能匹配出不同结果。所以如果你在网上看到一份特别老的SMARTS模式别直接抄先在RDKit里拿真实分子验证一遍。第二是解析失败的返回值。RDKit的MolFromSmarts遇到非法SMARTS时返回None同时往stderr丢一条warning而不是抛异常。这在批量处理时很容易被忽略某个模式解析失败返回了None你后面照常调用HasSubstructMatch就直接崩了。我自己写过不少类似的线上筛选脚本现在凡是涉及外部传入SMARTS的代码都会先加一层判空。第三是氢原子处理。Daylight规范里氢原子和氢原子数语义在某些Edge case下比较模糊RDKit把[H]显式氢、h隐式氢、H总氢数拆得很清楚。写到后面你就会发现RDKit的这套拆法其实更符合程序员直觉只是需要习惯。3. 从零手写实用SMARTS官能团、环系与药效团模式3.1 官能团三兄弟羰基、羧酸、酰胺光看语法不进实战记不住。我从最常见的官能团筛选开始展示怎么一步步把模式写准。先把三个经典模式摆出来# 酮羰基碳上连双键氧另外至少连一个碳 ketone Chem.MolFromSmarts([CX3](O)[#6]) # 羧酸羰基碳连OH carboxylic_acid Chem.MolFromSmarts([CX3](O)[OX2H1]) # 酰胺氮原子连羰基 amide Chem.MolFromSmarts([NX3][CX3](O))第一个需要解释的细节是为什么用[CX3]而不是[C]。羰基碳和氧之间是双键它的形式价态是4任何额外连接都会改变价态。写成[CX3]就把碳的显式连接数锁定为3这样醛、羧酸里的羰基碳都不会误入酮的匹配结果。这是SMARTS里最常用的价态锁定技巧。第二个细节是羧酸里的[OX2H1]。氧显式连接数为2、总氢数为1这正好是羧基OH的电子状态一个氢、一个连接到羰基碳。如果只写[OX2]酯的醚氧O连接数也是2就会被误匹配。所以氢数这个条件是区分羟基和醚键、氨基和酰胺键的关键。我是怎么发现这个问题的有一次我用[CX3](O)[OX2]去筛化合物库里的羧酸结果酯类化合物全进来了查了半天才意识到是氢数条件没写全。这类错误SMARTS本身没有语法问题解析不会报警只能靠测试集合来兜底后面第六节会细说。3.2 芳香环与杂环小写字母的规则芳香环是SMARTS里另一个高频场景。苯环的写法是c1ccccc1这里c是芳香碳数字1是环闭合标记。benzene Chem.MolFromSmarts(c1ccccc1) cyclohexane Chem.MolFromSmiles(C1CCCCC1) print(benzene.HasSubstructMatch(cyclohexane)) # False环己烷不是芳香环芳香杂环就得考虑氮的位置和氢数了。吡啶环和吡咯环是特别容易写错的两个例子pyridine Chem.MolFromSmarts(n1ccccc1) pyrrole Chem.MolFromSmarts([nH]1cccc1)吡啶的氮是芳香氮、不带氢直接写n就行。吡咯的氮虽然在芳香环里但它带一个氢如果不写[nH]RDKit的匹配引擎会认为吡咯环上那个氮是不存在的。这是因为吡咯的芳香性依赖氮上的孤对电子参与共轭氮的氢是明确存在于结构里的。我记得有次做激酶抑制剂的骨架分析要统计含有吡咯并嘧啶结构的分子。一开始用的模式是n1cccc1结果漏掉了大量含吡咯环的化合物补上[nH]之后才恢复正常。小写字母不够还要注意芳香原子的氢状态这是新手最容易漏的一环。3.3 带环境的复杂模式氢键供体与受体的经典写法药物化学里做药效团最常用的就是氢键供体HBD和氢键受体HBA定义。这两个定义用SMARTS写起来很有代表性因为要把环境也揉进去# 氢键供体带H的N/O/S以及带H的芳香氮 hbd Chem.MolFromSmarts( [$([N;!H0;v3,v41]),$([O,S;H1;v2]),$([n;H1;0])] ) # 氢键受体排除酚类环境后的O/S/N hba Chem.MolFromSmarts( [$([O,S;H1;v2;!$(*-*a)]),$([O,S;H0;v2]),$([O,S;-]),$([N;v3;!$(N-**)])] )这两条模式在RDKit的药效团代码里实际用过写法很经典。它们用到了三个技巧第一个是$()递归SMARTS它允许你在一对方括号里嵌入一个完整模式作为该原子的附加条件。[$([N;!H0;v3,v41])]的意思是这个氮原子满足氢数不等于0且价态为3或价态为4且带正电。注意这里同时用到了逗号的或和分号的与。第二个技巧是!$(*-*a)这种排除法。它的意思是不能连接到一个芳香原子上。为什么会出现在氢键受体里因为酚羟基的氧连着苯环它的配位环境和脂肪醇差很多在特定pH条件下不太适合当氢键受体所以先排除掉。这里的a匹配任意芳香原子。第三个技巧是把同一个原子拆成多条候选分支每条分支用一个递归条件描述一种环境。单个SMARTS写不出来既能是X状态又能是Y状态的复杂条件时分成多条并列分支再用中括号外部空格隔开这种写法往往能解决问题。4. 超越标准SMARTS递归查询与RDKit的扩展手段4.1 递归SMARTS$()的核心价值标准的原子表达式只能表达这个原子本身是什么样一旦需要表达这个原子周围的环境是什么样就要靠递归SMARTS。$()的语法是从Daylight规范里继承下来的允许你在一个原子条件内部嵌套另一个完整的SMARTS模式。它的意义有多大举个例子你要筛选非酰胺的伯胺也就是NH2不能直接连在羰基碳上# 错误写法只写伯胺会把酰胺的氮也筛进来 naive Chem.MolFromSmarts([NX3;H2]) # 正确写法排除羰基碳直接相连的情况 good Chem.MolFromSmarts([NX3;H2;!$(NCO)])!$(NCO)的意思是不存在这样的连接这个氮连到一个碳上该碳又双键连接一个氧。这个条件只靠原子属性是写不出来的因为你要约束的是键的远端环境。递归SMARTS相当于是把整个子结构匹配能力嵌到了原子查询里这也是SMARTS比正则表达式在化学问题上强大得多的根本原因。4.2 反向生成与自动扩展MolToSmarts和FindMCS很多人以为SMARTS只能人肉手写其实RDKit提供了反向工具能把一个具体分子片段自动转成SMARTSfrag Chem.MolFromSmiles(CC(O)Nc1ccccc1) smarts_str Chem.MolToSmarts(frag) print(smarts_str) # 会输出类似 [CH3][C]([O])[NH][c]1[c][c][c][c][c]1 的SMARTSMolToSmarts生成的模式通常比手写的更啰嗦因为它会把每个原子的连接度、氢数都显式带上。但它有个大用处当你手里有某个已知活性片段想快速得到一个可用的查询模式时先转成SMARTS再手动删掉多余条件比自己对着结构数原子快得多。更进阶的扩展场景是利用FindMCS从一堆分子里自动找公共子结构然后拿MCS的SMARTS当模式from rdkit.Chem import rdFMCS mol1 Chem.MolFromSmiles(CC(O)Nc1ccccc1) mol2 Chem.MolFromSmiles(CCC(O)Nc1ccccc1Cl) res rdFMCS.FindMCS([mol1, mol2], timeout20) print(res.smartsString) # 输出类似 C(O)Nc1ccccc1这条链路很有实用价值尤其是在做系列化合物分析时先跑一遍MCS拿到公共核心把它转成SMARTS然后到更大的化合物库里做子结构搜索找出所有包含该核心的分子。整个过程不需要手写一行模式却能得到一个性能可靠、语义正确的SMARTS查询。4.3 从子结构匹配到反应规则SMIRKSSMARTS的另一个扩展方向是反应SMARTS也就是SMIRKS。它把查询模式升级成了反应规则格式上依然沿用SMARTS语法只是用原子映射编号:1、:2来标记反应前后的对应关系from rdkit.Chem import rdChemReactions rxn rdChemReactions.ReactionFromSmarts( [C:1]([O:2])[O:3][C:4][C:1]([O:2])[O-:3].[C:4] ) reactants Chem.MolFromSmiles(CC(O)OCC) products rxn.RunReactants((reactants,)) print(len(products[0]))这个例子演示的是酯键断裂的反应羰基碳C:1、双键氧O:2、醚氧O:3、烷基碳C:4箭头右边表示断键后生成羧酸根和烷基正离子。原子映射编号是SMIRKS的灵魂没有它们反应引擎根本不知道断裂后哪个原子该去哪个产物。RDKit对SMIRKS的支持相当完整可以用来做产物枚举、代谢预测、反应数据库构建。我在实际项目中常拿它做组合化学的虚拟筛选写一个核心反应规则喂入一堆不同取代基的底物自动生成候选产物库然后再用SMARTS对新产物做子结构过滤。SMARTS负责筛SMIRKS负责生成配合起来非常顺。4.4 两阶段扩展SMARTS粗筛加Python细筛SMARTS再强大终究有表达不了的东西。比如你要筛选侧链含一个带正电荷的氮同时这个氮所在的取代基总碳数不超过5SMARTS的原子级约束写不了这种跨多个原子的算术条件。我的做法是两阶段策略rough Chem.MolFromSmarts([NX3;!H0]) # 第一阶段粗筛 for mol in molecule_list: matches mol.GetSubstructMatches(rough) for hit in matches: atom mol.GetAtomWithIdx(hit[0]) # 第二阶段Python细筛比如沿侧链数碳原子 # if count_carbons_along_chain(atom) 5: # accept.append(mol)粗筛用SMARTS保证速度细筛用Python保证灵活性。这样既不用写出天书一般的超级模式也不会因为正则式的局限漏掉真实命中。化学信息学项目做到后期几乎都会演变成SMARTS RDKit对象遍历 业务规则三层结构SMARTS只是第一层。5. 匹配引擎实战避坑芳香性、氢原子与性能5.1 芳香性判定差异导致的误匹配RDKit对芳香性的判断是整套匹配逻辑里最容易出意外的地方。我在第三节已经提过[C]和[c]的区别这里再深入说一个真实案例。有一次我要筛选含萘环的分子写的是c1ccc2ccccc2c1这个标准的萘SMARTS。所有含萘环的分子都能匹配到这个没问题。后来我拿它去跑一个包含几万个稠环化合物的库发现有一批茚满类化合物也被匹配到了。查了半天才发现茚满的六元环虽然名义上是饱和环但在RDKit的增强芳香性模型里某些融合环境下环会呈现部分芳香特征原先我以为反正不是芳香环的分子局部环确实被判定成了芳香。这件事给我的教训是不要凭直觉判断某段环是不是芳香要在代码里显式打印GetIsAromatic()验证。特别是做稠环体系时先画高亮图、先打印原子芳香标记再决定SMARTS该用大写还是小写。5.2 显式氢与隐式氢的纠缠RDKit的分子对象默认不保留氢原子氢是隐式的。这本来没什么问题但当你调用Chem.AddHs()把氢变成显式原子之后SMARTS的行为就会变。mol Chem.MolFromSmiles(CCO) pattern Chem.MolFromSmarts([#6][#6][#6]) print(mol.GetSubstructMatches(pattern)) # 正常脂肪醇能匹配到 mol_h Chem.AddHs(mol) # 此时模式 [*] 有可能匹配到显式氢原子导致额外命中如果你用了[*]或[a]这类宽泛匹配显式氢加入后氢原子本身也可能成为匹配的一部分因为*匹配任意原子。避免方法很简单要么在用AddHs之前完成子结构匹配要么在SMARTS里显式排除氢写成[!#1]。另外[H]这个写法在RDKit里有明确语义。如果你要查分子里是否存在某个特定位置的氢可以写[C][H]这种模式但要确保分子对象确实带显式氢否则永远匹配不到。5.3 大数据集上的性能调优经验SMARTS匹配的本质是子图同构算法最坏情况是指数复杂度。批量筛选上百万个分子时写得不讲究的SMARTS能把任务拖到不可接受。我的几条性能经验按性价比排序尽量用精确原子属性少用*。*会让匹配分支暴增能写[#6]就别写[*]能写[CX3]就别写[C]。优先匹配稀有原子。比如模式里含Br就先把Br匹配上再扩展周围原子匹配空间会小很多。用GetSubstructMatch代替GetSubstructMatches。只要判断是否有匹配前者找到第一个命中就返回速度快得多。大批量多模式场景用SubstructLibrary。它内部构建了分子指纹索引先用指纹做粗筛再做子图匹配整体吞吐量能提升一个数量级。解析SMARTS对象做复用。不要在循环里反复调Chem.MolFromSmarts那等于每次迭代都重新解析一遍模式把解析对象放到循环外面只对分子做匹配。from rdkit import Chem from rdkit.Chem import SubstructLibrary lib SubstructLibrary.SubstructLibrary() patterns [Chem.MolFromSmarts(c1ccccc1), Chem.MolFromSmarts([CX3]([OX1]))] for smi in huge_smiles_list: m Chem.MolFromSmiles(smi) if m is not None: lib.AddMol(m) for p in patterns: hits lib.GetMatches(p) # 处理命中结果这套组合拳救过我好几次。之前有个化合物库有两百多万个分子要跑二十几条PAINS类规则一开始直接嵌套循环跑预计要几个小时换成SubstructLibrary加模式复用之后十几分钟就出结果了。6. 验证与调试怎么确定你的SMARTS写得对6.1 用可视化确认匹配结果SMARTS写错不会报错只会默默返回错误结果。所以验证环节特别重要。我最快的验证方式是把匹配结果高亮画出来from rdkit.Chem import Draw from rdkit.Chem.Draw import IPythonConsole mol Chem.MolFromSmiles(CC(O)Oc1ccccc1C(O)O) pattern Chem.MolFromSmarts([CX3](O)[OX2H1]) matches mol.GetSubstructMatches(pattern) for hit in matches: img Draw.MolToImage(mol, highlightAtomslist(hit)) # Jupyter里直接显示img即可肉眼扫一遍高亮图比看一百行匹配日志更快。如果高亮部分不是你想查的那个片段说明SMARTS写偏了立刻就能发现。GetSubstructMatches返回的元组顺序不是固定的不要依赖顺序做后续的原子类型判断应该通过GetAtomWithIdx逐个取出来看属性。6.2 用正反例构建测试集所有需要长期维护的SMARTS模式我都建议配一个正反例测试集。做法很朴素准备若干分子标注预期是否命中然后写成断言test_cases { CC(O)O: True, # 乙酸应命中羧酸 CC(O)OC: True, # 乙酸甲酯含有羧酸酯片段 CCO: False, # 乙醇不应命中醛 CC(C)O: False, # 丙酮不含OH } pattern Chem.MolFromSmarts([CX3](O)[OX2H1]) for smi, expected in test_cases.items(): mol Chem.MolFromSmiles(smi) result mol.HasSubstructMatch(pattern) assert result expected, f{smi}: 期望{expected}实际{result}真实业务里很多模式是看起来对边界错误的。比如上面测试集里如果我把乙酸甲酯的预期标成False跑一下就能发现酯的醚氧并不会命中[OX2H1]因为酯氧没有氢。这种回归测试集在模式迭代时尤其有用——改了一条规则跑一遍历史用例马上知道有没有改坏。6.3 解析失败与空匹配的处理最后再强调一次MolFromSmarts返回None的问题。批量脚本里每来一个SMARTS就要判一次空def safe_smarts(smarts_str): pattern Chem.MolFromSmarts(smarts_str) if pattern is None: raise ValueError(f非法SMARTS: {smarts_str}) return pattern我还遇到过一种更隐蔽的情况SMARTS解析成功但某个分子怎么都匹配不上。这时候先别怀疑算法先确认两件事第一分子有没有成功构建MolFromSmiles如果返回了None后面所有匹配都会出问题第二分子的芳香性标记是否正常有时候从文件读进来的分子没做标准化芳香环还是单双键交替的状态小写字母的SMARTS就永远匹配不上。碰到这种问题我一般会用Chem.SanitizeMol(mol)强制走一遍标准化流程再重新匹配。这个操作能解决八成莫名其妙匹配不上的诡异现象。写SMARTS这件事说到底是用精确约束换精确结果。我做了两年化学信息学最深的体会是好的SMARTS模式不是一次写出来的而是一轮一轮拿真实数据磨出来的。每当你往模式里多加一个条件、多排除一种环境你对RDKit匹配引擎的理解就会深一层。刚开始我也会被[C]和[#6]、[nH]和n这些细枝末节气得头疼但踩过几次坑、建立起来正反例测试集之后SMARTS慢慢就从一个需要反复试错的工具变成了随手就能用的表达语言。最后再分享一个小技巧如果你在一个大项目里反复用同一组模式建议把它们集中放到一个patterns.py文件里每个模式配一句注释说明用途和使用条件。模式这种东西三个月后再看连自己都可能忘了当时为什么要加某个排除条件。留好注释和测试集等于给未来的自己省了一整天的排查时间。