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