Докинг
Докинг производился при помощи AutodockVina в GoogleColab, установленной согласно инструкции.
Лиганды (NAG с различными заместителями) были созданы при помощи редактирования smiles:
lig_smiles = 'CC(=O)NC1C(C(C(OC1O)CO)O)O'
mols = []
images = []
for rep in ['O', '[NH3+]', '[H]', 'C1=CC=C(C=C1)', 'C(=O)([O-])']:
lig_smiles_ = lig_smiles
lig_smiles_ = lig_smiles_.replace('CC(=O)N', f'{rep}C(=O)N')
m = oddt.toolkit.readstring('smi', lig_smiles_)
if not m.Mol.HasProp('3D'):
m.make3D(forcefield='mmff94', steps=150)
mols.append(m)
images.append((SVG(copy.deepcopy(m).write('svg'))))
display_svg(*images)
Белок (лизоцим 1lmp) был подготовлен к докингу, определён центр лиганда в структуре
import mdtraj as md
u = md.load('1lmp.pdb')
pdb = u.topology
prot = u.atom_slice(range(1, 999)) # граница белка найдена вручную по списку
prot.save_pdb('1lmp-clean.pdb', force_overwrite=True)
lig = u.atom_slice(range(1000, 1041))
lig_center = lig.xyz[0].mean(axis=0) * 10 # т.к. нм, а нужны А
prot = next(oddt.toolkit.readfile('pdb','1lmp-clean.pdb'))
Произведён докинг лигандов выше в кубическую область с центром в центре лиганда и размерами сторон 20, выдача по пять поз для каждого лиганда (num_modes), параметр, определящий "количество вычислительных усилий" (exhaustiveness) 8 (относительно низкий, может не найти оптимальную позу, которую нашла бы с большим exhaustiveness), сохранены будут только позы, отстающие от лучшей найденной не более, чем на 3 ккал/моль (energy_range=3).
dock_obj= oddt.docking.AutodockVina.autodock_vina(
protein=prot,
size=(20,20,20),
center=lig_center,
executable='/content/vina',
autocleanup=True,
num_modes=5)
print(" ".join(dock_obj.params))
> --center_x 42.507587 --center_y 43.94649 --center_z 28.304195 --size_x 20 --size_y 20 --size_z 20 --exhaustiveness 8 --num_modes 5 --energy_range 3
res = dock_obj.dock(mols,prot)
Сохраним структуры и выпишем характеристики поз:
res_sorted = sorted(res, key=lambda x: float(x.data['vina_affinity']))
for i,r in enumerate(res_sorted):
r.write(filename='r%s.pdb' % i, format='pdb', overwrite=True)
print(''.join(map(lambda x: f'{x: <12}', ['i', 'formula', 'affinity', 'rmsd_ub', 'resn', 'hbs', 'stack', 'phob'])))
for i,r in enumerate(res_sorted):
formula = Chem.rdMolDescriptors.CalcMolFormula(r.Mol)
hbs = oddt.interactions.hbonds(prot,r)
stack= oddt.interactions.pi_stacking(prot,r)
phob = oddt.interactions.hydrophobic_contacts(prot,r)
print(''.join(map(lambda x: f'{x: <12}', [i,
formula,
r.data['vina_affinity'],
r.data['vina_rmsd_ub'],
r.residues[0].name,
len(hbs[0]),
len(stack[0]),
len(phob[0])
])))
Результат
И табличка результатов (affinity — энергия, rmsd_ub — rmsd относительно лучшей позы для этого лиганда, hbs, stack и phob — водородные связи, стекинг и гидрофобные взаимодействия соответсвенно):
i formula affinity rmsd_ub resn hbs stack phob 0 C13H17NO6 -6.423 0 UNL 6 1 4 1 C13H17NO6 -6.415 4.049 UNL 9 0 4 2 C13H17NO6 -6.387 6.078 UNL 3 0 2 3 C13H17NO6 -6.375 7.407 UNL 3 1 5 4 C13H17NO6 -6.31 7.561 UNL 5 0 4 5 C7H13NO7 -5.819 0 UNL 10 0 0 6 C7H13NO7 -5.733 4.713 UNL 9 0 0 7 C7H13NO7 -5.671 5.623 UNL 12 0 0 8 C8H12NO8- -5.649 0 UNL 9 0 0 9 C8H12NO8- -5.605 5.203 UNL 6 0 0 10 C7H13NO7 -5.604 4.645 UNL 10 0 0 11 C7H13NO7 -5.554 3.586 UNL 4 0 0 12 C8H12NO8- -5.533 3.786 UNL 14 0 0 13 C7H13NO6 -5.438 0 UNL 14 0 0 14 C8H12NO8- -5.414 2.59 UNL 4 0 0 15 C8H12NO8- -5.288 3.003 UNL 12 0 0 16 C7H15N2O6+ -5.058 0 UNL 5 0 0 17 C7H13NO6 -4.969 3.572 UNL 10 0 0 18 C7H13NO6 -4.93 2.854 UNL 12 0 0 19 C7H15N2O6+ -4.901 5.808 UNL 10 0 0 20 C7H13NO6 -4.894 6.244 UNL 10 0 0 21 C7H15N2O6+ -4.809 5.844 UNL 10 0 0 22 C7H15N2O6+ -4.808 5.428 UNL 8 0 0 23 C7H13NO6 -4.671 5.982 UNL 4 0 0 24 C7H15N2O6+ -4.597 3.044 UNL 4 0 0
Первые пять строчек занимает лиганд с фенильным радикалом, вероятно, за счёт образования стекинга и гидрофобных контактов, в то время как остальные имеют лишь водородные связи (причём даже увеличение числа водородных связей вдвое не перебивает один стекинг+гидрофобику).
Поскольку автоматическое определение контактов может быть не вполне верным, сравним размеченные контакты с визуальным анализом
# hbonds according to oddt
protein ligand
id atype id atype
358 O.2 19 O.2
374 O.3 18 O.2
398 O.2 19 O.2
359 N.am 17 O.2
374 O.3 18 O.2
459 N.am 19 O.2
В целом найденные водородные связи действительно выглядят возможными кроме 359 N.am – 17 O.2, где у азота аспарагина водороды должны лежать в плоскости N-C-O и потому не могут донироваться под таким углом, а как акцептор он слаб в силу резонансного сопряжения с кислородом. Также из связей 358 O.2 — 19 O.2 и 398 O.2 – 19 O.2 единомоментно может существовать лишь одна, поскольку донорный водород у 19 O.2 лишь один. Ещё oddt по какой-то причине посчитал дважды связь 374 O.3 – 18 O.2 Таким образом, число возородных связей в описании oddt, вероятно, завышено.
В целом можно заметить, что в структурах отсутсвуют протоны, так что, видимо, в текущей версии они не добавляются автоматически и за этим нужно было проследить. Вероятно, это улучшило бы разметку взаимодействий и, возможно, улучшило бы результат докинга.
Методы из статьи про ODDT
В зависимости от того, какие дополнительные данные о лигандах у нас имеются — количественные (Kd, IC50) или качественные (активен/неактивен) можно было бы использовать регрессоры и классификаторы соответвенно для предсказания этих величин на наших лигандах после обучения на известных. Однако конкретно 1lmp в базе данных PDBbind нет, а более продвинутый поиск он позволяет только зарегистрированным пользователям, поэтому посмотреть, какие подходящие данные там есть и предложить конкретный метод не вышло.