In this task, we'll perform docking of ligands of some apparently random ligands (namely, hydroxybenzene, 1,4-Dihydroxybenzene and 3,5-Dihydroxybiphenyl), to the previously generated protein structure of 1lmpP00705.fasta.
In [1]:
from __future__ import division
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.Draw import IPythonConsole
import pmx # Модуль для манипулирования pdb
In [2]:
pdb=pmx.Model('../Task6/ligand_enabled.B99990001.pdb')
for r in pdb.residues[-15:]:
print r #посмотрим остатки
In [3]:
# there is a Model.fetch_residues() function, but it is woefully
# undocumented, and using magic numbers to split the protein is
# distasteful, so there.
protein_atoms = filter(lambda a: a.molecule.is_protein_residue(), pdb.atoms)
ligand_atoms = filter(lambda a: not a.molecule.is_protein_residue(), pdb.atoms)
newpdb = pmx.Model(atoms = protein_atoms)
lig = pmx.Model(atoms = ligand_atoms)
print(len(newpdb.residues), len(lig.residues))
Geometric center of ligand:
In [4]:
lig_center = sum((np.array(a.x) for a in lig.atoms))/len(lig.atoms)
lig_center
Out[4]:
In [5]:
newpdb.writePDB(fname="ligandless.pdb")
In [6]:
prot = oddt.toolkit.readfile('pdb','ligandless.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
# The flying duck?
# print('is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt)
print('MW of the mol: %f' % prot.molwt)
We want to take NAG and replace the first C (that is, CH3) in its' CC(=O)N- group with various other groups, such as:
- OH
- NH3+
- H
- Ph/C6H5
- COO-
(One possible) Smiles of NAG is CC(=O)NC1C(O)C(O)C(OC1O)CO, from which we construct smiles of alternate ligands:
In [7]:
smiles_common = "C(=O)NC1C(O)C(O)C(OC1O)CO"
smiles_variant = ["O", "[N+]", "", "c1ccccc1", "[O-]"]
smiles = [s+smiles_common for s in smiles_variant]
#smiles = ['OC(=O)NC1C(O)C(O)C(OC1O)CO',
# '[N+]C(=O)NC1C(O)C(O)C(OC1O)CO',
# 'C(=O)NC1C(O)C(O)C(OC1O)CO',
# 'c1ccccc1C(=O)NC1C(O)C(O)C(OC1O)CO',
# '[O-]C(=O)C(=O)NC1C(O)C(O)C(OC1O)CO']
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)
In [8]:
#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) # Опишите выдачу
In which autodock_vina parameters are described
center_x,y,zandsize_x,y,zdescribe dimensions and position of the docking box — fragment of space the ligand can be moved over in the process of docking;cpuquite probably limits cpu usage to one core;exhaustivenessparameter selects a point on the [computation time — probability of finding global minimum] tradeoff curve;num_modesdetermines the maximum number of conformations generated;energy_rangedefines maximum energy difference between best and worst displayed conformation, effectively filtering conformations with energy above a threshold.
In [9]:
res = dock_obj.dock(mols,prot)
In [25]:
res_sorted = sorted(res, key=lambda x: x.data['vina_affinity'], reverse=True)
#table = []
formulas = set(r.formula for r in res)
#for f in formulas:
# table.append(next(r for r in res_sorted if r.formula == f))
#table.sort(key=lambda x: x.data['vina_affinity'], reverse=True)
print("%4s%10s%9s%8s%8s") % ("#", "formula", "vina_aff", "rmsd", "name")
for i,r in enumerate(res_sorted):
hbs = oddt.interactions.hbonds(prot,r)
stack= oddt.interactions.pi_stacking(prot,r)
phob = oddt.interactions.hydrophobic_contacts(prot,r)
print "%4d%10s%9s%8s%8s" \
%(i,r.formula, r.data['vina_affinity'],
r.data['vina_rmsd_ub'], r.residues[0].name,
)
In [26]:
for i,r in enumerate(res_sorted):
hbs = oddt.interactions.hbonds(prot,r)
stack= oddt.interactions.pi_stacking(prot,r)
phob = oddt.interactions.hydrophobic_contacts(prot,r)
#print(hbs, stack, phob)
In [12]:
for i,r in enumerate(table):
r.write(filename='r%s.pdb' % i, format='pdb', overwrite=True)
Let's look at the best available ligand.
In [27]:
import nglview
import ipywidgets
from IPython.display import Image
In [59]:
view1 = nglview.show_structure_file("ligandless.pdb")
view1.add_structure(nglview.FileStructure('r0.pdb'))
#view1.clear_representations()
#view1[0].add_representation(repr_type="")
view1[1].add_representation(repr_type="licorice")
view1[1].center_view()
view1
In [61]:
view1.download_image("religanded.png")
Image("religanded.png")
Out[61]:
