Докинг низкомолекулярных лигандов в структуру белка¶
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.Chem.Draw import IPythonConsole
from rdkit import Chem
Подготовка белка¶
import mdtraj as md
u = md.load('seq_liganded.B99990001.pdb')
pdb = u.topology
for i,r in enumerate(pdb.atoms):
print(i,r)
#посмотрим на атомы
protein_end = 1256 # Последний атом белка
ligand_end = 1299 # последний атом лиганда
Сохраним белок без лиганда.
u2 = u.atom_slice(range(protein_end + 1))
u2.save_pdb('seq_liganded.B99990001_clean.pdb', force_overwrite=True)
Найдём геометрический центр лиганда.
lig = u.atom_slice(range(protein_end + 1, ligand_end))
lig_center = np.mean(lig.xyz, axis=1) * 10
lig_center = lig_center.reshape(3,)
lig.save_pdb('lig.pdb', force_overwrite=True)
print(lig_center)
(np.max(lig.xyz, axis=1) - np.min(lig.xyz, axis=1)) * 10
Подготовка белка для докинга¶
prot = next(oddt.toolkits.rdk.readfile('pdb','seq_liganded.B99990001_clean.pdb'))
print('is the first mol in 1lmp a protein?', prot.protein, ':) and MW of this mol is:', prot.molwt)
По умолчанию, наш белок воспринимается не как белок :(
Лиганды для докинга¶
smiles = ["CC(=O)NC1C(O)OC(CO)C(O)C1O", # Common 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_svg(*images)
Докинг¶
#create docking object
dock_obj= oddt.docking.AutodockVina.autodock_vina(protein=prot,
size=(20, 20, 20),
center=lig_center,
executable='./vina',
autocleanup=True,
num_modes=9,
prefix_dir='.')
print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params))
Выдача AutoDock Vina:
| Параметр | Суть |
|---|---|
| center_x | Центр коробки для докинга по x |
| center_y | Центр коробки для докинга по y |
| center_z | Центр коробки для докинга по z |
| size_x | Размер коробки по x |
| size_y | Размер коробки по y |
| size_z | Размер коробки по z |
| exhaustiveness | Долгота поиска глобального минимума (обычно эвристически как trade-off между числом атомов и временем работы) |
| num_modes | Наибольшее число моделей в выдаче |
| energy_range | Максимальное различие энергий между лучшей и худшей моделями в выдаче |
Запускаем докинг.
res = dock_obj.dock(mols, prot)
Результаты докинга¶
Посмотрим на результаты докинга в исходном виде.
for i,r in enumerate(res):
print("%4d%10s%8s%8s%8s" %(i, Chem.rdMolDescriptors.CalcMolFormula(r.Mol), r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name))
Видим, что лиганды сгруппированы по их типу, сортируются от лучшего к худшему внутри своего типа.
sorted_res = sorted([(i, Chem.rdMolDescriptors.CalcMolFormula(r.Mol), r.data['vina_affinity'], r.data['vina_rmsd_ub'], r.residues[0].name) for i, r in enumerate(res)],
key=lambda x: x[2],
reverse=True)
print('i\tformula\t\taff\trmsd\tname')
for r in sorted_res:
print(f"{r[0]}\t{r[1]}\t{r[2]}\t{r[3]}\t{r[4]}")
Гидрофобика победила всех: все Ph лучше всех остальных. Затем идут лиганды COO- и между ними обычный NAG.
Анализ докинга¶
Посмотрим на особенности взаимодействия лучших лигандов: все Ph, два COO- и один обычный.
best_res = [(i, res[i]) for i in list(range(36, 48)) + [0]]
print("\t".join(['No',
'formula',
'\taff',
'HB',
'HB_str',
'Stack',
'St_par',
'St_perp',
'Hphob']))
for i, r in best_res:
# hbonds
mol1_atoms, mol2_atoms, strict_bonds = oddt.interactions.hbonds(prot, r)
hbonds_total = len(mol1_atoms)
hbonds_strict = np.sum(strict_bonds)
# stacking
mol1_rings, mol2_rings, strict_parallel, strict_perp = oddt.interactions.pi_stacking(prot, r)
stack_total = len(mol1_rings)
stack_parallel_strict = np.sum(strict_parallel)
stack_perp_strict = np.sum(strict_perp)
mol1_hphob, mol2_hphob = oddt.interactions.hydrophobic_contacts(prot, r)
hphob_total = len(mol1_hphob)
print(i,
Chem.rdMolDescriptors.CalcMolFormula(r.Mol),
r.data['vina_affinity'],
hbonds_total,
hbonds_strict,
stack_total,
stack_parallel_strict,
stack_perp_strict,
hphob_total,
sep='\t')
Ph заместитель побеждает за счёт гидрофобики и некоторых водородных связей. Забавно, что у лиганда 0 есть гидрофобика, а у лигандов 46-47 её нет. Стекинг нигде не нашёлся.
Посмотрим визуально на лиганды 36, 37, 46, 47 и 0.
Визуализация¶
selected_vis = ((i, r) for i, r in best_res if i in [36, 37, 46, 47, 0])
for i, r in selected_vis:
r.write('pdb', f'lig_{i}.pdb', overwrite=True)
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd
Сделаем красиво.
cmd.reinitialize()
cmd.bg_color('white')
cmd.set('antialias', 2)
cmd.set('ray_trace_mode', 3)
Напишем функцию для подсчёта водородных связей (отсюда).
def find_hbonds(sele1, sele2, cutoff=3.5):
cmd.select('don', '(elem n,o and (neighbor hydro))')
cmd.select('acc, (elem o or (elem n and not (neighbor hydro)))')
sele1_hba = f'({sele1} and acc)'
sele2_hbd = f'({sele2} and don)'
sele1_hbd = f'({sele1} and don)'
sele2_hba = f'({sele2} and acc)'
cmd.dist('HBA', sele1_hba, sele2_hbd, cutoff)
cmd.dist('HBD', sele1_hbd, sele2_hba, cutoff)
amount = (cmd.count_atoms(f'{sele1_hba} within {cutoff} of {sele2_hbd}') +
cmd.count_atoms(f'{sele1_hbd} within {cutoff} of {sele2_hba}'))
cmd.delete('don')
cmd.delete('acc')
cmd.hide('labels', 'HBA')
cmd.hide('labels', 'HBD')
return amount
Загрузим модели.
cmd.load('seq_liganded.B99990001_clean.pdb', 'prot_model')
cmd.load('lig.pdb', 'lig_init')
cmd.load('lig_36.pdb', 'lig_36')
cmd.load('lig_37.pdb', 'lig_37')
cmd.load('lig_46.pdb', 'lig_46')
cmd.load('lig_47.pdb', 'lig_47')
cmd.load('lig_0.pdb', 'lig_0')
cmd.zoom()
Добавим водороды к белку и исходному лиганду.
cmd.h_add('prot_model')
cmd.h_add('lig_init')
Покрасим белок по гидрофобике (отсюда).
cmd.run('color_h.py')
cmd.do('color_h prot_model')
cmd.show('surface', 'prot_model')
Посмотрим на исходный лиганд и его водородные связи.
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_init')
cmd.orient('lig_init')
cmd.zoom('lig_init', 7)
cmd.rotate('x', 90)
cmd.rotate('y', 20)
find_hbonds('lig_init', 'prot_model')
cmd.png('lig_init.png', 1000, 1000, ray=1)
display(Image('lig_init.png'))
Вроде бы, всё нормально.
Посмотрим на лиганд 36 (Ph).
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_36')
cmd.delete('HBA or HBD')
find_hbonds('lig_36', 'prot_model')
0 водородных связей, хотя ранее было показано 5 (3).
cmd.png('lig_36.png', 1000, 1000, ray=1)
display(Image('lig_36.png'))
Лежит красиво. Заметим, фенольное кольцо лежит в менее гидрофобном кармане, чем пирановое кольцо, зато дальше от воды.
Посмотрим на лиганд 37 (Ph).
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_37')
cmd.delete('HBA, HBD')
find_hbonds('lig_37', 'prot_model')
Опять 0 водородных связей на расстоянии 3.5Å!
cmd.png('lig_37.png', 1000, 1000, ray=1)
display(Image('lig_37.png'))
Ну, выглядит нормально, но больше экспонировано на поверхность.
Посмотрим на лиганд 46 (COO-).
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_46')
cmd.delete('HBA or HBD')
find_hbonds('lig_46', 'prot_model')
3 водородные связи, а не 6, как показано ранее.
cmd.png('lig_46.png', 1000, 1000, ray=1)
display(Image('lig_46.png'))
Карбоксильная группа уехала внутрь белка, цепляясь там водородной связью.
Посмотрим на лиганд 47 (COO-).
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_47')
cmd.delete('HBA or HBD')
find_hbonds('lig_47', 'prot_model')
1 водородная связь, а не 6.
cmd.png('lig_47.png', 1000, 1000, ray=1)
display(Image('lig_47.png'))
У этого лиганда карбоксильная группа торчит наружу.
Посмотрим на лиганд 0 (NAG).
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_0')
cmd.delete('HBA or HBD')
find_hbonds('lig_0', 'prot_model')
Мда, 0 водородных связей (а было 2).
cmd.png('lig_0.png', 1000, 1000, ray=1)
display(Image('lig_0.png'))
Тоже болтается в кармане, говорят, что 2 пары атомов вовлечены в гидрофобные взаимодействия.
36, 37, 47 и 0 нависают полярным кислородом над TRP73 (такая гидрофобная площадка на переднем плане слева, которая поддерживает лиганды). 36 вообще красавец: Он цепляется OH-группой за TRP, далее фурановое кольцо лежит в гидрофобике, а фенил прячется в карман (не очень-то гидрофобный). 46 тоже хорошо проходится по примерно гидрофобному жёлобу, а его карбоксильная группа цепляется за водородные связи в кармане. Итог: 36 и 46 хорошо задочились, причём как аналитически, так и визуально.