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 import Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
import pmx # Модуль для манипулирования pdb
pdb=pmx.Model('seq_lig.B99990001.pdb')
for r in pdb.residues[-10:]:
print r # посмотрим остатки
Создадим объекты белок и лиганд
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-3]
Вычислим координаты центра лиганда
atom_coords = np.zeros((len(lig.atoms), 3))
for x, a in enumerate(lig.atoms):
atom_coords[x, :] = a.x
lig_center = np.mean(atom_coords, axis=0)
lig_center
newpdb.writePDB('myprot.pdb')
prot = oddt.toolkit.readfile('pdb','myprot.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
print 'is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt
Забавно, что оддт не считает наш белок белком :)
Создадим несколько производных NAG (замещаем один из метилов)
smiles = ['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)(c1ccccc1)',
'NC1C(C(C(OC1O)CO)O)O', 'NC1C(C(C(OC1O)CO)O)OC', '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.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)
#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)
Параметры Вины (по порядку): Первые 6 параметров задают бокс, внутри которого могут передвигаться атомы (в нашем случае только лиганда); exhaustiveness задает насколько долго мы ищем минимум (который ищется эвристически); num_modes задает сколько моделей должно быть в выдаче; energy_range задает масимальную разницу энергий (в ккал/моль) между лучшей и худшей моделью для выдачи.
# do it
res = dock_obj.dock(mols,prot)
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)
Мы посмотрим на самого лучшего для каждого из соединений
selected = [0, 9, 18, 27, 34, 43, 52]
for i,r in enumerate(res):
if i not in selected:
continue
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, hbonds_num, hbonds_strict, stack_num, stack_par, stack_per, hydroph_num
Видно, что гидрофобные контакты есть только у 27; больше всего водородных связей у 52, однако хороших водородных связей у него столько же, сколько и у 0, а у 9 их почти чтолько же. Стэкинг никто не образует.
Мы посмотрим на два самых лучших:
0 C7H13NO7 -6.4 0.000 LIG 27 C13H17NO6 -7.3 0.000 LIG
selected = [0, 27]
for i,r in enumerate(res):
if i not in selected:
continue
r.write(filename='r%s.pdb' % i, format='pdb')
display(Image('r0.png'))
display(Image('r27.png'))
Как видно из рисунков, самым афинным лигандом оказался NAG с фенилаланином, а вторым -- NAG с гидроксилом. Они оба находятся примерно в одной позе. Пимоль почему-то находит гораздо меньше водородных связей, чем оддт (там у 0 было аж 13). Ароматическое кольцо фенилаланина так же не было изуродовано. В общем и целом нормально все задочилось.