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

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

# organic
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

import pmx # pdb manipulation 
In [5]:
# protein preparation

pdb=pmx.Model('alignment_with_ligand.B99990001.pdb')
for r in pdb.residues[135:]:
    print r #show residues
<Molecule: id = 136 name = ASN chain_id =    natoms = 8>
<Molecule: id = 137 name = ARG chain_id =    natoms = 11>
<Molecule: id = 138 name = ASP chain_id =    natoms = 8>
<Molecule: id = 139 name = VAL chain_id =    natoms = 7>
<Molecule: id = 140 name = ARG chain_id =    natoms = 11>
<Molecule: id = 141 name = GLN chain_id =    natoms = 9>
<Molecule: id = 142 name = TYR chain_id =    natoms = 12>
<Molecule: id = 143 name = VAL chain_id =    natoms = 7>
<Molecule: id = 144 name = GLN chain_id =    natoms = 9>
<Molecule: id = 145 name = GLY chain_id =    natoms = 4>
<Molecule: id = 146 name = CYS chain_id =    natoms = 6>
<Molecule: id = 147 name = GLY chain_id =    natoms = 4>
<Molecule: id = 148 name = VAL chain_id =    natoms = 8>
<Molecule: id = 149 name = NAG chain_id =    natoms = 14>
<Molecule: id = 150 name = NAG chain_id =    natoms = 14>
<Molecule: id = 151 name = NDG chain_id =    natoms = 15>
In [6]:
# create protein and ligand objects
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
newpdb.writePDB(fname="protein_without_ligand.pdb")

ligand = pdb.copy()
for r in ligand.residues[:-3]:
    ligand.remove_residue(r)
In [7]:
# find geometric ligand center
center = np.mean([a.x for a in ligand.atoms],axis=0)
In [8]:
# prepare protein for docking

prot = oddt.toolkit.readfile('pdb','protein_without_ligand.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
Out[8]:
True

NAG содержит в себе CH3C(=O)NH группу

Создадим лиганды где метильный радикал этой группы будет заменён на:

OH

NH3+

H

Ph

COO-

In [9]:
# NAG SMILES notation: CC(=O)NC1C(C(C(OC1O)CO)O)O
# list ligands for docking
smiles = ['OC(=O)NC1C(C(C(OC1O)CO)O)O',
         '[NH3+]C(=O)NC1C(C(C(OC1O)CO)O)O',
         'C(=O)NC1C(C(C(OC1O)CO)O)O',
         'C1=CC=C(C=C1)C(=O)NC1C(C(C(OC1O)CO)O)O',
         '[O-]C(=O)C(=O)NC1C(C(C(OC1O)CO)O)O']
ligands = []
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()
        ligands.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 H O O O O O N O *****
***** - Open Babel Depiction H H H H H H H H O O O N 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 O O O O O N O *****
***** - Open Babel Depiction H H H H H O O O O O O N O HO *****
In [10]:
#create docking object

dock_obj= oddt.docking.AutodockVina.autodock_vina(protein=prot,size=(20,20,20),center=center, 
                                                  executable='/usr/bin/vina',autocleanup=True, num_modes=9)
print dock_obj.tmp_dir

# vina parameters:

# --center - coordinates of the center
# --size - size by dimensions
# --cpu - the number of CPUs to use
# --exhaustiveness - exhaustiveness of the global search (roughly proportional to time)
# --num_modes - maximum number of binding modes to generate
# energy_range - maximum energy difference between the best binding mode and the worst one displayed (kcal/mol)

print " ".join(dock_obj.params) 
/tmp/autodock_vina_D7ngZp
--center_x 54.2953488372 --center_y 38.7627906977 --center_z 28.0037209302 --size_x 20 --size_y 20 --size_z 20 --cpu 1 --exhaustiveness 8 --num_modes 9 --energy_range 3
In [11]:
# do docking

res = dock_obj.dock(ligands,prot)
In [12]:
import pandas as pd


# show results: affinity for every ligand

formula = []
affinity = []

for i,r in enumerate(res):
    formula.append(r.formula)
    affinity.append(r.data['vina_affinity'])
    
res_table = pd.DataFrame()
res_table['ligand'] = formula
res_table['affinity'] = affinity

# sort by affinity
res_table.sort_values('affinity', inplace=True)

res_table
Out[12]:
ligand affinity
26 C7H13NO6 -4.4
25 C7H13NO6 -4.5
24 C7H13NO6 -4.5
23 C7H13NO6 -4.7
17 C7H15N2O6 -4.8
22 C7H13NO6 -4.9
43 C8H13NO8 -4.9
16 C7H15N2O6 -4.9
44 C8H13NO8 -4.9
8 C7H13NO7 -4.9
42 C8H13NO8 -5.0
5 C7H13NO7 -5.0
21 C7H13NO6 -5.0
41 C8H13NO8 -5.0
14 C7H15N2O6 -5.0
15 C7H15N2O6 -5.0
7 C7H13NO7 -5.0
6 C7H13NO7 -5.0
40 C8H13NO8 -5.1
39 C8H13NO8 -5.2
4 C7H13NO7 -5.2
13 C7H15N2O6 -5.2
20 C7H13NO6 -5.2
19 C7H13NO6 -5.3
12 C7H15N2O6 -5.3
3 C7H13NO7 -5.3
2 C7H13NO7 -5.4
11 C7H15N2O6 -5.5
10 C7H15N2O6 -5.5
1 C7H13NO7 -5.6
38 C8H13NO8 -5.7
35 C13H17NO6 -5.7
9 C7H15N2O6 -5.7
18 C7H13NO6 -5.7
0 C7H13NO7 -5.7
33 C13H17NO6 -5.8
34 C13H17NO6 -5.8
37 C8H13NO8 -5.8
31 C13H17NO6 -5.9
30 C13H17NO6 -5.9
29 C13H17NO6 -5.9
32 C13H17NO6 -5.9
36 C8H13NO8 -6.1
28 C13H17NO6 -6.2
27 C13H17NO6 -6.9
In [14]:
# choose best ligand by its affinity to protein: C7H13NO6

print res_table.iloc[0].tolist()
['C7H13NO6', '-4.4']
In [19]:
# save pdbs for protein with best docked ligand
for i,r in enumerate(res):
    if r.formula == 'C7H13NO6':
        r.write(filename='docking-%s-%s.pdb' % (r.formula, i), format='pdb', overwrite=True)
In [1]:
# Launch PyMOL 
import __main__
__main__.pymol_argv = [ 'pymol', 'x' ]

import pymol
from pymol import cmd
pymol.finish_launching()

# Set image
import time
from IPython.display import Image
defaultImage = 'pymolimg.png'
def prepareImage(width=500, height=500, sleep=2, filename=defaultImage):
    cmd.ray(width, height)
    cmd.png(filename)
    time.sleep(sleep)
In [2]:
cmd.load('docking-C7H13NO6-26.pdb', 'ligand')
cmd.load('protein_without_ligand.pdb', 'protein')
In [3]:
cmd.show_as('sticks', 'ligand')
cmd.show_as('cartoon', 'protein')
cmd.color('white','protein')

cmd.center('ligand')
cmd.origin('ligand')
cmd.zoom('ligand', 4)
In [4]:
prepareImage()
Image(defaultImage)
Out[4]:
In [5]:
%%bash
jupyter nbconvert --to html hw7_filippova.ipynb