import numpy as np
import pandas as pd
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
1. Подготовка белка
pdb=pmx.Model('ligand.B99990001.pdb')
for r in pdb.residues[135:]:
print r
# создание объектов белок и лиганда
newpdb = pdb.copy()
for r in newpdb.residues[-3:]:
newpdb.remove_residue(r)
lig = pdb.copy()
del lig.residues[:-3]
#будем считать за центр среднее по всем координатам
coordinates = []
for a in lig.atoms:
coordinates.append(a.x)
center = np.mean(coordinates,axis=0)
print(center)
newpdb.writePDB('new_prot.pdb')
2. Подготовка белка для докинга
prot = oddt.toolkit.readfile('pdb','new_prot.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
3. Лиганды для докинга
NAG содержит в себе СH3C(=O)NH группу.
Заменяем метильный (CH3) радикал этой группы на :
1) OH -> O
2) NH3+ -> [N+]
3) H ->
4) Ph -> C1=CC=[C]C=C1
5) COO- -> [-O]C(=O)
SMILES для NAG 'CC(=O)NC1C(C(C(OC1O)CO)O)O': http://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/NAG
Значит, с точки зрения записи меняться будет самый левый углерод.
#http://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/NAG
smilesNAG = 'CC(=O)NC1C(C(C(OC1O)CO)O)O'
smiles = ['OC(=O)NC1C(C(C(OC1O)CO)O)O', #OH
'[NH3+]C(=O)NC1C(C(C(OC1O)CO)O)O', #NH3+
'C(=O)NC1C(C(C(OC1O)CO)O)O', #H
'C1=CC=[C]C=C1C(=O)NC1C(C(C(OC1O)CO)O)O', #PH
'[O-]C(=O)C(=O)NC1C(C(C(OC1O)CO)O)O'] #COO
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)
images.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*images)
4. Докинг
xx = center[0]
yy = center[1]
zz = center[2]
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('Center coordinates')
print " ".join(dock_obj.params[0:6])
print ("№ of points")
print " ".join(dock_obj.params[6:12])
print ("Number Of CPU Cores")
print " ".join(dock_obj.params[12:14])
print("Exhaustiveness parameter")
print " ".join(dock_obj.params[14:16])
print("Number of conformations generated by Autodock Vina")
print " ".join(dock_obj.params[16:18])
print("Energy range cutoff for Autodock Vina")
print " ".join(dock_obj.params[18:20])
res = dock_obj.dock(mols,prot)
5. Результаты докинга
formula = []
affinity = []
rsmd = []
name = []
for i,r in enumerate(res):
formula.append(r.formula)
affinity.append(r.data['vina_affinity'])
rsmd.append(r.data['vina_rmsd_ub'])
name.append(r.residues[0].name)
radical = ['OH']*9+['NH3+']*9+['H']*9+['Ph']*9+['COO']*9
#data = np.array()
table = pd.DataFrame([radical,formula,affinity,rsmd,name])
table = table.transpose()
table = table.rename(index=str, columns={0: 'Radical', 1: 'Formula',2:'Affinity',3:'RSMD',4:'Name'})
table.sort_values(by='Affinity')
6. Анализ докинга и сохранение результатов
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)
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb')
7. Визуализация
Критерий выбора для визуализации - наименьшее значение афинности по каждому радикалу.
import __main__
__main__.pymol_argv = ['pymol','-qc'] # Pymol: quiet and no GUI
from time import sleep
import time
import pymol
pymol.finish_launching()
defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=600, height=600, sleep=2, filename=defaultImage):
## To save the rendered image
pymol.cmd.ray(width, height)
pymol.cmd.png('/tmp/pymolimg.png')
time.sleep(sleep)
#OH
pymol.cmd.reinitialize()
pymol.cmd.load('new_prot.pdb')
pymol.cmd.load('r8.pdb')
pymol.cmd.show('sticks','r8')
pymol.cmd.show('cartoon','new_prot')
pymol.cmd.hide('lines','new_prot')
pymol.cmd.orient()
prepareImage()
Image(defaultImage)
#H
pymol.cmd.reinitialize()
pymol.cmd.load('new_prot.pdb')
pymol.cmd.load('r26.pdb')
pymol.cmd.show('sticks','r26')
pymol.cmd.show('cartoon','new_prot')
pymol.cmd.hide('lines','new_prot')
pymol.cmd.orient()
prepareImage()
Image(defaultImage)
#NH3
pymol.cmd.reinitialize()
pymol.cmd.load('new_prot.pdb')
pymol.cmd.load('r17.pdb')
pymol.cmd.show('sticks','r17')
pymol.cmd.show('cartoon','new_prot')
pymol.cmd.hide('lines','new_prot')
pymol.cmd.orient()
prepareImage()
Image(defaultImage)
#COO
pymol.cmd.reinitialize()
pymol.cmd.load('new_prot.pdb')
pymol.cmd.load('r43.pdb')
pymol.cmd.show('sticks','r43')
pymol.cmd.show('cartoon','new_prot')
pymol.cmd.hide('lines','new_prot')
pymol.cmd.orient()
prepareImage()
Image(defaultImage)
#Ph
pymol.cmd.reinitialize()
pymol.cmd.load('new_prot.pdb')
pymol.cmd.load('r35.pdb')
pymol.cmd.show('sticks','r35')
pymol.cmd.show('cartoon','new_prot')
pymol.cmd.hide('lines','new_prot')
pymol.cmd.orient()
prepareImage()
Image(defaultImage)
ODDT
Судя по всему, мы могли бы проделать весь Example 1: filtering, docking and re-scoring workflow, чтобы достичь нашей цели
1) score и rescore докинга с помощью RF-Score and NNScore
2) предсказывать афинность разными методами: случайным лесом, множественной линейной регрессией, PLS регрессией и нейронной сетью