Цель данного занятия ознакомится с возможностями докинга низкомолекулярного лиганда в структуру белка. Мы работаем с белком лизоцимом структуру которого построили на основе гомологичного моделирования на прошлом практикуме.
Загружаем модули:
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('P50717_with_ligand.B99990001.pdb')
for r in pdb.residues[135:]:
print r #посмотрим остатки
1-142 остаток принадлежит белку, 143-145 остатки являются лигандом.
# создаем объект белка newpdb и лиганда lig
newpdb = pdb.copy()
lig = pdb.copy()
del newpdb.residues[142:]
del lig.residues[:142]
# находим геометрический центр лиганда
def average_in_list(lst):
return reduce(lambda x, y: x+y, lst) / len(lst)
x=[]
y=[]
z=[]
for a in lig.atoms:
x.append(a.x[0])
y.append(a.x[1])
z.append(a.x[2])
centre_x=average_in_list(x)
centre_y=average_in_list(y)
centre_z=average_in_list(z)
# создаем pdb файл белка без лиганда
newpdb.writePDB("protein.pdb")
Подготовка белка для докинга:
prot = oddt.toolkit.readfile('pdb','protein.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
Мой белок не распознан как белок. Молекулярная масса 16865.9877.
Создаем лиганды для докинга:
smiles = ['c1cccc(O)c1', 'c1c(O)ccc(O)c1','c1(O)cc(c2ccccc2)cc(O)c1']
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=(centre_x, centre_y, centre_z),
executable='/usr/bin/vina',autocleanup=True, num_modes=20)
print dock_obj.tmp_dir
print " ".join(dock_obj.params)
# do it
res = dock_obj.dock(mols,prot)
Описываем результаты докинга:
print " ".join(["№", "formula", "affinity", "rmsd_ub", "name"])
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)
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 "%8s%8s%8s" %(len(hbs[1]),len(stack[1]), len(phob[1]))
N-Ацетилглюкозамин (NAG) содержит в себе СH3C(=O)NH группу. Создаем лиганды, где метильный радикал этой группы будет заменён на: OH, NH3+, H, Ph, COO-.
smiles_2 = ['CC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O', 'OC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O',
'[NH3+]C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O', 'C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O',
'C1=CC=C(C=C1)OC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O', '[O-]C(=O)C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O']
name=["C8H15NO6", "C7H13NO7", "C7H15N2O6", "C7H13NO6", "C13H17NO7", "C8H13NO8"]
mols_2= []
images_2 =[]
for i in range(len(smiles_2)):
m = oddt.toolkit.readstring('smi', smiles_2[i])
m.title=name[i]
if not m.OBMol.Has3D():
m.make3D(forcefield='mmff94', steps=150)
m.removeh()
m.OBMol.AddPolarHydrogens()
mols_2.append(m)
###with print m.OBMol.Has3D() was found that:
### deep copy needed to keep 3D , write svg make mols flat
images_2.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*images_2)
Запускаем докинг:
#create docking object
dock_obj_2= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=(centre_x, centre_y, centre_z),
executable='/usr/bin/vina',autocleanup=True, num_modes=20)
print dock_obj_2.tmp_dir
print " ".join(dock_obj_2.params)
# do it
res_2 = dock_obj_2.dock(mols_2,prot)
A=[]
for i,r in enumerate(res_2):
A.append([])
line=[i, r.formula, float(r.data['vina_affinity']), float(r.data['vina_rmsd_ub']), r.residues[0].name, len(oddt.interactions.hbonds(prot,r)[1]), len(oddt.interactions.pi_stacking(prot,r)[1]), len(oddt.interactions.hydrophobic_contacts(prot,r)[1])]
A[i].extend(line)
A = sorted(A, key=lambda k: k[2]) #сортируем по аффинности
print '\t'.join(["№", "formula", "affinity", "rmsd_ub", "name", "hbond", "stacking", "hydrophobic_contacts"])
for row in A:
print('\t'.join(list(map(str, row))))
Для всех модификаций N-ацетилглюкозамина афинность получилась невысокая. Самой лучшей находкой стала та, где метильный радикал был заменен на NH3+ - она имеет наименьшую энергию связывания и всего 6 водородных связей. Все остальные находки имеют большее количество связей, хотя и меньшую аффинность. Гидрофобные взаимодействия наблюдается лишь в случае метильного радикала и его замены на фенил.