RDKit|分子3D構象生成與優化

文章目錄

• 一、構象生成算法概述
• • 1.基于距離（distance-based）
  • 2.基于知識（knowledge-based）
• 二、代碼實作
• • 1.添加氫原子
  • 2.距離幾何算法生成3D結構
  • 3.距離幾何+ETKDG生成3D構象
  • 4.距離幾何+ETKDG生成多構象
  • 5.距離幾何+ETKDG+MMFF生成3D構象
  • 6.距離幾何+ETKDG+MMFF生成多構象
  • 7.多線程生成多構象
• 三、結語

一、構象生成算法概述

先來了解兩種rdkit生成3D構象的方法：

1.基于距離（distance-based）

傳統的距離幾何（Distance Geometry）法生成構象步驟：

• (1) 通過分子的連接配接表資訊和一套規則，生成分子的連接配接邊界矩陣（molecule’s distance bounds matrix ）
• (2) 使用三角形邊界平滑算法（triangle-bounds smoothing algorithm），對邊界矩陣進行平滑處理
• (3) 根據邊界矩陣，随機産生一個距離矩陣。
• (4) 把産生的距離矩陣映射到三維空間中，并為每個原子計算坐标。
• (5) 對計算的坐标結果使用力場和邊界矩陣進行粗略的優化。

這種方法雖然計算速度快，也能得到分子的三維坐标，但純粹基于距離計算會導緻部分結構的扭曲，也就是結果比較醜。想得到更好的結構，需要進行更細緻的力場優化，比如再使用rdkit的通用力場（Universal Force Field，UFF）進行處理。

2.基于知識（knowledge-based）

Riniker和Landrum開發了一種（Experimental-Torsion Basic Knowledge Distance Geometry，ETKDG）的方法。他們從晶體結構資料庫的小分子結構中總結了一些規則（扭轉角的傾向性，torsion angle preference，符合某種SMARTS的結構更傾向于生成某種固定的構象），并用這套規則來修正距離幾何算法産生的構象。如果使用了ETKDG，就可以不再使用力場進行優化了。

在rdkit中，目前已經将ETKDG作為預設的構象生成方法。使用代碼實作這種方法比上面的廢話都要簡單，幾行代碼就能搞定了。

>>> from rdkit import Chem
>>> from rdkit.Chem import AllChem
>>> from rdkit.Chem import Draw
>>> m = Chem.MolFromSmiles('c1ccccc1')
>>> m3d=Chem.AddHs(m)
>>> AllChem.EmbedMolecule(m3d, randomSeed=1)
>>> Draw.MolToImage(m3d, size=(250,250))
           

二、代碼實作

1.添加氫原子

• 加氫：Chem.AddHs()

  在rdkit中，分子在預設情況下是不顯示氫的，但氫原子對于真實的幾何構象計算有很大的影響，是以在計算3D構象前，需要使用Chem.AddHs()方法加上氫原子。

• 減氫：Chem.RemoveHs()

  相對應的，可以使用該方法把氫去除。

>>> m = Chem.MolFromSmiles('CC1=CC=C(C=C1)NC2=C3C(=C(C=C2)NC4=CC=C(C=C4)C)C(=O)C5=CC=CC=C5C3=O')
>>> m3d = Chem.AddHs(m)
>>> print(m.GetNumAtoms(), m3d.GetNumAtoms())
32 54
           

2.距離幾何算法生成3D結構

使用距離幾何算法初始化3D坐标。

• 生成3D構象：AllChem.EmbedMolecule(mol, randomSeed, clearConfs, useExpTorsionAnglePrefs, useBasicKnowledge, …)

  mol：傳入mol對象

  randomSeed：随機種子，友善結果重複

  clearConfs：清除已有構象，預設True

  useExpTorsionAnglePrefs和useBasicKnowledge兩個參數即控制是否使用ETKDG，預設都為True

>>> AllChem.EmbedMolecule(m3d, randomSeed=10, useExpTorsionAnglePrefs=False, useBasicKnowledge=False)
>>> Draw.MolToImage(m3d, size=(250,250))
           

• 檢視結果

  不知道如何生成圖檔的可以看這篇文章的第四部分。

  

Fig.1.距離幾何算法産生的結構

3.距離幾何+ETKDG生成3D構象

• 生成3D構象：AllChem.EmbedMolecule()

  先用距離幾何初始化3D坐标，再使用ETKDG算法優化

  參數同上，預設useExpTorsionAnglePrefs和useBasicKnowledge為True

>>> AllChem.EmbedMolecule(m3d, randomSeed=10)
>>> Draw.MolToImage(m3d, size=(250,250))
           

• 檢視結果

  比單獨使用距離算法确實順眼了點。

  

Fig.2.距離幾何+ETKDG産生的結構

4.距離幾何+ETKDG生成多構象

無論使用哪種方法，都可以對一個分子生成許多構象，實作也很簡單，無非使用不同的随機種子，多跑幾次距離幾何的計算過程就行了。

• 生成多個構象：AllChem.EmbedMultipleConfs(mol, numConfs, maxAttempts, randomSeed, clearConfs, pruneRmsThresh, useExpTorsionAnglePrefs, useBasicKnowledge, …)

  其中numConfs控制了生成構象的個數。

  mol：mol對象

  numConfs：生成的構象數量

  maxAttempts：嘗試生成構象的最多次數

  randomSeed：随機種子

  clearConfs：清除已有構象

  pruneRmsThresh：根據RMS進行合并

  ETKDG相關參數同上

>>> cids = AllChem.EmbedMultipleConfs(m3d, numConfs=10)
>>> print(len(cids))
10
           

• 擷取某個構象GetConformer(id)

  通過傳入id擷取指定構象

• 計算不同構象間的差異：AlignMolConformers(mol, RMSlist, …)

  對同一分子不同的構象先進行重疊排列，再計算RMS（root mean square）值

  RMSlist：用于接收RMS的清單，它由第一個構象與剩餘9個構象依次比對産生。

>>> rmslist = []
>>> AllChem.AlignMolConformers(m3d, RMSlist=rmslist)
>>> print(len(rmslist))
9
           

• 也可以計算任意兩個指定構象的RMS值GetConformerRMS(mol, confId1, confId2, atomIds, prealigned)

  confId1：第一個構象

  confId2：第二個構象

  atomIds：需要對比的原子，預設全部

  prealigned：構象是否已經對齊，預設False（沒有時，函數會自動将它們對齊）

>>> AllChem.GetConformerRMS(m2, 1, 9, prealigned=True)
1.5515373147254072
           

5.距離幾何+ETKDG+MMFF生成3D構象

對距離幾何産生的構象，進行ETKDG優化後，還可以繼續使用MMFF94等力場進行優化。不過需要注意的是，MMFF力場中的原子類型編碼采用了自身的芳香性模型，是以在使用MMFF相關方法後，分子的芳香屬性（aromaticity flags）會改變。

• 使用MMFF94進行優化構象：MMFFOptimizeMolecule(mol, mmffVariant, maxIters, …)

  mmffVariant：可選"MMFF94"或"MMFF94s"，預設MMFF94

  maxIters：最多疊代次數，預設200

>>> AllChem.EmbedMolecule(m3d, randomSeed=10)
>>> AllChem.MMFFOptimizeMolecule(m3d)
>>> Draw.MolToImage(m3d, size=(250,250))
           

6.距離幾何+ETKDG+MMFF生成多構象

其實當使用ETKDG算法生成3D構象時，通常無需再使用MMFF94力場優化。如果執意要這麼做，有一個友善的函數可以調用

• 使用MMFF94優化多個構象：MMFFOptimizeMoleculeConfs()

  傳回的結果是元組組成的清單，每個元組表示每個構象的收斂值和能量(not_converged, energy)。如果not_converged是0，則收斂，接近最穩态，否則沒有到達最穩态。

>>> res = AllChem.MMFFOptimizeMoleculeConfs(m3d)
>>> res
[(0, 102.4284038486388),
 (1, 102.42249902079311),...]
           

7.多線程生成多構象

可以通過numThreads參數進行設定。

• EmbedMultipleConfs(mol, numThreads)

• MMFFOptimizeMoleculeConfs(mol, numThreads)

  numThreads：預設為1，表示單程序執行。設定為0表示将會使用本機最大線程數執行。

>>> cids = AllChem.EmbedMultipleConfs(m3d, numThreads=0)
>>> res = AllChem.MMFFOptimizeMoleculeConfs(m3d, numThreads=0)
           

三、結語

最後需要注意的是，構象生成是一項複雜且漫長的工作。rdkit生成的3D結構并不能取代真實的構象，它僅僅為需要3D結構的場景提供了一種快速實作的方法。不過，rdkit提供的ETKDG方法還是可以滿足大多數3D需求場景的。

本文參考自rdkit官方文檔。

代碼及源檔案在這裡。