In [1]:
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 rdkit.Chem.Lipinski as Lipinksi
import pubchempy as pcp
import nglview as nv
from rdkit.Chem.Draw import SimilarityMaps
from IPython.display import SVG
from itertools import chain
Наш герой -- ибупрофен!¶
Цель практикума -- предложить аналог ибупрофена, если бы он был продуктом клик-химии
In [67]:
ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)
In [14]:
def mol_draw(smiles):
mol=Chem.MolFromSmiles(smiles)
AllChem.Compute2DCoords(mol)
display(mol)
Молекула маленькая, "легкая", за счет кабоксильной группы может образовывать водородную связь как донор или акцептор. Молекула ибупрофена удовлетворяет правилу 5 Липински:
In [15]:
print(f"Количество атомов-доноров водородной связи: {Lipinksi.NumHDonors(ibu)} <5")
print(f"Количество атомов-акцепторов водородной связи: {Lipinksi.NumHAcceptors(ibu)} <10")
print(f"Молекулярная масса: {Lipinksi.rdMolDescriptors.CalcExactMolWt(ibu):.2f} Да <500 Да")
print(f"log P: {Lipinksi.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0]:.2f} <5")
Количество атомов-доноров водородной связи: 1 <5 Количество атомов-акцепторов водородной связи: 1 <10 Молекулярная масса: 206.13 Да <500 Да log P: 3.07 <5
Зададим функцию, которая оценивает удовлетворяет ли молекула правилу 5:
In [68]:
def lipinsky_istrue(mol):
a = (Lipinksi.NumHDonors(mol)<5) and \
(Lipinksi.NumHAcceptors(mol) < 10) and \
(Lipinksi.rdMolDescriptors.CalcExactMolWt(mol)) <500 and \
(Lipinksi.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0] <5)
if a:
return True
return False
lipinsky_istrue(ibu) # ибупрофен удовлетворяет правилу Липински
Out[68]:
True
Проэмулируем клик-химию на ибупрофене. Для этого ибупрофену нужен алкин. Если у нас уже есть ибупрофен, то проще всего провести реакцию по карбоксильной группе (с пропаргиламином или пропаргиловым спиртом например). Продукт реакции ибупрофена с пропаргиловым спиртом:
In [16]:
mol_draw('CC(C)Cc1ccc(cc1)C(C)C(=O)OCC#C')
Добавим азид:
In [33]:
mol_draw('CC(C)Cc1ccc(cc1)C(C)C(=O)OCC1CNN=N1')
In [45]:
r_smarts = 'CC(C)Cc1ccc(cc1)C(C)C(=O)OCC#C>>CC(C)Cc1ccc(cc1)C(C)C(=O)OCC1CNN=N1'
r = AllChem.ReactionFromSmarts(r_smarts, useSmiles=True)
r
Out[45]:
Ну, примерно такой продукт клик-химии должен получиться, если у азида нет никаких радикалов
In [34]:
def add_azid(azid):
return f'CC(C)Cc1ccc(cc1)C(C)C(=O)OCC1CN({azid})N=N1'
Поиск азидов в базе pubchempy¶
Вот такие паттерны мы будем искать:
In [3]:
groups = ['N=[N+]=[N-]', '[N-]=[N+]=[N-]', '[N-][N+]#[N]', '[N]#[N+][N-]', '[N-][N]#[N+]', '[N+]#[N][N-]']
In [ ]:
l = []
per_page = 10**2 # снизила тк выдавал ошибку типо сервер перегружен
for smiles in groups:
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
)
with open("./output_pcp.txt", "a") as file:
print(*[i["ConnectivitySMILES"] for i in a], sep="\n", file=file)
l.extend(a)
except Exception as e:
print(f"{e}")
# Такая например
# PUGREST.ServerError: GetOperationStatus fault: Server, message: NCBI C++ Exception
print(f"Retrieved page {i+1} of {smiles} search")
Откроем полученые smiles:
In [89]:
smiles = []
with open('./output_pcp.txt','r') as file:
smiles_raw = [i.strip() for i in file.readlines()]
l_filter = lambda smile: (len(smile) < 30) and ('.' not in smile)
smiles = list(filter(l_filter, smiles_raw))
В итоге из-за ошибок и снижения числа запросов скачалось около 5к smiles. Это немного, но поскольку цели учебные, не буду перегружать сервер
In [90]:
len(smiles)
Out[90]:
5587
In [91]:
smiles[0]
Out[91]:
'C(COCCOCCOCC=O)N=[N+]=N'
Получение клик-продуктов и фильтрация¶
Теперь получим продукты клик-химии. Для этого я решила воспользоваться возможностью rdkit: генерировать продукт по smarts-шаблону (https://github.com/rdkit/rdkit-tutorials/blob/master/notebooks/003_SMARTS_ReactionsExamples.ipynb)
In [120]:
copper_click_smarts = '''
[CH:1]#[C:2].[N:3]=[N+:4]=[N-:5]>>
[c:1]1[c:2][n-0:3][n-0:4][n-0:5]1.
[c:2]1[c:1][n-0:3][n-0:4][n-0:5]1
'''
copper_click = AllChem.ReactionFromSmarts(copper_click_smarts)
copper_click
Out[120]:
In [126]:
template
Out[126]:
In [149]:
result = []
template = Chem.MolFromSmiles("CC(C)Cc1ccc(cc1)C(C)C(=O)OCC#C")
num_to_check = 1500
for smi in smiles[:num_to_check]:
for group in groups:
if group in smi:
#newsmi = smi.replace(group, template)
# R = smi.replace(group, "")
# newsmi = f"CC(C)Cc1ccc(cc1)C(C)C(=O)OCC1CN{R}N=N1"
try:
radical = Chem.MolFromSmiles(smi)
products_tuples = copper_click.RunReactants((template, radical))
products = list(chain.from_iterable(products_tuples))
[Chem.SanitizeMol(prod) for prod in products]
mol = products[0] #возьму только 1 вариант...
if lipinsky_istrue(mol):
result.append(mol)
break
except:
print(f"FAIL. Smiles {smi}")
continue
[18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted [18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted [18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted
FAIL. Smiles COCl(N=[N+]=[N-])O FAIL. Smiles COCl(C#N)N=[N+]=[N-] FAIL. Smiles CO[Cl-]N=[N+]=[N-] FAIL. Smiles COCl(N=[N+]=[N-])O FAIL. Smiles COCl(C#N)N=[N+]=[N-] FAIL. Smiles CO[Cl-]N=[N+]=[N-]
[18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted [18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted [18:01:36] Explicit valence for atom # 2 Cl, 3, is greater than permitted
Какие-то беды в записях соединений с хлором
Визуализация¶
Вот такие молекулы получились:
In [151]:
Draw.MolsToGridImage(result[:12],useSVG=True, molsPerRow=3, subImgSize=(250, 200))
Out[151]:
Больше всего впечатляет соединение азид - мышьяк - азид.
In [170]:
fp = SimilarityMaps.GetMorganFingerprint(result[0], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, result[4], SimilarityMaps.GetMorganFingerprint)
В 3D:
In [171]:
m3d=Chem.AddHs(result[4])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
nv.show_rdkit(m3d)
NGLWidget()
In [172]:
display(m3d)
