Будем работать с белком лизоцимом, структуру которого построили на основе гомологичного моделирования на прошлом практикуме: LYS_MERLU. Используем пакеты Autodock Vina и oddt. Необходимо определить место докинга, удалить лиганд, запустить докинг и провести анализ.
import numpy as np
import copy
# Отображение структур
import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image
# Open Drug Discovery Toolkit
import oddt
import oddt.docking
import oddt.interactions
# Органика
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
import pmx # Модуль для манипулирования pdb
Подготовка белка:
pdb=pmx.Model('MERLU-lig.B99990001.pdb')
for r in pdb.residues[100:]:
print r #посмотрим концевые остатки
# создание объектов белок и лиганд
newpdb = pdb.copy()
lig = pdb.copy()
del newpdb.residues[122:]
del lig.residues[:122]
# находим координаты центра лиганда
coords= []
x = []
y = []
z = []
for a in lig.atoms:
# print a.x
coords.append(a.x) #добываем координаты
for c in coords:
x.append(c[0])
y.append(c[1])
z.append(c[2])
x_center = np.mean(x)
y_center = np.mean(y)
z_center = np.mean(z)
print x_center,y_center,z_center
newpdb.writePDB("onlyprot.pdb")
Подготовка белка для докинга:
prot = oddt.toolkit.readfile('pdb','onlyprot.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
print 'is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt
Нам говорят, что первая молекула в файле onlyprot.pdb - не белок. Но на самом деле в файле белок и все хорошо.
Лиганды для докинга.
NAG (N-ацетил-D-глюкозамин) содержит в себе СH3C(=O)NH группу. Будем проводить докинг с ним самим и с лигандами на его основе, где метильный радикал этой группы заменён на OH, NH3+, H, Ph, COO-.
smiles = ['CC(=O)NC1C(C(C(OC1O)CO)O)O', 'OC(=O)NC1C(C(C(OC1O)CO)O)O', '[NH3+]C(=O)NC1C(C(C(OC1O)CO)O)O',
'C(=O)NC1C(C(C(OC1O)CO)O)O', 'C1=CC=CC=C1C(=O)NC1C(C(C(OC1O)CO)O)O', "[O-]C(=O)C(=O)NC1C(C(C(OC1O)CO)O)O"]
mols= []
images =[]
for s in smiles:
m = oddt.toolkit.readstring('smi', s)
if not m.OBMol.Has3D():
m.make3D(forcefield='mmff94', steps=150)
m.removeh()
m.OBMol.AddPolarHydrogens()
mols.append(m)
###with print m.OBMol.Has3D() was found that:
### deep copy needed to keep 3D , write svg make mols flat
images.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*images)
Запускаем докинг.
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=[x_center,y_center,z_center],
executable='/usr/bin/vina',autocleanup=True, num_modes=20)
print dock_obj.tmp_dir
print " ".join(dock_obj.params)
Выдача:
координаты центра лиганда
размеры ячейки (с этим центром), в которой будет идти докинг
используем 1 ядро для работы
точность (полнота расчетов)
маскимальное число положений лиганда
максимальная разница энергии лучшего положения с худшим.
# do it
res = dock_obj.dock(mols,prot)
Результаты докинга:
for i,r in enumerate(res):
print "%4d%10s%8s%8s%8s" %(i,r.formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name)
Анализ докинга:
print "i ligand hbonds stacking hf contacts"
for i,r in enumerate(res):
hbs = oddt.interactions.hbonds(prot,r)
stack= oddt.interactions.pi_stacking(prot,r)
phob = oddt.interactions.hydrophobic_contacts(prot,r)
print "%s %8s%8s%8s%8s" % (i,r.formula,len(hbs[1]), len(stack[1]), len(phob[1]))
Посмотрим глазами:
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb')
Получили 53 pdb-файла с всеми положениями всех лигандов.
Место для всех лигандов, найденное докингом, отличается от первоначального (зеленая молекула это все три лиганда):
from IPython.display import Image
Image(filename='all.png')
В таблице показаны заместители от лучшего к худшему. Сортировали по vina_affinity, в случае одинаковых значений сравнивали по водородным связям.
№ | group | ligand | vina_affinity | h_bonds |
9 | OH | C7H13NO7 | -6.9 | 11 |
45 | COO- | C8H13NO8 | -6.7 | 7 |
0 | none | C8H15NO6 | -6.6 | 5 |
18 | NH3+ | C7H16N2O6 | -6.5 | 6 |
27 | H | C7H13NO6 | -6.5 | 10 |
36 | Ph | C13H17NO6 | -5.5 | 10 |