Задание №7

Загрузим необходимые пакеты.

In [1]:
import sys
stdout = sys.stdout

import numpy as np
import copy
import base64

# Отображение структур
import IPython.display
import ipywidgets
from IPython.display import display, display_svg, SVG, Image, HTML

# 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

import time
import pymol
from pymol import cmd
pymol.finish_launching()

sys.stdout = stdout

/usr/local/lib/python2.7/dist-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)

 Adjusting settings to improve performance for Intel cards.

Реализуем вспомогательные функции.

In [2]:
def prepareImage(ind, width=500, height=500, sleep=5):
    cmd.ray(width, height)
    cmd.png('pymolimg%d.png' % ind)
    time.sleep(sleep)

def generate_html(title, ind):
    image_file = 'pymolimg%s.png' % ind
    html_code = '<h1>%s</h1><br><img src=\"data:image/jpeg;base64,%s\"/>'
    img_encoding = ''
    with open(image_file, 'rb') as f:
        img_encoding = f.read().encode('base64')
    
    return html_code % (title, img_encoding)

Посмотрим на файл с белком и удалим из него лиганд.

In [3]:
# Загрузим белок
pdb=pmx.Model('model_ligand.B99990001.pdb')

for r in pdb.residues[120:]:
    print r

<Molecule: id = 121 name = GLN chain_id =    natoms = 9>
<Molecule: id = 122 name = ALA chain_id =    natoms = 5>
<Molecule: id = 123 name = TRP chain_id =    natoms = 14>
<Molecule: id = 124 name = ILE chain_id =    natoms = 8>
<Molecule: id = 125 name = ARG chain_id =    natoms = 11>
<Molecule: id = 126 name = GLY chain_id =    natoms = 4>
<Molecule: id = 127 name = CYS chain_id =    natoms = 6>
<Molecule: id = 128 name = ARG chain_id =    natoms = 11>
<Molecule: id = 129 name = LEU chain_id =    natoms = 9>
<Molecule: id = 130 name = NAG chain_id =    natoms = 14>
<Molecule: id = 131 name = NAG chain_id =    natoms = 14>
<Molecule: id = 132 name = NDG chain_id =    natoms = 15>
In [4]:
# Создание объектов белка и лиганда
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
    newpdb.remove_residue(r)
newpdb.writePDB('model.pdb')

lig = pdb.copy()
del lig.residues[:-3]
ligand_center = sum(map(lambda atom: np.array(atom.x), lig.atoms)) / len(lig.atoms)

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

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

prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
Out[5]:
True

Выпишем смайлсы лигандов для каждого из приведенных в задании заместителей.

In [6]:
smiles = ['OCC1OC(O)C(NC(O)=O)C(O)C1O', '[NH3+]C(=O)NC1C(O)OC(CO)C(O)C1O', '[H]C(=O)NC1C(O)OC(CO)C(O)C1O', 
          'OCC1OC(O)C(NC(=O)C2=CC=CC=C2)C(O)C1O', 'OCC1OC(O)C(NC(=O)C([O-])=O)C(O)C1O']
labels = ['OH', 'NH3+', 'H', 'Ph', 'COO-']
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 O O N O O O O H H H H H H
OBDepict ***** N O N O O O O O H H H H H H H
OBDepict ***** O N O O O O O H H H H H
OBDepict ***** O O O N O O O H H H H H
OBDepict ***** O O O N O O HO O O H H H H H

Проведем докинг всех построенных лигандов в структуру белка.

In [7]:
dock_obj= oddt.docking.AutodockVina.autodock_vina(
    protein=prot,size=(20,20,20),center=ligand_center,
    executable='autodockvina/bin/vina',autocleanup=True, num_modes=9)
In [8]:
cmd.load('model.pdb')
cmd.do('select mod, model model')
cmd.show_as('cartoon', 'mod')
affinities = []
for ind, molecule in enumerate(mols):
    res = dock_obj.dock([molecule], prot)
    res[0].write(filename='result.pdb', format='pdb', overwrite=True)
    affinities.append(float(res[0].data['vina_affinity']))
    cmd.load('result.pdb')
    cmd.do('select res, model result')
    cmd.do('select neighbourhood, mod within 7 of res')
    cmd.do('zoom neighbourhood')
    cmd.do('show_as sticks, res')
    prepareImage(ind)
    display(HTML(generate_html(labels[ind], ind)))
    cmd.do('delete res')
    cmd.do('delete result')
    cmd.do('delete neighbourhood')
    time.sleep(5)

OH


NH3+


H


Ph


COO-


In [9]:
alts = sorted(zip(labels, affinities), key=lambda el: el[1])
print "%10s%10s" % ('Ligand', 'Affinity')
for lab, aff in alts:
    print "%10s%10s" % (lab, aff)

    Ligand  Affinity
        Ph      -6.6
      COO-      -6.0
      NH3+      -5.9
        OH      -5.4
         H      -5.3

Видно, что лучшим заместителем является фенильная группа.