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

Загружаем необходимые модули:

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

Нарисуем ибупрофен:

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 Lipinksy
print (Lipinksy.NumHDonors(ibu))
print (Lipinksy.NumHAcceptors(ibu))
print (Lipinksy.rdMolDescriptors.CalcExactMolWt(ibu))
print (Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0])

1
1
206.130679816
3.073200000000001

Допустим, реакция подет по механизму азид-алкинового циклоприсоединения. Продуктом такой клик-химии в реакции азида с модифицированным ибупрофеном может быть соединение следующего строения:

prod=Chem.MolFromSmiles('N2N=NC=C2C(C)CC1=CC=C(C=C1)C(C)C(=O)O')
display(prod)

print (Lipinksy.NumHDonors(prod))
print (Lipinksy.NumHAcceptors(prod))
print (Lipinksy.rdMolDescriptors.CalcExactMolWt(prod))
print (Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(prod)[0])

2
3
259.132076784
2.339

Загрузим с сайта PubChem smiles всех соединений, содержащих азидную группу. При этом воспользуемся следующими фильтрами: MW<300, HBD<10, HBA<15. В полученном файле отфильтруем соединения, содержащие нековалентные связи:

infile = open('pubchem_output.txt', 'r')
lines = infile.read().split("\n")[:-1]
smiles = [] 
for s in lines:
    s = s.split()
    if len(s[1]) < 30 and'.' not in s[1]:
        smiles.append(s[1])

В итоге получим список со smiles из 98624 соединений. Заменим в них азидную группу на модифицированный ибупрофен, отфильтруем результаты с использованием правила Лепински. Таким образом получим список из 87624 соединений -- возможных продуктов нашей клик-химии.

products=[]
for smi in smiles:   
    if ('N=[N+]=[N-]' or 'N=[N+]=[NH]') in smi:
        newsmi=smi.replace('N=[N+]=[N-]','N2N=NC=C2C(C)CC1=CC=C(C=C1)C(C)C(=O)O')
    else:
        continue
    try:
        newmol=Chem.MolFromSmiles(newsmi)
        if Lipinksy.NumHDonors(newmol) <= 5 and Lipinksy.NumHAcceptors(newmol) <= 10 and Lipinksy.rdMolDescriptors.CalcExactMolWt(newmol) < 500 and Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(newmol)[0] <=5:
            products.append(newmol)     
    except:
        pass

RDKit ERROR: [23:25:25] Explicit valence for atom # 4 Cl, 3, is greater than permitted
RDKit ERROR: [23:25:31] Explicit valence for atom # 1 Cl, 2, is greater than permitted

Как выглядят некоторые из полученных продуктов:

Chem.Draw.MolsToGridImage( products[:6], molsPerRow=3, subImgSize=(400, 400))

Построим Similiraty Map ибупрофена с пятым веществом из полученного массива продуктов:

from rdkit.Chem.Draw import SimilarityMaps
fp = SimilarityMaps.GetMorganFingerprint(products[4], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, products[4], SimilarityMaps.GetMorganFingerprint)

Построим 3D структуру этого продукта:

m3d=Chem.AddHs(products[4])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
m3d

# не работает
import nglview as nv
nv.show_rdkit(m3d)

Widget Javascript not detected.  It may not be installed properly. Did you enable the widgetsnbextension? If not, then run "jupyter nbextension enable --py --sys-prefix widgetsnbextension"

Chem.rdmolfiles.MolToPDBFile(m3d, 'prod4.pdb')

from multiprocessing import Process

def launch_pymol():
    !pymol  -R &> /dev/null
Process(target=launch_pymol).start()

import xmlrpc.client as xmlrpclib
cmd = xmlrpclib.ServerProxy('http://localhost:9123')

import time
cmd.load('prod4.pdb')
cmd.do('''
bg_color white
orient
mset 1x200
frame 1
mview store
frame 100
turn y, 180
mview store
''')

cmd.do('''
python
cmd.mstop()
for x in range(0, 200, 3):
    cmd.mpng('mv/', x + 1, x + 1, mode=1, width=720, height=480)
python end
''')

import imageio
import os
filenames = sorted(os.listdir('mv'))[1:]
with imageio.get_writer('movie.gif', mode='I', fps=15) as writer:   
    for filename in filenames:
        filepath = os.path.join('mv', filename)
        image = imageio.imread(filepath)
        writer.append_data(image)

from IPython.display import Image, display
display(Image(filename='movie.gif', format='png'))