Задание 3¶
Предварительные ласки.
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
from tqdm import tqdm, trange
import rdkit.Chem.Lipinski as Lipinksy
def lipinsky(smiles):
return (Lipinksy.NumHDonors(smiles),
Lipinksy.NumHAcceptors(smiles),
Lipinksy.rdMolDescriptors.CalcExactMolWt(smiles),
Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(smiles)[0])
def is_lipinsky(smiles, ideal=(5, 10, 500, 5)):
return all(i < j for i, j in zip(lipinsky(smiles), ideal))
Создадим молекулу ибупрофена из SMILES и нарисуем её.
ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)
is_lipinsky(ibu)
Ибупрофен удовлетворяет правилу 5 Липински.
Модификация ибупрофена для клик-химии¶
Популярная реакция в клик-химии это азид-алкиновое циклоприсоединение.
display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/6/66/CuAAC.svg/768px-CuAAC.svg.png'))
Так как мы ищем азиды в PubChem, то нам нужно получить производное ибупрофена с ацетиленовой группой. Один из способов — это нуклеофильное замещение алкинидов по следующему механизму:
display(Image('https://upload.wikimedia.org/wikipedia/commons/f/fd/Alkynespreparation-3.png'))
Следовательно, нам надо ввести галоген в молекулу ибупрофена. С помощью реакции Гелля-Фольгарда-Зелинского можно галогенировать Cα-атом карбоновых кислот.
display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/4/4d/HVZReaction.png/600px-HVZReaction.png'))
Стоит отметить, что Cα-атом в ибупрофене хирален. Однако в медицине используется рацемат, поэтому нам не надо будет думать о стереспецифичности синтезов.
Посмотрим на хлор-производное ибупрофена.
cl_ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(Cl)(C)C(=O)O')
AllChem.Compute2DCoords(cl_ibu)
display(cl_ibu)
Затем на алкиновое производное.
al_ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C#C)(C)C(=O)O')
AllChem.Compute2DCoords(al_ibu)
display(al_ibu)
И на клик-производное.
al_ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C1=CN=NN1)(C)C(=O)O')
AllChem.Compute2DCoords(al_ibu)
display(al_ibu)
Теперь перепишем SMILES, чтобы производное было удобно использовать в качестве шаблона для модификцаии азидов.
template = Chem.MolFromSmiles('N1N=NC(=C1)C(C)(C(=O)O)C1=CC=C(C=C1)CC(C)C')
AllChem.Compute2DCoords(template)
display(template)
Азиды в PubChem¶
Теперь найдём все соединения, содержащие азид, в PubChem.
import pubchempy as pcp
compounds = []
per_page = 10**5
for smiles in tqdm(["N=N=N", "NN#N", "N=[N+]=[N-]"]):
for i in trange(200):
try:
a = pcp.get_properties(properties="CanonicalSMILES",
identifier=smiles, namespace="smiles",
searchtype="substructure",
RingsNotEmbedded=False,
listkey_count=per_page,
listkey_start=i * per_page)
except Exception as e:
tqdm.write(e)
tqdm.write("Retrieved page {} of {} search".format(i+1, smiles))
compounds.extend(a)
Всё очень плохо с API, поэтому будем искать ручками. PubChem CID для азида 33558.
Скачали из PubChem (substructure search, download as SMILES).
with open('azide_search.txt', 'r') as infile:
compounds = [smiles
for smiles in map(lambda i: i.strip().split('\t')[1], infile)
if len(smiles) < 30 and len(smiles) > 1 and '.' not in smiles]
len(compounds)
Теперь заменим азиды на производное ибупрофена.
Азид может быть терминальным и нетерминальным, в некоторых SMILES Это учитывается.
azides = [Chem.MolFromSmiles('N=[N+]=[N-]'),
Chem.MolFromSmiles('N=[N+]=[NH]')]
template
broken_count = 0
unlinked_count = 0
modified_compounds = set()
for compound in compounds:
mol = Chem.MolFromSmiles(compound)
if mol is None: # some SMILES are broken
broken_count += 1
continue
modified = False
for azide in azides:
if mol.HasSubstructMatch(azide):
mol = Chem.rdmolops.ReplaceSubstructs(mol=mol,
query=azide,
replacement=template,
replaceAll=True)[0]
modified = True
else:
continue
if modified:
modified_compound = Chem.MolToSmiles(mol)
if '.' not in modified_compound: # some molecules break up upon modification
modified_compounds.add(modified_compound)
else:
unlinked_count += 1
print(broken_count, unlinked_count, len(modified_compounds))
Осталось только 10% молекул от исходных азидов. Видимо, как-то плохо искали. Теперь мы отфильтруем полученные молекулы по Липински.
final_compounds = [mol
for mol in map(lambda i: Chem.MolFromSmiles(i), modified_compounds)
if mol is not None and is_lipinsky(mol)]
len(final_compounds)
Потеряли все молекулы. Видимо, плохо искали в PubChem. Возьмём файл от однокурсников.
with open('azids.csv', 'r') as infile:
compounds = [smiles
for smiles in map(lambda i: i.strip().split(',')[2], infile)
if len(smiles) < 30 and len(smiles) > 1 and '.' not in smiles]
len(compounds)
Молекул в 2 раза меньше, чем было в первый раз.
broken_count = 0
unlinked_count = 0
modified_compounds = set()
for compound in compounds:
mol = Chem.MolFromSmiles(compound)
if mol is None: # some SMILES are broken
broken_count += 1
continue
modified = False
for azide in azides:
if mol.HasSubstructMatch(azide):
mol = Chem.rdmolops.ReplaceSubstructs(mol=mol,
query=azide,
replacement=template,
replaceAll=True)[0]
modified = True
else:
continue
if modified:
modified_compound = Chem.MolToSmiles(mol)
if '.' not in modified_compound: # some molecules break up upon modification
modified_compounds.add(modified_compound)
else:
unlinked_count += 1
print(broken_count, unlinked_count, len(modified_compounds))
final_compounds = [mol
for mol in map(lambda i: Chem.MolFromSmiles(i), modified_compounds)
if mol is not None and is_lipinsky(mol)]
len(final_compounds)
display(Draw.MolsToGridImage(final_compounds, maxMols=9, molsPerRow=3, subImgSize=(200, 200)))
Всё очень плохо. Видимо, слишком умный способ ломает SMILES и ничего не получается. Попробуем по-простому.
final_compounds = []
for compound in compounds:
if "N=[N+]=[N-]" in compound:
compound = compound.replace("N=[N+]=[N-]", 'N1N=NC(=C1)C(C)(C(=O)O)C1=CC=C(C=C1)CC(C)C')
else:
continue
try:
mol = Chem.MolFromSmiles(compound)
if is_lipinsky(mol):
final_compounds.append(mol)
except:
pass
len(final_compounds)
Зато почти все сохранились после модификации.
display(Draw.MolsToGridImage(final_compounds, maxMols=9, molsPerRow=3, subImgSize=(400, 400)))
Построим карты схожестей производных с ибупрофеном.
from rdkit.Chem.Draw import SimilarityMaps
similarity_data = [SimilarityMaps.GetSimilarityMapForFingerprint(ibu, mol, SimilarityMaps.GetMorganFingerprint)
for mol in final_compounds[:5]]
list(zip(*similarity_data))[1]
Веса схожести небольшие, кажется. Странно, что часть исходной молекулы ибупрофена покрашена розовым, а не зелёным, мол, не очень похожа на исходную молекулу.
3D-изображение молекул¶
m3d = Chem.AddHs(final_compounds[0])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d, maxIters=500, nonBondedThresh=200)
display(m3d)
Жизнь потрепала молекулу.
import nglview as nv
nv.show_rdkit(m3d)
Прикольно получилось. А главное, удалось завести nglview в jupyter lab. Но в html оно не работает((
view = nv.show_rdkit(m3d)
view
view.render_image()
view._display_image()
