基于RDKit绘制化学反应

导入库

from rdkit import RDConfig
import unittest
import random
from rdkit import Chem
from rdkit.Chem import Draw, AllChem
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit import Geometry
%matplotlib inline
from numpy.polynomial.polynomial import polyfit
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
from IPython.display import SVG, display
import seaborn as sns; sns.set(color_codes=True)

定义反应

rxn = AllChem.ReactionFromSmarts('[CH3:1][C:2](=[O:3])[OH:4].[CH3:5][NH2:6]>CC(O)C.[Pt]>[CH3:1][C:2](=[O:3])[NH:6][CH3:5].[OH2:4]',useSmiles=True)
d = Draw.MolDraw2DSVG(900, 300)
d.DrawReaction(rxn)
d.FinishDrawing()

绘制反应



定义绘制反应并高亮

rxn = AllChem.ReactionFromSmarts('[CH3:1][C:2](=[O:3])[OH:4].[CH3:5][NH2:6]>CC(O)C.[Pt]>[CH3:1][C:2](=[O:3])[NH:6][CH3:5].[OH2:4]',useSmiles=True)
colors=[(0.3,0.7,0.9),(0.9,0.7,0.9),(0.6,0.9,0.3),(0.9,0.9,0.1)]
d = Draw.MolDraw2DSVG(900, 300)
d.DrawReaction(rxn,highlightByReactant=True,highlightColorsReactants=colors)
d.FinishDrawing()

svg = d.GetDrawingText()
svg2 = svg.replace('svg:','')
svg3 = SVG(svg2)
display(svg3)



 

参考

1.https://www.rdkit.org/docs/source/rdkit.Chem.rdChemReactions.html

2. https://nodepit.com/node/org.rdkit.knime.nodes.onecomponentreaction2.RDKitOneComponentReactionNodeFactory

3. https://www.kesci.com/home/project/5c7685191ce0af002b556cc5

 

 

Mol对象的组成

由于Mol对象是分子，因此它们自然是由原子组成的。分子是通过原子间键的形成而形成的。

 

导入库

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
print('rdkit version: {}'.format(rdBase.rdkitVersion))


rdkit version: 2019.09.3


载入数据

suppl = Chem.SDMolSupplier('sdf_20200114172835.sdf')
mols = [x for x in suppl if x is not None]
len(mols)

 

Atom对象


Mol.GetAtoms()
Mol.GetAtomWithIdx(idx) 


 

for mol in mols[:5]:
    print(mol.GetProp('PRODUCT_NAME'))
    for a in mol.GetAtoms():
         if a.GetSymbol() == 'C' and str(a.GetHybridization()) == 'SP3':
            print('index for sp3 carbon: {}'.format(a.GetIdx()))
    print('###')


Methyl Acetylsalicylate

index for sp3 carbon: 9

index for sp3 carbon: 12

###

Methyl o-Anisate

index for sp3 carbon: 9

index for sp3 carbon: 11

###

Dimethyl 4-Acetoxyisophthalate

index for sp3 carbon: 12

index for sp3 carbon: 13

index for sp3 carbon: 16

###

Ethyl Acetylsalicylate

index for sp3 carbon: 10

index for sp3 carbon: 11

index for sp3 carbon: 13

###

(+)-Bicuculline

index for sp3 carbon: 2

index for sp3 carbon: 4

index for sp3 carbon: 5

index for sp3 carbon: 8

index for sp3 carbon: 10

index for sp3 carbon: 16

index for sp3 carbon: 25

###


 

Bond对象


Mol.GetBonds()

Mol.GetBondWithIdx(idx)

Mol.GetBondBetweenAtoms()


 

for mol in mols[-3:]:
    print(mol.GetProp('PRODUCT_NAME'))
    for b in mol.GetBonds():
         if b.GetIsAromatic():
                print('bond between {}-{} is aromatic.'.
             format(b.GetBeginAtomIdx(),b.GetEndAtomIdx()))
    print('###')


2-Ethylhexyl Salicylate

bond between 3-4 is aromatic.

bond between 3-8 is aromatic.

bond between 4-5 is aromatic.

bond between 5-6 is aromatic.

bond between 6-7 is aromatic.

bond between 7-8 is aromatic.

###

Propyl Salicylate

bond between 3-4 is aromatic.

bond between 3-8 is aromatic.

bond between 4-5 is aromatic.

bond between 5-6 is aromatic.

bond between 6-7 is aromatic.

bond between 7-8 is aromatic.

###

3,3,5-Trimethylcyclohexyl Salicylate (cis- and trans- mixture)

bond between 0-1 is aromatic.

bond between 1-2 is aromatic.

bond between 2-3 is aromatic.

bond between 3-4 is aromatic.

bond between 4-5 is aromatic.

bond between 0-5 is aromatic.

###


 

属性


GetProp(name)
SetProp(name, value)
GetPropNames()


可以将称为属性的任意值添加到Mol对象。使用SetProp进行设置，并使用GetProp作为参考。如果从头开始在原始SDF中设置属性，则在创建对象的同时设置它们的值。

 

for p in mols[0].GetPropNames():
     print('{}: {}'.format(p, mols[0].GetProp(p)))


PRODUCT_NUMBER: A0114

PRODUCT_NAME: Methyl Acetylsalicylate

MOLECULAR_FORMULA: C10H10O4

MOLECULAR_WEIGHT: 194.19

CAS_NUMBER: 580-02-9

MDL_NUMBER: MFCD00014978


 

分子的2D结构图


1) 一个分子绘制

a. Chem.Draw.MolToImage(mol)

b. Chem.Draw.MolToFile(mol, file_name)


2) 多个分子的绘制

a. Chem.Draw.MolsToImage(mols)

b. Chem.Draw.MolsToGridImage(mols, molsPerRow, subImgSize, legends)


处理的所有分子都具有原子坐标，但是，如果未初始化它们，则需要使用AllChem.Compute2DCoordinate计算。

Draw.MolToImage(mols[0])



 

Draw.MolsToGridImage(mols[:9], molsPerRow=3, subImgSize=(300,200),
                     legends=[x.GetProp('PRODUCT_NUMBER') for x in mols[:9]])



基于模板结构并行化结构

由于化合物的排列在网格图像中不匹配，因此一眼就很难理解类似的结构。因此，让我们以水杨酸为模板旋转分子。在Chem.AllChem中使用GenerateDepictionMatching2DStructure。所要做的就是为模板创建一个Mol对象并计算2D坐标。

import pubchempy as pcp
tmp = pcp.get_compounds('salicylic acid', 'name')
tmp = tmp[0]
tmp_smiles = tmp.canonical_smiles
template = Chem.MolFromSmiles(tmp_smiles)
AllChem.Compute2DCoords(template)

for mol in mols:
    if mol.HasSubstructMatch(template):
        AllChem.GenerateDepictionMatching2DStructure(mol, template)

Draw.MolsToGridImage(mols[:9], molsPerRow=3, subImgSize=(300,200),
                    legends=[x.GetProp('PRODUCT_NUMBER') for x in mols[:9]])



 

 亚结构的处理

 

aspirin = pcp.get_compounds('aspirin', 'name')
aspirin = aspirin[0]
aspirin_sm = aspirin.canonical_smiles
aspirin_mol = Chem.MolFromSmiles(aspirin_sm)
AllChem.Compute2DCoords(aspirin_mol)

match = []
for mol in mols:
    if mol.HasSubstructMatch(aspirin_mol):
        match.append(mol)
print(len(match)) # 6
Draw.MolsToGridImage(match, molsPerRow=3, subImgSize=(300,200),
                     legends=[x.GetProp('PRODUCT_NAME') for x in match])



 

参考

1. http://www.rdkit.org/docs/index.html

2. http://www.rdkit.org/docs/api-docs.html

基于Python和RDKit对化合物数据进行预处理。

 

环境

MolVS是专门用于化合物预处理的库。

rdkit 2020.03
	
molvs 0.1.1
 

化合物（Compound）预处理

RDKit：SanitizeMol

Kekure的形成，化合价的确认，芳香性的设定，结合等。

参考：http://rdkit.org/docs/source/rdkit.Chem.rdmolops.html

 

MolVS : Normarize

参考：https://molvs.readthedocs.io/en/latest/guide/standardize.html

进行了一系列转换，以修复常见错误并标准化特征组。

from rdkit import Chem
from molvs.normalize import Normalizer, Normalization

old_smiles = "[Na]OC(=O)c1ccc(C[S+2]([O-])([O-]))cc1"
print("PREV:" + old_smiles)
old_mol = Chem.MolFromSmiles(old_smiles)
normalizer = Normalizer(normalizations=[Normalization('Sulfone to S(=O)(=O)', '[S+2:1]([O-:2])([O-:3])>>[S+0:1](=[O-0:2])(=[O-0:3])')])
new_mol = normalizer.normalize(old_mol)
new_smiles = Chem.MolToSmiles(new_mol)
print("NEW:" + new_smiles)

 以上，选择性地执行在“S（＝ O）（＝ O）”中定义的归一化处理。

结果如下，硫原子和氧原子的电荷发生了变化。

如果生成不带参数的规范化器，则将执行MolVS中预先定义的所有规范化过程。


PREV:[Na]OC(=O)c1ccc(C[S+2]([O-])([O-]))cc1
NEW:O=C(O[Na])c1ccc(C[S](=O)=O)cc1


 

MolVS : TautomerCanonicalizer

参考：https://molvs.readthedocs.io/en/latest/guide/tautomer.html

互变异构体似乎是一组通过氢原子的移动而易于彼此交换的分子。

from rdkit import Chem
from molvs.tautomer import TAUTOMER_TRANSFORMS, TAUTOMER_SCORES, MAX_TAUTOMERS, TautomerCanonicalizer, TautomerEnumerator, TautomerTransform


tautomerCanonicalizer = TautomerCanonicalizer((
    TautomerTransform('1,7 aromatic heteroatom H shift r', '[#7,S,O,Se,Te,CX4;!H0]-[#6,#7X2]=[#6]-[#6,#7X2]=[#6,#7X2]-[#6,#7X2]=[NX2,S,O,Se,Te]'),
    ))

mol = Chem.MolFromSmiles("O=C1CC=CC=C1")
print("prev:" + Chem.MolToSmiles(mol))
mol2 = tautomerCanonicalizer.canonicalize(mol)
print("after: "+ Chem.MolToSmiles(mol2))


prev:O=C1C=CC=CC1
after: Oc1ccccc1


 

MolVS : LargestFragmentChooser

参考：https://molvs.readthedocs.io/en/latest/api.html#molvs-fragment

当包含多个分子时，它返回最大的分子。

from rdkit import Chem
from molvs.fragment import LargestFragmentChooser

flagmentChooser1 = LargestFragmentChooser()
old_smiles = "O=S(=O)(Cc1[nH]c(-c2ccc(Cl)s2)c[s+]1)c1cccs1.[Br-]"
print("prev:" + old_smiles)
mol = Chem.MolFromSmiles(old_smiles)
mol2 = flagmentChooser1(mol)
print("after:" + Chem.MolToSmiles(mol2))


prev:O=S(=O)(Cc1[nH]c(-c2ccc(Cl)s2)c[s+]1)c1cccs1.[Br-]
after:O=S(=O)(Cc1[nH]c(-c2ccc(Cl)s2)c[s+]1)c1cccs1


 

MolVS : Uncharger

参考：https://molvs.readthedocs.io/en/latest/api.html#molvs-charge

试图中和分子上的离子化酸和碱。

from molvs.charge import Reionizer, Uncharger

uncharger = Uncharger()
mol = Chem.MolFromSmiles("c1cccc[nH+]1")
print("prev:" + Chem.MolToSmiles(mol))
mol2 = uncharger(mol)
print("after:" + Chem.MolToSmiles(mol2))


prev:c1cc[nH+]cc1
after:c1ccncc1


 

不夸张，我用过最香的rdkit安装

引用：https://anaconda.org/rdkit/rdkit