rdkit 處理2D、3D分子

Smiles 可以看成分子的1D形式，分子的平面結構可以看成分子的2D形式。該算法能夠減少分子中原子在平面内的碰撞，使得繪制的分子更加清晰。

文章目錄

• 一、引入所需庫
• 二、處理2D分子
• • 2.1 計算分子的2D坐标函數解析
  • 2.2 生成一個分子的描述函數解析
  • 2.3 計算分子的2D坐标 結果坐标存儲在分子的每個原子上
• 三、處理3D分子
• • 3.1 産生3D構象
  • 3.2 産生多個3D構象
• 四、 計算構象的RMS值
• • 4.1 計算其他構象和第一個構象的RMS值
  • 4.2 計算指定兩個構象的RMS值
  • 4.3 MMFF立場對構象進行優化
• 五、儲存分子對象
• • 5.1使用pickling機制進行儲存。
  • 5.2 儲存為普通二進制檔案

一、引入所需庫

#! /usr/bin/python
# coding: utf-8
# rdkit 處理2D、3D分子


from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw

import pickle

           

二、處理2D分子

如果有多個分子共享一個骨架，我們希望繪圖的時候能夠保證在這些分子中的骨架方向是一緻的。 首先把骨架提取成繪圖模闆，然後在繪圖上添加基團。

我們可以看到這4個結構都含有吡咯并嘧啶（c1nccc2n1ccc2）的子結構， 但是在圖檔上他們的子結構取向不一緻，不便于比較他們的結構。 我們可以采用模闆繪制法，固定公共子結構的取向。

smis = [
    'COC1=C(C=CC(=C1)NS(=O)(=O)C)C2=CN=CN3C2=CC=C3',
    'C1=CC2=C(C(=C1)C3=CN=CN4C3=CC=C4)ON=C2C5=CC=C(C=C5)F',
    'COC(=O)C1=CC2=CC=CN2C=N1',
    'C1=C2C=C(N=CN2C(=C1)Cl)C(=O)O',
]
template = Chem.MolFromSmiles('c1nccc2n1ccc2')
           

2.1 計算分子的2D坐标函數解析

Compute2DCoords（
    （Mol）mol   # 目标分子
    [,（bool）canonOrient = True #以規範的方式定向分子
    [,（bool）clearConfs = True # 如果為真，則該分子上的所有現有構象将被清除
    [,（dict）coordMap = {}  # 映射原子ID的字典-> Point2D對象與應該鎖定其位置的原子的起始坐标。
    [,（int）nFlipsPerSample = 0  # 可旋轉債券的數量 一次随機翻轉。
    [,（int）nSample = 0   # 可旋轉鍵的随機采樣數。
    [,（int ）sampleSeed = 0  # 随機抽樣過程的種子。
    [,（bool）permuteDeg4Nodes = False  # 允許鍵的排列為4度
    [,（float）bondLength = -1.0   # 更改描述的預設鍵長
    [,（bool）forceRDKit = False 
    ] ] ] ] ] ] ] ] ] 
)
           

2.2 生成一個分子的描述函數解析

GenerateDepictionMatching2DStructure(
    （Mol）mol, # mol-要排列的分子，它會回來一個單一的conformer。
    （Mol）reference  # 與參考原子對齊的分子；這應該有一個描述。
    [,（int）confId = -1  # confId-（可選）要使用reference的參考構象的idPattern-（可選）要用于的查詢分子生成分子和參考之間的原子映射。
    [,（AtomPairsParameters）refPatt = None 
    [,（bool）acceptFailure = False  # （可選）如果為True，則将生成标準描述  對于沒有與參考比對的子結構的分子；如果為False，則抛出DepictException。
    [,（bool）forceRDKit = False 
    ] ] ] ] 
)
           

2.3 計算分子的2D坐标 結果坐标存儲在分子的每個原子上

AllChem.Compute2DCoords(template)
mols = []
for smi in smis:
    mol = Chem.MolFromSmiles(smi)
    # 生成一個分子的描述，其中一部分 分子被限制為具有與參考相同的坐标。
    AllChem.GenerateDepictionMatching2DStructure(mol, template)
    mols.append(mol)

# 基于分子檔案輸出分子結構
img = Draw.MolsToGridImage(
    mols,   # mol對象
    molsPerRow=4,
    subImgSize=(200,200),
    legends=['' for x in mols]
)
img.save('/Users/zeoy/st/drug_development/st_rdcit/img/mol9.jpg')
           

分子結構

三、處理3D分子

RDKit可以使用兩種不同的方法為分子生成構象。原始方法使用距離幾何。請注意，此過程産生的構象往往很難看。應使用力場對其進行清理。可以在RDKit中使用通用力場（UFF）的實作來完成。還有一種Riniker和Landrum方法的實作，該方法使用劍橋幾何結構資料庫（CSD）中的扭轉角首選項在已經使用距離幾何圖形生成它們之後校正矯正器。使用這種方法，無需使用最小化步驟來清理結構。

從RDKit的2018.09版本開始，ETKDG是預設的構象生成方法。

3.1 産生3D構象

預設情況下，RDKit分子在圖中沒有明确存在H原子，但是H原子對于獲得逼真的幾何形狀很重要，是以通常在産生3D構象前，需要為分子添加H原子。

mol = Chem.MolFromSmiles('c1nccc2n1ccc2')
m2 = Chem.AddHs(mol) # 加氫原子
AllChem.EmbedMolecule(m2) # 2D->3D化
AllChem.MMFFOptimizeMolecule(m2)   # 使用MMFF94最小化RDKit生成的構象
m3 = Chem.RemoveHs(m2) # 删除氫原子
           

3.2 産生多個3D構象

AllChem.EmbedMultipleConfs(m2, numConfs=10) # 為m2分子産生10個構象，儲存在m2分子中。
m2.GetConformer(1)  # 通路指定構象
m2.GetConformers()  # 擷取構象

# 産生多個構象比較耗時，支援多線程加速。
cids = AllChem.EmbedMultipleConfs(m2,numConfs=10, numThreads=4)
print('構象個數=',len(cids))
           

四、 計算構象的RMS值

4.1 計算其他構象和第一個構象的RMS值

rmslist = []
AllChem.AlignMolConformers(m2, RMSlist=rmslist)
print('第一個構想者與所有其他構想者之間的RMS值')
print(rmslist)
# [0.031617482029182735, 0.027100784678672763, 0.03211558617566241, 0.03876593076230512, 0.03981170998684297, 0.023608843156566642, 0.023619922568325093, 0.037621957801700204, 0.0532861818792842]
           

4.2 計算指定兩個構象的RMS值

rms = AllChem.GetConformerRMS(m2, 1, 9, prealigned=True)
print('計算指定兩個構象的RMS值=', rms) # 計算指定兩個構象的RMS值= 0.049340647238401286
           

4.3 MMFF立場對構象進行優化

res = AllChem.MMFFOptimizeMoleculeConfs(m2)
print(res)
# [(0, 44.63060837612333), (0, 44.6306083325595), (0, 44.63060829914013), (0, 44.63060834780363), (0, 44.63060829995848), (0, 44.63060839000017), (0, 44.630608353905146), (0, 44.6306086577808), (0, 44.630608361427775), (0, 44.63060856588376)]
# 構象優化比較耗時，支援多線程加速
res = AllChem.MMFFOptimizeMoleculeConfs(m2, numThreads=10)
           

五、儲存分子對象

RDKit内置了2種儲存分子對象的方法。

5.1使用pickling機制進行儲存。

pkl = pickle.dumps(mol)
print(pkl)
m2 = pickle.loads(pkl)
mol2 = Chem.MolToSmiles(m2)
print(mol2)  # c1nccc2n1ccc2
           

5.2 儲存為普通二進制檔案

rdkit的pickling檔案非常緊湊，從pickling檔案中加載分子比從mol檔案和SMILES字元串開很多。

将經常使用的分子，儲存為picklingz是一個很好的辦法。

二進制檔案的大小和分子大小有關系，而pickle檔案和分子大小關系不明顯。推薦用pickle儲存檔案。

binStr = mol.ToBinary()
m2 = Chem.Mol(binStr)
mol2 = Chem.MolToSmiles(M2)