Суть задания заключается в предсказани свзывания фрагмента антитела (VHH domain) c панкреатической альфа-амилазой с использованием программы ZDOCK.
Для успешной работы zrank и zdock из pdb файлов надо удалить строчки MODEL и TER Установлено, что:
программе zrank очень нужна пустая строка до начала записей об атомах.
полезно убрать именно вторую строчку из копии zdock.out (zdock.out.cp) и использовать её для работы zrank.
Укажем путь к ZDOCK:
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
export PATH=${PATH}:/home/preps/golovin/progs/bin
Скопируем в рабочую директорию файлы amilase.pdb,camelid.pdb и uniCHARMM из:
%%bash
wget http://kodomo.cmm.msu.ru/~golovin/zdock/uniCHARMM
wget http://kodomo.cmm.msu.ru/~golovin/zdock/amylase.pdb
wget http://kodomo.cmm.msu.ru/~golovin/zdock/camelid.pdb
Добавим водороды к обоим pdb-файлам, исполльзуя pdb2gmx из GROMACS:
! gmx pdb2gmx -f amylase.pdb -o amylase_h.pdb -p -ff gromos53a6 -ignh -water none
! gmx pdb2gmx -f camelid.pdb -o camelid_h.pdb -p -ff gromos53a6 -ignh -water none
Параметры запуска:
Утилита mark_sur для каждого атома определяет ASA (accesible surface area - поверхность, доступная для растворителя), которая обычно определяется "качением" молекулу воды радиуса 1.4 А по срезам структуры. Если атом имеет ASA больше 1 А, он отмечается утилитой mark_sur как поверхностный. Проведём препроцессинг файлов pdb:
Usage: mark_sur in_pdb_file out_pdb_file
Mark the surface residues of a PDB file
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
export PATH=${PATH}:/home/domain/data/prog/zdock/zdock3.0.2_linux_x64
mark_sur amylase_h.pdb amylase_m.pdb > amylase.txt
mark_sur camelid_h.pdb camelid_m.pdb > camelid.txt
# Формат файлов (последние 10 строк 1-го файла и перые 10 строк 2-го файла)
with open("amylase.txt") as rf1, open("camelid.txt") as rf2:
f1 = rf1.readlines()
f2 = rf2.readlines()
for i in f1[-10:]:
print(i.strip())
print("\n")
for j in f2[:10]:
print(j.strip())
Посмотрим, как работать с ZDOCK:
Usage: zdock [options] -R [receptor.pdb] -L [ligand.pdb]
Standard PDB format files must be processed by mark_sur to add atom radius, charge and atom type. If you know that some atoms are not in the binding site, you can block them by changing their atom type (column 55-56) to 19.
Available options:
-o filename : set the output filename (default to zdock.out)
-N int : set the number of output predictions, (default to 2000)
-S int : set the seed for the randomization of the starting PDBs (default to no randomization)
-D : set rotational sampling to be dense, used only if zdock output will be further processed by more detailed binding energy calculations (default to none)
-F : fix the receptor, preventing its rotation and switching with ligand during execution.
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
export PATH=${PATH}:/home/domain/data/prog/zdock/zdock3.0.2_linux_x64
zdock -N 10 -R camelid_m.pdb -L amylase_m.pdb
! cat zdock.out
Вырезаем строчки, содержащие MODEL и TER, в файлах camelid_m.pdb и amylase_m.pdb и сохраняем в camelid_m.pdb и amylase_m.pdb:
with open("amylase_m.pdb") as wf1, open("camelid_m.pdb") as wf2:
w1 = wf1.readlines()
w2 = wf2.readlines()
with open("amylase_new.pdb", "w+") as wwf1, open("camelid_new.pdb", "w+") as wwf2:
for l1 in w1:
if l1.find("MODEL") == -1 and l1.find("TER") == -1:
wwf1.write(l1)
for l2 in w2:
if l2.find("MODEL") == -1 and l1.find("TER") == -1:
wwf2.write(l1)
Создаем файл-копию zdock.out.cp с заменой camelid_m.pdb и amylase_m.pdb на сamelid_new.pdb и amylase_new.pdb и удаленной второй строчкой:
with open("zdock.out") as rf, open("zdock.out.cp", "w+") as wf:
new_z = []
for z in rf:
z.replace("amylase_m.pdb", "amylase_new.pdb")
z.replace("camelid_m.pdb", "camelid_new.pdb")
new_z.append(z)
z = new_z[:1] + new_z[2:]
for n in z:
wf.write(n)
! cat zdock.out.cp
Проведем предварительный анализ результатов докинга с помощью утилиты zrank:
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
export PATH=${PATH}:/home/domain/data/prog/zdock/zrank_linux_64bit
zrank zdock.out.cp 1 10
sort -n -k2 zdock.out.cp.zr.out | head
%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/preps/golovin/progs/lib
export PATH=${PATH}:/home/domain/data/prog/zdock/zdock3.0.2_linux_x64
ln -s /home/domain/data/prog/zdock/zdock3.0.2_linux_x64/create_lig create_lig
create.pl zdock.out
Скачаем структуру 1KXT из банка PDB:
%%bash
wget https://files.rcsb.org/view/1KXT.pdb
Визуально сравним результаты с известной структурой 1kxt (модели расположены в порядке выдачи zrank):
from IPython.display import display,Image
a = Image(filename='vhh8.png')
b = Image(filename='vhh10.png')
c = Image(filename='vhh1.png')
d = Image(filename='vhh6.png')
e = Image(filename='vhh4.png')
f = Image(filename='vhh3.png')
g = Image(filename='vhh9.png')
h = Image(filename='vhh5.png')
i = Image(filename='vhh7.png')
j = Image(filename='vhh2.png')
display(a, b, c, d, e, f, g, h, i, j)
Посчитаем RMSD с помощью rms пакета GROMACS, чтобы найти модель наиболее близкую к результату РСА:
gmx rms -f complex.i.pdb -s 1KXT.pdb -o rmsd_i, i = 1, ..., 10
RMSD считалось по группе "Protein", т. е. 1
Structure | RMSD |
---|---|
1 | 4.7041268 |
2 | 4.3613758 |
3 | 4.7204866 |
4 | 4.3944025 |
5 | 4.6301455 |
6 | 4.4350629 |
7 | 4.4898782 |
8 | 4.6887312 |
9 | 4.3714395 |
10 | 4.8258491 |
Согласно данным RMSD наилучшей структурой является 2, однако по результатам zrank лидирует 8, а 2, наоборот, занимает последнее место в топ-10.
Визуализация: PyMOLWiki
Макромолекулярный докинг: ZDOCK