Хемоинформатика

Импортировать необходимые модули

from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, Draw
import numpy as np
from IPython.display import display, Image, display_svg

Создать молекулу ибупрофена

ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)

Посчитать параметры для правила Лепински

import rdkit.Chem.Lipinski as Lip
print (Lip.NumHDonors(ibu))
print (Lip.NumHAcceptors(ibu))
print (Lip.rdMolDescriptors.CalcExactMolWt(ibu))
print (Lip.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0])

1
1
206.130679816
3.073200000000001

Создать шаблон для замены азидного радикала

template='N1N=NC(=C1)CC2=CC=C(C=C2)C(C)C(=O)O'
display(Chem.MolFromSmiles(template))
print (Lip.NumHDonors(Chem.MolFromSmiles(template)))
print (Lip.NumHAcceptors(Chem.MolFromSmiles(template)))
print (Lip.rdMolDescriptors.CalcExactMolWt(Chem.MolFromSmiles(template)))
print (Lip.rdMolDescriptors.CalcCrippenDescriptors(Chem.MolFromSmiles(template))[0])

<rdkit.Chem.rdchem.Mol at 0x18c866c93f0>

2
3
231.100776656
1.5836000000000001

Скачать с PubChem список азидов в формате SMILES (вручную), прочитать из файла

strings=np.genfromtxt('azides.txt',dtype=np.str)
az=[]
for line in strings:
    if len(line[1]) < 30 and not '.' in line[1]:
        az.append( line[1])
print (len(az))

1032

Заменить азидную группу в каждом элементе списка на ибупрофен (как если бы произошла реакция Хьюсгена с аналогом ибупрофена, имеющим терминальную алкиновую группу), проверить на соответствие правилу 5

new=[]
for a in az:  
    if 'N=[N+]=[N-]' in a:
        b=a.replace('N=[N+]=[N-]',template)
    else:
        continue
    try:
        newmol=Chem.MolFromSmiles(b)
        AllChem.Compute2DCoords(newmol)
        if ((Lip.NumHDonors(newmol)<=5)and (Lip.NumHAcceptors(newmol)<=10) and 
            (Lip.rdMolDescriptors.CalcExactMolWt(newmol)<500) and 
            (Lip.rdMolDescriptors.CalcCrippenDescriptors(newmol)[0]<=5)):             
            new.append(newmol)  
    except:
        pass
print (len(new))

783

Вывести все полученные молекулы на экран

display_svg(Draw.MolsToGridImage(new,useSVG=True, molsPerRow=3, subImgSize=(300, 200),maxMols=1000))

Показать 6 элемент списка вместе с ибупрофеном

display_svg(Draw.MolsToGridImage([ibu,new[5]],useSVG=True, molsPerRow=2, subImgSize=(400, 200)))

Построить Similiraty Map ибупрофена и 6 элемента списка

from rdkit.Chem.Draw import SimilarityMaps
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, new[5], 
                                                               SimilarityMaps.GetMorganFingerprint)
print(maxweight)

0.2085415720127215

Построить трехмерную структуру молекулы, показать

m3d=Chem.AddHs(new[5])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
display(m3d)

import nglview as nv
nv.show_rdkit(m3d)

Nglview не работает, поэтому использовался drawMol3D. Для его использования необходимо предварительно скачать py3Dmol (conda install -c conda-forge py3Dmol)

import py3Dmol
from rdkit.Chem.Draw import IPythonConsole
IPythonConsole.drawMol3D(m3d)

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol