Учебный сайт Ксении Березиной

Назад к семестру

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

Будем работать с белком лизоцимом, структуру которого построили на основе гомологичного моделирования на прошлом практикуме: LYS_MERLU. Используем пакеты 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

# 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-lig.B99990001.pdb')

for r in pdb.residues[100:]:
    print r #посмотрим концевые остатки
<Molecule: id = 101 name = LYS chain_id =    natoms = 9>
<Molecule: id = 102 name = LYS chain_id =    natoms = 9>
<Molecule: id = 103 name = GLY chain_id =    natoms = 4>
<Molecule: id = 104 name = SER chain_id =    natoms = 6>
<Molecule: id = 105 name = THR chain_id =    natoms = 7>
<Molecule: id = 106 name = ILE chain_id =    natoms = 8>
<Molecule: id = 107 name = GLY chain_id =    natoms = 4>
<Molecule: id = 108 name = TYR chain_id =    natoms = 12>
<Molecule: id = 109 name = TRP chain_id =    natoms = 14>
<Molecule: id = 110 name = ASN chain_id =    natoms = 8>
<Molecule: id = 111 name = ARG chain_id =    natoms = 11>
<Molecule: id = 112 name = LEU chain_id =    natoms = 8>
<Molecule: id = 113 name = GLN chain_id =    natoms = 9>
<Molecule: id = 114 name = LYS chain_id =    natoms = 9>
<Molecule: id = 115 name = ILE chain_id =    natoms = 8>
<Molecule: id = 116 name = SER chain_id =    natoms = 6>
<Molecule: id = 117 name = GLY chain_id =    natoms = 4>
<Molecule: id = 118 name = CYS chain_id =    natoms = 6>
<Molecule: id = 119 name = HIS chain_id =    natoms = 10>
<Molecule: id = 120 name = GLY chain_id =    natoms = 4>
<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()
lig = pdb.copy()
del newpdb.residues[122:]
del lig.residues[:122]
In [11]:
# находим координаты центра лиганда
coords= []
x = []
y = []
z = []
for a in lig.atoms:
#    print a.x
    coords.append(a.x) #добываем координаты

for c in coords:
    x.append(c[0])
    y.append(c[1])
    z.append(c[2])
    
x_center = np.mean(x)
y_center = np.mean(y)
z_center = np.mean(z)
print x_center,y_center,z_center
34.8182094943 37.8708286894 18.6351919505
In [10]:
newpdb.writePDB("onlyprot.pdb")

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

In [13]:
prot = oddt.toolkit.readfile('pdb','onlyprot.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: 14002.98796

Нам говорят, что первая молекула в файле onlyprot.pdb - не белок. Но на самом деле в файле белок и все хорошо.

Лиганды для докинга.
NAG (N-ацетил-D-глюкозамин) содержит в себе СH3C(=O)NH группу. Будем проводить докинг с ним самим и с лигандами на его основе, где метильный радикал этой группы заменён на OH, NH3+, H, Ph, COO-.

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

Запускаем докинг.

In [15]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=[x_center,y_center,z_center],
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params)
/tmp/autodock_vina_opWpBA
--center_x 34.8182094943 --center_y 37.8708286894 --center_z 18.6351919505 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

Выдача:
координаты центра лиганда
размеры ячейки (с этим центром), в которой будет идти докинг
используем 1 ядро для работы
точность (полнота расчетов)
маскимальное число положений лиганда
максимальная разница энергии лучшего положения с худшим.

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

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

In [17]:
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.6   0.000     UNL
   1  C8H15NO6    -6.4   4.477     UNL
   2  C8H15NO6    -6.3   4.484     UNL
   3  C8H15NO6    -6.3   4.770     UNL
   4  C8H15NO6    -6.2   3.181     UNL
   5  C8H15NO6    -6.1   2.433     UNL
   6  C8H15NO6    -6.1   3.186     UNL
   7  C8H15NO6    -6.1   5.160     UNL
   8  C8H15NO6    -6.1   3.085     UNL
   9  C7H13NO7    -6.9   0.000     UNL
  10  C7H13NO7    -6.8   4.818     UNL
  11  C7H13NO7    -6.7   2.989     UNL
  12  C7H13NO7    -6.6   5.052     UNL
  13  C7H13NO7    -6.5   4.520     UNL
  14  C7H13NO7    -6.4   3.743     UNL
  15  C7H13NO7    -6.4   4.224     UNL
  16  C7H13NO7    -6.3   3.108     UNL
  17  C7H13NO7    -6.3   3.170     UNL
  18 C7H16N2O6    -6.5   0.000     UNL
  19 C7H16N2O6    -6.4   2.910     UNL
  20 C7H16N2O6    -6.4   4.330     UNL
  21 C7H16N2O6    -6.3   4.958     UNL
  22 C7H16N2O6    -6.2   3.587     UNL
  23 C7H16N2O6    -6.2   4.760     UNL
  24 C7H16N2O6    -6.2   3.172     UNL
  25 C7H16N2O6    -6.1   4.799     UNL
  26 C7H16N2O6    -6.1   3.565     UNL
  27  C7H13NO6    -6.5   0.000     UNL
  28  C7H13NO6    -6.5   3.895     UNL
  29  C7H13NO6    -6.4   3.946     UNL
  30  C7H13NO6    -6.2   4.620     UNL
  31  C7H13NO6    -6.2   3.909     UNL
  32  C7H13NO6    -6.2   4.243     UNL
  33  C7H13NO6    -6.2   3.729     UNL
  34  C7H13NO6    -6.1   4.602     UNL
  35  C7H13NO6    -6.0   2.683     UNL
  36 C13H17NO6    -5.5   0.000     UNL
  37 C13H17NO6    -5.4   3.027     UNL
  38 C13H17NO6    -5.2   6.520     UNL
  39 C13H17NO6    -4.5   2.125     UNL
  40 C13H17NO6    -4.4   6.595     UNL
  41 C13H17NO6    -4.2   4.398     UNL
  42 C13H17NO6    -3.3   5.864     UNL
  43 C13H17NO6    -2.9   3.702     UNL
  44 C13H17NO6    -2.8   3.892     UNL
  45  C8H13NO8    -6.7   0.000     UNL
  46  C8H13NO8    -6.6   5.036     UNL
  47  C8H13NO8    -6.6   2.012     UNL
  48  C8H13NO8    -6.5   5.543     UNL
  49  C8H13NO8    -6.3   5.835     UNL
  50  C8H13NO8    -6.2   2.933     UNL
  51  C8H13NO8    -5.9   3.439     UNL
  52  C8H13NO8    -5.9   4.678     UNL
  53  C8H13NO8    -5.8   4.229     UNL

Анализ докинга:

In [28]:
print "i ligand     hbonds    stacking    hf contacts"
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 "%s %8s%8s%8s%8s" % (i,r.formula,len(hbs[1]), len(stack[1]), len(phob[1]))
i ligand     hbonds    stacking    hf contacts
0 C8H15NO6       5       0       2
1 C8H15NO6       6       0       0
2 C8H15NO6       4       0       3
3 C8H15NO6       6       0       3
4 C8H15NO6       5       0       0
5 C8H15NO6       5       0       3
6 C8H15NO6       6       0       0
7 C8H15NO6       4       0       1
8 C8H15NO6       8       0       3
9 C7H13NO7      11       0       0
10 C7H13NO7      11       0       0
11 C7H13NO7       6       0       0
12 C7H13NO7      11       0       0
13 C7H13NO7      10       0       0
14 C7H13NO7       7       0       0
15 C7H13NO7       6       0       0
16 C7H13NO7       8       0       0
17 C7H13NO7       4       0       0
18 C7H16N2O6       6       0       0
19 C7H16N2O6      10       0       0
20 C7H16N2O6       5       0       0
21 C7H16N2O6       8       0       0
22 C7H16N2O6       7       0       0
23 C7H16N2O6       8       0       0
24 C7H16N2O6       9       0       0
25 C7H16N2O6      10       0       0
26 C7H16N2O6       4       0       0
27 C7H13NO6      10       0       0
28 C7H13NO6       5       0       0
29 C7H13NO6       6       0       0
30 C7H13NO6       8       0       0
31 C7H13NO6       7       0       0
32 C7H13NO6       7       0       0
33 C7H13NO6       5       0       0
34 C7H13NO6       7       0       0
35 C7H13NO6       5       0       0
36 C13H17NO6      10       0      10
37 C13H17NO6       7       0      11
38 C13H17NO6       9       0       7
39 C13H17NO6      10       0       8
40 C13H17NO6       6       0       9
41 C13H17NO6       5       0      14
42 C13H17NO6      11       0      16
43 C13H17NO6      11       0      12
44 C13H17NO6       8       0       9
45 C8H13NO8       7       0       0
46 C8H13NO8       9       0       0
47 C8H13NO8      11       0       0
48 C8H13NO8      11       0       0
49 C8H13NO8       6       0       0
50 C8H13NO8       8       0       0
51 C8H13NO8       7       0       0
52 C8H13NO8       9       0       0
53 C8H13NO8       9       0       0

Посмотрим глазами:

In [29]:
for i,r in enumerate(res):
    r.write(filename='r%s.pdb' % i, format='pdb')

Получили 53 pdb-файла с всеми положениями всех лигандов.

Место для всех лигандов, найденное докингом, отличается от первоначального (зеленая молекула это все три лиганда):

In [30]:
from IPython.display import Image
Image(filename='all.png')
Out[30]:

В таблице показаны заместители от лучшего к худшему. Сортировали по vina_affinity, в случае одинаковых значений сравнивали по водородным связям.

groupligandvina_affinityh_bonds
9OHC7H13NO7-6.911
45COO-C8H13NO8-6.77
0noneC8H15NO6-6.65
18NH3+C7H16N2O6-6.56
27HC7H13NO6-6.510
36PhC13H17NO6-5.510

Назад к семестру