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

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 import Chem
from rdkit.Chem.Draw import IPythonConsole
In [2]:
import mdtraj as md

Загрузим файл из прошлого практикума

In [3]:
u = md.load('sequence.B99990006.pdb')
pdb = u.topology

Сохраним файл без лиганда

In [4]:
u2 = u.atom_slice(range(1096))
u2.save_pdb('sequence.B99990006-clean.pdb', force_overwrite=True)

Найдём геометрический центр лиганда

In [5]:
lig = u.atom_slice(range(1096,1138))
lig_center = np.mean(lig.xyz[0], axis=0) * 10
lig_center
Out[5]:
array([49.79834 , 38.615234, 29.907879], dtype=float32)

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

In [6]:
prot = next(oddt.toolkit.readfile('pdb','sequence.B99990006-clean.pdb'))
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: 15662.003999999964
In [7]:
smiles = ["CC(=O)NC1C(O)OC(CO)C(O)C1O", # NAG
          "OC(=O)NC1C(O)OC(CO)C(O)C1O", # OH
          "[NH3+]C(=O)NC1C(O)OC(CO)C(O)C1O", # NH3+
          "C(=O)NC1C(O)OC(CO)C(O)C1O", # H
          "C=1(C=CC=CC1)C(=O)NC1C(O)OC(CO)C(O)C1O", # Ph
          "C(C([O-])=O)(=O)NC1C(O)OC(CO)C(O)C1O"] # COO-
mols = []
images =[]

for s in smiles:    
    m = oddt.toolkit.readstring('smi', s)
    # m.addh()
    m.make3D(forcefield='mmff94', steps=150)
    mols.append(m)
    images.append((SVG(copy.deepcopy(m).write('svg'))))
    
display_svg(*images)

RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()
RDKit WARNING: [11:37:00] Molecule does not have explicit Hs. Consider calling AddHs()

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

In [8]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(protein=prot,
                                                  size=(20, 20, 20),
                                                  center=lig_center)

print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params))

/tmp/autodock_vina_cacqv7lp
--center_x 49.79834 --center_y 38.615234 --center_z 29.907879 --size_x 20 --size_y 20 --size_z 20 --exhaustiveness 8 --num_modes 9 --energy_range 3

center_x, center_y, center_z - координаты центра области для докинга; size_x, size_y, size_z - размер области для докинга; exhaustiveness - то, насколько "дотошно" будет происходить поиск. Иными словами параметр, который указывает, как долго будет происходить поиск минимума энергии; num_mode - число вариантов, которое найдёт программа; energy_range - в каком диапазоне от лучшего варианта будут находиться остальные или размер отсечки по энергии. Варианты хуже не будут рассматриваться.

In [9]:
res = dock_obj.dock(mols, prot)

Анализ результатов докинга

Посчитаем определённое количество метрик для каждого лиганда и отсортируем по affinity:

In [16]:
data = []

for i,r in enumerate(res):
    data.append([
        Chem.rdMolDescriptors.CalcMolFormula(r.Mol),
        r.data['vina_affinity'],
        r.data['vina_rmsd_ub'],
        len(oddt.interactions.hbonds(prot,r)[0]),
        len(oddt.interactions.pi_stacking(prot,r)[0]),
        len(oddt.interactions.hydrophobic_contacts(prot,r)[0])
    ])

data.sort(key=lambda x:x[0])

import pandas as pd

table = pd.DataFrame(data, columns=['Formula', 'Affinity', 'RMSD', 'H bonds', 'Stacking', 'Hydrophobic contacts'])
table
Out[16]:
Formula Affinity RMSD H bonds Stacking Hydrophobic contacts
0 C13H17NO6 -6.6 0.000 10 0 5
1 C13H17NO6 -6.6 3.628 3 1 10
2 C13H17NO6 -6.3 2.943 5 0 10
3 C13H17NO6 -6.2 3.153 2 1 18
4 C13H17NO6 -6.1 3.707 2 1 13
5 C13H17NO6 -6.0 7.148 3 0 3
6 C13H17NO6 -5.6 6.957 1 0 3
7 C13H17NO6 -5.6 6.877 2 0 0
8 C13H17NO6 -5.5 7.561 6 1 13
9 C7H13NO6 -4.6 0.000 6 0 0
10 C7H13NO6 -4.4 4.344 5 0 0
11 C7H13NO6 -4.4 5.015 4 0 0
12 C7H13NO6 -4.4 7.333 4 0 0
13 C7H13NO6 -4.4 4.393 3 0 0
14 C7H13NO6 -4.3 4.829 6 0 0
15 C7H13NO6 -4.3 4.491 5 0 0
16 C7H13NO6 -4.3 4.061 4 0 0
17 C7H13NO6 -4.1 3.567 5 0 0
18 C7H13NO7 -5.2 0.000 4 0 0
19 C7H13NO7 -5.1 3.007 4 0 0
20 C7H13NO7 -5.0 3.898 7 0 0
21 C7H13NO7 -4.9 2.300 2 0 0
22 C7H13NO7 -4.7 5.080 3 0 0
23 C7H13NO7 -4.6 5.583 3 0 0
24 C7H13NO7 -4.5 3.062 5 0 0
25 C7H13NO7 -4.5 3.386 3 0 0
26 C7H13NO7 -4.5 3.488 4 0 0
27 C7H15N2O6+ -5.2 0.000 8 0 0
28 C7H15N2O6+ -4.9 2.579 4 0 0
29 C7H15N2O6+ -4.8 5.501 4 0 0
30 C7H15N2O6+ -4.7 4.186 3 0 0
31 C7H15N2O6+ -4.5 5.548 2 0 0
32 C7H15N2O6+ -4.5 5.407 2 0 0
33 C7H15N2O6+ -4.5 5.850 6 0 0
34 C7H15N2O6+ -4.4 3.363 2 0 0
35 C7H15N2O6+ -4.4 4.057 3 0 0
36 C8H12NO8- -5.8 0.000 4 0 0
37 C8H12NO8- -5.5 3.018 5 0 0
38 C8H12NO8- -5.4 3.042 4 0 0
39 C8H12NO8- -5.3 5.842 3 0 0
40 C8H12NO8- -5.3 2.257 6 0 0
41 C8H12NO8- -5.2 6.254 2 0 0
42 C8H12NO8- -5.1 6.047 3 0 0
43 C8H12NO8- -5.0 3.679 3 0 0
44 C8H12NO8- -4.9 6.005 4 0 0
45 C8H15NO6 -5.1 0.000 1 0 1
46 C8H15NO6 -4.9 5.758 3 0 1
47 C8H15NO6 -4.9 1.454 1 0 4
48 C8H15NO6 -4.8 5.235 1 0 0
49 C8H15NO6 -4.8 4.472 5 0 0
50 C8H15NO6 -4.5 5.451 6 0 2
51 C8H15NO6 -4.5 4.330 5 0 1
52 C8H15NO6 -4.4 6.128 1 0 1
53 C8H15NO6 -4.3 4.000 1 0 0

Лиганд с наивысшей аффинностью — C13H17NO6. Он также образует много водородных связей.