Structural bionformatics. Homework №9¶
Elen Tevanyan
Зная, что мне предстоит копировать все файлы из папки на Ломоносов, все материалы я внесла в отдельную папку в моей директории. Т.к. ноутбук создан внутри него, дополнительно перезадавать директорию не нужно - все записывается ко мне в подпапку "hwtev".
Часть команд запускались отсюда; часть команд - из терминала, они указаны в текстовых ячейках
Подготовка к моделированию¶
import __main__
__main__.pymol_argv = ['pymol','-qc'] # Pymol: quiet and no GUI
from time import sleep
import time
import pymol
import IPython.display
from IPython.display import display,display_svg,SVG,Image, HTML
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
pymol.finish_launching()
defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=500, height=500, sleep=2, filename=defaultImage):
## To save the rendered image
pymol.cmd.ray(width, height)
pymol.cmd.png('/tmp/pymolimg.png')
time.sleep(sleep)
Скачиваем файлы, указанные в задании:
https://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac/SelfAss
%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
%%bash
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
pymol.cmd.reinitialize()
pymol.cmd.load('dppc.pdb')
pymol.cmd.zoom('dppc',5)
prepareImage()
Image(defaultImage)
pymol.cmd.reinitialize()
pymol.cmd.load('b_64.pdb')
pymol.cmd.zoom('b_64',5)
prepareImage()
Image(defaultImage)
Был один липид, создали ячейку с 64, что и видим на двух картинках выше.
В файле btop заменила число 34 на 64 и делаю небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды.
%%bash
editconf -f b_64.gro -o b_ec -d 0.5
Запуск отсюда - довольно глупая затея, поэтому команды ниже запустила в терминале, где изменения силы хорошо видны.$$F_{max}^0 = 4.37970e+05$$ $$F_{max}^{fin} = 6.16887e+02$$
grompp -f em -c b_ec -p b -o b_em -maxwarn 2 mdrun -deffnm b_em -v
Следующую серию команд запускала также в терминале.
genbox -cp b_em -p b -cs spc216 -o b_s
grompp -f em -c b_s -p b -o b_empr -maxwarn 1
mdrun -deffnm b_empr -v
grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
Визуализируем pdb-файл с добавленной водой.
pymol.cmd.reinitialize()
pymol.cmd.load('b_s.pdb')
pymol.cmd.zoom('b_s',5)
prepareImage()
Image(defaultImage)
pymol.cmd.reinitialize()
pymol.cmd.load('b_pr.pdb')
pymol.cmd.zoom('b_pr',5)
prepareImage()
Image(defaultImage)
На глаз не особо много что изменилось, но мы "тресли" воду, и можно заметить, что вода более равномерно и плотно расположилась внутри ячейки.
Моделирование¶
Перенос на суперкомпьютер
cd
scp -r ./hse/tevanyan/hwtev lom:_scratch/chem
Тестовое моделирование
ssh lom
cd _scratch/chem/hwtev
cp /home/users/golovin/progs/share/gromacs/top/residuetypes.dat .
grompp -f md -c b_pr -p b -o b_md -maxwarn 1
sbatch -n 4 -e error.log -o output.log -t 5 -p test impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm b_md -v
Основное моделирование : номер задачи 1653801
sbatch -N1 --ntasks-per-node=2 -e error-gpu.log -o output.log -t 350 -p gpu impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm b_md -v
Копирование результатов :
cd ./hse/tevanyan/hwtev
scp -r lom:_scratch/chem/hwtev ./
Анализ результатов¶
Визуализация¶
cd ./hse/tevanyan/hwtev
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
Обнаружила, что первой командой выдается неадектваный результат, очень хаотичный итог получается, слои "смешиваются", поэтому запустила вторую команду.
pymol.cmd.reinitialize()
pymol.cmd.load('b_pbc_1.pdb')
prepareImage()
Image(defaultImage)
Скачала файл к себе, записала видео с помощью десктопного PyMol и утилиты ffmpeg. Т.к. у меня образовались кадры всех моделирований, могла спокойно их рассмотреть.
import io
from base64 import b64encode
from IPython.core.display import HTML, display
video = io.open('bilip2.mp4', 'r+b').read()
video_encoded = b64encode(video)
html_code = '<video controls alt="PyMol Movie" src="data:video/mp4;base64,{}" type="video/mp4">'.format(video_encoded)
display(HTML(html_code))
Это финальная модель.
Image('mov0101.png')
С 55 красиво собирались как гидрофобные хвосты, так и гидрофильные головки - ну, на мой вкус, конечно же. Это bilayer in water t= 27000.00000
Image('mov0055.png')
Площадь липида¶
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
В файле есть 24 некрасивые строчки. Планировала удалить командой sed, но она неудачно прошла для файла (выдавала совершенно пустой), поэтому стираю вручную, гружу файл и сортирую на 4 вектора: время, х-координата, y-координата, z-координата
b = np.loadtxt('box_1.xvg')
t = b[:,0]
x = b[:,1]
y = b[:,2]
z = b[:,3]
Определим нормаль к поверхности белка. Наблюдения из файла box_1.xvg - вариации координаты Х почти нет, в то время как Y и Z изменяются => Х будет нормалью к поверхности.
Соответственно, площадь считаем как произведение по осям Y и Z и нормируем на 32 (число липидов в слое)
S_yz = y*z/32
plt.plot(t,S_yz, label = 'Square of Y and Z')
plt.title('Square and Time Dependency')
plt.legend()
Мы видим явно нисходящий тренд: со временем площадь липида уменьшается. Кажется, что это разумно - в процессе моделирования структура компактно собирается, приходя в оптимальное состояние. Также отмечу, что на прошлом шаге я отметила время t = 27000 как момент образования бислоя. Судя по графику, я почти права.
Изменение гидрофобной и гидрофильной поверхности¶
g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
Опять неизящное изменение файла руками (удаление первых 22 строк), чтобы из нулевого столбца вытащить время, из первого - инфорацию о гидрофобных, из второго - о гидрофильных.
s = np.loadtxt('sas_b.xvg')
t_1 = s[:,0]
phob = s[:,1]
phil = s[:,2]
plt.plot(t_1,phob, label = 'Hydrophobic surface')
plt.plot(t_1,phil, label = 'Hydrophilic surface', color ='red')
plt.title('Surface Changes within Time')
plt.legend()
И здесь тоже видим нисходящий тренд: размер поверхностей уменьшается. Уменьшается энергия системы, что приводит к образованию фосфолипидного бислоя.
Мера порядка¶
%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
# xaxis label "Atom"
# yaxis label "S"
ord_end = np.loadtxt('ord_end.xvg')
atom = ord_end[:,0]
end = ord_end[:,3]
ord_start = np.loadtxt('ord_start.xvg')
start = ord_start[:,3]
plt.plot(atom,end, label = 'End')
plt.plot(atom,start, label = 'Start', color ='red')
plt.title('Order Parameter')
plt.legend()
В начале траектории все значения выше; к концу - ниже. Это значит, что молекулы в начале подвижнее, а в конце - менее, что логично, т.к. они образовали фосфолипидный бислой.
import subprocess
converted = subprocess.call(["jupyter-nbconvert", '--to', 'html',"SB. HW9. Tevanyan.ipynb"], shell=False)
if converted==0:
print 'Your ipynb-file was sucessfully converted!'
else:
print 'Smth went wrong. for instance, check the filename...'
show_link = '<a href="%s" target="_blank">You can download it here!</a>'%('SB. HW9. Tevanyan.html')
display(HTML(show_link))