Докинг низкомолекулярных лигандов в структуру белка¶
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.Draw import IPythonConsole
import pandas as pd
import mdtraj as md
u = md.load('lys_lig.B99990001.pdb')
pdb = u.topology
for i,r in enumerate(pdb.atoms):
print(i,r, end=("\t" if i % 8 else "\n"))
#посмотрим на атомы
# сохраним pdb без лиганда
u2 = u.atom_slice(range(1166))
u2.save_pdb('lys_lig-clean.pdb', force_overwrite=True)
#определим геометрический центр лиганда
lig = u.atom_slice(range(1166, 1209))
center = np.mean(lig.xyz, axis=1) * 10 # В ангстремы
center = center.reshape(3,)
center
prot = next(oddt.toolkits.rdk.readfile('pdb', 'lys_lig-clean.pdb'))
print('Is the first mol in model is protein?',prot.protein,':) and MW of this mol is:', prot.molwt )
В смысле, не белок?
smiles = ["CC(=O)NC1C(O)OC(CO)C(O)C1O", # NAG
"OC(=O)NC1C(O)OC(CO)C(O)C1O", # -OH
"[NH3+]C(=O)NC1C(O)OC(CO)C(O)C1O", # -NH3+
"C(=O)NC1C(O)OC(CO)C(O)C1O", # -H
"C=1(C=CC=CC1)C(=O)NC1C(O)OC(CO)C(O)C1O", # -Ph
"C(C([O-])=O)(=O)NC1C(O)OC(CO)C(O)C1O" # -COO-
]
mols= []
images =[]
for s in smiles:
m = oddt.toolkits.rdk.readstring('smi', s)
m.make3D(forcefield='mmff94', steps=150)
mols.append(m)
print(Chem.rdMolDescriptors.CalcMolFormula(m.Mol))
images.append((SVG(copy.deepcopy(m).write('svg'))))
display(*images)
Докинг
center
# create docking object
# центром докинга дадим геометрический центр лиганда,
# найденный нами ранее - приблизительное положение сайта связывания
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,size=(20,20,20),center=center,
executable='C:\\Program Files (x86)\\The Scripps Research Institute\\Vina\\vina.exe',autocleanup=True, num_modes=9)
print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params))
#C:\\Program Files (x86)\\The Scripps Research Institute\\Vina\\vina.exe
center_x/y/z - геометрический центр лиганда
size_x/y/z - размер области, в которой будет производиться поиск
exhaustiveness = параметр, который говорит программе, насколько тщательно производить поиск
num_modes - количество конформаций лиганда в активном центре, которые будут найдены программой
energy_range - максимально допустимый разброс энергии для разных моделей одного и того же лиганда
res = dock_obj.dock(mols,prot)
# отсортируем результаты по аффинности (параметр vina_affinity):
res_sorted = sorted(res, key = lambda x : float(x.data['vina_affinity']))
#что вышло
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))
Анализ докинга
hbtotal = []
hbstrict = []
phob = []
formulas = []
for i,r in enumerate(res_sorted):
formula = Chem.rdMolDescriptors.CalcMolFormula(r.Mol)
formulas.append(formula)
int1, int2, strict = oddt.interactions.hbonds(prot,r)
hbtotal.append(len(int1))
hbstrict.append(strict.sum())
ph1, ph2 = oddt.interactions.hydrophobic_contacts(prot,r)
phob.append(len(ph1))
sort_table = pd.DataFrame({'Formula': formulas,
'Affinity, kcal/mol':[r.data['vina_affinity'] for r in res_sorted],
'Total number of h-bonds':hbtotal, 'Strict number of h-bonds':hbstrict, 'Hydrophobic bonds':phob,
'RMSD':[r.data['vina_rmsd_ub'] for r in res_sorted]})
sort_table[:]
Гидрофобные взаимодействия выявились у лигандов с ароматическими группами в составе. Самая лучшшая афииность у лиганда C13H17NO6. Неожиданно мало водородных связей у лучших хитов.
Визуализируем: посмотрим лучший (номер 36) и худший (номер 35) лиганд, а также лиганд с наибольшим количеством гидрофобных связей (номер 38).
display(Image('lig_36.png'))
Ароматическая группа лучшего лиганда находится в гидрофобном кармане белка (валином и CH2 группой аспарагина). К тому же, дополнительно позицию лиганжа укрепляют 2 водородные связи.
display(Image('lig_35.png'))
Худший лиганд. Всего одна водородная связь. Расположен в том же кармане, что и предыдущий лиганд. Позиция не усилена дополнительно гидрофобными взаимодействиями или стекингом.
display(Image('lig_38.png'))
Лиганд с наибольшим количеством гидрофобных связей. Дополнительно имеется 5 водородных связей. Высокая афинность. Странно, что по энергии связывания лиганд проигрывает лучшему хиту.