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 Lipinski
Небольшая функция для удобства
def chemdraw(smiles):
mol = Chem.MolFromSmiles(smiles)
AllChem.Compute2DCoords(mol)
display(mol)
Нарисуем ибупрофен
chemdraw('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
Можно предложить множество позиций для введения алкильного радикала для осуществления затем комбинаторного синтеза с помощью клик-химии. Ниже я привел три примера.
chemdraw('CC(C)C(C#C)C1=CC=C(C=C1)C(C)C(=O)O')
chemdraw('CC(C)CC1=CC=C(C=C1)C(C(C#C))C(=O)O')
chemdraw('CC(C)CC1=CC=C(C(C#C)=C1)C(C)C(=O)O')
Теперь нарисуем второго участника реакции — какое-либо вещество с азидом N=[N+]=[N-]
chemdraw('C1C(COC1=O)CN=[N+]=[N-]')
Перестроим SMILES модификаций ибупрофена, чтобы мимикрировать результат клик-реакции. Атом азота, к которому присоединяется радикал при азиде, должен стоять первым для удобства дальнейшей манипуляции.
template_1 = 'N1N=NC(=C1)C(C(C)(C))C1=CC=C(C=C1)C(C)C(=O)O'
chemdraw(template_1)
template_2 = 'N1N=NC(=C1)CC(C(=O)O)C1=CC=C(C=C1)CC(C)C'
chemdraw(template_2)
template_3 = 'N1N=NC(=C1)C1=CC(CC(C)C)=CC=C1(C(C)C(=O)O)'
chemdraw(template_3)
Сделав так, мы можем "осуществить" клие-реакцию простой функцией:
def click(ibu,azide):
return azide.replace('N=[N+]=[N-]',ibu)
Протестируем на уже знакомом нам азиде
test = click(template_1, 'C1C(COC1=O)CN=[N+]=[N-]')
chemdraw(test)
Путем структурного поиска в PubChem были найдены вещества, содержащие азид, и их SMILES сохранены в файл. Отберем только небольшие азид-содержащие вещества:
smiles = []
with open('pubchem.txt','r') as pchem:
for line in pchem:
tokens = line.split()
smile = tokens[1]
if len(smile) < 30 and '.' not in smile: # we exclude '.': '.' means 2+ separated ions/molecules
smiles.append(smile)
Сколько веществ осталось в выборке?
len(smiles)
Проведем клик-реакцию для всех азидов против трех алкил-содержащих ибупрофенов и отсеим только те производные, что отвечают правилу 5-ти Липински (я использовал его расширенную формулировку, включающую еще число способных вращаться связей и индекс рефракции. Есть и еще более расширенные варианты)
hits = np.zeros((3*len(smiles),8), dtype=np.float16)
allgood = []
templates = (template_1, template_2, template_3)
good_k = 0
for smile in smiles:
if 'N=[N+]=[N-]' in smile:
for i in range(3):
template = templates[i]
hit_smile = click(template,smile)
try:
hit = Chem.MolFromSmiles(hit_smile)
hit_hbdon = Lipinski.NumHDonors(hit)
hit_hbacc = Lipinski.NumHAcceptors(hit)
hit_mw = Lipinski.rdMolDescriptors.CalcExactMolWt(hit)
hit_nrotb = Lipinski.NumRotatableBonds(hit)
hit_cripp = Lipinski.rdMolDescriptors.CalcCrippenDescriptors(hit)
hit_logp = hit_cripp[0]
hit_refr = hit_cripp[1]
# Lipinsky rule of 5 (extended with rotatable bonds and refraction coefficient)
if hit_hbdon <= 5 and hit_hbacc <= 10 and hit_mw <= 500 and hit_nrotb <= 10 and hit_logp <= 5 and (abs(hit_refr - 85) <= 45):
allgood.append(hit)
hits[good_k,:] = [good_k,i,hit_hbdon,
hit_hbacc,hit_mw,hit_nrotb,
hit_logp,hit_refr]
good_k += 1
except:
continue
hits = hits[:good_k,:]
len(allgood) == len(hits)
Возьмем только лучшие 50 по параметру logP
top50idx = []
for c,hit in enumerate(sorted(hits, key = lambda x: x[6])):
if c < 50:
print('%d\tlogP=%.3f' % (c,hit[6]))
idx = int(hit[0])
top50idx.append(idx)
mol = allgood[idx]
AllChem.Compute2DCoords(mol)
display(mol)
Попадаются дубликаты, например, 0, 2, 3. Видимо, это следствие избыточности исходной выборки из PubChem (либо хиральность, не отражаемая в 2D). Любопытно, что для этих веществ вычисленные значения logP различны.
def to3d(m):
m3d=Chem.AddHs(m)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200)
return m3d
top20mols = [allgood[x] for x in top50idx[:20]]
Chem.Draw.MolsToGridImage(top20mols, molsPerRow=4, subImgSize=(200,150))
top20_3d = list(map(to3d,top20mols))
import nglview
import ipywidgets
Некоторые варианты в 3D
nglview.show_rdkit(top20_3d[0])
nglview.show_rdkit(top20_3d[7])
nglview.show_rdkit(top20_3d[17])