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

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 import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

import pmx # Модуль для манипулирования pdb 

Подготовка белка

In [2]:
pdb=pmx.Model('seq_lig.B99990001.pdb')

for r in pdb.residues[-10:]:
    print r # посмотрим остатки
<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 [3]:
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-3]

Вычислим координаты центра лиганда

In [5]:
atom_coords = np.zeros((len(lig.atoms), 3))
for x, a in enumerate(lig.atoms):
   atom_coords[x, :] = a.x
lig_center = np.mean(atom_coords, axis=0)
lig_center
Out[5]:
array([42.4899741 , 34.20196614, 18.99605876])
In [6]:
newpdb.writePDB('myprot.pdb')

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

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

Забавно, что оддт не считает наш белок белком :)

Лиганды для докинга

Создадим несколько производных NAG (замещаем один из метилов)

In [8]:
smiles = ['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)(c1ccccc1)',
          'NC1C(C(C(OC1O)CO)O)O', 'NC1C(C(C(OC1O)CO)O)OC', '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)
OBDepict ***** O N O O O O O O H H H H H H
OBDepict ***** O N O O O O O N H H H H H H H
OBDepict ***** O N O O O O O H H H H H
OBDepict ***** O N O O O O O H H H H H
OBDepict ***** N O O O O O H H H H H H
OBDepict ***** N O O O O O CH 3 H H H H H
OBDepict ***** O N O O O O O O O H H H H H H

Докинг

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

print dock_obj.tmp_dir
print " ".join(dock_obj.params)
/tmp/autodock_vina_4Grt9L
--center_x 42.489974103585624 --center_y 34.20196613545814 --center_z 18.99605876494022 --size_x 20 --size_y 20 --size_z 20 --exhaustiveness 8 --num_modes 9 --energy_range 3

Параметры Вины (по порядку): Первые 6 параметров задают бокс, внутри которого могут передвигаться атомы (в нашем случае только лиганда); exhaustiveness задает насколько долго мы ищем минимум (который ищется эвристически); num_modes задает сколько моделей должно быть в выдаче; energy_range задает масимальную разницу энергий (в ккал/моль) между лучшей и худшей моделью для выдачи.

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

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

In [11]:
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  C7H13NO7    -6.4   0.000     LIG
   1  C7H13NO7    -6.3   5.110     LIG
   2  C7H13NO7    -6.2   3.347     LIG
   3  C7H13NO7    -5.6   4.180     LIG
   4  C7H13NO7    -5.6   2.764     LIG
   5  C7H13NO7    -5.6   4.810     LIG
   6  C7H13NO7    -5.2   5.521     LIG
   7  C7H13NO7    -5.2   5.024     LIG
   8  C7H13NO7    -5.1   4.939     LIG
   9 C7H14N2O6    -6.2   0.000     LIG
  10 C7H14N2O6    -6.0   4.823     LIG
  11 C7H14N2O6    -5.7   2.765     LIG
  12 C7H14N2O6    -5.6   4.885     LIG
  13 C7H14N2O6    -5.6   3.145     LIG
  14 C7H14N2O6    -5.4   5.592     LIG
  15 C7H14N2O6    -5.4   2.857     LIG
  16 C7H14N2O6    -5.4   5.594     LIG
  17 C7H14N2O6    -5.1   5.194     LIG
  18  C7H13NO6    -6.0   0.000     LIG
  19  C7H13NO6    -5.9   4.839     LIG
  20  C7H13NO6    -5.6   4.899     LIG
  21  C7H13NO6    -5.4   4.585     LIG
  22  C7H13NO6    -5.0   3.507     LIG
  23  C7H13NO6    -4.9   2.103     LIG
  24  C7H13NO6    -4.8   4.836     LIG
  25  C7H13NO6    -4.7   5.576     LIG
  26  C7H13NO6    -4.7   4.835     LIG
  27 C13H17NO6    -7.3   0.000     LIG
  28 C13H17NO6    -5.9   1.708     LIG
  29 C13H17NO6    -5.7   2.636     LIG
  30 C13H17NO6    -5.5   2.692     LIG
  31 C13H17NO6    -5.3   6.981     LIG
  32 C13H17NO6    -4.5   3.936     LIG
  33 C13H17NO6    -4.5   4.366     LIG
  34  C6H13NO5    -5.6   0.000     LIG
  35  C6H13NO5    -5.4   3.716     LIG
  36  C6H13NO5    -5.1   2.463     LIG
  37  C6H13NO5    -5.0   3.028     LIG
  38  C6H13NO5    -4.9   3.218     LIG
  39  C6H13NO5    -4.7   3.269     LIG
  40  C6H13NO5    -4.5   2.210     LIG
  41  C6H13NO5    -4.4   3.069     LIG
  42  C6H13NO5    -4.3   4.357     LIG
  43  C7H15NO5    -5.3   0.000     LIG
  44  C7H15NO5    -5.2   3.709     LIG
  45  C7H15NO5    -5.2   4.030     LIG
  46  C7H15NO5    -4.9   4.097     LIG
  47  C7H15NO5    -4.8   3.813     LIG
  48  C7H15NO5    -4.8   4.078     LIG
  49  C7H15NO5    -4.7   4.311     LIG
  50  C7H15NO5    -4.7   3.337     LIG
  51  C7H15NO5    -4.5   4.420     LIG
  52  C8H13NO8    -6.1   0.000     LIG
  53  C8H13NO8    -5.9   5.755     LIG
  54  C8H13NO8    -5.8   4.304     LIG
  55  C8H13NO8    -5.5   5.759     LIG
  56  C8H13NO8    -5.2   5.477     LIG
  57  C8H13NO8    -5.2   2.824     LIG
  58  C8H13NO8    -5.1   5.604     LIG
  59  C8H13NO8    -5.0   2.293     LIG
  60  C8H13NO8    -4.6   5.500     LIG

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

Мы посмотрим на самого лучшего для каждого из соединений

In [20]:
selected = [0, 9, 18, 27, 34, 43, 52]
for i,r in enumerate(res):
    if i not in selected:
        continue
    
    mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
    hbonds_num = len(mol1)
    hbonds_strict = strict.sum()
    
    r1, r2, strict_par, strict_per = oddt.interactions.pi_stacking(prot,r)
    stack_num = len(r1)
    stack_par = strict_par.sum()
    stack_per = strict_per.sum()
    
    mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
    hydroph_num = len(mol1)
    
    print i, hbonds_num, hbonds_strict, stack_num, stack_par, stack_per, hydroph_num
0 13 13 0 0 0 0
9 13 12 0 0 0 0
18 12 10 0 0 0 0
27 11 7 0 0 0 6
34 11 7 0 0 0 0
43 8 8 0 0 0 0
52 17 13 0 0 0 0

Видно, что гидрофобные контакты есть только у 27; больше всего водородных связей у 52, однако хороших водородных связей у него столько же, сколько и у 0, а у 9 их почти чтолько же. Стэкинг никто не образует.

Визуализация

Мы посмотрим на два самых лучших:

 0  C7H13NO7    -6.4   0.000     LIG
27 C13H17NO6    -7.3   0.000     LIG
In [21]:
selected = [0, 27]
for i,r in enumerate(res):
    if i not in selected:
        continue

    r.write(filename='r%s.pdb' % i, format='pdb')
In [2]:
display(Image('r0.png'))
display(Image('r27.png'))

Как видно из рисунков, самым афинным лигандом оказался NAG с фенилаланином, а вторым -- NAG с гидроксилом. Они оба находятся примерно в одной позе. Пимоль почему-то находит гораздо меньше водородных связей, чем оддт (там у 0 было аж 13). Ароматическое кольцо фенилаланина так же не было изуродовано. В общем и целом нормально все задочилось.