• 一、描述符計算模塊
  • 1.rdkit.Chem.Lipinski模塊
  • 2.rdkit.Chem.Descriptors模塊
  • 3.rdkit.ML.Descriptors.MoleculeDescriptors模塊
• 二、原子描述符可視化
  • 1.原子partial charge可視化
  • 2.原子logP可視化

一、描述符計算模塊

1.rdkit.Chem.Lipinski模塊

rdkit中提供了許多描述符的計算方法，可用于分子篩選、成藥性評估等。以lipinski類藥的相關規(guī)則為例，可以通過rdkit.Chem.Lipinski模塊進行計算。常用的一些性質舉例：

• 氫鍵受體數(shù)NumHAcceptors
• 氫鍵供體數(shù)NumHDonors
• 可旋轉鍵數(shù)NumRotatableBonds
• 脂肪環(huán)數(shù)量NumAliphaticRings
• 芳香環(huán)數(shù)量NumAromaticRings
• SP3雜化碳原子比例FractionCSP3
• ……

>>> from rdkit import Chem
>>> from rdkit.Chem import Lipinski
>>> mol = Chem.MolFromSmiles('c1ccccc1C(=O)O')
>>> Lipinski.NumHAcceptors(mol)
1

• 將描述符名稱(key)和值(val)添加到分子屬性中：m.SetProp(key, val)
• 獲取分子屬性：m.GetProp(key)

>>> Ha = Lipinski.NumHAcceptors(mol)
>>> mol.SetProp('Ha', '%s'%Ha)
>>> mol.GetProp('Ha')
'1'

2.rdkit.Chem.Descriptors模塊

大部分的描述符都可以通過rdkit.Chem.Descriptors模塊進行計算。該模塊也包含了Lipinski的描述符。常用的一些描述符舉例：

• 分子量MolWt
• 脂水分配系數(shù)MolLogP
• 拓撲極表面積TPSA
• 以計算TPSA為例：Descriptors.TPSA()

>>> from rdkit.Chem import Descriptors
>>> print(Descriptors.TPSA(mol), Descriptors.MolLogP(mol), Descriptors.MolWt(mol))
37.3 1.3848 122.12299999999998

3.rdkit.ML.Descriptors.MoleculeDescriptors模塊

該模塊可以批量計算描述符。

• 先指定一個列表des_list，包含所要計算的描述符名稱，
• 使用MolecularDescriptorCalculator創(chuàng)建一個計算描述符的對象，傳入要計算的des_list
• 調用對象的CalcDescriptors方法，傳入要計算的mol對象，得到所需的描述符

>>> from rdkit.ML.Descriptors import MoleculeDescriptors
>>> des_list = ['MolWt', 'NumHAcceptors', 'NumHDonors', 'MolLogP', 'NumRotatableBonds']
>>> calculator = MoleculeDescriptors.MolecularDescriptorCalculator(des_list)
>>> calculator.CalcDescriptors(mol)
(122.12299999999998, 1, 1, 1.3848, 1)

• 查看各種描述符的含義：GetDescriptorSummaries()

>>> calculator.GetDescriptorSummaries()
['The average molecular weight of the molecule',
 'Number of Hydrogen Bond Acceptors',
 'Number of Hydrogen Bond Donors',
 'Wildman-Crippen LogP value',
 'Number of Rotatable Bonds']

• 獲取所有描述符：Descriptors._descList

>>> des_list = [x[0] for x in Descriptors._descList]
>>> len(des_list)
200

• 將calculator保存起來：SaveState(filename)
• 再次調用該calculator：pickle.load()

>>> import pickle
>>> calculator.SaveState('data/descriptor_calculator')
>>> with open('data/descriptor_calculator', 'rb') as f:
>>>     calc = pickle.load(f)
>>> calc.CalcDescriptors(mol)
(122.12299999999998, 1, 1, 1.3848, 1)

二、原子描述符可視化

可以用相似性地圖來查看每個原子對描述符的貢獻，更多相似性地圖的應用可以查看這篇文章的第二部分。

1.原子partial charge可視化

計算partial charge要復雜一點。計算出的partial charge存儲在每個原子的屬性中，可以通過GetDoubleProp（浮點數(shù)）或GetProp（字符串）來獲取。

partial charge可以表示電子的分布。分子中的化學鍵是由分布在相連原子周圍的電子對組成的。但由于原子的電負性不同，成鍵電子并不是均勻分布。電負性大的原子吸電子能力強，成鍵的電子對會更偏向該原子，導致該原子帶有負電的性質。相對應的，電負性小的原子則帶有正電性質。Partial charge衡量了電子偏向的程度。

• 先計算partial charge：AllChem.ComputeGasteigerCharges()
• 獲取某一個原子：GetAtomWithIdx()
• 獲取原子的partial charge：GetDoubleProp('_GasteigerCharge')

>>> from rdkit.Chem import AllChem
>>> mol = Chem.MolFromSmiles('COc1cccc2cc(C(=O)NCCCCN3CCN(c4cccc5nccnc54)CC3)oc21')
>>> AllChem.ComputeGasteigerCharges(mol)
>>> atom = mol.GetAtomWithIdx(0)
>>> atom.GetDoubleProp('_GasteigerCharge')
0.07771844728655561

• 獲取每個原子的partial charge，放到contribs中

>>> contribs = [round(mol.GetAtomWithIdx(i).GetDoubleProp('_GasteigerCharge'), 2) for i in range(mol.GetNumAtoms())]
>>> print(contribs)
[0.08, -0.49, 0.16, -0.02, -0.06, ...]

• 根據(jù)給定的權重，生成分子權重圖：GetSimilarityMapFromWeights(mol, weights, colorMap, contourLines, ...)
  mol：要繪制的mol對象
  weights：權重
  colorMap：matplotlib中的色系
  contourLines：等高線數(shù)量

>>> from rdkit.Chem.Draw import SimilarityMaps
>>> fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, contribs, contourLines=10)

1

2.原子logP可視化

LogP表示脂水分配系數(shù)，該值認為與細胞通透性有一定相關性。rdkit中提供的Descriptors.MolLogP()方法可以粗略計算logP值，該方法首先做了一個原子分類系統(tǒng)，根據(jù)原子及其相連原子的不同而進行分類，再對化學性質相似、logP貢獻相似的類別做合并，最終得到了68種精確的原子類別和4種通配類別，并用SMARTS表示。計算時，對一個分子中所有原子進行分類，再乘以每一類的權重并加和，最終得到LogP值。該方法在9920個分子的訓練集上的r2為0.918，標準差為0.667。此外摩爾折射率（molar refractivity，MR）也可以通過這種方法計算得到。

• 計算每個原子的logP和MR值：rdMolDescriptors._CalcCrippenContribs
  返回結果是每個原子logP和MR元組的列表

>>> from rdkit.Chem import rdMolDescriptors
>>> contribs = rdMolDescriptors._CalcCrippenContribs(mol)
>>> contribs[:3]
[(-0.2035, 2.753), (-0.4195, 1.182), (0.5437, 3.853)]

• 生成分子權重圖：GetSimilarityMapFromWeights()

>>> fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, [x for x,y in contribs])

2

本文參考自rdkit官方文檔。
代碼及源文件在這里。