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


Цель данного занятия ознакомится с возможностями докинга низкомолекулярного лиганда в структуру белка. В этом занятии был использован пакет Autodock Vina и oddt. Надо определить место докинга, удалить лигнад, запустить докинг и провести анализ.

Подключение модулей

In [2]:
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 [3]:
pdb=pmx.Model('WithLigand.B99990001.pdb')

for r in pdb.residues[135:]:
    print r #посмотрим остатки
<Molecule: id = 136 name = ILE chain_id =    natoms = 8>
<Molecule: id = 137 name = SER chain_id =    natoms = 6>
<Molecule: id = 138 name = ASP chain_id =    natoms = 8>
<Molecule: id = 139 name = CYS chain_id =    natoms = 7>
<Molecule: id = 140 name = NAG chain_id =    natoms = 14>
<Molecule: id = 141 name = NAG chain_id =    natoms = 14>
<Molecule: id = 142 name = NDG chain_id =    natoms = 15>
In [5]:
# создание объектов белок и лиганда
newpdb = pdb.copy()
lig = pdb.copy()
del newpdb.residues[139:]
del lig.residues[:139]
In [6]:
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])
cen_x=average_in_list(x)
cen_y=average_in_list(y)
cen_z=average_in_list(z)
    # найдите геометрический центр лиганда
In [7]:
newpdb.writePDB("hyace.pdb")

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

In [8]:
prot = oddt.toolkit.readfile('pdb','hyace.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: 16500.6646
In [9]:
smiles = ['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=CC=C1C(=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']
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 O H H H H H O O CH 3 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 H H H O O N O 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 O O O O N O *****
***** - Open Babel Depiction H H H H H O O O O O O N O OH *****

Докинг

In [10]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=(cen_x, cen_y, cen_z),
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
/tmp/autodock_vina_GiMYXW
--center_x 40.1232095238 --center_y 35.7878415584 --center_z 21.0768008658 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

/tmp/autodock_vina_q6iT_Q --center_x 40.1232095238 --center_y 35.7878415584 --center_z 21.0768008658 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

In [12]:
res = dock_obj.dock(mols,prot)
In [13]:
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)
   0  C8H15NO6     6.7   0.000     UNL
   1  C8H15NO6     6.8   1.922     UNL
   2  C8H15NO6     7.2   5.372     UNL
   3  C8H15NO6     8.4   4.407     UNL
   4  C8H15NO6     8.7   3.672     UNL
   5  C8H15NO6     8.7   3.213     UNL
   6  C8H15NO6     8.8   7.742     UNL
   7  C8H15NO6     8.9   3.960     UNL
   8  C8H15NO6     9.4   2.810     UNL
   9  C7H13NO7    -3.2   0.000     UNL
  10  C7H13NO7    -2.9   3.024     UNL
  11 C7H16N2O6    -3.2   0.000     UNL
  12 C7H16N2O6    -3.0   2.989     UNL
  13 C7H16N2O6    -2.9   3.876     UNL
  14 C7H16N2O6    -2.1   3.993     UNL
  15 C7H16N2O6    -1.9   4.190     UNL
  16 C7H16N2O6    -0.3   4.775     UNL
  17  C7H13NO6     4.7   0.000     UNL
  18  C7H13NO6     5.1   1.467     UNL
  19  C7H13NO6     5.7   4.312     UNL
  20  C7H13NO6     6.5   4.333     UNL
  21  C7H13NO6     6.6   5.199     UNL
  22  C7H13NO6     7.2   1.538     UNL
  23  C7H13NO6     7.3   4.446     UNL
  24 C13H17NO6     3.2   0.000     UNL
  25  C8H13NO8     7.8   0.000     UNL
  26  C8H13NO8     7.8   5.621     UNL
  27  C8H13NO8     7.9   2.732     UNL
  28  C8H13NO8     8.5   5.475     UNL

0 C6H6O -3.0 0.000 UNL 1 C6H6O -3.0 1.829 UNL 2 C6H6O -3.0 4.389 UNL 3 C6H6O -2.8 3.005 UNL 4 C6H6O -2.7 3.747 UNL 5 C6H6O -2.7 2.933 UNL 6 C6H6O -2.7 2.755 UNL 7 C6H6O -2.6 3.665 UNL 8 C6H6O -2.6 4.117 UNL 9 C6H6O2 -3.0 0.000 UNL 10 C6H6O2 -3.0 3.656 UNL 11 C6H6O2 -3.0 3.231 UNL 12 C6H6O2 -3.0 4.154 UNL 13 C6H6O2 -2.9 4.191 UNL 14 C6H6O2 -2.9 4.397 UNL 15 C6H6O2 -2.7 2.282 UNL 16 C6H6O2 -2.7 3.675 UNL 17 C6H6O2 -2.1 2.743 UNL 18 C12H10O2 -0.1 0.000 UNL 19 C12H10O2 -0.0 2.188 UNL 20 C12H10O2 2.5 5.750 UNL 21 C12H10O2 2.8 5.403 UNL

In [14]:
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]))
       4       0       6
       6       0       3
       3       0       5
       3       0       9
       4       0       5
       5       0       2
       4       0       6
       5       0       6
       2       0       5
       4       0       0
       2       0       0
       5       0       0
       2       0       0
       4       0       0
       1       0       0
       1       0       0
       0       0       0
       5       0       0
       6       0       0
       4       0       0
       5       0       0
       6       0       0
       4       0       0
       3       0       0
       1       0       4
       6       0       0
       8       0       0
       6       0       0
       6       0       0
/home/preps/golovin/miniconda2/lib/python2.7/site-packages/oddt/interactions.py:103: RuntimeWarning: invalid value encountered in greater
  strict = (((angle1 > (base_angle / a_neighbors_num - tolerance)) |
/home/preps/golovin/miniconda2/lib/python2.7/site-packages/oddt/interactions.py:105: RuntimeWarning: invalid value encountered in greater
  ((angle2 > (base_angle / d_neighbors_num - tolerance)) | np.isnan(angle2))).all(axis=-1)

1 0 8 1 0 8 2 0 6 1 0 9 2 1 9 0 0 7 0 0 7 1 1 6 0 0 5 1 0 6 1 0 6 1 0 6 2 0 5 2 0 3 2 0 4 2 0 2 2 1 8 0 0 4 2 4 38 2 4 37 3 2 37 2 2 36

Отсортировали от лучшего заместителя к худшему, результат представлен в таблице.
Из таблицы видно, что лучшим заместителем оказался C7H16N2O6, то есть тот лиганд, где метильная группа была заменена на NH3+. Также хочется отметить, что в соединении в любым типом лигандов вообще отсутствует стэкинг и очень мало гидрофобных взаимодействий.

In [ ]: