Хемоинформатика¶
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
# Draw molecules from SMILES
def chemdraw(smiles):
mol = Chem.MolFromSmiles(smiles)
AllChem.Compute2DCoords(mol)
display(mol)
ibu = 'CC(C)CC1=CC=C(C=C1)C(C)C(=O)O'
chemdraw(ibu)
Посчитаем параметры для "правила Липинского", согласно которому соединение, чтобы «быть похожим» на лекарство, должно:
• Иметь менее 5 атомов-доноров водородной связи;
• Обладать молекулярным весом менее 500;
• Иметь липофильность (log P — коэффициент распределения вещества на границе раздела вода-октанол) менее 5;
• Иметь суммарно не более 10 атомов азота и кислорода (грубая оценка количества акцепторов водородной связи)
В связи с тем, что числа, указанные в вышеперечисленных пунктах, делятся нацело на 5, данное правило также называется "правилом пяти".
import rdkit.Chem.Lipinski as Lipinksy
ibu_mol = Chem.MolFromSmiles(ibu)
print Lipinksy.NumHDonors(ibu_mol)
print Lipinksy.NumHAcceptors(ibu_mol)
print Lipinksy.rdMolDescriptors.CalcExactMolWt(ibu_mol)
print Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu_mol)[0]
Одна из наиболее используемых клик-реакций — азид-алкиновое циклоприсоединение (реакция между азидами и алкинами с образованием 1,2,3-триазолов) с использованием медного катализатора (CuAAC: Cu-catalyzed azide-alkyne cycloaddition). Каталитический вариант реакции не протекает синхронно, а имеет постадийный механизм, который протекает через промежуточное образование ацетиленидов меди. По этой причине высокую реакционную способность в данной реакции демонстрируют лишь терминальные алкины. Поэтому для эмуляции продукта Click Chemistry SMILES ибупрофена нужно изменять так, чтобы в молекуле обязательна присутствовала концевая тройная связь.
# SMILES модифицированного ибупрофена
ibu1 = 'C#CCC1=CC=C(C=C1)C(C)C(=O)O'
ibu2 = 'CC(C)CC1=CC=C(C=C1)C(C#C)C(=O)O'
ibu3 = 'CC(C#C)CC1=CC=C(C=C1)C(C)C(=O)O'
chemdraw(ibu1)
chemdraw(ibu2)
chemdraw(ibu3)
Так как существуют три варианта модификации ибцпрофена (три места присоединения тройной связи), то и рекция азид-алкинового циклоприсоедениение может идти по трем положениям.
# Продукты реакции азид-алкиновое циклоприсоединение
tr1 = 'N1N=NC(=C1)C1=CC(CC(C)C)=CC=C1(C(C)C(=O)O)'
tr2 = 'N1N=NC(=C1)C(C(C)(C))C1=CC=C(C=C1)C(C)C(=O)O'
tr3 = 'N1N=NC(=C1)CC(C(=O)O)C1=CC=C(C=C1)CC(C)C'
chemdraw(tr1)
chemdraw(tr2)
chemdraw(tr3)
Cписок азидов можно скачать вручную или с помощью https://pubchempy.readthedocs.io/en/latest/:
# Загрузим скаченные данные и отфильтруем (SMILES-нотации должны содержать менее 30 символов и не иметь разрывов)
strings = np.genfromtxt('pubchem.txt', dtype = np.str)
smiles = []
for line in strings:
if len(line[1]) < 30 and not '.' in line[1]:
smiles.append(line[1])
print (len(smiles))
Коэффициент распределения в системе октанол-вода (logP) характеризует способность транспорта лекарственных веществ через клеточные мембраны, определяя их абсорбцию и распределение в различных системах организма. Отрицательное значение logP означает, что соединение имеет более высокую близость к водной фазе (более гидрофильно); когда logP = 0, то соединение одинаково распределено между липидной и водной фазами; положительное значение logP означает более высокую концентрацию соединения в липидной фазе (более липофильно). Поэтому отфильтруем модификации по logP (более водорастворимые соединения будут лучше выводиться из организма) и покажем в конце лучшую.
# Пстроим новые молекулы, заменив в отфильтрованных радикалах азидную группу на модифцированный ибупрофен, для первых 1500 соединений
new_smi = []
N = 0
all_mol = 1500
templates = (tr1, tr2, tr3)
min_mol = Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu_mol)[0]
min_ind = 0
# Новую молекулу лучше создавать в try из-за битых Smiles
for smi in smiles[:all_mol]:
if 'N=[N+]=[N-]' in smi:
for t in range(3):
template = templates[t]
#new_smi.append(smi.replace('N=[N+]=[N-]', template))
newsmi = smi.replace('N=[N+]=[N-]', template)
else:
continue
# Если новая молекула удолтворяет правилу пяти Липинского, сохраним в массив и покажем
#for n_smi in new_smi:
try:
new_mol=Chem.MolFromSmiles(newsmi)
if (Lipinksy.NumHDonors(new_mol) <= 5 and
Lipinksy.NumHAcceptors(new_mol) <= 10 and
Lipinksy.rdMolDescriptors.CalcExactMolWt(new_mol) <=500 and
Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(new_mol)[0] <= min_mol):
new_smi.append(new_mol)
# Display all molecules
#AllChem.Compute2DCoords(new_mol)
#display(new_mol)
# Display index of each molecule
#print N
N = N + 1
if Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(new_mol)[0] < min_mol:
min = Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(new_mol)[0]
min_ind = N - 1
except:
pass
print(' ')
print ('The amount of good structures: %s' % (len(new_smi)))
print ('The best hit is number %i with the minimal value of logP %i' % (min_ind, min_mol))
display(new_smi[min_ind])
from IPython.display import SVG
img = Chem.Draw.MolsToGridImage(new_smi[-15:], molsPerRow=3, subImgSize=(300, 300))
img
Similiraty Map¶
Построим Similiraty Map ибупрофена с пятым веществом из массива
from rdkit.Chem.Draw import SimilarityMaps
somemol = new_smi[4]
fp = SimilarityMaps.GetMorganFingerprint(somemol, fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu_mol, somemol, SimilarityMaps.GetMorganFingerprint)
3D-структура¶
m3d = Chem.AddHs(somemol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d, maxIters=500, nonBondedThresh=200)
m3d
# К сожалению, данная команда не сохраняет изображение в html-странице
import nglview as nv
view = nv.show_rdkit(m3d)
view
Создадим pdb-файл молекулы
pdb_mol = Chem.MolToPDBBlock(m3d)
print(pdb_mol)
Image(filename='pdb_mol.png')
Ссылки¶
- Хемоинформатика: RDKit, Open Babel
https://github.com/rdkit/rdkit
http://openbabel.org/wiki/Main_Page