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.Draw import IPythonConsole
import pmx # Модуль для манипулирования pdb
In [3]:
best_model_pdb_path = "red2.BL00010001.pdb"
In [52]:
pdb=pmx.Model(best_model_pdb_path)
for r in pdb.residues[135:]:
print r #посмотрим остатки
In [54]:
newpdb = pdb.copy()
for i in range(151, len(newpdb.residues)):
newpdb.remove_residue(newpdb.residues[151])
In [56]:
lig = pdb.copy()
for i in range(0, 151):
lig.remove_residue(lig.residues[0])
In [58]:
mean_center = np.mean(np.array([a.x for a in lig.atoms]), axis = 0)
median_center = np.median(np.array([a.x for a in lig.atoms]), axis = 0)
In [59]:
print mean_center
print median_center
Будем использовать mean center
In [60]:
out_pdb_file_name = "lys_whout_ligand.pdb"
newpdb.writePDB(out_pdb_file_name, title="Modelled Lysozyme with deleted ligand")
In [61]:
prot = oddt.toolkit.readfile('pdb',out_pdb_file_name).next()
In [62]:
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
Out[62]:
In [63]:
print 'is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt
In [66]:
prot.protein = True
In [67]:
smiles = ['c1cccc(O)c1', 'c1c(O)ccc(O)c1','c1(O)cc(c2ccccc2)cc(O)c1']
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 [68]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=mean_center,
executable='/usr/bin/vina',autocleanup=True, num_modes=20)
print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
In [74]:
res = dock_obj.dock(mols, prot)
In [75]:
for i,r in enumerate(res):
print "%4d%10s%8s%8s%8s" %(i,r.formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name)
In [84]:
?oddt.interactions.hydrophobic_contacts
In [177]:
?oddt.interactions.hbonds
In [179]:
for i,r in enumerate(res):
mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
mol1, mol2, strict_parallel, strict_perpendicular= oddt.interactions.pi_stacking(prot,r)
mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
Создадим лиганды где метильный радикал группы СH3C(=O)NH будет заменён на :
1.OH
2.NH3+
3.H
4.Ph
5.COO-
In [161]:
nag_mod_smiles = ["CC(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O", \
"C(O)(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O",
"C([NH3+])(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O",
"C([H])(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O",
"C(C(=O)[O-])(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O",
"C(c1ccccc1)(=O)N[C@@H]1[C@H]([C@@H]([C@H](O[C@H]1O)CO)O)O"]
names = ["NAG", "OH_mod", "NH3_mod", "H_mod", "COO_mod", "Ph_mod"]
In [235]:
nag_mols= []
nag_images =[]
for s, name in zip(nag_mod_smiles, names):
m = oddt.toolkit.readstring('smi', s)
if not m.OBMol.Has3D():
m.make3D(forcefield='mmff94', steps=150)
m.removeh()
m.OBMol.AddPolarHydrogens()
m.name = name
nag_mols.append(m)
###with print m.OBMol.Has3D() was found that:
### deep copy needed to keep 3D , write svg make mols flat
nag_images.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*nag_images)
In [249]:
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=mean_center,
executable='/usr/bin/vina',autocleanup=False, num_modes=20)
print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
In [250]:
nag_res = dock_obj.dock(nag_mols, prot)
In [255]:
nag_res_sorted = sorted(nag_res, key = lambda x : float(x.data['vina_affinity']))
In [256]:
print "Place\t\tFormula\tAffinity\tRMSD\tName\tHb_tot\tHb_strict\tPi_tot\tPi_par\tPi_str_perp\tHydphob"
for i,r in enumerate(nag_res_sorted):
mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
hbonds_num_total = len(mol1)
hbonds_num_strict = strict.sum()
mol1, mol2, strict_parallel, strict_perpendicular= oddt.interactions.pi_stacking(prot,r)
pi_stack_num_total = len(mol1)
pi_stack_num_strict_par = strict_parallel.sum()
pi_stack_num_strict_perp = strict_perpendicular.sum()
mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
hydrophob_num = len(mol1)
print "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}".format(i,r.formula, r.data['vina_affinity'],\
r.data['vina_rmsd_ub'], r.residues[0].name,
hbonds_num_total, hbonds_num_strict,\
pi_stack_num_total, pi_stack_num_strict_par, pi_stack_num_strict_perp,\
hydrophob_num)
Кажется странным целых 17 водородных связей. Посмотрим на их характер
In [257]:
mol1, mol2, strict = oddt.interactions.hbonds(prot,nag_res_sorted[1])
In [258]:
for x,y in zip(mol1, mol2):
print "{}:{}-->{}:{}".format(x[9], x[5],y[9] or "LIG", y[5])
Видим, что часть из них "дублируется", т.е у них одинаковый акцептор/донор
In [270]:
for ind, res in enumerate(nag_res_sorted):
res.write(filename='r{}.pdb'.format(ind), format='pdb')
Посмотрим на этот лиганд в структуре белка
In [323]:
from IPython.display import Image
Image(filename='1.png')
Out[323]:
В принципе, расположен он удачно
Отсортируем более по-умному За каждую strict водородную связб будем давать 7.5 очков, strict pi-связь - 2, гидрофобную - 5, за не strict - половину соответствующей энергии.
In [289]:
def score_func(prot, r):
mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
hbonds_num_total = len(mol1)
hbonds_num_strict = strict.sum()
mol1, mol2, strict_parallel, strict_perpendicular= oddt.interactions.pi_stacking(prot,r)
pi_stack_num_total = len(mol1)
pi_stack_num_strict_par = strict_parallel.sum()
pi_stack_num_strict_perp = strict_perpendicular.sum()
mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
hydrophob_num = len(mol1)
score = hbonds_num_strict * 7.5 + (hbonds_num_total - hbonds_num_strict) * 0.5 +\
pi_stack_num_strict_par * 2 + pi_stack_num_strict_perp * 2 +\
(pi_stack_num_total - pi_stack_num_strict_par - pi_stack_num_strict_perp) * 1+\
hydrophob_num * 5
return score
In [292]:
nag_res_cl_sorted = sorted(nag_res, key = lambda x : -score_func(prot, x))
In [322]:
print "Place Formula Score Hb_tot Hb_str Pi_tot Pi_par Pi_perp Hydphob"
for i,r in enumerate(nag_res_cl_sorted):
mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
hbonds_num_total = len(mol1)
hbonds_num_strict = strict.sum()
mol1, mol2, strict_parallel, strict_perpendicular= oddt.interactions.pi_stacking(prot,r)
pi_stack_num_total = len(mol1)
pi_stack_num_strict_par = strict_parallel.sum()
pi_stack_num_strict_perp = strict_perpendicular.sum()
mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
hydrophob_num = len(mol1)
print "{:5} {:10} {:5} {:5}{:5}{:7}{:7}{:7}{:7}".format(i,r.formula, \
score_func(prot, r),
hbonds_num_total, hbonds_num_strict,\
pi_stack_num_total, pi_stack_num_strict_par, pi_stack_num_strict_perp,\
hydrophob_num)
Видим, что ситуация не особо изменилась в данном случае (для топ-5, к примеру)