Докинг лизоцима¶
Проведем импорт библиотек для работы
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 mdtraj as md
u = md.load('lys_bommo.pdb')
pdb = u.topology
for i,r in enumerate(pdb.atoms):
if str(r).startswith('NAG') or str(r).startswith("HOH") or str(r).startswith("NDG"):
print(i, r)
Сохраним новый файл без лиганда и проверим правомерность наших манипуляций
u2 = u.atom_slice(range(1096))
u2.save_pdb('lys-clean.pdb', force_overwrite=True)
m = md.load('lys-clean.pdb')
pdb_n = m.topology
for it,ra in enumerate(pdb_n.atoms):
print(it,ra)
#нахождение геометрического центра лиганда
lig = u.atom_slice(range(1096,1123)) #список номера атомов лиганда
for v in lig.xyz:
lig_center = np.mean(v, axis=0)*10
lig_center
lig.save_pdb('lig.pdb', force_overwrite=True)
Подготовим белок и лиганды для докинга
prot = next(oddt.toolkit.readfile('pdb','lys-clean.pdb'))
print('is it the first mol in 1lmp is 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='/usr/bin/vina',
autocleanup=True,
num_modes=9,
prefix_dir='.')
print(dock_obj.tmp_dir)
print(" ".join(dock_obj.params))
res = dock_obj.dock(mols, prot)
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]}")
Выше всего афинность у комплекса белка с фениловым производным гликозида, затем идет производная с дополнительной карбоксильной группой и лучший из комплексов с исходным лигандом.
import pandas as pd
hbtotal = []
hbstrict = []
phob = []
stack = []
nice_table= pd.DataFrame(columns=['Formula', 'Affinity, kcal/mol', 'Total number of h-bonds', 'Strict number of h-bonds',
'Stacking total', 'Strict parallel stacking', 'Strict perpendicular stacking', 'Hydrophobic',
'RMSD'])
for i,r in enumerate(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)
# phobia
mol1_hphob, mol2_hphob = oddt.interactions.hydrophobic_contacts(prot, r)
hphob_total = len(mol1_hphob)
nice_table=nice_table.append({'Formula': Chem.rdMolDescriptors.CalcMolFormula(r.Mol),
'Affinity, kcal/mol': r.data['vina_affinity'],
'Total number of h-bonds': hbonds_total,
'Strict number of h-bonds': hbonds_strict,
'Stacking total': stack_total,
'Strict parallel stacking': stack_parallel_strict,
'Strict perpendicular stacking': stack_perp_strict,
'Hydrophobic': hphob_total,
'RMSD': r.data['vina_rmsd_ub']}, ignore_index=True)
nice_table.sort_values(by=['Affinity, kcal/mol'], ascending=False)[:]
Стэкинг-взаимодействия не были выявлены ни у одной модели. Фенол-производное сахара обладало высокой афинностью за счет большого количества гидрофобных взаимодействий и водородных связей. Для дальнейшей визуализации были выбраны модели под номером 36, 37, 46, 45, 0, 9, 18
for i in [36, 37, 46, 45, 0, 9, 18]:
res[i].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)
cmd.load('lys-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_45.pdb', 'lig_45')
cmd.load('lig_18.pdb', 'lig_18')
cmd.load('lig_9.pdb', 'lig_9')
cmd.load('lig_0.pdb', 'lig_0')
cmd.zoom()
Добавим в программу функцию для подсчета водородных связей (Источник pymolwiki) и дополним лиганд и исходный белок водородами.
cmd.h_add('prot_model')
cmd.h_add('lig_init')
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.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 (фенольный довесок). Водородной связи из докинга pymol не нашел. Похоже, что гидрофобно взаимодействует дальний от нас бок фенольного кольца
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_36')
cmd.delete('HBA or HBD')
find_hbonds('lig_36', 'prot_model')
cmd.png('lig_36.png', 1000, 1000, ray=1)
display(Image('lig_36.png'))
В модели 37 было найдено только 3 водородные связи, гидрофобным фенольным кольцом он тоже обращен в сторону гидрофобного участка белка.
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_37')
cmd.delete('HBA or HBD')
find_hbonds('lig_37', 'prot_model')
cmd.png('lig_37.png', 1000, 1000, ray=1)
display(Image('lig_37.png'))
В модели 46 pymol детектирует две водородные связи из трех предполагаемых в докинге
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_46')
cmd.delete('HBA or HBD')
find_hbonds('lig_46', 'prot_model')
cmd.png('lig_46.png', 1000, 1000, ray=1)
display(Image('lig_46.png'))
В модели 45 видна одна из шести водородных связей, причем от оксолатного остатка.
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_45')
cmd.delete('HBA or HBD')
find_hbonds('lig_45', 'prot_model')
cmd.png('lig_45.png', 1000, 1000, ray=1)
display(Image('lig_45.png'))
Оригинальный лиганд прошел докинг, потеряв две водородные связи
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_0')
cmd.delete('HBA or HBD')
find_hbonds('lig_0', 'prot_model')
cmd.png('lig_0.png', 1000, 1000, ray=1)
display(Image('lig_0.png'))
Модель 9 с гидроксигруппой имеет 3 из пяти водородных связей. Стоит отметить что полярная гидроксигруппа развернута из кармана белка в направление сольвента.
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_9')
cmd.delete('HBA or HBD')
find_hbonds('lig_9', 'prot_model')
cmd.png('lig_9.png', 1000, 1000, ray=1)
display(Image('lig_9.png'))
Модель 18 - с аминогруппой. Из пяти контактов pymol видит три.
cmd.hide('everything', 'lig*')
cmd.show('sticks', 'lig_18')
cmd.delete('HBA or HBD')
find_hbonds('lig_18', 'prot_model')
cmd.png('lig_18.png', 1000, 1000, ray=1)
display(Image('lig_18.png'))
В общем, pymol недооценивает найденные в ходе докинга водородные связи, а сам доккинг не всегда способен выдать нам исходную структуру. В белке, похоже, есть четко гидрофобный карман, что повышает афинность фенольных производных. Может стоит проверить еще и с циклогексановым довеском?