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

В этом практикуме мы будем проводить докинг разных лигандов в структуру лизоцима жабы, которая была определа в прошлом практикуме. Делать это будем с помощью пакетов Autodock Vina и oddt. Для работы возьмем третью модель BUFGA.pdb (у нее была скор 1074).

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

Отделим лиганд от аминокислотных остатков белка: найдем в пдб-файле остатки лиганда и разделим информацию на 2 пдб-файла: только про белок и только про лиганд.

In [3]:
for r in pdb.residues[135:]:
    print r #посмотрим остатки
<Molecule: id = 136 name = ASN chain_id =    natoms = 8>
<Molecule: id = 137 name = LEU chain_id =    natoms = 8>
<Molecule: id = 138 name = SER chain_id =    natoms = 6>
<Molecule: id = 139 name = GLN chain_id =    natoms = 9>
<Molecule: id = 140 name = TRP chain_id =    natoms = 14>
<Molecule: id = 141 name = THR chain_id =    natoms = 7>
<Molecule: id = 142 name = GLN chain_id =    natoms = 9>
<Molecule: id = 143 name = GLY chain_id =    natoms = 4>
<Molecule: id = 144 name = CYS chain_id =    natoms = 6>
<Molecule: id = 145 name = LYS chain_id =    natoms = 9>
<Molecule: id = 146 name = LEU chain_id =    natoms = 9>
<Molecule: id = 147 name = NAG chain_id =    natoms = 14>
<Molecule: id = 148 name = NAG chain_id =    natoms = 14>
<Molecule: id = 149 name = NDG chain_id =    natoms = 15>
In [4]:
prot = pdb.copy()
for r in prot.residues[-3:]:
    prot.remove_residue(r)
lig = pdb.copy()
for r in lig.residues[:-3]:
    lig.remove_residue(r)
In [5]:
prot.writePDB('BUFGA_prot.pdb')
lig.writePDB('lig.pdb')

Находим геометрический центр лиганда:

In [6]:
coords = np.zeros((len(lig.atoms),3))
for c,a in enumerate(lig.atoms):
    coords[c,:] = a.x
geom_center = np.mean(coords, axis = 0)
geom_center
Out[6]:
array([ 40.72974419,  44.23083721,  27.9565814 ])

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

In [7]:
prot = oddt.toolkit.readfile('pdb','BUFGA_prot.pdb').next()

prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()


print 'is it the first mol in BUFGA_prot is protein?',prot.protein,':) and MW of this mol is:', prot.molwt 
is it the first mol in BUFGA_prot is protein? False :) and MW of this mol is: 16599.6016

Странно, что oddt считает, что в BUFGA_prot.pdb не белок, потому что ничего другого там нет.

In [12]:
! less BUFGA_prot.pdb
TITLE    PMX MODEL
MODEL    1
ATOM      1  N   SER     1      33.781   0.019  30.825  1.00 54.36
ATOM      2  CA  SER     1      34.268   1.155  31.636  1.00 54.36
ATOM      3  CB  SER     1      33.075   1.889  32.274  1.00 54.36
ATOM      4  OG  SER     1      33.499   3.095  32.892  1.00 54.36
ATOM      5  C   SER     1      35.168   0.633  32.707  1.00 54.36
ATOM      6  O   SER     1      35.926  -0.311  32.485  1.00 54.36
ATOM      7  N   GLY     2      35.094   1.228  33.911  1.00 73.67
ATOM      8  CA  GLY     2      35.946   0.798  34.977  1.00 73.67
ATOM      9  C   GLY     2      37.179   1.637  34.936  1.00 73.67
ATOM     10  O   GLY     2      37.222   2.667  34.264  1.00 73.67
ATOM     11  N   LYS     3      38.219   1.207  35.676  1.00  6.50
ATOM     12  CA  LYS     3      39.428   1.972  35.729  1.00  6.50
ATOM     13  CB  LYS     3      40.464   1.410  36.725  1.00  6.50
ATOM     14  CG  LYS     3      41.671   2.320  36.997  1.00  6.50
ATOM     15  CD  LYS     3      42.527   2.631  35.767  1.00  6.50
ATOM     16  CE  LYS     3      43.770   3.469  36.075  1.00  6.50
ATOM     17  NZ  LYS     3      44.922   2.583  36.359  1.00  6.50
ATOM     18  C   LYS     3      40.030   1.953  34.365  1.00  6.50
ATOM     19  O   LYS     3      40.153   0.905  33.734  1.00  6.50
ATOM     20  N   TYR     4      40.415   3.148  33.883  1.00 99.31
ATOM     21  CA  TYR     4      41.021   3.299  32.595  1.00 99.31
BUFGA_prot.pdb lines 1-23/1169 2%

Затем создаем набор лигандов для докинга на основе известного NAG.

In [8]:
smiles = ['O=C(NC1C(O)C(O)C(OC1O)CO)O', 'O=C(NC1C(O)C(O)C(OC1O)CO)[NH3+]',
          'O=C(NC1C(O)C(O)C(OC1O)CO)', 'O=C(NC1C(O)C(O)C(OC1O)CO)(c1ccccc1)',
          'NC1C(C(C(OC1O)CO)O)O', 'NC1C(C(C(OC1O)CO)O)OC', 
          'O=C(NC1C(O)C(O)C(OC1O)CO)C(=O)O', 'O=C(NC1C(C)C(C)C(OC1O)CC)C', 'O=C(NC1C(O)C(O)C(OC1O)CO)C']
In [9]:
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 O O O O O O O N *****
***** - Open Babel Depiction O H H H H H H H H N O O O O O N *****
***** - Open Babel Depiction O H H H H H O O O O O N *****
***** - Open Babel Depiction H H H H H O O O O O O N *****
***** - Open Babel Depiction O H H H H H H O O N O O *****
***** - Open Babel Depiction O H H H H H CH 3 O O N O O *****
***** - Open Babel Depiction O H H H H H H O O O O O O O N *****
***** - Open Babel Depiction H H H 3 C CH 3 O O O CH 3 H 3 C N *****
***** - Open Babel Depiction H H H H H CH 3 O O O O O O N *****

Сделала NAG (последняя структура), несколько его модификаций с заменой метильного радикала СH3C(=O)NH группы, а также лиганд с тремя метильными радикалами вместо ОН в сахаре (предпоследняя структура).

Дальше проводим докинг всех этих лигандов в жабий лизоцим. Центром докинга сделаем найденный геометрический центр лиганда.

In [10]:
# create docking object
dock_object = 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_object.tmp_dir
print " ".join(dock_object.params)
/tmp/autodock_vina_DF7SA2
--center_x 40.729744186 --center_y 44.2308372093 --center_z 27.9565813953 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3
In [11]:
# do it
res = dock_object.dock(mols, prot)

Визуализация результатов докинга

In [12]:
import pandas as pd
In [15]:
# отсортируем результаты по аффинности (параметр vina_affinity):
res_sorted = sorted(res, key = lambda x : float(x.data['vina_affinity']))
In [35]:
for i,r in enumerate(res):
    print "%4d%10s%8s%8s" %(i, r.formula, r.data['vina_affinity'],  r.data['vina_rmsd_ub'])
   0  C7H13NO7    -5.9   0.000
   1  C7H13NO7    -5.8   4.879
   2  C7H13NO7    -5.6   3.230
   3  C7H13NO7    -5.6   5.030
   4  C7H13NO7    -5.4   6.497
   5  C7H13NO7    -5.4   1.914
   6  C7H13NO7    -5.3   4.565
   7  C7H13NO7    -5.1   4.639
   8  C7H13NO7    -5.1   3.829
   9 C7H16N2O6    -6.0   0.000
  10 C7H16N2O6    -5.9   4.900
  11 C7H16N2O6    -5.6   3.192
  12 C7H16N2O6    -5.6   4.533
  13 C7H16N2O6    -5.3   5.240
  14 C7H16N2O6    -5.2   5.058
  15 C7H16N2O6    -5.1   4.633
  16 C7H16N2O6    -5.1   2.870
  17 C7H16N2O6    -5.0   5.236
  18  C7H13NO6    -5.7   0.000
  19  C7H13NO6    -5.5   5.087
  20  C7H13NO6    -5.4   4.626
  21  C7H13NO6    -5.4   3.688
  22  C7H13NO6    -5.2   4.933
  23  C7H13NO6    -5.1   4.929
  24  C7H13NO6    -4.9   2.137
  25  C7H13NO6    -4.8   7.907
  26  C7H13NO6    -4.7   6.141
  27 C13H17NO6    -7.0   0.000
  28 C13H17NO6    -6.8  11.724
  29 C13H17NO6    -6.7   8.006
  30 C13H17NO6    -6.5   7.654
  31 C13H17NO6    -6.4   7.275
  32 C13H17NO6    -6.3  11.134
  33 C13H17NO6    -6.3   9.970
  34 C13H17NO6    -6.3  13.318
  35 C13H17NO6    -6.2   1.788
  36  C6H13NO5    -5.0   0.000
  37  C6H13NO5    -5.0   2.821
  38  C6H13NO5    -5.0   3.363
  39  C6H13NO5    -4.9   4.327
  40  C6H13NO5    -4.9   3.002
  41  C6H13NO5    -4.8   3.803
  42  C6H13NO5    -4.8   3.688
  43  C6H13NO5    -4.8   3.722
  44  C6H13NO5    -4.6   2.683
  45  C7H15NO5    -5.1   0.000
  46  C7H15NO5    -5.1   4.647
  47  C7H15NO5    -5.0   4.448
  48  C7H15NO5    -4.6   3.810
  49  C7H15NO5    -4.5   4.667
  50  C7H15NO5    -4.5   3.623
  51  C7H15NO5    -4.3   3.287
  52  C7H15NO5    -4.3   3.036
  53  C7H15NO5    -4.0   4.905
  54  C8H13NO8    -6.5   0.000
  55  C8H13NO8    -6.1   6.355
  56  C8H13NO8    -5.9   6.489
  57  C8H13NO8    -5.4   6.878
  58  C8H13NO8    -5.4   4.288
  59  C8H13NO8    -5.3   6.334
  60  C8H13NO8    -5.3   6.518
  61  C8H13NO8    -5.1   7.496
  62  C8H13NO8    -5.1  10.200
  63 C11H21NO3    -6.6   0.000
  64 C11H21NO3    -5.8   3.471
  65 C11H21NO3    -5.5   5.381
  66 C11H21NO3    -5.4   6.783
  67 C11H21NO3    -5.3   6.748
  68 C11H21NO3    -5.2   3.025
  69 C11H21NO3    -5.2   6.328
  70 C11H21NO3    -5.1   6.157
  71 C11H21NO3    -5.0   7.774
  72  C8H15NO6    -6.4   0.000
  73  C8H15NO6    -5.8   3.645
  74  C8H15NO6    -5.3   6.557
  75  C8H15NO6    -5.3   4.866
  76  C8H15NO6    -5.2   3.144
  77  C8H15NO6    -5.0   5.252
  78  C8H15NO6    -5.0   6.663
  79  C8H15NO6    -5.0   2.929
  80  C8H15NO6    -5.0   3.920
In [45]:
hbtotal = []
hbstrict = []
phob = []
for i,r in enumerate(res_sorted):
    int1, int2, strict = oddt.interactions.hbonds(prot,r)
    hbtotal.append(len(int1))
    hbstrict.append(strict.sum())
    ph1, ph2 = oddt.interactions.hydrophobic_contacts(prot,r)
    phob.append(len(ph1))
sort_table = pd.DataFrame({'Formula':[r.formula for r in res_sorted],
                   'Affinity, kcal/mol':[r.data['vina_affinity'] for r in res_sorted],
                   'Total number of h-bonds':hbtotal, 'Strict number of h-bonds':hbstrict, 'Hydrophobic':phob, 
                    'RMSD':[r.data['vina_rmsd_ub'] for r in res_sorted]})
In [46]:
sort_table[:]
Out[46]:
Affinity, kcal/mol Formula Hydrophobic RMSD Strict number of h-bonds Total number of h-bonds
0 -7.0 C13H17NO6 2 0.000 8 8
1 -6.8 C13H17NO6 3 11.724 7 9
2 -6.7 C13H17NO6 8 8.006 9 11
3 -6.6 C11H21NO3 5 0.000 6 6
4 -6.5 C13H17NO6 7 7.654 6 7
5 -6.5 C8H13NO8 0 0.000 10 12
6 -6.4 C13H17NO6 7 7.275 7 9
7 -6.4 C8H15NO6 3 0.000 10 10
8 -6.3 C13H17NO6 10 11.134 4 5
9 -6.3 C13H17NO6 1 9.970 4 4
10 -6.3 C13H17NO6 12 13.318 8 8
11 -6.2 C13H17NO6 5 1.788 5 7
12 -6.1 C8H13NO8 0 6.355 11 12
13 -6.0 C7H16N2O6 0 0.000 7 8
14 -5.9 C7H13NO7 0 0.000 11 13
15 -5.9 C7H16N2O6 0 4.900 8 12
16 -5.9 C8H13NO8 0 6.489 10 11
17 -5.8 C7H13NO7 0 4.879 9 9
18 -5.8 C11H21NO3 3 3.471 3 3
19 -5.8 C8H15NO6 1 3.645 8 9
20 -5.7 C7H13NO6 0 0.000 8 8
21 -5.6 C7H13NO7 0 3.230 10 11
22 -5.6 C7H13NO7 0 5.030 11 12
23 -5.6 C7H16N2O6 0 3.192 11 11
24 -5.6 C7H16N2O6 0 4.533 7 10
25 -5.5 C7H13NO6 0 5.087 6 7
26 -5.5 C11H21NO3 6 5.381 2 2
27 -5.4 C7H13NO7 0 6.497 10 10
28 -5.4 C7H13NO7 0 1.914 7 7
29 -5.4 C7H13NO6 0 4.626 7 9
... ... ... ... ... ... ...
51 -5.1 C7H15NO5 0 0.000 5 5
52 -5.1 C7H15NO5 0 4.647 10 10
53 -5.1 C8H13NO8 0 7.496 7 7
54 -5.1 C8H13NO8 0 10.200 5 7
55 -5.1 C11H21NO3 3 6.157 4 4
56 -5.0 C7H16N2O6 0 5.236 2 4
57 -5.0 C6H13NO5 0 0.000 6 9
58 -5.0 C6H13NO5 0 2.821 7 7
59 -5.0 C6H13NO5 0 3.363 10 10
60 -5.0 C7H15NO5 0 4.448 10 10
61 -5.0 C11H21NO3 8 7.774 1 1
62 -5.0 C8H15NO6 0 5.252 2 3
63 -5.0 C8H15NO6 3 6.663 4 4
64 -5.0 C8H15NO6 1 2.929 3 6
65 -5.0 C8H15NO6 3 3.920 6 7
66 -4.9 C7H13NO6 0 2.137 6 7
67 -4.9 C6H13NO5 0 4.327 7 8
68 -4.9 C6H13NO5 0 3.002 7 8
69 -4.8 C7H13NO6 0 7.907 4 4
70 -4.8 C6H13NO5 0 3.803 9 9
71 -4.8 C6H13NO5 0 3.688 5 5
72 -4.8 C6H13NO5 0 3.722 7 7
73 -4.7 C7H13NO6 0 6.141 5 6
74 -4.6 C6H13NO5 0 2.683 5 7
75 -4.6 C7H15NO5 0 3.810 8 8
76 -4.5 C7H15NO5 0 4.667 8 10
77 -4.5 C7H15NO5 0 3.623 2 4
78 -4.3 C7H15NO5 0 3.287 3 3
79 -4.3 C7H15NO5 0 3.036 4 6
80 -4.0 C7H15NO5 0 4.905 4 4

81 rows × 6 columns

Самым высокоаффинным лигандом оказался C13H17NO6 (r0.pdb) - это лиганд с фенильной группой вместо метильного радикала. Посмотрим в PyMol его, а также NAG (C8H15NO6, r7.pdb) и еще парочку из представленных выше (C11H21NO3 r55.pdb и C7H15NO5 r80.pdb, который является наименее аффинным, если судить по таблице).

In [37]:
for i in (0, 7, 55, 80):
    res_sorted[i].write(filename='res%s.pdb' % str(i), format='pdb')

C13H17NO6 (r0.pdb)
NAG, C8H15NO6 (r7.pdb)
C7H15NO (r55.pdb)
C7H15NO5 (r80.pdb)

Рис. 1. Придуманные лиганды в структуре белка: C13H17NO6 (r0), C8H15NO6 (NAG, r7), C11H21NO3 (r55) и C7H15NO5 (r80). Желтым показаны водородные связи между лигандом и аминокислотными остатками белка, имеющие длину менее или равную 3.5 ангстрем.

Глюкоза у всех четырех лигандов находится в конформации "кресло". Первый лиганд хорошо ложится в "карман" белка и, вероятно, стабилизируется стэкинговым взаимодействием бензольного кольца с аминокислотными остатками белка, поэтому у него высокая аффинность. Видно, что нативный лиганд NAG взаимодействует с белком большим количеством водородных связей, правда, в одном месте кислород ОН-группы образует аж 3 водородные связи - что-то многовато даже для кислорода. У последнего лиганда азот NH2-группы тоже почему-то образует 2 водородные связи с белком.