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

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

В занятии был использован белок с прошлого занятия - лизоцим из тутового шелкопряда. Последовательность действий состояла в следующем: оределить место докинга, удалить лигнад, запустить докинг и провести анализ.

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 

Посмотрим на то, какие молекулы являются лигандами - это последние три

pdb=pmx.Model('new.B99990001.pdb')

for r in pdb.residues[135:]:
    print r #посмотрим остатки

<Molecule: id = 136 name = SER chain_id =    natoms = 6>
<Molecule: id = 137 name = CYS chain_id =    natoms = 7>
<Molecule: id = 138 name = NAG chain_id =    natoms = 14>
<Molecule: id = 139 name = NAG chain_id =    natoms = 14>
<Molecule: id = 140 name = NDG chain_id =    natoms = 15>

Создадим белок без лиганда и лиганд для докинга

# создание объектов белок и лиганда
newpdb = pdb.copy()
del newpdb.residues[-3:]
lig = pdb.copy()
del lig.residues[:-3]

Теперь наша задача - найти каким местом лиганда его пихать в белок. Поэтому найдём геометрический центр лиганда - среднеарифметическое по каждой его координате, но сначала найдём где искать эти координаты.

print lig

< Model: nchain = 1 nres = 3 natom = 1139 >

print lig.atoms[0]

ATOM      1  N   MET     1      36.652   2.235  54.647  1.00 40.90

print lig.atoms[0].x

[36.652, 2.235, 54.647]

Собственно ищем геометрический центр

x_list = []
y_list = []
z_list = []
for i in lig.atoms:
    x_list.append(i.x[0])
    y_list.append(i.x[1])
    z_list.append(i.x[2])

ave_x=sum(x_list)/float(len(x_list))
ave_y=sum(y_list)/float(len(y_list))
ave_z=sum(z_list)/float(len(z_list))

Запишем теперь белок без лиганда в новый файл

newpdb.writePDB("no_lig_lysbommo.pdb")

prot = oddt.toolkit.readfile('pdb','no_lig_lysbommo.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: 16289.32596

Создадим лиганды для докинга. Создадим лиганды где метильный радикал СH3C(=O)NH группы NAG будет заменён на: OH, NH3+, H, Ph, COO-

smiles = ['O=C(NC1C(O)C(O)C(OC1O)CO)C', '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)C1=CC=CC=C1','O=C(NC1C(O)C(O)C(OC1O)CO)C(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)

#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=[ave_x,ave_y,ave_z],
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу

/tmp/autodock_vina_GVnBwf
--center_x 45.6327989464 --center_y 31.3904969271 --center_z 24.939953468 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

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

# do it
res = dock_obj.dock(mols,prot)

Посмотрим на лиганды, они не подряд по аффинности.

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    -1.7   0.000     UNL
   1  C8H15NO6    -0.8   4.392     UNL
   2  C8H15NO6    -0.1  14.116     UNL
   3  C8H15NO6     0.0   4.818     UNL
   4  C8H15NO6     0.3  14.285     UNL
   5  C8H15NO6     0.6   3.557     UNL
   6  C8H15NO6     0.9   6.066     UNL
   7  C8H15NO6     1.0  15.263     UNL
   8  C8H15NO6     1.1  17.056     UNL
   9  C7H13NO7    -1.9   0.000     UNL
  10  C7H13NO7    -0.7  18.727     UNL
  11  C7H13NO7    -0.3   2.776     UNL
  12  C7H13NO7    -0.1   7.855     UNL
  13  C7H13NO7     0.3   6.864     UNL
  14  C7H13NO7     0.4  18.011     UNL
  15  C7H13NO7     0.8  17.041     UNL
  16  C7H13NO7     1.0   9.374     UNL
  17 C7H15N2O6    -1.9   0.000     UNL
  18 C7H15N2O6    -1.9  16.507     UNL
  19 C7H15N2O6    -0.5   8.390     UNL
  20 C7H15N2O6    -0.0   6.627     UNL
  21 C7H15N2O6     0.7   5.047     UNL
  22 C7H15N2O6     0.8  13.190     UNL
  23 C7H15N2O6     0.9   7.824     UNL
  24 C7H15N2O6     1.0  15.013     UNL
  25  C7H13NO6    -3.6   0.000     UNL
  26  C7H13NO6    -1.2   2.588     UNL
  27  C7H13NO6    -1.1   4.215     UNL
  28  C7H13NO6    -1.0   4.464     UNL
  29  C7H13NO6    -0.6   9.341     UNL
  30 C13H17NO6    -0.9   0.000     UNL
  31 C13H17NO6    -0.0  13.505     UNL
  32 C13H17NO6     1.5   2.226     UNL
  33 C13H17NO6     1.8   7.148     UNL
  34 C13H17NO6     1.9  13.496     UNL
  35  C8H15NO8     0.5   0.000     UNL
  36  C8H15NO8     1.1  16.668     UNL
  37  C8H15NO8     1.2  14.491     UNL
  38  C8H15NO8     1.5  13.565     UNL
  39  C8H15NO8     1.6  13.222     UNL
  40  C8H15NO8     1.8  15.992     UNL
  41  C8H15NO8     1.8   2.890     UNL
  42  C8H15NO8     2.2  13.835     UNL
  43  C8H15NO8     2.3   5.876     UNL

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

А теперь отсортируем лиганды по аффинности

sort_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)
    sort_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(sort_res, key=lambda resi: resi[2])

for j in sort_res:
    print "%4d%10s%8s%8s%8s%4d%4d%4d" %(j[0],j[1],j[2],j[3],j[4],j[5],j[6],j[7])

  25  C7H13NO6    -3.6     0.0     UNL   9   0   0
   9  C7H13NO7    -1.9     0.0     UNL  13   0   0
  17 C7H15N2O6    -1.9     0.0     UNL   5   0   0
  18 C7H15N2O6    -1.9  16.507     UNL   5   0   0
   0  C8H15NO6    -1.7     0.0     UNL   4   0   1
  26  C7H13NO6    -1.2   2.588     UNL  10   0   0
  27  C7H13NO6    -1.1   4.215     UNL  11   0   0
  28  C7H13NO6    -1.0   4.464     UNL  10   0   0
  30 C13H17NO6    -0.9     0.0     UNL   5   0  17
   1  C8H15NO6    -0.8   4.392     UNL   3   0   0
  10  C7H13NO7    -0.7  18.727     UNL   2   0   0
  29  C7H13NO6    -0.6   9.341     UNL   3   0   0
  19 C7H15N2O6    -0.5    8.39     UNL   7   0   0
  11  C7H13NO7    -0.3   2.776     UNL  10   0   0
   2  C8H15NO6    -0.1  14.116     UNL  12   0   3
  12  C7H13NO7    -0.1   7.855     UNL  12   0   0
   3  C8H15NO6     0.0   4.818     UNL   1   0   0
  20 C7H15N2O6    -0.0   6.627     UNL   9   0   0
  31 C13H17NO6    -0.0  13.505     UNL   9   0   8
   4  C8H15NO6     0.3  14.285     UNL  16   0   4
  13  C7H13NO7     0.3   6.864     UNL  16   0   0
  14  C7H13NO7     0.4  18.011     UNL   3   0   0
  35  C8H15NO8     0.5     0.0     UNL   5   0   0
   5  C8H15NO6     0.6   3.557     UNL   3   0   3
  21 C7H15N2O6     0.7   5.047     UNL   6   0   0
  15  C7H13NO7     0.8  17.041     UNL   6   0   0
  22 C7H15N2O6     0.8   13.19     UNL  10   0   0
   6  C8H15NO6     0.9   6.066     UNL   7   0   4
  23 C7H15N2O6     0.9   7.824     UNL   9   0   0
   7  C8H15NO6     1.0  15.263     UNL   3   0   5
  16  C7H13NO7     1.0   9.374     UNL   4   0   0
  24 C7H15N2O6     1.0  15.013     UNL   2   0   0
   8  C8H15NO6     1.1  17.056     UNL   9   0   4
  36  C8H15NO8     1.1  16.668     UNL  18   0   0
  37  C8H15NO8     1.2  14.491     UNL   4   0   0
  32 C13H17NO6     1.5   2.226     UNL   3   0  14
  38  C8H15NO8     1.5  13.565     UNL  15   0   0
  39  C8H15NO8     1.6  13.222     UNL  10   0   0
  33 C13H17NO6     1.8   7.148     UNL  21   1  24
  40  C8H15NO8     1.8  15.992     UNL   9   0   0
  41  C8H15NO8     1.8    2.89     UNL   4   0   0
  34 C13H17NO6     1.9  13.496     UNL   8   0  10
  42  C8H15NO8     2.2  13.835     UNL   9   0   0
  43  C8H15NO8     2.3   5.876     UNL   8   0   0

Наилучшей оказалась структура 25 C7H13NO6 с энергией -3.6. Она содержит 9 водородных связей, и не содержит стэкингов и гидрофобных контактов.