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

Цель данного занятия - ознакомиться с возможностями докинга низкомолекулярного лиганда в структуру белка.

Для этого использовали пакет Autodock Vina и oddt.

Работа проводилась с белком лизоцимом восточного моллюска Meretrix lusoria.

Цель – определить место докинга, удалить лиганд, запустить докинг и провести анализ.

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-файл
pdb=pmx.Model('new_seq.B99990001.pdb')
#смотрим на остатки
for r in pdb.residues[:125]:
    print r
<Molecule: id = 1 name = PHE chain_id =    natoms = 11>
<Molecule: id = 2 name = ALA chain_id =    natoms = 5>
<Molecule: id = 3 name = GLY chain_id =    natoms = 4>
<Molecule: id = 4 name = GLY chain_id =    natoms = 4>
<Molecule: id = 5 name = ILE chain_id =    natoms = 8>
<Molecule: id = 6 name = VAL chain_id =    natoms = 7>
<Molecule: id = 7 name = SER chain_id =    natoms = 6>
<Molecule: id = 8 name = GLN chain_id =    natoms = 9>
<Molecule: id = 9 name = ARG chain_id =    natoms = 11>
<Molecule: id = 10 name = CYS chain_id =    natoms = 6>
<Molecule: id = 11 name = LEU chain_id =    natoms = 8>
<Molecule: id = 12 name = SER chain_id =    natoms = 6>
<Molecule: id = 13 name = CYS chain_id =    natoms = 6>
<Molecule: id = 14 name = ILE chain_id =    natoms = 8>
<Molecule: id = 15 name = CYS chain_id =    natoms = 6>
<Molecule: id = 16 name = LYS chain_id =    natoms = 9>
<Molecule: id = 17 name = MET chain_id =    natoms = 8>
<Molecule: id = 18 name = GLU chain_id =    natoms = 9>
<Molecule: id = 19 name = SER chain_id =    natoms = 6>
<Molecule: id = 20 name = GLY chain_id =    natoms = 4>
<Molecule: id = 21 name = CYS chain_id =    natoms = 6>
<Molecule: id = 22 name = ARG chain_id =    natoms = 11>
<Molecule: id = 23 name = ASN chain_id =    natoms = 8>
<Molecule: id = 24 name = VAL chain_id =    natoms = 7>
<Molecule: id = 25 name = GLY chain_id =    natoms = 4>
<Molecule: id = 26 name = CYS chain_id =    natoms = 6>
<Molecule: id = 27 name = LYS chain_id =    natoms = 9>
<Molecule: id = 28 name = MET chain_id =    natoms = 8>
<Molecule: id = 29 name = ASP chain_id =    natoms = 8>
<Molecule: id = 30 name = MET chain_id =    natoms = 8>
<Molecule: id = 31 name = GLY chain_id =    natoms = 4>
<Molecule: id = 32 name = SER chain_id =    natoms = 6>
<Molecule: id = 33 name = LEU chain_id =    natoms = 8>
<Molecule: id = 34 name = SER chain_id =    natoms = 6>
<Molecule: id = 35 name = CYS chain_id =    natoms = 6>
<Molecule: id = 36 name = GLY chain_id =    natoms = 4>
<Molecule: id = 37 name = TYR chain_id =    natoms = 12>
<Molecule: id = 38 name = PHE chain_id =    natoms = 11>
<Molecule: id = 39 name = GLN chain_id =    natoms = 9>
<Molecule: id = 40 name = ILE chain_id =    natoms = 8>
<Molecule: id = 41 name = LYS chain_id =    natoms = 9>
<Molecule: id = 42 name = GLU chain_id =    natoms = 9>
<Molecule: id = 43 name = ALA chain_id =    natoms = 5>
<Molecule: id = 44 name = TYR chain_id =    natoms = 12>
<Molecule: id = 45 name = TRP chain_id =    natoms = 14>
<Molecule: id = 46 name = ILE chain_id =    natoms = 8>
<Molecule: id = 47 name = ASP chain_id =    natoms = 8>
<Molecule: id = 48 name = CYS chain_id =    natoms = 6>
<Molecule: id = 49 name = GLY chain_id =    natoms = 4>
<Molecule: id = 50 name = ARG chain_id =    natoms = 11>
<Molecule: id = 51 name = PRO chain_id =    natoms = 7>
<Molecule: id = 52 name = GLY chain_id =    natoms = 4>
<Molecule: id = 53 name = SER chain_id =    natoms = 6>
<Molecule: id = 54 name = SER chain_id =    natoms = 6>
<Molecule: id = 55 name = TRP chain_id =    natoms = 14>
<Molecule: id = 56 name = LYS chain_id =    natoms = 9>
<Molecule: id = 57 name = SER chain_id =    natoms = 6>
<Molecule: id = 58 name = CYS chain_id =    natoms = 6>
<Molecule: id = 59 name = ALA chain_id =    natoms = 5>
<Molecule: id = 60 name = ALA chain_id =    natoms = 5>
<Molecule: id = 61 name = SER chain_id =    natoms = 6>
<Molecule: id = 62 name = SER chain_id =    natoms = 6>
<Molecule: id = 63 name = TYR chain_id =    natoms = 12>
<Molecule: id = 64 name = CYS chain_id =    natoms = 6>
<Molecule: id = 65 name = ALA chain_id =    natoms = 5>
<Molecule: id = 66 name = SER chain_id =    natoms = 6>
<Molecule: id = 67 name = LEU chain_id =    natoms = 8>
<Molecule: id = 68 name = CYS chain_id =    natoms = 6>
<Molecule: id = 69 name = VAL chain_id =    natoms = 7>
<Molecule: id = 70 name = GLN chain_id =    natoms = 9>
<Molecule: id = 71 name = ASN chain_id =    natoms = 8>
<Molecule: id = 72 name = TYR chain_id =    natoms = 12>
<Molecule: id = 73 name = MET chain_id =    natoms = 8>
<Molecule: id = 74 name = LYS chain_id =    natoms = 9>
<Molecule: id = 75 name = ARG chain_id =    natoms = 11>
<Molecule: id = 76 name = TYR chain_id =    natoms = 12>
<Molecule: id = 77 name = ALA chain_id =    natoms = 5>
<Molecule: id = 78 name = LYS chain_id =    natoms = 9>
<Molecule: id = 79 name = TRP chain_id =    natoms = 14>
<Molecule: id = 80 name = ALA chain_id =    natoms = 5>
<Molecule: id = 81 name = GLY chain_id =    natoms = 4>
<Molecule: id = 82 name = CYS chain_id =    natoms = 6>
<Molecule: id = 83 name = PRO chain_id =    natoms = 7>
<Molecule: id = 84 name = LEU chain_id =    natoms = 8>
<Molecule: id = 85 name = ARG chain_id =    natoms = 11>
<Molecule: id = 86 name = CYS chain_id =    natoms = 6>
<Molecule: id = 87 name = GLU chain_id =    natoms = 9>
<Molecule: id = 88 name = GLY chain_id =    natoms = 4>
<Molecule: id = 89 name = PHE chain_id =    natoms = 11>
<Molecule: id = 90 name = ALA chain_id =    natoms = 5>
<Molecule: id = 91 name = ARG chain_id =    natoms = 11>
<Molecule: id = 92 name = GLU chain_id =    natoms = 9>
<Molecule: id = 93 name = HIS chain_id =    natoms = 10>
<Molecule: id = 94 name = ASN chain_id =    natoms = 8>
<Molecule: id = 95 name = GLY chain_id =    natoms = 4>
<Molecule: id = 96 name = GLY chain_id =    natoms = 4>
<Molecule: id = 97 name = PRO chain_id =    natoms = 7>
<Molecule: id = 98 name = ARG chain_id =    natoms = 11>
<Molecule: id = 99 name = GLY chain_id =    natoms = 4>
<Molecule: id = 100 name = CYS chain_id =    natoms = 6>
<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 [3]:
# создаём отдельные объекты для белка и для лиганда
newpdb = pdb.copy()
for r in newpdb.residues[122:]:
    newpdb.remove_residue(r)
lig = pdb.copy()
for r in lig.residues[:122]:
    lig.remove_residue(r)
In [4]:
# Поиск геометрического центра лиганда
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 x_mean, y_mean, z_mean
[36.689, 47.293, 29.716]
[36.672, 48.845, 29.982]
[35.598, 49.239, 30.974]
[34.179, 48.547, 30.724]
[34.404, 47.047, 30.511]
[33.216, 46.355, 30.036]
[38.258, 49.984, 29.171]
[39.464, 50.596, 29.909]
[37.76, 49.221, 30.158]
[35.188, 50.634, 30.8]
[33.23, 48.95, 31.724]
[35.383, 46.765, 29.494]
[32.564, 46.986, 28.986]
[37.999, 50.322, 27.991]
[39.575, 43.993, 26.494]
[38.019, 43.705, 26.527]
[37.294, 44.947, 27.033]
[37.973, 45.65, 28.297]
[39.474, 45.739, 28.148]
[40.159, 46.084, 29.416]
[37.036, 41.852, 25.66]
[36.7, 41.623, 24.184]
[37.512, 43.114, 25.714]
[35.959, 44.725, 27.5]
[37.354, 46.965, 28.352]
[39.97, 44.403, 27.861]
[39.771, 45.374, 30.563]
[36.826, 40.862, 26.385]
[43.626, 41.355, 25.057]
[43.549, 42.888, 25.204]
[42.335, 43.282, 26.047]
[41.006, 42.702, 25.544]
[41.035, 41.268, 24.979]
[40.232, 41.314, 23.71]
[45.412, 44.036, 26.199]
[45.348, 45.104, 25.068]
[42.356, 40.854, 24.514]
[42.25, 44.702, 26.357]
[40.067, 42.675, 26.662]
[40.611, 42.393, 22.802]
[46.489, 44.136, 26.928]
[44.457, 43.063, 26.447]
[43.847, 40.961, 26.147]
38.9964186047 45.1291395349 27.5343023256

Т. е. координаты геометрического центра лиганда - 38.9964186047 45.1291395349 27.5343023256. Эта информация нам пригодится далее для запуска докинга.

In [42]:
# Сохраняем белок без лиганда в отдельный pdb-файл
newpdb.writePDB("prot_without_lig.pdb")
In [5]:
# Готовим белок для докинга
prot = oddt.toolkit.readfile('pdb','prot_without_lig.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 

Молекулярная масса нашего белка - 13375.39512

Далее подготовим лиганды для докинга (NAG и его модификации, в которых вместо метильного радикала СH3C(=O)NH группы - OH, NH3+, H, Ph, COO-)

In [6]:
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','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 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 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 [7]:
#Создание объекта для докинга
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_1Yn35e
--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 (y,z): геометрический центр лиганда
  • ---size_x (y,z): размер области, в которой будет производиться поиск
  • ---cpu: сколько процессоров будет использоваться для выполнения нашей задачи
  • --exhaustiveness: параметр, который говорит программе, насколько тщательно производить поиск
  • --num_modes: количество конформаций лиганда в активном центре, которые будут найдены программой
  • --energy_range: максимальная разница по энергии между лучшей и худшей конформацией из найденных
In [39]:
#запуск докинга
res = dock_obj.dock(mols,prot)
In [40]:
# проанализируем результаты докинга
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 [41]:
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])
  45 C13H17NO6    -7.6     0.0     UNL   8   1   4
  36  C8H13NO8    -6.7     0.0     UNL  10   0   0
  46 C13H17NO6    -6.7    7.81     UNL   6   0   7
  47 C13H17NO6    -6.7    3.02     UNL   3   0   1
   0  C8H15NO6    -6.5     0.0     UNL  10   0   2
   9  C7H13NO7    -6.5     0.0     UNL  11   0   0
  48 C13H17NO6    -6.5    8.13     UNL   5   0   2
  27 C7H16N2O6    -6.4     0.0     UNL  10   0   0
  28 C7H16N2O6    -6.4   2.847     UNL  10   0   0
   1  C8H15NO6    -6.3   2.723     UNL   8   0   0
  10  C7H13NO7    -6.3   2.888     UNL  11   0   0
  29 C7H16N2O6    -6.2   4.353     UNL   6   0   0
  30 C7H16N2O6    -6.2   4.854     UNL   8   0   0
  37  C8H13NO8    -6.2   5.737     UNL  10   0   0
  38  C8H13NO8    -6.2   2.413     UNL  10   0   0
  11  C7H13NO7    -6.1    4.01     UNL  11   0   0
  12  C7H13NO7    -6.1   3.841     UNL  10   0   0
  18  C7H13NO6    -6.1     0.0     UNL  11   0   0
  31 C7H16N2O6    -6.1   2.606     UNL   9   0   0
  39  C8H13NO8    -6.1   4.304     UNL   9   0   0
  49 C13H17NO6    -6.1   7.435     UNL   7   0   3
  50 C13H17NO6    -6.1    4.23     UNL   6   0   4
   2  C8H15NO6    -6.0   4.035     UNL   7   0   2
  13  C7H13NO7    -6.0   4.691     UNL  12   0   0
  14  C7H13NO7    -6.0    3.58     UNL   8   0   0
  15  C7H13NO7    -6.0   5.002     UNL  10   0   0
  19  C7H13NO6    -6.0   4.128     UNL   9   0   0
  32 C7H16N2O6    -6.0   2.337     UNL   7   0   0
   3  C8H15NO6    -5.9   2.403     UNL   6   0   0
  20  C7H13NO6    -5.8   2.661     UNL   9   0   0
  33 C7H16N2O6    -5.8   3.102     UNL   7   0   0
  21  C7H13NO6    -5.7   4.271     UNL   7   0   0
  22  C7H13NO6    -5.7   4.036     UNL   5   0   0
  40  C8H13NO8    -5.7   7.531     UNL  10   0   0
   4  C8H15NO6    -5.6    4.97     UNL   5   0   1
   5  C8H15NO6    -5.6    3.68     UNL   6   0   0
  16  C7H13NO7    -5.6   3.655     UNL  10   0   0
  23  C7H13NO6    -5.6   2.341     UNL   7   0   0
  24  C7H13NO6    -5.6   4.081     UNL   9   0   0
  34 C7H16N2O6    -5.6   3.045     UNL   7   0   0
  51 C13H17NO6    -5.6   8.155     UNL   7   0   4
  52 C13H17NO6    -5.6   7.906     UNL   5   0   1
   6  C8H15NO6    -5.5    4.05     UNL   6   0   0
  17  C7H13NO7    -5.5   3.356     UNL   9   0   0
  25  C7H13NO6    -5.5   3.249     UNL   7   0   0
  35 C7H16N2O6    -5.5   8.397     UNL   9   0   0
  41  C8H13NO8    -5.5   7.867     UNL  10   0   0
  26  C7H13NO6    -5.4    5.07     UNL   7   0   0
  53 C13H17NO6    -5.4   7.811     UNL   8   0   0
  42  C8H13NO8    -5.2   8.773     UNL   7   0   0
   7  C8H15NO6    -5.1  10.061     UNL  10   0   0
   8  C8H15NO6    -5.1    8.53     UNL   8   0   0
  43  C8H13NO8    -5.1   9.068     UNL   8   0   0
  44  C8H13NO8    -5.0   9.692     UNL   6   0   0

Наилучшая структура – структура пд номером 45, на месте заместителя – фенил-радикал.

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

Исходные лиганды структуры выделены зелёным цветом. Сиреневым цветов выделен лиганд из наихудшей структуры (44), синим (36) и красным (45) - из двух наилучших. Обе лучшие структуры попадают глубоко в карман белка в отличие от исходных лигандов.