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

import pmx
In [2]:
from os import system
In [3]:
def chemdraw(smiles):
    mol = MolFromSmiles(smiles)
    AllChem.Compute2DCoords(mol)
    display(mol)

Для работы возьмем модель с прошлого практикума

In [4]:
pdb=pmx.Model('../pr10/QUERY.B99990005.pdb')

Локализуем лиганд

In [5]:
for r in pdb.residues[-5:]:
    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>

Сделаем новые пдб - только лиганд и только модель

In [6]:
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
lig = pdb.copy()
for r in lig.residues[:-3]:
    lig.remove_residue(r)
In [7]:
lig.writePDB('ligand.pdb')
newpdb.writePDB('base.pdb')

Найдем геометрический центр лиганда - приближение положения сайта связывания

In [8]:
coords = np.zeros((len(lig.atoms),3))
for c,a in enumerate(lig.atoms):
    coords[c,:] = a.x
geom_center = np.mean(coords, axis = 0)
print geom_center

[ 41.87662791  45.44525581  29.17669767]

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

In [9]:
prot = oddt.toolkit.readfile('pdb','base.pdb').next()

prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
Out[9]:
True
In [10]:
prot.molwt
Out[10]:
15661.733120000134

Придумаем субстраты

In [11]:
NAG = "CC(=O)NC1C(C(C(OC1O)CO)O)O"
In [12]:
chemdraw(NAG)

Замещаем группу -NH-C(CH3)=O

In [13]:
NAG_OH = "OC1C(C(C(OC1O)CO)O)O"
NAG_N = "[NH3+]C1C(C(C(OC1O)CO)O)O"
NAG_H = "C1C(C(C(OC1O)CO)O)O"
NAG_Ph = "C1=CC=C(C=C1)C1C(C(C(OC1O)CO)O)O"

Посмотрим, что это за вещества, и подготовим их структуры к докингу

In [14]:
mols = []
for nag in (NAG_OH, NAG_N, NAG_H, NAG_Ph):
    chemdraw(nag)
    m = oddt.toolkit.readstring('smi', nag)
    if not m.OBMol.Has3D(): 
        m.make3D(forcefield='mmff94', steps=1000)
        m.removeh()
        m.OBMol.AddPolarHydrogens()

    mols.append(m)

Докинг

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

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

/tmp/autodock_vina_pEh42g
--center_x 41.876627907 --center_y 45.445255814 --center_z 29.1766976744 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 20 --energy_range 3
In [16]:
res = dock_obj.dock(mols,prot)
In [25]:
display(Image('all.png'))

Запишем пдб-файлы и сохраним в таблицу аффинность

In [18]:
results = np.zeros((len(res),3),dtype=object)
for i,r in enumerate(res):    
    r.write(filename='r%s.pdb' % i, format='pdb')
    results[i] = [i, r.formula, float(r.data['vina_affinity'])]
In [19]:
top3 = sorted(results, key=lambda x: x[2])[:3]
In [20]:
top3
Out[20]:
[array([27, 'C12H16O5', -6.3], dtype=object),
 array([28, 'C12H16O5', -6.2], dtype=object),
 array([29, 'C12H16O5', -6.2], dtype=object)]
In [26]:
display(Image('27.png'))
In [27]:
display(Image('28.png'))
In [28]:
display(Image('29.png'))

В топ-3 вошли позиции модификации с фенилом, что не удивительно. Подобная модификация мимикрирует еще одно моносахаридное кольцо по форме. Также 29 образует пи-пи стекинг

Посмотрим лучших по остальным трем веществам

In [23]:
res_nag_h = results[results[:,1] == 'C6H12O5',:]
res_nag_oh = results[results[:,1] == 'C6H12O6',:]
res_nag_nh3 = results[results[:,1] == 'C6H14NO5',:]
In [24]:
print res_nag_h[res_nag_h[:,2] == np.min(res_nag_h[:,2]),0][0]
print res_nag_oh[res_nag_oh[:,2] == np.min(res_nag_oh[:,2]),0][0]
print res_nag_nh3[res_nag_nh3[:,2] == np.min(res_nag_nh3[:,2]),0][0]

18
0
9
In [29]:
display(Image('18.png'))
In [30]:
display(Image('0.png'))
In [31]:
display(Image('9.png'))