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

#Загрузим нужные пакеты
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
import os

#Создаем объект rdkit из SMILES-формулы ибупрофена, рисуем.
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

#Добавляем тройную связь к молекуле ибупрофена, чтобы она могла вступать в реакцию азид-алкинового циклоприсоединения с азидами.
ibu_alkine=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(CC#C)C(=O)O')
AllChem.Compute2DCoords(ibu_alkine)
display(ibu_alkine)

#Шаблон для симуляции реакции циклоприсоединения.
template=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(CC1=CN=NN1)C(=O)O')
AllChem.Compute2DCoords(template)
display(template)

#Меняем нумерацию в формуле SMILES шаблона, чтобы получать формулу SMILES продукта реакции можно заменой азидной группы в SMILES формуле азида на SMILES формулу шаблона.
template=Chem.MolFromSmiles('N1N=NC=C(CC(C(=O)O)C2=CC=C(CC(C)C)C=C2)1')
AllChem.Compute2DCoords(template)
display(template)
template1 = '(N1N=NC=C(CC(C(=O)O)C2=CC=C(CC(C)C)C=C2)1)'

#Используем PubChemPy для поиска азидов в базе данных.
import pubchempy as pcp

structures = []
per_page = 10**5

for smiles in ["N=N=N", "NN#N", "N#NN"]:
    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:
            break
#         print("Retrieved page {} of {} search".format(i+1, smiles))
        structures.extend(a)

len(structures)

2390

#Соберем SMILES в список:
compounds = {str(i['CID']).strip():i['CanonicalSMILES'].strip() for i in structures}
smiles = list(compounds.values())

#Найдем 25 структур, которые после присоединения к ним модифицированного ибупрфена будут удовлетворять 5 условиям Лепински:
res = []
for smi in smiles:
    if len(res) > 25:
        break
    if '=N[N+]#N' in smi:
        newsmi=smi.replace('=N[N+]#N','N1C=C(N=N1)C(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):
            AllChem.Compute2DCoords(newmol)
            res.append(newmol)
            display(newmol)
    except:
        pass

#Отобразим полученные структуры в более удобном виде.
from IPython.display import display
display(Draw.MolsToGridImage(res, useSVG=True, molsPerRow=3, subImgSize=(300, 200)))

#Построим SimilarityMap ибупрофена с пятым веществом из нашего массива:
from rdkit.Chem.Draw import SimilarityMaps
ibu = Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
fp = SimilarityMaps.GetMorganFingerprint(res[5], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, res[5], SimilarityMaps.GetMorganFingerprint)

#Построим 3D структуру нашего лиганда:
m3d = Chem.AddHs(res[5])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
display(m3d)