import numpy as np
import copy
#Отображение структур
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
pdb=pmx.Model('with_legand.B99990001.pdb')#/home/shad/hse/kononkova/
for r in pdb.residues[135:]:
print r
# создание объектов белок и лиганда
newpdb = pdb.copy()
for r in newpdb.residues[137:]:
newpdb.remove_residue(r)
lig = pdb.copy()
for r in lig.residues[0:137]:
lig.remove_residue(r)
for r in lig.residues[0:]:
print(r)
Находим геометрический центр лиганда
tmp=[]
for a in lig.atoms:
tmp.append(a.x)
tmp=np.array(tmp)
cent=tmp.mean(axis=0)
newpdb.writePDB('lonely_protein.pdb')
prot = oddt.toolkit.readfile('pdb','lonely_protein.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
Программа считает, что это не белок. Ну и что. Заменим метильный радикал на другие (OH,NH3+,H,Ph, COO-) и проведем докинг.
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=[C]C=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.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=cent,
executable='/usr/bin/vina',autocleanup=True, num_modes=9)
print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
" ".join(dock_obj.params):
center - геометрический центр лиганда, size - размер docking box - пространство, которым будет ограничиваться докинг, cpu - количество процессоров, exhaustiveness - "полнота охвата" поиска (пропорциональна времени поиска) или количество попыток (чтобы найти наилучшую конформацию) в одном докинг-эксперименте, num_modes - количество конформаций c энергиями, energy_range - максимальная разница в "очках" для конформаций: если для какой-то конформации отклонение от наилучшего score окажется больше этой величины (по дефолту - 3 ккал/моль), то эта конформация исключается.
res = dock_obj.dock(mols,prot)
print 'description: №, formula, affinity, rmsd, name, hbonds, pi_stacking, hydrophobic_contacts'
for i,r in enumerate(res):
hbs = oddt.interactions.hbonds(prot,r)# водородные связи
stack= oddt.interactions.pi_stacking(prot,r)# стэкинг-взаимодействия
phob = oddt.interactions.hydrophobic_contacts(prot,r)# гидрофобные контакты
print "%4d%10s%8s%8s%8s%8s%8s%8s" %(i,r.formula, r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name,len(hbs[1]),len(stack[1]), len(phob[1]))
По результатам докинга лучшим лигандом оказался лиганд с заместителем - фенильной группой. У него самая лучшая афинность (самая малая энергия образования комплекса белок-лиганд) и больше всего гидрофобных контактов. Заметим, что замещение метильной группы на фенильную приближает лиганд к первоначальному лиганду из трех остатков, для которого считали центр.
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb')
Белок с "родным" лигандом
display(Image('native.png'))
Все лиганды (по одной конформации на лиганд (с разными заместителями))
display(Image('ligands.png'))
3 конформации для лучшего заместителя
display(Image('3best.png'))
Лиганд с лучшим заместителем
display(Image('best.png'))
Лиганд с разными заместителями метильной группы от лучшего к худшему, средняя афинность для 9 конформаций, ккал/моль. Афинность здесь пропорциональна dG. Победила фенильная группа.
from IPython.display import HTML, display
import tabulate
table = [["ligand", "mean affinity"],["C13H17NO6",-6.12],["C7H15N2O6",-5.66],["C7H13NO7",-5.53],["C7H13NO6",-5.42],["C8H13NO8",-5.35],["C8H15NO6",-5.26]]
display(HTML(tabulate.tabulate(table, tablefmt='html')))
По данным из статьи с помощью oddt можно протестировать разные scoring functions для наших молекул, можно попробовать предсказать афинность с помощью разных моделей (random forests, support vector machines, neural networks and multiple linear regression). Также в статье говорится о том, что oddt дает возможность