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 import AllChem
from rdkit.Chem import MolFromSmiles
from rdkit.Chem.Draw import IPythonConsole
import pandas as pd
from os import system
def chemdraw(smiles):
mol = MolFromSmiles(smiles)
AllChem.Compute2DCoords(mol)
display(mol)
import mdtraj as md
u = md.load('homolog.pdb')
pdb = u.topology
for i,r in enumerate(pdb.atoms):
print(i,r)
#посмотрим на атомы
# сохраним pdb без лиганда
u2 = u.atom_slice(range(1096))
u2.save_pdb('1lmp-clean.pdb', force_overwrite=True)
lig = u.atom_slice(range(1096,1124))
#lig= lig.topology
for v in lig.xyz:
lig_center = np.mean(v, axis=0)
lig_center
prot = next(oddt.toolkit.readfile('pdb','1lmp-clean.pdb'))
### В новой версии oddt всё сделает за вас
#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 )
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)',
'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)
###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)
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=[54.88893 , 43.46129 , 27.558858],
executable='/usr/bin/vina',autocleanup=True, num_modes=9)
print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params)) #
Первые 3 параметра - координаты центра, далее 3 параметра - размер ячейки, exhaustiveness - насколько долго мы ищем минимум ( ищется эвристически), num_modes - число моделей которое должно быть в выдаче, energy_range - масимальная разница энергий (в ккал/моль) между лучшей и худшей моделью для выдачи
res = dock_obj.dock(mols,prot)
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))
Посмотрим на резульаты докинга!
for i,r in enumerate(res):
r.write(filename='r%s.pdb' % i, format='pdb',overwrite=True)
best = np.zeros((len(res),3),dtype=object)
for i,r in enumerate(res):
formula = Chem.rdMolDescriptors.CalcMolFormula(r.Mol)
best[i] = [i, formula, float(r.data['vina_affinity'])]
top = sorted(best, key=lambda x: x[2])[:2] #посмотрим например на 28 позу
top
Лучшая поза докинга
top = sorted(best, key=lambda x: x[2])[:3]
top
Топ 3 лучших позы пришлись на лиганд с фенилом. Данный радикал среди остальных наиболее похож на моносахаридное кольцо, а также формирует пи-пи стекинг.
hb_total = []
hb_strict = []
phob_all = []
pi_par=[]
pi_perp=[]
pi=[]
for i,r in enumerate(res):
# mol2 - акцептор водородной связи, mol2 - донор, 'strict' связь (прошла все отсечки)
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)
hb_total.append(hbonds_num)
hb_strict.append(hbonds_strict)
phob_all.append(hydroph_num)
pi_par.append(stack_par)
pi_perp.append(stack_per)
pi.append(stack_num)
sort_table = pd.DataFrame({
'Formula':[Chem.rdMolDescriptors.CalcMolFormula(r.Mol) for r in res],
'Affinity (kcal/mol)':[r.data['vina_affinity'] for r in res],
'Total number of h-bonds':hbtotal,
'Strict number of h-bonds':hb_strict,
'Pi-stacking':pi,
'Pi-stacking (parallel)':pi_par,
'Pi-stacking(perp)':pi_perp,
'Hydrophobic':phob_all,
})
sort_aff = sort_table.sort_values(by=['Affinity (kcal/mol)'],ascending=False) #отсортируем по афинности
sort_aff
Посмотрим также на лучшие позы для других заместителей
C8H13NO8 По данным в таблице - 2 водородные связи, которые удовлетворяют всем условиям. Посмотрев, на окружение 3.5 ангстрема и потенциальные связи, данные связи были обнаружены.
C7H13NO7. Для данного заместителя было найдено 3 связи, удовлеворяющие критериям, всего 4. Все связи были найдены.
C7H15N2O6+ Тут в принципе 4 водородных связей и правда есть. Но вот еще 4 - непонятно. Ну и сам лиганд явно сильно закрутило.
C7H13NO6 А в данном случае все 5 водородных связей есть!
