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 tabulate import tabulate
# Органика
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
import pmx # Модуль для манипулирования pdb
Подготовим белок:
pdb=pmx.Model('P08905_with_ligand.B99990001.pdb')
for r in pdb.residues[135:]:
print r #посмотрим остатки
# создание объектов белок и лиганда
X=3
newpdb = pdb.copy()
for r in newpdb.residues[-X:]:
newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-X]
Все успешно разделилось:
for i in lig.residues:
print i
Найдем геометрический центр лиганда:
s=np.zeros(3)
for a in lig.atoms:
s=s+a.x
centr=s/len(lig.atoms)
print centr # найдите геометрический центр лиганда
#пишем белок без лиганда
newpdb.writePDB('newpdb.pdb')
Подготовка белка для докинга:
prot = oddt.toolkit.readfile('pdb','newpdb.pdb').next()
prot.OBMol.AddPolarHydrogens()
prot.OBMol.AutomaticPartialCharge()
mols=[]
print 'is it the first mol in 1lmp is protein?',prot.protein,':) and MW of this mol is:', prot.molwt
Функция отрисовки смайлов:
def drawSmile(smiles,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)
smiles = ['c1cccc(O)c1', 'c1c(O)ccc(O)c1','c1(O)cc(c2ccccc2)cc(O)c1']
drawSmile(smiles,mols)
Функция для докинга:
def docking(centr,prot,mols):
#create docking object
xx=centr[0]
yy=centr[1]
zz=centr[2]
#установлен лимит в 1-9 моделей
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=[xx,yy,zz],
executable='/usr/bin/vina',autocleanup=True, num_modes=9)
print dock_obj.tmp_dir
print " ".join(dock_obj.params) # Опишите выдачу
# do it
res = dock_obj.dock(mols,prot)
#анализ докинга
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)
return res
Функция для результатов докинга:
Визуализация
def results(res,n):
param=len(res)/n
d_affinity = [[] for x in xrange(n)]
d_rmsd=[[] for x in xrange(n)]
names=[]
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)
d_affinity[i / param].append(r.data['vina_affinity'])
d_rmsd[i / param].append(r.data['vina_rmsd_ub'])
names.append(r.formula)
setlist=set(names)
# на выходе - данные по афинности, rmsd и названия молекул
return d_affinity,d_rmsd,list(setlist)
def visualization(res,name):
for i,r in enumerate(res):
r.write(filename=name % i, format='pdb',overwrite=True)
#pymol .pdb r*pdb
Функция для создания красивого вывода таблиц:
def tablichka(affin,rmsd,names):
row_format ="{:>15}" * (len(names) + 1)
rows=list(range(0,len(affin)*len(affin[0])))
print ("AFFINITY")
print tabulate(map(list, zip(*affin)),headers=names)
print ("===================================")
print ("RMSD")
print tabulate(map(list, zip(*rmsd)),headers=names)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,3)
visualization(res,'r%s.pdb')
tablichka(affin,rmsd,names)
Чем меньше афинность, тем лучше. Поэтому видно, что для C12H10O2 1 вариант в таблице со значением афинности -6.1 и rmsd 0 будет наилучшим. Изобразим белок и лиганд:
Image(filename="Screenshot_1.png")
Создадим лиганды где метильный радикал СH3C(=O)NH группы будет заменен на следующие группы:
smile=(['C1(C(C(C(C(O1)CO[H])O[H])O[H])O[H])[OH]'])
mols=[]
drawSmile(smile,mols)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,1)
tablichka(affin,rmsd,names)
Первый вариант - лучший.
visualization(res,'r_oh%s.pdb')
Пристраиваем белок:
Image("oh.png")
smile=(['C1(C(C(C(C(O1)CO[H])O[H])O[H])O[H])[NH3+]'])
mols=[]
drawSmile(smile,mols)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,1)
tablichka(affin,rmsd,names)
visualization(res,"r_nh3%s.pdb")
Image('nh3.png')
smile=(['C1(C(C(C(C(O1)CO[H])O[H])O[H])O[H])[H]'])
mols=[]
drawSmile(smile,mols)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,1)
tablichka(affin,rmsd,names)
visualization(res,"r_h%s.pdb")
Image('h.png')
smile=(['C1(C(C(C(C(O1)O[H])O[H])O[H])C2CCCCC2)O[H]'])
mols=[]
drawSmile(smile,mols)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,1)
tablichka(affin,rmsd,names)
visualization(res,"r_ph%s.pdb")
Image("ph.png")
smile=(['C1(C(C(C(C(O1)O[H])O[H])O[H])C(=O)O)O[H]'])
mols=[]
drawSmile(smile,mols)
res=docking(centr,prot,mols)
affin,rmsd,names=results(res,1)
tablichka(affin,rmsd,names)
visualization(res,"r_coo%s.pdb")
Image("coo.png")
Из статьи , например , можно было бы использовать autodock vina для докинга.