Докинг низкомолекулярных лигандов в структуру белка

В этом задании будем осуществлять докинг в структуру, модель которой делали в прошлом практикуме.

In [3]:
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 

Для работы возьмем из прошлого практикума по гомологичному моделированию структуру с наибольшим скором.

In [4]:
pdb=pmx.Model('MERLU.B99990002.pdb')

for r in pdb.residues[120:]:
    print r #посмотрим остатки, чтобы найти лиганд
<Molecule: id = 121 name = VAL chain_id =    natoms = 7>
<Molecule: id = 122 name = GLN chain_id =    natoms = 10>
<Molecule: id = 123 name = NAG chain_id =    natoms = 14>
<Molecule: id = 124 name = NAG chain_id =    natoms = 14>
<Molecule: id = 125 name = NDG chain_id =    natoms = 15>
In [5]:
# разделяем на отдельные объекты белок и лиганд
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-3]
In [6]:
# ищем геометрический центр лиганда
x = []
y = []
z = []
for a in lig.atoms:
    coord = a.x 
    x.append(coord[0])
    y.append(coord[1])
    z.append(coord[2])
    
geom_center = [np.mean(x), np.mean(y), np.mean(z)]
geom_center
Out[6]:
[36.95716511867905, 38.747262125902985, 19.413063983488129]
In [7]:
newpdb.writePDB("myprot.pdb")
lig.writePDB("lig.pdb")

Подготовка белка для докинга

In [10]:
prot = oddt.toolkit.readfile('pdb','myprot.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
is it the first mol in 1lmp is protein? False :) and MW of this mol is: 13377.411

Странно, но белок не распознается как белок.

Лиганды для докинга

In [11]:
smiles = ['[NH3+]C(=O)NC1C(C(C(OC1O)CO)O)O', 'C(=O)NC1C(C(C(OC1O)CO)O)O',
          '[OH]C(=O)NC1C(C(C(OC1O)CO)O)O', 'C1=CC=C(C=C1)C(=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)
***** - Open Babel Depiction H H H H H H H H O O O N O O N O *****
***** - Open Babel Depiction O H H H H H O O O O N O *****
***** - Open Babel Depiction O H H H H H H O O O O O N O *****
***** - Open Babel Depiction H H H H H O O O O O N O *****
***** - Open Babel Depiction H H H H H O O O O O O N O HO *****

Докинг

In [12]:
#create docking object
#в качестве центра укажем геометрический центр лиганда"
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=geom_center,
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params)
/tmp/autodock_vina_qdVhmz
--center_x 36.9571651187 --center_y 38.7472621259 --center_z 19.4130639835 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

В выдаче указаны координаты геометрического центра, размер объекта и технические характеристики.

In [13]:
# do it
res = dock_obj.dock(mols,prot)

Результаты докинга

In [16]:
import pandas as pd
In [35]:
hbs_total = []
hbs_strict = []
stack = []
phob = []
formulas = []
aff = []
rmsd = []
for i,r in enumerate(res):
    hbs_total.append(len(oddt.interactions.hbonds(prot,r)[0]))
    hbs_strict.append(oddt.interactions.hbonds(prot,r)[2].sum())
    stack.append(len(oddt.interactions.pi_stacking(prot,r)[0]))
    phob.append(len(oddt.interactions.hydrophobic_contacts(prot,r)[0]))
    formulas.append(r.formula)
    aff.append(r.data['vina_affinity'])
    rmsd.append(r.data['vina_rmsd_ub'])
    
restable = pd.DataFrame({"Ligand_formulae": formulas, "Affinity": aff, "RMSD": rmsd, "H-bonds (total)": hbs_total, 
                         "H-bonds (strict)": hbs_strict, "Stacking": stack, "Hydrophobic_contacts": phob})
    
In [39]:
restable.sort_values(by = ['Affinity'], ascending=False)
Out[39]:
Affinity H-bonds (strict) H-bonds (total) Hydrophobic_contacts Ligand_formulae RMSD Stacking
27 -7.4 6 6 10 C13H17NO6 0.000 0
28 -7.2 9 10 5 C13H17NO6 7.440 1
29 -7.0 2 4 2 C13H17NO6 6.983 1
30 -7.0 2 5 5 C13H17NO6 2.396 0
31 -6.9 6 6 6 C13H17NO6 4.559 0
32 -6.9 2 4 4 C13H17NO6 7.209 1
36 -6.9 7 11 0 C8H13NO8 0.000 0
37 -6.8 8 8 0 C8H13NO8 5.017 0
33 -6.7 9 10 1 C13H17NO6 6.230 0
38 -6.7 10 11 0 C8H13NO8 4.146 0
18 -6.6 9 15 0 C7H13NO7 0.000 0
1 -6.6 5 8 0 C7H15N2O6 3.078 0
0 -6.6 7 10 0 C7H15N2O6 0.000 0
35 -6.5 5 5 5 C13H17NO6 3.693 0
34 -6.5 2 2 7 C13H17NO6 2.204 0
39 -6.4 5 7 0 C8H13NO8 3.377 0
40 -6.4 13 14 0 C8H13NO8 4.891 0
42 -6.3 5 5 0 C8H13NO8 4.285 0
41 -6.3 11 12 0 C8H13NO8 5.497 0
19 -6.3 6 9 0 C7H13NO7 3.008 0
3 -6.2 6 6 0 C7H15N2O6 4.558 0
2 -6.2 1 3 0 C7H15N2O6 5.035 0
9 -6.2 6 7 0 C7H13NO6 0.000 0
43 -6.1 3 6 0 C8H13NO8 3.736 0
22 -6.1 6 7 0 C7H13NO7 4.580 0
44 -6.1 5 9 0 C8H13NO8 5.573 0
21 -6.1 10 12 0 C7H13NO7 3.783 0
20 -6.1 4 6 0 C7H13NO7 4.857 0
4 -6.1 5 6 0 C7H15N2O6 4.924 0
5 -6.1 7 10 0 C7H15N2O6 4.571 0
6 -6.1 5 5 0 C7H15N2O6 2.564 0
7 -6.1 2 3 0 C7H15N2O6 5.056 0
10 -6.1 9 12 0 C7H13NO6 4.725 0
11 -6.0 7 7 0 C7H13NO6 3.816 0
24 -6.0 7 7 0 C7H13NO7 4.894 0
23 -6.0 9 10 0 C7H13NO7 2.229 0
8 -5.9 6 6 0 C7H15N2O6 2.461 0
25 -5.9 5 7 0 C7H13NO7 5.061 0
12 -5.9 8 11 0 C7H13NO6 2.744 0
26 -5.9 13 13 0 C7H13NO7 4.660 0
14 -5.8 6 8 0 C7H13NO6 2.270 0
15 -5.8 7 8 0 C7H13NO6 3.485 0
13 -5.8 6 6 0 C7H13NO6 4.036 0
16 -5.7 8 8 0 C7H13NO6 3.763 0
17 -5.6 7 7 0 C7H13NO6 4.497 0

Топ-5 находок принадлежат лиганду с фенильным радикалом. Помимо него в десятку лучших попало карбокси-производное NAG. Интересно, что две лучшие находки добиваются хорошей аффинности разными способами: первая - в большей степени за счет гидрофобных взаимодействий, а вторая - водородных связей и стекинга. Посмотрим на них поближе.

In [40]:
for i in (27, 28):
    res[i].write(filename='%s.pdb' % str(i), format='pdb')
Рис. 1a. В этом взаимодействии не удалось найти ни одной из 6 заявленных водородных связей. Гидрофобные остатки (желтые) слишком далеко от лиганда, чтобы образовывать гидрофобные контакты. Рис. 1b. Нашлось 3 из 9 водородных связей, одна из них с тирозином, который и, может быть, образует стэкинг с лигандом.