In [9]:
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 import Chem
from rdkit.Chem.Draw import IPythonConsole
In [12]:
# Надо убрать старый лигнад из белка
import mdtraj as md
u = md.load('with_lig.B99990001.pdb')
pdb = u.topology
for i,r in enumerate(pdb.atoms):
print(i,r)
#посмотрим на атомы
In [15]:
# сохраним pdb без лиганда
u2 = u.atom_slice(range(982))
u2.save_pdb('glam.pdb', force_overwrite=True)
In [16]:
lig = u.atom_slice(range(982,997))
gcenter = [lig.xyz[:,:,0].mean()*10,lig.xyz[:,:,1].mean()*10,lig.xyz[:,:,2].mean()*10]
gcenter
Out[16]:
Подготовка белка для докинга
In [17]:
prot = next(oddt.toolkit.readfile('pdb','glam.pdb'))
print('is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt )
In [20]:
smiles = ['O=C(NC1C(O)C(O)C(OC1O)CO)C', 'O=C(NC1C(O)C(O)C(OC1O)CO)O','O=C(NC1C(O)C(O)C(OC1O)CO)[NH3+]','O=C(NC1C(O)C(O)C(OC1O)CO)','O=C(NC1C(O)C(O)C(OC1O)CO)C1=CC=CC=C1','O=C(NC1C(O)C(O)C(OC1O)CO)C([O-])=O']
mols= []
images =[]
for s in smiles:
m = oddt.toolkit.readstring('smi', s)
if not m.Mol.HasProp('3D'):
m.make3D(forcefield='mmff94', steps=150)
mols.append(m)
images.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*images)
Докинг
In [31]:
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=gcenter,
executable='bin/vina',autocleanup=True, num_modes=9)
print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params))
Сначала - координаты геометрического центра, затем - объем, в котором могут перемещаться атомы лигандов при докинге, время поиска минимума энергии, количество моделей, разница в энергиях лучшей и худшей модели.
In [33]:
res = dock_obj.dock(mols,prot)
Результаты докинга
In [35]:
for i,r in enumerate(res):
### Узнаем дескриптор "формула"
formula = Chem.rdMolDescriptors.CalcMolFormula(r.Mol)
print ("%4d%10s%8s%8s%8s" %(i, formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name))
In [36]:
import pandas as pd
In [40]:
results_df = pd.DataFrame([[Chem.rdMolDescriptors.CalcMolFormula(r.Mol), r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name] for r in res],
columns=['formula', 'affinity', 'rmsd', 'name'])
results_df.sort_values(by=['affinity'], ascending=False)
Out[40]:
In [44]:
analyze = [0, 8, 16, 22, 30, 34, 46,52]
print("i\tformula\thbonds_num\thbonds_strict\tstack_num\tstack_par\tstack_per\thydroph_num")
for i,r in enumerate(res):
if i not in analyze:
continue
formula = Chem.rdMolDescriptors.CalcMolFormula(r.Mol)
mol1, mol2, strict = oddt.interactions.hbonds(prot,r)
hbonds_num = len(mol1)
hbonds_strict = strict.sum()
r1, r2, strict_par, strict_per = oddt.interactions.pi_stacking(prot,r)
stack_num = len(r1)
stack_par = strict_par.sum()
stack_per = strict_per.sum()
mol1, mol2 = oddt.interactions.hydrophobic_contacts(prot,r)
hydroph_num = len(mol1)
print(i,"\t",formula,"\t",hbonds_num,"\t","\t",hbonds_strict,"\t","\t",stack_num,"\t","\t",stack_par,"\t","\t",stack_per,"\t","\t",hydroph_num)
Визуализация
In [45]:
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb')
In [46]:
from IPython.display import Image
In [48]:
display(Image('r0.png'))
In [50]:
display(Image('r10.png'))
In [ ]: