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

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

In [1]:
import numpy as np
import copy

import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image

import oddt
import oddt.docking
import oddt.interactions

from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

import pmx

Подготовим белок

In [4]:
pdb=pmx.Model('lys2_antmy.B99990001.pdb')

for r in pdb.residues[:123]:
    print r
<Molecule: id = 1 name = LYS chain_id =    natoms = 9>
<Molecule: id = 2 name = ARG chain_id =    natoms = 11>
<Molecule: id = 3 name = PHE chain_id =    natoms = 11>
<Molecule: id = 4 name = THR chain_id =    natoms = 7>
<Molecule: id = 5 name = ARG chain_id =    natoms = 11>
<Molecule: id = 6 name = CYS chain_id =    natoms = 6>
<Molecule: id = 7 name = GLY chain_id =    natoms = 4>
<Molecule: id = 8 name = LEU chain_id =    natoms = 8>
<Molecule: id = 9 name = VAL chain_id =    natoms = 7>
<Molecule: id = 10 name = ASN chain_id =    natoms = 8>
<Molecule: id = 11 name = GLU chain_id =    natoms = 9>
<Molecule: id = 12 name = LEU chain_id =    natoms = 8>
<Molecule: id = 13 name = ARG chain_id =    natoms = 11>
<Molecule: id = 14 name = LYS chain_id =    natoms = 9>
<Molecule: id = 15 name = GLN chain_id =    natoms = 9>
<Molecule: id = 16 name = GLY chain_id =    natoms = 4>
<Molecule: id = 17 name = PHE chain_id =    natoms = 11>
<Molecule: id = 18 name = ASP chain_id =    natoms = 8>
<Molecule: id = 19 name = GLU chain_id =    natoms = 9>
<Molecule: id = 20 name = ASN chain_id =    natoms = 8>
<Molecule: id = 21 name = LEU chain_id =    natoms = 8>
<Molecule: id = 22 name = MET chain_id =    natoms = 8>
<Molecule: id = 23 name = ARG chain_id =    natoms = 11>
<Molecule: id = 24 name = ASP chain_id =    natoms = 8>
<Molecule: id = 25 name = TRP chain_id =    natoms = 14>
<Molecule: id = 26 name = VAL chain_id =    natoms = 7>
<Molecule: id = 27 name = CYS chain_id =    natoms = 6>
<Molecule: id = 28 name = LEU chain_id =    natoms = 8>
<Molecule: id = 29 name = VAL chain_id =    natoms = 7>
<Molecule: id = 30 name = GLU chain_id =    natoms = 9>
<Molecule: id = 31 name = ASN chain_id =    natoms = 8>
<Molecule: id = 32 name = GLU chain_id =    natoms = 9>
<Molecule: id = 33 name = SER chain_id =    natoms = 6>
<Molecule: id = 34 name = ALA chain_id =    natoms = 5>
<Molecule: id = 35 name = ARG chain_id =    natoms = 11>
<Molecule: id = 36 name = TYR chain_id =    natoms = 12>
<Molecule: id = 37 name = THR chain_id =    natoms = 7>
<Molecule: id = 38 name = ASP chain_id =    natoms = 8>
<Molecule: id = 39 name = LYS chain_id =    natoms = 9>
<Molecule: id = 40 name = ILE chain_id =    natoms = 8>
<Molecule: id = 41 name = ALA chain_id =    natoms = 5>
<Molecule: id = 42 name = ASN chain_id =    natoms = 8>
<Molecule: id = 43 name = VAL chain_id =    natoms = 7>
<Molecule: id = 44 name = ASN chain_id =    natoms = 8>
<Molecule: id = 45 name = LYS chain_id =    natoms = 9>
<Molecule: id = 46 name = ASN chain_id =    natoms = 8>
<Molecule: id = 47 name = GLY chain_id =    natoms = 4>
<Molecule: id = 48 name = SER chain_id =    natoms = 6>
<Molecule: id = 49 name = ARG chain_id =    natoms = 11>
<Molecule: id = 50 name = ASP chain_id =    natoms = 8>
<Molecule: id = 51 name = TYR chain_id =    natoms = 12>
<Molecule: id = 52 name = GLY chain_id =    natoms = 4>
<Molecule: id = 53 name = LEU chain_id =    natoms = 8>
<Molecule: id = 54 name = PHE chain_id =    natoms = 11>
<Molecule: id = 55 name = GLN chain_id =    natoms = 9>
<Molecule: id = 56 name = ILE chain_id =    natoms = 8>
<Molecule: id = 57 name = ASN chain_id =    natoms = 8>
<Molecule: id = 58 name = ASP chain_id =    natoms = 8>
<Molecule: id = 59 name = LYS chain_id =    natoms = 9>
<Molecule: id = 60 name = TYR chain_id =    natoms = 12>
<Molecule: id = 61 name = TRP chain_id =    natoms = 14>
<Molecule: id = 62 name = CYS chain_id =    natoms = 6>
<Molecule: id = 63 name = SER chain_id =    natoms = 6>
<Molecule: id = 64 name = LYS chain_id =    natoms = 9>
<Molecule: id = 65 name = GLY chain_id =    natoms = 4>
<Molecule: id = 66 name = SER chain_id =    natoms = 6>
<Molecule: id = 67 name = THR chain_id =    natoms = 7>
<Molecule: id = 68 name = PRO chain_id =    natoms = 7>
<Molecule: id = 69 name = GLY chain_id =    natoms = 4>
<Molecule: id = 70 name = LYS chain_id =    natoms = 9>
<Molecule: id = 71 name = ASP chain_id =    natoms = 8>
<Molecule: id = 72 name = CYS chain_id =    natoms = 6>
<Molecule: id = 73 name = ASN chain_id =    natoms = 8>
<Molecule: id = 74 name = VAL chain_id =    natoms = 7>
<Molecule: id = 75 name = THR chain_id =    natoms = 7>
<Molecule: id = 76 name = CYS chain_id =    natoms = 6>
<Molecule: id = 77 name = SER chain_id =    natoms = 6>
<Molecule: id = 78 name = GLN chain_id =    natoms = 9>
<Molecule: id = 79 name = LEU chain_id =    natoms = 8>
<Molecule: id = 80 name = LEU chain_id =    natoms = 8>
<Molecule: id = 81 name = THR chain_id =    natoms = 7>
<Molecule: id = 82 name = ASP chain_id =    natoms = 8>
<Molecule: id = 83 name = ASP chain_id =    natoms = 8>
<Molecule: id = 84 name = ILE chain_id =    natoms = 8>
<Molecule: id = 85 name = THR chain_id =    natoms = 7>
<Molecule: id = 86 name = VAL chain_id =    natoms = 7>
<Molecule: id = 87 name = ALA chain_id =    natoms = 5>
<Molecule: id = 88 name = SER chain_id =    natoms = 6>
<Molecule: id = 89 name = THR chain_id =    natoms = 7>
<Molecule: id = 90 name = CYS chain_id =    natoms = 6>
<Molecule: id = 91 name = ALA chain_id =    natoms = 5>
<Molecule: id = 92 name = LYS chain_id =    natoms = 9>
<Molecule: id = 93 name = LYS chain_id =    natoms = 9>
<Molecule: id = 94 name = ILE chain_id =    natoms = 8>
<Molecule: id = 95 name = TYR chain_id =    natoms = 12>
<Molecule: id = 96 name = LYS chain_id =    natoms = 9>
<Molecule: id = 97 name = ARG chain_id =    natoms = 11>
<Molecule: id = 98 name = THR chain_id =    natoms = 7>
<Molecule: id = 99 name = LYS chain_id =    natoms = 9>
<Molecule: id = 100 name = PHE chain_id =    natoms = 11>
<Molecule: id = 101 name = ASP chain_id =    natoms = 8>
<Molecule: id = 102 name = ALA chain_id =    natoms = 5>
<Molecule: id = 103 name = TRP chain_id =    natoms = 14>
<Molecule: id = 104 name = SER chain_id =    natoms = 6>
<Molecule: id = 105 name = GLY chain_id =    natoms = 4>
<Molecule: id = 106 name = TRP chain_id =    natoms = 14>
<Molecule: id = 107 name = ASP chain_id =    natoms = 8>
<Molecule: id = 108 name = ASN chain_id =    natoms = 8>
<Molecule: id = 109 name = HIS chain_id =    natoms = 10>
<Molecule: id = 110 name = CYS chain_id =    natoms = 6>
<Molecule: id = 111 name = ASN chain_id =    natoms = 8>
<Molecule: id = 112 name = HIS chain_id =    natoms = 10>
<Molecule: id = 113 name = SER chain_id =    natoms = 6>
<Molecule: id = 114 name = ASN chain_id =    natoms = 8>
<Molecule: id = 115 name = PRO chain_id =    natoms = 7>
<Molecule: id = 116 name = ASP chain_id =    natoms = 8>
<Molecule: id = 117 name = ILE chain_id =    natoms = 8>
<Molecule: id = 118 name = SER chain_id =    natoms = 6>
<Molecule: id = 119 name = SER chain_id =    natoms = 6>
<Molecule: id = 120 name = CYS chain_id =    natoms = 7>
<Molecule: id = 121 name = NAG chain_id =    natoms = 14>
<Molecule: id = 122 name = NAG chain_id =    natoms = 14>
<Molecule: id = 123 name = NDG chain_id =    natoms = 15>
In [5]:
# создаём отдельные объекты для белка и для лиганда
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
lig = pdb.copy()
for r in lig.residues[:-3]:
    lig.remove_residue(r)
In [7]:
# Ищем геометрический центр лиганда
count=0
x=0
y=0
z=0
for a in lig.atoms:
    x=x+a.x[0]
    y=y+a.x[1]
    z=z+a.x[2]    
    print a.x
    count=count+1
x_mean=x/count
y_mean=y/count
z_mean=z/count
print "Координаты геометрического центра лиганда:"
print x_mean, y_mean, z_mean
[38.606, 46.391, 31.413]
[38.518, 47.93, 31.743]
[37.49, 48.219, 32.819]
[36.093, 47.46, 32.619]
[36.395, 45.978, 32.343]
[35.222, 45.233, 31.922]
[39.973, 49.194, 30.877]
[41.195, 49.842, 31.56]
[39.585, 48.367, 31.861]
[36.987, 49.595, 32.722]
[35.195, 47.76, 33.701]
[37.313, 45.802, 31.247]
[34.465, 45.859, 30.944]
[39.616, 49.562, 29.733]
[41.478, 43.421, 27.877]
[39.947, 43.035, 27.982]
[39.183, 44.211, 28.582]
[39.888, 44.894, 29.84]
[41.369, 45.067, 29.63]
[42.104, 45.402, 30.864]
[39.012, 41.175, 27.08]
[38.618, 40.992, 25.615]
[39.425, 42.457, 27.166]
[37.888, 43.884, 29.106]
[39.201, 46.165, 29.99]
[41.92, 43.78, 29.247]
[41.807, 44.621, 31.999]
[38.903, 40.142, 27.764]
[45.577, 41.034, 26.17]
[45.448, 42.554, 26.386]
[44.245, 42.869, 27.278]
[42.923, 42.254, 26.8]
[42.988, 40.842, 26.176]
[42.137, 40.922, 24.935]
[47.293, 43.735, 27.382]
[47.16, 44.837, 26.281]
[44.308, 40.506, 25.647]
[44.112, 44.266, 27.659]
[42.027, 42.143, 27.946]
[42.445, 42.053, 24.065]
[48.389, 43.844, 28.084]
[46.387, 42.711, 27.608]
[45.846, 40.604, 27.235]
Координаты геометрического центра лиганда:
40.8995581395 44.4560930233 29.0208837209
In [8]:
#Сохраняем белок без лиганда в отдельный pdb-файл
newpdb.writePDB("lys_without_lig.pdb")
In [22]:
# Готовим белок для докинга
prot = oddt.toolkit.readfile('pdb','lys_without_lig.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
print 'is it the first mol in 1lmp is protein?',prot.protein
print 'MW of this mol is:', prot.molwt 
is it the first mol in 1lmp is protein? False
MW of this mol is: 13735.20358

Теперь подготовим лиганды для докинга. Создадим 5 лигандов, где метильный радикал группы СH3C(=O)NH из NAG будет заменён на:

  • OH
  • NH3+
  • H
  • Ph
  • COO-
In [23]:
smiles = ['OC(=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',
          '[NH3+]C(=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',
          'C2=CC=C(C=C2)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 H O 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 H H H O O N O O O N O *****
***** - Open Babel Depiction H H H H H O O O O O O N O OH *****
***** - Open Babel Depiction O H H H H H O O O O N O *****

Теперь для каждого из этих лигандов проведем докинг.

In [16]:
# Создаём объект для докинга
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=[38.9964186047, 45.1291395349, 27.5343023256],
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params)
/tmp/autodock_vina_xyg1fq
--center_x 38.9964186047 --center_y 45.1291395349 --center_z 27.5343023256 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

Результат выдачи предыдущего дейсвия это параметры, которые используются для запуска докинга:

  • center_x, center_y, center_z - центр куба, куда будет помещена молекула;
  • size_x, size_y, size_z - размер куба;
  • exhaustiveness - параметр Autodock Vina;
  • num_modes - число конформаций, генерируемых Autodock Vina;
  • energy_range - cut-off по энергии.
In [17]:
# запускаем сам докинг
res = dock_obj.dock(mols,prot)
In [20]:
# проанализируем результаты докинга
from operator import itemgetter
dock_res=[]
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)#гидрофобные контакты
    dock_res.append([i,r.formula, float(r.data['vina_affinity']),  float(r.data['vina_rmsd_ub']), 
                     r.residues[0].name, len(hbs[1]), len(stack[1]), len(phob[1])])
sort_res = sorted(dock_res, key=itemgetter(2))
In [21]:
for resi in sort_res:
    print "%4d%10s%8s%8s%8s%4d%4d%4d" %(resi[0],resi[1],resi[2],resi[3],resi[4],resi[5],resi[6],resi[7])
  36 C13H17NO6    -7.1     0.0     UNL  15   0   9
  37 C13H17NO6    -6.8   2.225     UNL   8   0   6
  38 C13H17NO6    -6.4   4.763     UNL   3   1  15
  39 C13H17NO6    -6.3   2.402     UNL   9   1   6
  40 C13H17NO6    -6.2   6.354     UNL   9   0   3
   0  C7H13NO7    -6.0     0.0     UNL  15   0   0
  18 C7H15N2O6    -6.0     0.0     UNL  14   0   0
  19 C7H15N2O6    -5.9   2.607     UNL  12   0   0
  27  C8H13NO8    -5.9     0.0     UNL  13   0   0
  28  C8H13NO8    -5.9   3.158     UNL  16   0   0
  41 C13H17NO6    -5.9   2.415     UNL  10   0   3
   1  C7H13NO7    -5.8   2.649     UNL  12   0   0
  29  C8H13NO8    -5.8    2.94     UNL  12   0   0
  42 C13H17NO6    -5.8   4.168     UNL   3   1  11
  30  C8H13NO8    -5.7   5.596     UNL  16   0   0
  31  C8H13NO8    -5.7   2.808     UNL  11   0   0
  32  C8H13NO8    -5.7   6.005     UNL  17   0   0
  43 C13H17NO6    -5.7  11.441     UNL   3   1  13
   2  C7H13NO7    -5.6   3.106     UNL  12   0   0
   9  C7H13NO6    -5.6     0.0     UNL  13   0   0
  33  C8H13NO8    -5.6   5.254     UNL  16   0   0
  44 C13H17NO6    -5.6     5.5     UNL   4   0   3
   3  C7H13NO7    -5.5   5.087     UNL  14   0   0
  20 C7H15N2O6    -5.5   3.166     UNL  10   0   0
   4  C7H13NO7    -5.4   5.599     UNL  16   0   0
  10  C7H13NO6    -5.4   2.709     UNL  13   0   0
  21 C7H15N2O6    -5.4   5.544     UNL  12   0   0
  22 C7H15N2O6    -5.4   5.152     UNL  12   0   0
  11  C7H13NO6    -5.3   3.249     UNL  14   0   0
  34  C8H13NO8    -5.3   4.959     UNL  15   0   0
   5  C7H13NO7    -5.2   2.261     UNL  13   0   0
   6  C7H13NO7    -5.2   5.234     UNL  17   0   0
  35  C8H13NO8    -5.2   6.254     UNL  15   0   0
   7  C7H13NO7    -5.1   4.741     UNL  17   0   0
   8  C7H13NO7    -5.1    3.53     UNL  13   0   0
  12  C7H13NO6    -5.1   5.105     UNL  13   0   0
  13  C7H13NO6    -5.1   2.242     UNL  13   0   0
  14  C7H13NO6    -5.0    4.23     UNL  13   0   0
  15  C7H13NO6    -5.0   4.416     UNL  13   0   0
  23 C7H15N2O6    -5.0   4.648     UNL   9   0   0
  16  C7H13NO6    -4.9   3.968     UNL  13   0   0
  17  C7H13NO6    -4.9   4.918     UNL  10   0   0
  24 C7H15N2O6    -4.9   5.431     UNL   4   0   0
  25 C7H15N2O6    -4.7   4.813     UNL  13   0   0
  26 C7H15N2O6    -4.7   5.008     UNL   6   0   0

Наилучшей оказалась структура под номером 36 (C13H17NO6), которая содержит фенил-группу в лиганде, а наихудшей - 26 структура (C7H15N2O6), содержащая группу NH3 в лиганде. Теперь получим трехмерные структуры всех полученных структур.

In [24]:
#запишем результаты в файлы pdb
for i,r in enumerate(res):
    r.write(filename='r%s.pdb' %i, format='pdb')
In [25]:
from IPython.display import Image
Image(filename='lys_docking.png')
Out[25]:

Рис.1. Белок лизоцим с лигандом. Черный цветом покрашен исходный лиганд, розовым - лиганд из наилучшей структуры по результатам докинга (36), зеленым - лиганд из худшей структуры по результатам докинга (26).
Все три лиганда располагаются в одном кармане и достаточно близко к друг другу, но лиганд из лучшей структуры "сидит" наиболее глубоко в этом кармане.