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

Цель данного занятия ознакомится с возможностями докинга низкомолекулярного лиганда в структуру белка.
Объект - белок лизоцим, структуру которого постронеа на основе гомологичного моделирования

In [4]:
#Модели
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 [6]:
pdb=pmx.Model('WithLigand.B99990001.pdb')

for r in pdb.residues[135:]:
    print r 
<Molecule: id = 136 name = PHE chain_id =    natoms = 11>
<Molecule: id = 137 name = ASN chain_id =    natoms = 8>
<Molecule: id = 138 name = LYS chain_id =    natoms = 9>
<Molecule: id = 139 name = PRO chain_id =    natoms = 7>
<Molecule: id = 140 name = VAL chain_id =    natoms = 7>
<Molecule: id = 141 name = SER chain_id =    natoms = 6>
<Molecule: id = 142 name = SER chain_id =    natoms = 6>
<Molecule: id = 143 name = ASN chain_id =    natoms = 8>
<Molecule: id = 144 name = SER chain_id =    natoms = 6>
<Molecule: id = 145 name = ASN chain_id =    natoms = 8>
<Molecule: id = 146 name = GLN chain_id =    natoms = 9>
<Molecule: id = 147 name = ASN chain_id =    natoms = 8>
<Molecule: id = 148 name = ASN chain_id =    natoms = 8>
<Molecule: id = 149 name = GLN chain_id =    natoms = 9>
<Molecule: id = 150 name = THR chain_id =    natoms = 7>
<Molecule: id = 151 name = GLY chain_id =    natoms = 4>
<Molecule: id = 152 name = GLY chain_id =    natoms = 4>
<Molecule: id = 153 name = MET chain_id =    natoms = 8>
<Molecule: id = 154 name = ILE chain_id =    natoms = 8>
<Molecule: id = 155 name = LYS chain_id =    natoms = 9>
<Molecule: id = 156 name = MET chain_id =    natoms = 8>
<Molecule: id = 157 name = TYR chain_id =    natoms = 12>
<Molecule: id = 158 name = LEU chain_id =    natoms = 8>
<Molecule: id = 159 name = ILE chain_id =    natoms = 8>
<Molecule: id = 160 name = ILE chain_id =    natoms = 8>
<Molecule: id = 161 name = GLY chain_id =    natoms = 4>
<Molecule: id = 162 name = LEU chain_id =    natoms = 8>
<Molecule: id = 163 name = ASP chain_id =    natoms = 8>
<Molecule: id = 164 name = ASN chain_id =    natoms = 8>
<Molecule: id = 165 name = SER chain_id =    natoms = 6>
<Molecule: id = 166 name = GLY chain_id =    natoms = 4>
<Molecule: id = 167 name = LYS chain_id =    natoms = 9>
<Molecule: id = 168 name = ALA chain_id =    natoms = 5>
<Molecule: id = 169 name = LYS chain_id =    natoms = 9>
<Molecule: id = 170 name = HIS chain_id =    natoms = 10>
<Molecule: id = 171 name = TRP chain_id =    natoms = 14>
<Molecule: id = 172 name = TYR chain_id =    natoms = 12>
<Molecule: id = 173 name = VAL chain_id =    natoms = 7>
<Molecule: id = 174 name = SER chain_id =    natoms = 6>
<Molecule: id = 175 name = ASP chain_id =    natoms = 8>
<Molecule: id = 176 name = GLY chain_id =    natoms = 4>
<Molecule: id = 177 name = VAL chain_id =    natoms = 7>
<Molecule: id = 178 name = SER chain_id =    natoms = 6>
<Molecule: id = 179 name = VAL chain_id =    natoms = 7>
<Molecule: id = 180 name = ARG chain_id =    natoms = 11>
<Molecule: id = 181 name = HIS chain_id =    natoms = 10>
<Molecule: id = 182 name = VAL chain_id =    natoms = 7>
<Molecule: id = 183 name = ARG chain_id =    natoms = 11>
<Molecule: id = 184 name = THR chain_id =    natoms = 7>
<Molecule: id = 185 name = ILE chain_id =    natoms = 8>
<Molecule: id = 186 name = ARG chain_id =    natoms = 11>
<Molecule: id = 187 name = MET chain_id =    natoms = 8>
<Molecule: id = 188 name = LEU chain_id =    natoms = 8>
<Molecule: id = 189 name = GLU chain_id =    natoms = 9>
<Molecule: id = 190 name = ASN chain_id =    natoms = 8>
<Molecule: id = 191 name = TYR chain_id =    natoms = 12>
<Molecule: id = 192 name = GLN chain_id =    natoms = 9>
<Molecule: id = 193 name = ASN chain_id =    natoms = 8>
<Molecule: id = 194 name = LYS chain_id =    natoms = 9>
<Molecule: id = 195 name = TRP chain_id =    natoms = 14>
<Molecule: id = 196 name = ALA chain_id =    natoms = 5>
<Molecule: id = 197 name = LYS chain_id =    natoms = 9>
<Molecule: id = 198 name = LEU chain_id =    natoms = 8>
<Molecule: id = 199 name = ASN chain_id =    natoms = 8>
<Molecule: id = 200 name = LEU chain_id =    natoms = 8>
<Molecule: id = 201 name = PRO chain_id =    natoms = 7>
<Molecule: id = 202 name = VAL chain_id =    natoms = 7>
<Molecule: id = 203 name = ASP chain_id =    natoms = 8>
<Molecule: id = 204 name = THR chain_id =    natoms = 7>
<Molecule: id = 205 name = MET chain_id =    natoms = 8>
<Molecule: id = 206 name = PHE chain_id =    natoms = 11>
<Molecule: id = 207 name = ILE chain_id =    natoms = 8>
<Molecule: id = 208 name = ALA chain_id =    natoms = 5>
<Molecule: id = 209 name = GLU chain_id =    natoms = 9>
<Molecule: id = 210 name = ILE chain_id =    natoms = 8>
<Molecule: id = 211 name = GLU chain_id =    natoms = 9>
<Molecule: id = 212 name = ALA chain_id =    natoms = 5>
<Molecule: id = 213 name = GLU chain_id =    natoms = 9>
<Molecule: id = 214 name = PHE chain_id =    natoms = 11>
<Molecule: id = 215 name = GLY chain_id =    natoms = 4>
<Molecule: id = 216 name = ARG chain_id =    natoms = 11>
<Molecule: id = 217 name = LYS chain_id =    natoms = 9>
<Molecule: id = 218 name = ILE chain_id =    natoms = 8>
<Molecule: id = 219 name = ASP chain_id =    natoms = 8>
<Molecule: id = 220 name = MET chain_id =    natoms = 8>
<Molecule: id = 221 name = ALA chain_id =    natoms = 5>
<Molecule: id = 222 name = SER chain_id =    natoms = 6>
<Molecule: id = 223 name = GLY chain_id =    natoms = 4>
<Molecule: id = 224 name = GLU chain_id =    natoms = 9>
<Molecule: id = 225 name = VAL chain_id =    natoms = 7>
<Molecule: id = 226 name = LYS chain_id =    natoms = 10>
<Molecule: id = 227 name = NAG chain_id =    natoms = 14>
<Molecule: id = 228 name = NAG chain_id =    natoms = 14>
<Molecule: id = 229 name = NDG chain_id =    natoms = 15>
In [7]:
# создание объектов белок и лиганда
newpdb = pdb.copy()
lig = pdb.copy()
del newpdb.residues[226:]
del lig.residues[:226]
In [8]:
# был найден геометрический центр лиганда
x_c = 0
y_c = 0
z_c = 0

for a in lig.atoms:
    x_c = x_c+a.x[0]
    y_c = y_c+a.x[1]
    z_c = z_c+a.x[2]
    
x_c = x_c / len(lig.atoms)
y_c = y_c / len(lig.atoms)
z_c = z_c / len(lig.atoms)

print  "center ["+str(x_c)+", "+str(x_c)+", "+str(x_c)+"]" 
center [25.3029178982, 25.3029178982, 25.3029178982]
In [9]:
newpdb.writePDB("ProteinDelResid.pdb")

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

In [10]:
prot = oddt.toolkit.readfile('pdb','ProteinDelResid.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: 25999.2393

Создание лигандов, где метильный радикал этой группы будет заменён на [OH,NH3+,H,Ph,COO-]

In [12]:
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem.Draw import IPythonConsole 
from rdkit.Chem import Draw
import numpy as np
from IPython.display import display,Image

ibu=Chem.MolFromSmiles('CC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu)
display(ibu)

ibu1=Chem.MolFromSmiles('OC(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu1)
display(ibu1)

ibu1=Chem.MolFromSmiles('C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu1)
display(ibu1)

ibu1=Chem.MolFromSmiles('[NH3+]C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu1)
display(ibu1)

ibu1=Chem.MolFromSmiles('[O-]C(=O)C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu1)
display(ibu1)

ibu1=Chem.MolFromSmiles('c2ccccc2C(=O)N[C@H]1[C@H](O)O[C@H](CO)[C@@H](O)[C@@H]1O')
AllChem.Compute2DCoords(ibu1)
display(ibu1)
In [14]:
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','c2ccccc2C(=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 [17]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=[25.3029178982, 25.3029178982, 25.3029178982],
    executable='/usr/bin/vina',autocleanup=True, num_modes=20)

print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
/tmp/autodock_vina_XOZTHS
--center_x 25.3029178982 --center_y 25.3029178982 --center_z 25.3029178982 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3

Докинг и его результаты

In [18]:
# do it
res = dock_obj.dock(mols,prot)
In [19]:
out= open("results.txt","w")
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)
    out.write("%4d%10s%8s%8s%8s\n" %(i,r.formula, r.data['vina_affinity'],  r.data['vina_rmsd_ub'], r.residues[0].name))
   0  C8H15NO6    -2.9   0.000     UNL
   1  C8H15NO6    -2.5   5.344     UNL
   2  C8H15NO6    -2.4   4.475     UNL
   3  C8H15NO6    -2.4   3.122     UNL
   4  C8H15NO6    -2.4   3.261     UNL
   5  C8H15NO6    -2.4   2.723     UNL
   6  C8H15NO6    -2.3  10.227     UNL
   7  C8H15NO6    -2.3   4.729     UNL
   8  C8H15NO6    -2.2   5.486     UNL
   9  C7H13NO7    -2.9   0.000     UNL
  10  C7H13NO7    -2.6   5.151     UNL
  11  C7H13NO7    -2.6   3.056     UNL
  12  C7H13NO7    -2.5   5.136     UNL
  13  C7H13NO7    -2.4   5.760     UNL
  14  C7H13NO7    -2.3   2.697     UNL
  15  C7H13NO7    -2.3   4.546     UNL
  16  C7H13NO7    -2.3   4.912     UNL
  17  C7H13NO7    -2.3   4.872     UNL
  18  C7H13NO6    -2.5   0.000     UNL
  19  C7H13NO6    -2.5   2.353     UNL
  20  C7H13NO6    -2.4   4.226     UNL
  21  C7H13NO6    -2.4   4.418     UNL
  22  C7H13NO6    -2.4   2.324     UNL
  23  C7H13NO6    -2.3   4.519     UNL
  24  C7H13NO6    -2.2   4.640     UNL
  25  C7H13NO6    -2.1   4.524     UNL
  26  C7H13NO6    -2.1   3.810     UNL
  27 C7H15N2O6    -3.0   0.000     UNL
  28 C7H15N2O6    -2.6   5.127     UNL
  29 C7H15N2O6    -2.6   3.269     UNL
  30 C7H15N2O6    -2.6   2.830     UNL
  31 C7H15N2O6    -2.5   2.464     UNL
  32 C7H15N2O6    -2.4   2.887     UNL
  33 C7H15N2O6    -2.4   3.048     UNL
  34 C7H15N2O6    -2.4   6.310     UNL
  35 C7H15N2O6    -2.4   4.451     UNL
  36  C8H13NO8    -2.8   0.000     UNL
  37  C8H13NO8    -2.8   3.873     UNL
  38  C8H13NO8    -2.7   5.130     UNL
  39  C8H13NO8    -2.7   1.889     UNL
  40  C8H13NO8    -2.5   4.800     UNL
  41  C8H13NO8    -2.5   3.971     UNL
  42  C8H13NO8    -2.4   5.554     UNL
  43  C8H13NO8    -2.4   5.221     UNL
  44  C8H13NO8    -2.3   5.388     UNL
  45 C13H17NO6    -3.4   0.000     UNL
  46 C13H17NO6    -3.2   3.281     UNL
  47 C13H17NO6    -3.2   5.374     UNL
  48 C13H17NO6    -3.1   4.956     UNL
  49 C13H17NO6    -3.1   6.142     UNL
  50 C13H17NO6    -2.9   4.746     UNL
  51 C13H17NO6    -2.9   5.245     UNL
  52 C13H17NO6    -2.8   5.473     UNL
  53 C13H17NO6    -2.8   8.450     UNL
In [49]:
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 "%8s%8s%8s" %(len(hbs[0]),len(stack[0]),len(phob[0]))
#распечатать с результатами докинга
       2       0       0
       1       0       1
       1       0       0
       1       0       1
       3       0       0
       1       0       3
       2       0       2
       1       0       3
       1       0       0
       3       0       0
       2       0       0
       3       0       0
       1       0       0
       2       0       0
       3       0       0
       1       0       0
       1       0       0
       1       0       0
       2       0       0
       2       0       0
       2       0       0
       1       0       0
       2       0       0
       1       0       0
       1       0       0
       1       0       0
       1       0       0
       2       0       0
       1       0       0
       2       0       0
       2       0       0
       2       0       0
       1       0       0
       1       0       0
       0       0       0
       1       0       0
       5       0       0
       3       0       0
       2       0       0
       1       0       0
       2       0       0
       2       0       0
       1       0       0
       5       0       0
       2       0       0
       1       0       3
       3       0       5
       1       0      11
       1       0       5
       1       0      10
       1       0      11
       3       0       8
       1       0       5
       1       0       4

Отсортировали от лучшего заместителя к худшему, результат представлен в таблице.

45    C13H17NO6    -3.4    0        UNL
46    C13H17NO6    -3.3    5.162    UNL
47    C13H17NO6    -3.1    5.253    UNL
48    C13H17NO6    -3.1    2.288    UNL
49    C13H17NO6    -3.1    4.802    UNL
00    С8H15NO6    -2.9    0        UNL
09    C7H13NO7    -2.9    0        UNL
27    C7H15N2O6    -2.9    0        UNL
50    C13H17NO6    -2.9    5.453    UNL
36    C8H13NO8    -2.8    0        UNL
51    C13H17NO6    -2.8    9.345    UNL
52    C13H17NO6    -2.8    3.049    UNL
53    C13H17NO6    -2.8    8.249    UNL
18    C7H13NO6    -2.7    0        UNL
37    C8H13NO8    -2.7    3.224    UNL
01    C8H15NO6    -2.6    5.28    UNL
10    C7H13NO7    -2.6    5.136    UNL
11    C7H13NO7    -2.6    3.203    UNL
19    C7H13NO6    -2.6    3.588    UNL
20    C7H13NO6    -2.6    4.748    UNL
38    C8H13NO8    -2.6    4.939    UNL
39    C8H13NO8    -2.6    5.084    UNL
40    C8H13NO8    -2.6    6.975    UNL
02    C8H15NO6    -2.5    2.823    UNL
03    C8H15NO6    -2.5    5.231    UNL
04    C8H15NO6    -2.5    5.357    UNL
05    C8H15NO6    -2.5    3.209    UNL
06    C8H15NO6    -2.5    4.75    UNL
12    C7H13NO7    -2.5    3.039    UNL
13    C7H13NO7    -2.5    2.271    UNL
21    C7H13NO6    -2.5    2.498    UNL
28    C7H15N2O6    -2.5    5.169    UNL
41    C8H13NO8    -2.5    5.924    UNL
42    C8H13NO8    -2.5    3.506    UNL
43    C8H13NO8    -2.5    2.213    UNL
44    C8H13NO8    -2.5    3.844    UNL
07    C8H15NO6    -2.4    4.467    UNL
14    C7H13NO7    -2.4    1.946    UNL
29    C7H15N2O6    -2.4    3.211    UNL
30    C7H15N2O6    -2.4    1.879    UNL
08    C8H15NO6    -2.3    2.376    UNL
15    C7H13NO7    -2.3    4.97    UNL
16    C7H13NO7    -2.3    2.861    UNL
17    C7H13NO7    -2.3    4.424    UNL
22    C7H13NO6    -2.3    4.571    UNL
23    C7H13NO6    -2.3    2.481    UNL
24    C7H13NO6    -2.3    4.998    UNL
31    C7H15N2O6    -2.3    4.56    UNL
32    C7H15N2O6    -2.3    4.507    UNL
25    C7H13NO6    -2.2    4.572    UNL
33    C7H15N2O6    -2.2    4.787    UNL
34    C7H15N2O6    -2.2    4.266    UNL
35    C7H15N2O6    -2.2    2.119    UNL
26    C7H13NO6    -2.1    3.564    UNL

Лучшим заместителем оказался C13H17NO6 с лигандом - Ph. Ни в каком из соединений нет стакинг взаимодействий.

In [24]:
for i,r in enumerate(res):
    r.write(filename='r%s.pdb' % i, format='pdb')
In [25]:
# на рисунке представлен лучший лиганд
import nglview
import ipywidgets
w1 = nglview.show_structure_file('r45.pdb')
w1