Цель данного занятия - ознакомиться с возможностями докинга низкомолекулярного лиганда в структуру белка. Белок - предсказанный в прошлом практикуме с помощью гомологичного моделирования. Лиганд - измененный NAG, где метильный радикал СH3C(=O)NH группы заменен на:
- OH
- NH3+
- H
- Ph
- COO- </ul>
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('LYZ_lig.B99990001.pdb')
for r in pdb.residues[:137]:
print r #посмотрим остатки
# создание объектов белок и лиганда
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-3]
for r in newpdb.residues:
print r #посмотрим остатки
print
for r in lig.residues:
print r #посмотрим остатки
help(pmx.Atom)
x_list = []
x_coord = []
y_coord = []
z_coord = []
for a in lig.atoms:
x_list.append(a.x)# найдите геометрический центр лиганда
for coord in x_list:
x_coord.append(coord[0])
y_coord.append(coord[1])
z_coord.append(coord[2])
center = (np.mean(x_coord),np.mean(y_coord),np.mean(z_coord))
print center
newpdb.writePDB('protein_no_lig.pdb')
Подготовим белок для докинга. Пишет, что в PDB файле первая молекула - не белок, но с массой 15274. Но если посмотреть собственно в файл, то там координаты аминокислот. Не верим выдаче, но файл проверить стоит.
prot = oddt.toolkit.readfile('pdb','/tmp/prot.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
print 'is the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt
Подготовим лиганд
#smiles = ['c1cccc(O)c1', 'c1c(O)ccc(O)c1','c1(O)cc(c2ccccc2)cc(O)c1']
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=center,
executable='/usr/bin/vina',autocleanup=True, num_modes=20)
print dock_obj.tmp_dir #это папка докинга
print " ".join(dock_obj.params)
Вначале идут координаты центра лиганда, затем размеры grid-box, 1 ядро, точность расчетов, количество положений лиганда в выдаче, максимальная разница энергии с наилучшей находкой.
res = dock_obj.dock(mols,prot)
Результаты докинга¶
Покажем заместители от лучшего по энергии к худшему:
#print '#\tformula\tenergy\trmsd\tprot_hbonds\tlig_hbonds\tprot_stack\tlig_stack\tprot_phob\tlig_phob'
res_energy = {}
res_dict = {}
print '#\tenergy\tformula\t\trmsd\thbonds\tstack\tphob'
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)
res_energy[i] = r.data['vina_affinity']
res_dict[i] = [r.formula, r.data['vina_rmsd_ub'],str(len(hbs[0])),str(len(stack[0])),str(len(phob[0]))]
# print '\t'.join([str(i),r.formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'],str(len(hbs[0])),str(len(stack[0])),str(len(phob[0]))])
# print '\t'.join([str(i),r.formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'],str(len(hbs[0])),str(len(hbs[1])),str(len(stack[0])),str(len(stack[1])),str(len(phob[0])),str(len(phob[1]))])
for k, v in sorted(res_energy.items(), key=lambda x:x[1],reverse=True):
print '\t'.join([str(k)]+[str(v)]+res_dict[k])
Визуализация:
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb', overwrite=True)
Лиганд №9 (жёлтый),с наименьшей энергией, располагается приблизительно в том же месте, что и NAG (синий). А вот №35 (зелёный), с наибольшей энергией располагается совсем не там и слишком близко расположен с остову белка.
from IPython.display import Image
Image(filename='result.png')
