Домашнее задание 3

На сайте PubChem найдем все радикалы c азидом для Click Chemistry и скачаем их SMILES нотации.

Формула ибупрофена с сайта pubchem: CC(C)CC1=CC=C(C=C1)C(C)C(=O)O Формула изопропила: C(3)H(8)O или CH(3)-CH-CH(3), а ацителена C(2)H_(2), где есть тройная связь между атомами углерода (с вики).

Тогда измененный ибупрофен можно получить следующей заменой: C#CCC1=CC=C(C=C1)C(C)C(=O)O , где # - обозначает тройную связь между C.

In [75]:
import rdkit
import pubchempy
import re
from IPython.display import Image
import nglview as nv
import numpy as np
from rdkit.Chem import rdmolfiles
from rdkit import Chem
from rdkit.Chem import Draw
import rdkit.Chem.Lipinski as Lipinski
from rdkit.Chem import AllChem
from IPython.display import display
from rdkit.Chem.Draw import IPythonConsole 

#from rdkit import Chem
#import pubchempy as pcp
#from rdkit.Chem import Draw, rdMolDescriptors, Descriptors, Crippen
#import itertools
In [12]:
ibuprofen=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
ibuprofen_modified=Chem.MolFromSmiles('C#CCC1=CC=C(C=C1)C(C)C(=O)O')
print ('ibuprofen:')
display(ibuprofen)
print ('ibuprofen_modified:')
display(ibuprofen_modified)

ibuprofen:

ibuprofen_modified:

Правила Лепински с wiki:

  • No more than 5 hydrogen bond donors (the total number of nitrogen–hydrogen and oxygen–hydrogen bonds)
  • No more than 10 hydrogen bond acceptors (all nitrogen or oxygen atoms)
  • A molecular mass less than 500 daltons
  • An octanol-water partition coefficient[5] log P not greater than 5
In [14]:
def Check_Lipinski(compound):
    if ((Lipinski.rdMolDescriptors.CalcExactMolWt(compound)<500) & (Lipinski.NumHDonors(compound) <= 5) & (Lipinski.NumHAcceptors(compound) <= 10) & (Lipinski.rdMolDescriptors.CalcCrippenDescriptors(compound)[0]<=5)):
        return True
    else:
        return False

Проверим, что ибупрофен и измененный ибупрофен удовлетворяют условиям Липински

In [15]:
print(Check_Lipinski(ibuprofen),Check_Lipinski(ibuprofen_modified))

(True, True)

Шаблон:

In [16]:
template="N2C=C(N=N2)CC1=CC=C(C=C1)C(C)C(=O)O"
display(Chem.MolFromSmiles(template))
In [22]:
#possible azids
smiles=['N=[N+]=[N-]','[N-][N+]#N','[N-]N=[N+]']
dat={}
molecules={}
for i,smile in enumerate(smiles):
    print(i)
    dat[i] = pubchempy.get_properties( 
    properties="CanonicalSMILES", identifier=smile, namespace="smiles", 
    searchtype="substructure", RingsNotEmbedded=True, as_dataframe=True )
    molecules[i] = np.array(dat[i]).flatten()

0
1
2
In [23]:
l=['C:/hw3/file1.txt','C:/hw3/file2.txt','C:/hw3/file3.txt']
for j in range(3):
    f=open(l[j],'w')
    for i in range(np.size(molecules[j])):
        f.write(molecules[j][i])
        f.write('\n')
    f.close()

Возможные виды написания азидов: N=[N+]=[N-],[N-][N+]#N,[N-]N=[N+] И насколько я понял необходимо, чтобы они 'торчали с одного из концов' соединения. Поэтому найденные соединения для каждого из написания азидов я сохранил в отдельные файлы и отфильтровал их следующей функцией:

In [24]:
def Filtration(filepath,type):
    answer=[]
    smiles_pattern=['N=\[N\+\]=\[N\-\]','\[N\-\]\[N\+\]\#N','\[N-\]N=\[N\+\]']
    smiles=['N=[N+]=[N-]','[N-][N+]#N','[N-]N=[N+]']
    p=re.compile(smiles_pattern[type])
    with open(filepath) as f:
        for line in f:
            line=line.strip()
            index=list(p.finditer(line))
            if (len(index)==1):
                if ((index[0].start()==0)|(index[0].start()==(len(line)-len(smiles[type]))))&(len(line) < 30) & ('.' not in line):
                    answer.append(line)
    return answer
In [25]:
list1=Filtration(l[0],0)
list2=Filtration(l[1],1)
list3=Filtration(l[2],2)

print (len(list1),len(list2),len(list3))

(8591, 466, 0)
In [26]:
def buildMolecules(smiles,azid):
    newMolecules=[]
    for smile in smiles:
        if azid in smile:
            new_smile = smile.replace(azid,template)
            try:
                new_mol = Chem.MolFromSmiles(new_smile)
                newMolecules.append(new_mol) if Check_Lipinski(new_mol) else None
            except:
                pass
        else:
            continue
    print('Number of molecules: ')
    print (len(newMolecules))
    return newMolecules

Построим теперь из этих smiles новые молекулы: третий лист не берем - в него ни один азид не попал.

In [27]:
newMolList1 = buildMolecules(list1,smiles[0])
newMolList2=buildMolecules(list2,smiles[1])

Number of molecules: 
8457
Number of molecules: 
455

Отобразим некоторые из молекул обоих списков:

In [28]:
for i in range (5):
    display(newMolList1[i])
In [29]:
for i in range (5):
    display(newMolList2[i])

Теперь попробуем представить какую-нибудь молекулу в 3D:

In [77]:
m = Chem.AddHs(newMolList2[1])
_ = AllChem.EmbedMultipleConfs(m, useExpTorsionAnglePrefs=True, useBasicKnowledge=True)
view = nv.show_rdkit(m)
view

Вроде бы писал все как в документации, но вывести изображение не получилось. Поэтому решил перевести smile в pdb формат и отобразить в pymol молекулу.

In [73]:
Chem.MolToSmiles(newMolList2[1])
display(newMolList2[1])
In [76]:
Image('C:/hw3/screen.png')
Out[76]: