1. 從一張畫(huà)滿(mǎn)熒光筆的分子圖說(shuō)起做化學(xué)信息學(xué)的人遲早會(huì)碰到一個(gè)繞不開(kāi)的需求手里有一堆分子想快速知道某個(gè)片段到底出現(xiàn)在哪些結(jié)構(gòu)里而且最好能一眼看出來(lái)它長(zhǎng)在哪個(gè)位置。我第一次遇到這個(gè)場(chǎng)景是要從一批候選化合物里篩出所有含苯并咪唑骨架的分子當(dāng)時(shí)傻乎乎地一個(gè)個(gè)用肉眼比對(duì)幾十個(gè)分子下來(lái)眼睛都快瞎了還漏掉了好幾個(gè)。后來(lái)才意識(shí)到這種活兒本來(lái)就該交給 RDKit 去做——它不僅能幫你把匹配的子結(jié)構(gòu)找出來(lái)還能直接在分子圖上用高亮標(biāo)出來(lái)省時(shí)省力還不容易出錯(cuò)。這篇內(nèi)容就是圍繞“搜索并高亮展示分子中匹配的子結(jié)構(gòu)”這件事展開(kāi)的。核心工具是 RDKit一個(gè)在化學(xué)信息學(xué)領(lǐng)域被廣泛使用的開(kāi)源工具包。它提供了子結(jié)構(gòu)搜索Substructure Search和高亮繪制Highlight Rendering兩大能力前者負(fù)責(zé)判斷“這個(gè)片段在不在分子里”后者負(fù)責(zé)把“在哪兒”可視化出來(lái)。適合的讀者包括做虛擬篩選的藥物化學(xué)研究者、處理化合物庫(kù)的數(shù)據(jù)科學(xué)家、以及任何需要批量分析分子結(jié)構(gòu)的人。哪怕你之前只用 RDKit 畫(huà)過(guò)簡(jiǎn)單的分子式這篇內(nèi)容也能讓你把子結(jié)構(gòu)搜索和高亮這套組合拳打順。我下面會(huì)從最基礎(chǔ)的匹配原理講起然后一步步拆到代碼層面把 SMARTS 的寫(xiě)法、匹配函數(shù)的區(qū)別、高亮參數(shù)的調(diào)優(yōu)、以及批量處理時(shí)的性能坑都過(guò)一遍。中間會(huì)穿插我自己踩過(guò)的坑和實(shí)測(cè)有效的技巧盡量讓你少走彎路。2. 子結(jié)構(gòu)搜索的底層邏輯RDKit 到底在比什么2.1 圖匹配而非字符串匹配很多人第一次聽(tīng)到“子結(jié)構(gòu)搜索”下意識(shí)會(huì)以為是在做字符串包含判斷就像 Python 里abc in abcdef那樣。但分子結(jié)構(gòu)的本質(zhì)是圖Graph原子是節(jié)點(diǎn)化學(xué)鍵是邊。RDKit 做子結(jié)構(gòu)搜索實(shí)際上是在做**子圖同構(gòu)Subgraph Isomorphism**的判斷——問(wèn)的是“查詢(xún)片段這個(gè)圖能不能嵌入到目標(biāo)分子這個(gè)圖里”。這個(gè)區(qū)別很關(guān)鍵。字符串匹配是線性的、位置敏感的而圖匹配要考慮原子的連接關(guān)系、環(huán)的閉合方式、以及化學(xué)環(huán)境。舉個(gè)例子查詢(xún)片段是“苯環(huán)”目標(biāo)分子里有一個(gè)苯環(huán)但苯環(huán)上還連著其他基團(tuán)。字符串匹配可能會(huì)因?yàn)椤岸嗔俗址倍〉珗D匹配只關(guān)心苯環(huán)這個(gè)子圖是否存在不關(guān)心它外面還掛了什么。這就是為什么子結(jié)構(gòu)搜索必須用圖算法而不能用簡(jiǎn)單的文本比對(duì)。RDKit 底層用的是 VF2 算法的變體來(lái)做子圖同構(gòu)判斷。VF2 是一種經(jīng)典的回溯搜索算法通過(guò)逐步擴(kuò)展匹配節(jié)點(diǎn)、剪枝不滿(mǎn)足條件的候選最終找到所有可能的嵌入方式。你不需要自己實(shí)現(xiàn)這個(gè)算法但理解它的存在有助于你明白為什么有些查詢(xún)會(huì)返回多個(gè)匹配結(jié)果因?yàn)橥粋€(gè)片段可能在分子里有多種嵌入方式以及為什么某些復(fù)雜查詢(xún)會(huì)變慢因?yàn)榛厮菟阉鞯臓顟B(tài)空間可能很大。2.2 SMARTS描述子結(jié)構(gòu)的專(zhuān)用語(yǔ)言要讓 RDKit 知道你想找什么片段你得用一種它能讀懂的語(yǔ)言來(lái)描述。這個(gè)語(yǔ)言就是SMARTSSMiles ARbitrary Target Specification。SMARTS 是 SMILES 的擴(kuò)展SMILES 用來(lái)描述一個(gè)具體分子而 SMARTS 用來(lái)描述一類(lèi)結(jié)構(gòu)模式。舉個(gè)最直觀的例子。苯的 SMILES 是c1ccccc1這是一個(gè)具體的分子。如果你想表達(dá)“任意一個(gè)苯環(huán)”SMARTS 可以寫(xiě)成c1ccccc1但更通用的寫(xiě)法是[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1其中[#6]表示任意碳原子:表示芳香鍵。當(dāng)然實(shí)際使用中很少有人寫(xiě)這么啰嗦通常直接用c1ccccc1就夠了因?yàn)樾?xiě)c在 SMARTS 里就代表芳香碳。SMARTS 的強(qiáng)大之處在于它支持各種通配符和邏輯操作。比如[#6]匹配任意碳[C,N]匹配碳或氮[!O]匹配非氧原子*匹配任意原子。還有針對(duì)化學(xué)環(huán)境的修飾符比如[C;H3]表示恰好連了三個(gè)氫的碳[C;R]表示處于環(huán)上的碳。這些語(yǔ)法讓你能精確地描述“我要找的是什么”。我個(gè)人的經(jīng)驗(yàn)是剛開(kāi)始用 SMARTS 時(shí)最容易犯的錯(cuò)是把 SMILES 當(dāng) SMARTS 用。比如你想找羧基寫(xiě)了C(O)O這在很多情況下確實(shí)能匹配到但它也會(huì)匹配到酯基里的C(O)O部分因?yàn)轷セ舶@個(gè)子圖。如果你只想找游離的羧酸得寫(xiě)成[CX3](O)[OX2H1]明確指定羥基上的氧連了一個(gè)氫。這種細(xì)節(jié)在批量篩選時(shí)特別重要寫(xiě)錯(cuò)一個(gè)修飾符結(jié)果可能差出幾百個(gè)分子。2.3 匹配函數(shù)的家族HasSubstructMatch 與 GetSubstructMatchesRDKit 提供了兩個(gè)層次的子結(jié)構(gòu)搜索函數(shù)用途不同返回值也不同。HasSubstructMatch返回一個(gè)布爾值只告訴你“有沒(méi)有匹配”不告訴你“匹配在哪里”。這個(gè)函數(shù)適合做快速篩選比如你有一個(gè)十萬(wàn)級(jí)的化合物庫(kù)想先粗篩出所有含某個(gè)片段的分子用這個(gè)函數(shù)就夠了速度快內(nèi)存占用小。GetSubstructMatches返回一個(gè)元組的元組每個(gè)內(nèi)層元組是一組原子索引表示一次匹配中查詢(xún)片段的各個(gè)原子分別對(duì)應(yīng)目標(biāo)分子的哪些原子。這個(gè)函數(shù)適合需要精確定位的場(chǎng)景比如你要在高亮圖上標(biāo)出匹配位置就必須知道具體是哪些原子參與了匹配。還有一個(gè)GetSubstructMatch注意沒(méi)有復(fù)數(shù) s只返回第一個(gè)匹配的原子索引元組。如果你只關(guān)心“找到一個(gè)就行”用這個(gè)比GetSubstructMatches更省資源因?yàn)樗业降谝粋€(gè)匹配就停了不會(huì)繼續(xù)搜索所有可能的嵌入方式。選擇哪個(gè)函數(shù)取決于你的下游需求。如果只是統(tǒng)計(jì)有多少分子含有某片段HasSubstructMatch足夠如果要畫(huà)高亮圖必須用GetSubstructMatches或GetSubstructMatch。我見(jiàn)過(guò)有人用HasSubstructMatch篩完之后又對(duì)每個(gè)陽(yáng)性分子調(diào)GetSubstructMatches來(lái)畫(huà)圖這其實(shí)多了一次搜索不如直接調(diào)后者然后判斷返回是否為空。3. 高亮繪制的參數(shù)調(diào)優(yōu)讓匹配位置一眼可見(jiàn)3.1 高亮原理原子索引到顏色的映射RDKit 的繪制模塊rdkit.Chem.Draw提供了Draw.MolToImage和Draw.MolToFile等函數(shù)它們都接受一個(gè)highlightAtoms參數(shù)。這個(gè)參數(shù)是一個(gè)原子索引的序列RDKit 會(huì)把這些原子以及它們之間的鍵用高亮顏色畫(huà)出來(lái)。關(guān)鍵點(diǎn)在于高亮的是原子鍵的高亮是附帶的。當(dāng)你傳入一組原子索引時(shí)RDKit 會(huì)檢查這些原子之間的鍵如果一條鍵的兩個(gè)端點(diǎn)都在高亮集合里這條鍵也會(huì)被高亮。所以如果你只高亮了一個(gè)原子它周?chē)逆I不會(huì)變色如果你高亮了整個(gè)苯環(huán)的六個(gè)碳環(huán)上的六條鍵也會(huì)跟著高亮。這個(gè)機(jī)制意味著GetSubstructMatches返回的原子索引元組可以直接喂給highlightAtoms不需要額外轉(zhuǎn)換。比如匹配結(jié)果返回(0, 1, 2, 3, 4, 5)你直接把這個(gè)元組傳給highlightAtoms苯環(huán)就會(huì)被完整高亮。3.2 顏色與樣式的自定義默認(rèn)的高亮顏色是黃色說(shuō)實(shí)話(huà)在白色背景上對(duì)比度一般尤其是當(dāng)分子本身有顏色或者打印出來(lái)的時(shí)候黃色高亮經(jīng)?？床磺?。RDKit 允許你通過(guò)highlightColor參數(shù)自定義顏色傳一個(gè) RGB 元組就行。我通常用(0.8, 0.2, 0.2)這種偏紅的顏色在白底上非常醒目。如果你需要同時(shí)高亮多個(gè)不同的子結(jié)構(gòu)比如一個(gè)分子里既有苯環(huán)又有羧基你想用不同顏色區(qū)分那就得用highlightAtomColors參數(shù)。這個(gè)參數(shù)接受一個(gè)字典鍵是原子索引值是 RGB 元組。你可以把第一個(gè)匹配的原子染成紅色第二個(gè)匹配的染成藍(lán)色以此類(lèi)推。還有一個(gè)容易被忽略的參數(shù)是highlightBonds它允許你單獨(dú)指定要高亮的鍵。通常情況下你不需要用它因?yàn)樵痈吡習(xí)詣?dòng)帶動(dòng)鍵高亮。但在某些特殊場(chǎng)景下比如你想高亮一條鍵但不想高亮它的端點(diǎn)原子就得手動(dòng)指定。不過(guò)這種需求很少見(jiàn)大部分時(shí)候用原子高亮就夠了。3.3 圖例與標(biāo)注的添加畫(huà)高亮圖時(shí)經(jīng)常需要在圖旁邊加個(gè)圖例說(shuō)明“紅色高亮的是苯環(huán)藍(lán)色高亮的是羧基”。RDKit 的Draw.MolToImage支持legend參數(shù)可以傳一個(gè)字符串作為圖例。如果你用Draw.MolsToGridImage畫(huà)多個(gè)分子的網(wǎng)格圖legends參數(shù)可以傳一個(gè)字符串列表每個(gè)分子對(duì)應(yīng)一個(gè)圖例。圖例的字體大小和位置 RDKit 沒(méi)有直接暴露參數(shù)來(lái)調(diào)但你可以通過(guò)Draw.DrawingOptions來(lái)全局調(diào)整。比如Draw.DrawingOptions.legendFontSize可以改圖例字號(hào)Draw.DrawingOptions.legendFraction可以改圖例占整張圖的高度比例。這些參數(shù)在批量出圖時(shí)很有用尤其是當(dāng)分子圖很小、圖例文字?jǐn)D在一起的時(shí)候。我自己的習(xí)慣是如果只是自己看不加圖例也行高亮顏色本身就說(shuō)明了問(wèn)題。但如果是給別人看或者放進(jìn)報(bào)告里圖例一定要加而且要寫(xiě)清楚 SMARTS 表達(dá)式這樣別人能復(fù)現(xiàn)你的篩選條件。4. 從單分子到批量完整實(shí)操鏈路4.1 環(huán)境準(zhǔn)備與依賴(lài)安裝RDKit 的安裝方式取決于你的操作系統(tǒng)和 Python 環(huán)境。最省事的方式是用 condaconda install -c conda-forge rdkit如果你用 pip也可以裝pip install rdkit但要注意pip 版本的 RDKit 在某些平臺(tái)上可能缺少一些可選依賴(lài)比如 Cairo 繪圖后端。如果你發(fā)現(xiàn)畫(huà)出來(lái)的圖是空白的或者報(bào)錯(cuò)說(shuō)找不到繪圖后端大概率是這個(gè)問(wèn)題。conda 版本通常更完整推薦優(yōu)先用 conda。安裝完之后驗(yàn)證一下from rdkit import Chem from rdkit.Chem import Draw print(Chem.__version__ if hasattr(Chem, __version__) else ok)如果沒(méi)報(bào)錯(cuò)說(shuō)明環(huán)境沒(méi)問(wèn)題。我建議在虛擬環(huán)境里裝 RDKit因?yàn)樗蕾?lài)的庫(kù)比較多跟其他科學(xué)計(jì)算庫(kù)混在一起容易出沖突。4.2 單分子搜索與高亮的完整代碼下面這段代碼是我平時(shí)最常用的模板從 SMILES 讀入分子搜索子結(jié)構(gòu)然后畫(huà)高亮圖from rdkit import Chem from rdkit.Chem import Draw # 目標(biāo)分子 mol Chem.MolFromSmiles(c1ccc2c(c1)nc(n2)c3ccccc3) # 查詢(xún)片段苯并咪唑骨架 query Chem.MolFromSmarts(c1ccc2c(c1)nc(n2)) # 搜索所有匹配 matches mol.GetSubstructMatches(query) # 畫(huà)高亮圖 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是一個(gè)元組的元組。如果分子里只有一個(gè)苯并咪唑matches就只有一個(gè)元素如果有多個(gè)就會(huì)有多個(gè)。我通常只高亮第一個(gè)匹配因?yàn)槎鄠€(gè)匹配同時(shí)高亮?xí)寛D變得很亂。如果你確實(shí)需要展示所有匹配可以用highlightAtoms傳一個(gè)扁平化的列表把所有匹配的原子都放進(jìn)去但顏色就沒(méi)法區(qū)分了。4.3 批量處理與網(wǎng)格圖輸出當(dāng)你有一批分子要處理時(shí)逐個(gè)畫(huà)圖再拼起來(lái)太麻煩。RDKit 的Draw.MolsToGridImage可以直接畫(huà)網(wǎng)格圖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)這里有個(gè)細(xì)節(jié)highlightAtomLists參數(shù)接受一個(gè)列表每個(gè)元素對(duì)應(yīng)一個(gè)分子的高亮原子索引。對(duì)于沒(méi)有匹配的分子傳空元組()就行RDKit 會(huì)正常畫(huà)出分子但不加高亮。legends參數(shù)同理每個(gè)分子對(duì)應(yīng)一個(gè)圖例字符串。實(shí)測(cè)下來(lái)MolsToGridImage在分子數(shù)量超過(guò) 50 個(gè)時(shí)生成的圖片會(huì)非常大加載和渲染都會(huì)變慢。如果你要處理幾百個(gè)分子建議分批出圖每批 20 到 30 個(gè)或者直接用MolsToFile逐個(gè)保存成單獨(dú)的圖片文件。5. 那些文檔里不會(huì)寫(xiě)的踩坑經(jīng)驗(yàn)5.1 芳香性感知導(dǎo)致的匹配失敗這是我最常遇到的坑沒(méi)有之一。RDKit 在讀入分子時(shí)會(huì)做芳香性感知Aromaticity Perception決定哪些環(huán)是芳香環(huán)。如果你的 SMARTS 里用了小寫(xiě)字母如c表示芳香碳但目標(biāo)分子在 RDKit 眼里那個(gè)環(huán)不是芳香環(huán)匹配就會(huì)失敗。舉個(gè)例子吡啶的 SMILES 是c1ccncc1RDKit 會(huì)正確識(shí)別為芳香環(huán)。但如果你從某些文件格式讀入分子或者分子經(jīng)過(guò)了某些化學(xué)修飾RDKit 可能沒(méi)有正確感知芳香性導(dǎo)致c1ccncc1這個(gè) SMARTS 匹配不上。解決辦法是用Chem.SanitizeMol強(qiáng)制做一次完整的化學(xué)感知或者在 SMARTS 里用大寫(xiě)字母加芳香性修飾符比如[C;R]表示環(huán)上的碳不依賴(lài)芳香性判斷。我個(gè)人的習(xí)慣是在批量處理前先對(duì)每個(gè)分子調(diào)一次Chem.SanitizeMol確保 RDKit 對(duì)分子的理解是一致的。雖然這會(huì)稍微增加處理時(shí)間但能避免很多莫名其妙的匹配失敗。5.2 匹配數(shù)量爆炸與性能陷阱有些 SMARTS 查詢(xún)會(huì)返回大量的匹配結(jié)果尤其是當(dāng)查詢(xún)片段很小或者很通用的時(shí)候。比如你用[#6]去匹配一個(gè)含 50 個(gè)碳的分子會(huì)返回 50 個(gè)匹配每個(gè)匹配就是一個(gè)碳原子。如果你用GetSubstructMatches并且沒(méi)有設(shè)置uniquifyFalseRDKit 會(huì)返回所有可能的嵌入方式數(shù)量可能是指數(shù)級(jí)的。GetSubstructMatches有一個(gè)uniquify參數(shù)默認(rèn)是True表示對(duì)結(jié)果去重。去重的規(guī)則是如果兩組匹配的原子集合相同只保留一組。但即使去重之后匹配數(shù)量仍然可能很大。還有一個(gè)maxMatches參數(shù)可以限制返回的最大匹配數(shù)默認(rèn)是 1000。如果你只關(guān)心“有沒(méi)有匹配”或者“前幾個(gè)匹配”把maxMatches設(shè)小一點(diǎn)能顯著提速。另外GetSubstructMatches還有一個(gè)useChirality參數(shù)默認(rèn)是False。如果你需要區(qū)分手性比如只想匹配特定構(gòu)型的片段得把它設(shè)成True。但開(kāi)啟手性匹配會(huì)變慢因?yàn)樗惴ㄐ枰~外比較原子的立體化學(xué)信息。5.3 高亮圖的原子索引錯(cuò)位這個(gè)問(wèn)題比較隱蔽通常出現(xiàn)在你從文件讀入分子之后。RDKit 在讀入分子時(shí)原子索引的順序可能跟你預(yù)期的不同。比如你從 SDF 文件讀入一個(gè)分子SDF 里原子的順序是 A、B、C但 RDKit 內(nèi)部可能會(huì)重新排列原子索引。如果你用外部工具生成的原子索引去高亮就會(huì)標(biāo)錯(cuò)位置。解決辦法是始終用 RDKit 自己返回的索引。GetSubstructMatches返回的索引就是 RDKit 內(nèi)部的原子索引直接傳給highlightAtoms不會(huì)錯(cuò)位。如果你需要跟外部數(shù)據(jù)對(duì)照先用mol.GetAtomWithIdx(i).GetSymbol()確認(rèn)一下索引對(duì)應(yīng)的原子是不是你想要的。還有一個(gè)相關(guān)的問(wèn)題是當(dāng)你對(duì)分子做了Chem.AddHs或者Chem.RemoveHs之后原子索引會(huì)變。加氫會(huì)插入新的氫原子索引全部后移去氫會(huì)刪除氫原子索引前移。如果你在高亮之前做了這些操作一定要重新獲取匹配結(jié)果不能用之前的索引。5.4 繪圖后端的兼容性問(wèn)題RDKit 的繪圖功能依賴(lài)后端默認(rèn)在 Jupyter 里用IPython.display顯示在腳本里用PIL或Cairo保存文件。如果你在無(wú)圖形界面的服務(wù)器上跑代碼可能會(huì)遇到ImportError或者OSError提示找不到顯示設(shè)備。解決辦法是明確指定用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是純文件輸出的不依賴(lài)圖形界面在服務(wù)器上跑很穩(wěn)。MolDraw2DSVG則是輸出 SVG 矢量圖適合放進(jìn)網(wǎng)頁(yè)或者需要無(wú)損縮放的場(chǎng)景。我一般在服務(wù)器上批量出圖時(shí)都用MolDraw2DCairo從來(lái)沒(méi)出過(guò)問(wèn)題。6. 進(jìn)階玩法讓子結(jié)構(gòu)搜索更貼合實(shí)際需求6.1 多片段聯(lián)合搜索有時(shí)候你要找的不是單個(gè)片段而是“同時(shí)含有 A 和 B 兩個(gè)片段”的分子。最直接的做法是分別搜索然后取交集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: # 兩個(gè)片段都存在 all_atoms list(match_a[0]) list(match_b[0]) img Draw.MolToImage(mol, highlightAtomsall_atoms)但這樣高亮出來(lái)兩個(gè)片段是同一種顏色分不清哪個(gè)是哪個(gè)。更好的做法是用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) # 藍(lán)色 img Draw.MolToImage(mol, highlightAtomslist(colors.keys()), highlightAtomColorscolors)這樣苯環(huán)是紅色羧基是藍(lán)色一目了然。如果你有更多片段可以繼續(xù)加顏色但建議不要超過(guò)三種否則圖會(huì)變得很花。6.2 用 SMARTS 做官能團(tuán)篩選子結(jié)構(gòu)搜索最常見(jiàn)的應(yīng)用場(chǎng)景之一就是官能團(tuán)篩選。RDKit 官方文檔里有一份常用的 SMARTS 列表我摘幾個(gè)我經(jīng)常用的官能團(tuán)SMARTS說(shuō)明羧酸[CX3](O)[OX2H1]注意要指定羥基氫酯基[CX3](O)[OX2][#6]氧另一端連碳酰胺[CX3](O)[NX3]氮可以是伯胺、仲胺、叔胺磺酸基[SX4](O)(O)[OX2H1]硫連兩個(gè)雙鍵氧和一個(gè)羥基硝基[NX3](O)O氮連兩個(gè)雙鍵氧芳香環(huán)c1ccccc1苯環(huán)更通用的是c1ccccc1或[#6;R]1:[#6;R]:[#6;R]:[#6;R]:[#6;R]:[#6;R]:1這些 SMARTS 我實(shí)測(cè)過(guò)大部分但要注意不同的 RDKit 版本對(duì)芳香性的處理可能有細(xì)微差異同一個(gè) SMARTS 在不同版本上匹配結(jié)果可能不完全一致。如果你要做嚴(yán)格的篩選建議在目標(biāo)版本上先做一輪驗(yàn)證。6.3 匹配結(jié)果的可視化統(tǒng)計(jì)當(dāng)你處理完一批分子后通常需要統(tǒng)計(jì)一下匹配情況。比如“這批分子里有多少含羧酸多少含酰胺”。這個(gè)用 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)如果你想把統(tǒng)計(jì)結(jié)果可視化可以用 matplotlib 畫(huà)個(gè)柱狀圖。但要注意RDKit 本身不提供統(tǒng)計(jì)圖表功能得靠外部庫(kù)。我一般用 matplotlib 或者 seaborn簡(jiǎn)單畫(huà)個(gè) bar chart 就夠了。還有一個(gè)技巧是把匹配結(jié)果導(dǎo)出成 CSV方便后續(xù)在 Excel 里做透視表。RDKit 的PandasTools模塊可以直接把分子和屬性轉(zhuǎn)成 DataFrame配合 pandas 做統(tǒng)計(jì)分析非常方便。7. 我個(gè)人的幾條實(shí)操建議第一SMARTS 寫(xiě)完一定要驗(yàn)證。我習(xí)慣拿幾個(gè)已知含目標(biāo)片段的分子和幾個(gè)已知不含的分子做正負(fù)對(duì)照確認(rèn) SMARTS 的匹配行為符合預(yù)期。這一步花不了幾分鐘但能避免后面批量篩選時(shí)出現(xiàn)系統(tǒng)性錯(cuò)誤。第二高亮圖不要追求一次畫(huà)完所有匹配。如果一個(gè)分子里有多個(gè)匹配全部高亮?xí)寛D變得難以辨認(rèn)。我通常只高亮第一個(gè)匹配或者用不同顏色高亮前兩三個(gè)剩下的在文字里說(shuō)明。第三批量處理時(shí)注意內(nèi)存。RDKit 的分子對(duì)象雖然不算大但如果你一次加載幾萬(wàn)個(gè)分子并全部保留在內(nèi)存里還是會(huì)吃不少資源。我一般用生成器逐行讀 SDF 或 SMILES 文件處理完一個(gè)就釋放一個(gè)避免內(nèi)存堆積。第四繪圖尺寸要跟用途匹配。如果是放進(jìn)論文或報(bào)告用size(600, 450)或更大保證打印清晰如果只是在終端里快速看一眼size(300, 200)就夠了。SVG 格式適合需要縮放的場(chǎng)景PNG 適合直接嵌入文檔。第五保留中間結(jié)果。子結(jié)構(gòu)搜索的結(jié)果匹配的原子索引建議保存下來(lái)不要每次重新算。尤其是當(dāng)你的 SMARTS 比較復(fù)雜、分子比較大的時(shí)候重新搜索可能很耗時(shí)。把匹配結(jié)果存成 JSON 或 pickle下次直接讀能省不少時(shí)間。這套流程我在多個(gè)項(xiàng)目里反復(fù)用過(guò)從幾十個(gè)分子的手工驗(yàn)證到幾萬(wàn)個(gè)分子的批量篩選基本都能覆蓋。RDKit 的子結(jié)構(gòu)搜索和高亮繪制雖然只是它眾多功能中的一小塊但把這一塊吃透能解決化學(xué)信息學(xué)里相當(dāng)一部分日常需求。