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

Цель занятия используя пакет моудлей RDkit предложить аналог ибупрофена :

Hа сайте PubChem найти все радикалы c азидом для Click Chemistry и скачать их SMILES нотации или использовать pubchempy. Найти формулу ибупрофена и предложить способ изменения его SMILES для эмуляции продукта Click Chemistry. Заменить в найденых радикалах азидную группу на модифцированный ибупрофен. Превратить новые SMILES в объекты-молекулы. Отобрать те молекулы, которые удовлетворяют правилу пяти Lipinski.

Установим всякое, а потом будет импортировать:

conda install -c rdkit rdkit conda install -c jirinovo pubchempy conda install -c conda-forge matplotlib conda install nglview -c conda-forge

In [1]:
import numpy as np
from IPython.display import display, Image

import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [2]:
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 rdkit.Chem.Lipinski as Lipinksy
import pandas as pd
import nglview as nv

Загрузим SMILES ибупрофена и отрисуем его:

In [3]:
smiles_string = 'CC(C)CC1=CC=C(C=C1)C(C)C(=O)O'
ibu = Chem.MolFromSmiles(smiles_string)
AllChem.Compute2DCoords(ibu)
display(ibu)

Функция для проверки соответствия правилам Липински:

In [4]:
def Lipinsky_checking(comp):
    a1 = Lipinksy.NumHDonors(comp) <= 5 
    a2 = Lipinksy.NumHAcceptors(comp) <= 10
    a3 = Lipinksy.rdMolDescriptors.CalcExactMolWt(comp) < 500
    a4 = Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0] <= 5
    return  a1&a2&a3&a4

Проверим ибупрофен.

In [5]:
Lipinsky_checking(ibu)
Out[5]:
True

Модификация ибупрофена для Click Chemistry¶

Хотим сделать производное ибупрофена, выглядящее как продукт азид-алкинового циклоприсоединения.

In [6]:
display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/6/66/CuAAC.svg/768px-CuAAC.svg.png'))
In [7]:
ibu_for_azids = 'C1NN=NC=1C(C(O)=O)C(C=C1)=CC=C1CC(C)C'
click = Chem.MolFromSmiles(ibu_for_azids)
AllChem.Compute2DCoords(click)
display(click)
In [8]:
Lipinsky_checking(click)
Out[8]:
True

Загрузка азидов¶

In [9]:
#Можно дополнить код из задания и скачать азиды самостоятетельно, но это очень долго.
compounds = []
per_page = 10**5
for smiles in ["N=N=N", "NN#N",]:
    for i in range(200):
        try:
            a = pcp.get_properties(
              properties="CanonicalSMILES", 
              identifier=smiles, namespace="smiles", 
              searchtype="substructure",
              RingsNotEmbedded=True,
              listkey_count=per_page, listkey_start=i*per_page
            )
        except:
            print("Page {} returns 500, which is"
                 "our signal to stop.".format(i+1))
            break
        print("Retrieved page {} of {} search".format(i+1, smiles))
        compounds.extend(a)
In [10]:
#А можно взять список у тех, кто это уже сделал.
azids = pd.read_csv("azids.csv")
print(azids.shape)
azids.head()
Out[10]:
Unnamed: 0 CID CanonicalSMILES
0 0 145926031 C1=CC(=CC=C1C2=NOC(=C2)C3=NC(=CS3)C4=CC=NC=C4)...
1 1 145927246 C(CCl)N=[N+]=N
2 2 145927073 COC(=O)C1CC(CN1C(=O)OCC2=CC=CC=C2)N=[N+]=N
3 3 145926982 CN1CCN(CC1)C2=C(C=C(C=C2)N=[N+]=N)[N+](=O)[O-]
4 4 145926828 CSC1=NC=C(C(=N1)CN=[N+]=N)Br
In [11]:
#Отфилтруем нерабочее.
work_azids =[i for i in azids.CanonicalSMILES if len(i) < 30 and len(i) > 1 and '.' not in i]
len(work_azids), work_azids[0]
Out[11]:
(15441, 'C(CCl)N=[N+]=N')

Замена азидных групп в найденных молекулах на модифицированный ибупрофен.¶

Заменяем азидную группу на модифицированный ибупрофен, тех, что удовлетворяют правилам Липински, записываем в продукты.

In [12]:
products = []
for azid in work_azids:
    if "N=[N+]=[N-]" in azid:
        azid_ibu = azid.replace("N=[N+]=[N-]", ibu_for_azids)
        try: 
            newmol = Chem.MolFromSmiles(azid_ibu)
            if Lipinsky_checking(newmol):
                products.append(newmol)
        except Exception:
            pass
        
len(products)
Out[12]:
14555
In [13]:
for_pic = [products[i] for i in  np.random.randint(0, len(products)-1, size=(1,12)).ravel()]
Chem.Draw.MolsToGridImage(for_pic,useSVG=True, molsPerRow=3, subImgSize=(300, 300))
Out[13]:

Остатки присоединены не к тому атому, не похоже на формулу из википедии, надо переписать SMILES.

In [14]:
ibu_for_azids_re = 'N1N=NC(=C1)C(C(O)=O)C2=CC=C(C=C2)CC(C)C'
click_re = Chem.MolFromSmiles(ibu_for_azids_re)
AllChem.Compute2DCoords(click_re)
display(click_re)

То же соединение, но запись стартует с того атома, к которому я хочу присоединить остаток азида.

In [15]:
products = []
for azid in work_azids:
    if "N=[N+]=[N-]" in azid:
        azid_ibu = azid.replace("N=[N+]=[N-]", ibu_for_azids_re)
        try: 
            newmol = Chem.MolFromSmiles(azid_ibu)
            if Lipinsky_checking(newmol):
                products.append(newmol)
        except Exception:
            pass
        
len(products)
Out[15]:
14572
In [16]:
for_pic = [products[i] for i in  np.random.randint(0, len(products)-1, size=(1,12)).ravel()]
Chem.Draw.MolsToGridImage(for_pic,useSVG=True, molsPerRow=3, subImgSize=(300, 300))
Out[16]:
In [ ]:
Построим карту схожести пятого соединения и ибупрофена.
In [66]:
from rdkit.Chem.Draw import SimilarityMaps

fp = SimilarityMaps.GetMorganFingerprint(products[4], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, products[4],
                                                               SimilarityMaps.GetMorganFingerprint)
In [68]:
maxweight
Out[68]:
0.1983367983367983

Программа окрасила весь остаток ибупрофена в зелёный, т.е. распознала его.

А теперь в 3D.

In [28]:
m3d = Chem.AddHs(products[4])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200)
Out[28]:
1
In [19]:
import nglview as nv
nv.show_rdkit(m3d)
NGLWidget()

Работает! Вы в html этого не увидете, поэтому у меня есть кое-что получше чем 3D-структура, скриншоты 3D-структуры:

In [29]:
display(Image('one.jpg'))
In [30]:
display(Image('two.jpg'))
In [ ]: