Virtual Structural Biology Group » Учебные курсы » Структурная биоинформатика в ВШЭ

h1. Метадинамика

*Создайте рабочую директорию.*

* Координаты белка возьмите из записи PDB:1L2Y .

* файл настроек для минимизации энергии "em.mdp":https://gitlab.sirius-web.org/golovin/protein-design/-/raw/master/pracs/p5/em.mdp 

* файл настроек для "утряски" воды "pr.mdp":https://gitlab.sirius-web.org/golovin/protein-design/-/raw/master/pracs/p5/pr.mdp

* файл настроек для молекулярной динамики "md.mdp":https://gitlab.sirius-web.org/golovin/protein-design/-/raw/master/pracs/p5/md.mdp

* Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs.

<pre><code class="bash">

gmx pdb2gmx -f p.pdb -o p -p  -ff amber99sb -water tip3p

</code></pre>

* Делаем небольшой отступ в ячейке от белка.

<pre><code class="bash">

gmx editconf -f p.gro -o p_ec -d 1.5

</code></pre>

* Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты в молекуле.

<pre><code class="bash">

gmx grompp -f em -c p_ec -p  -o p_em -maxwarn 1

gmx mdrun -deffnm dna_p -v

</code></pre>

*Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Знаесите начальное и конечное значение максимальной силы.

* Добавим в ячейку молекулы воды

<pre><code class="bash">

gmx genbox -cp p_em -p  -cs  -o p_s

</code></pre>

* Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp обратите внимание на информацию о заряде системы.

<pre><code class="bash">

gmx grompp -f em -p  -c p_s -o p_s

gmx genion -s p_s -o p_si -p  -neutral -conc 0.15

</code></pre>

* Проведём "утряску" воды:

<pre><code class="bash">

gmx grompp -f pr -c p_si -p  -o p_pr -maxwarn 1

gmx mdrun -deffnm p_pr -v

</code></pre>

* Cравните визуально изменения в системах p_pr.gro и p_si.gro . Занесите наблюдения в отчёт. 

* Запускаем тестовое моделирование молекулярной динамики.

<pre><code class="bash">

gmx grompp -f md -c p_pr -p  -o p_md -maxwarn 1

gmx mdrun  -deffnm p_md -v -nsteps 1000000 -ntomp 4

</code></pre>

Просмотреть ход счёта можно в файле p_md.log

h3. Метадинамика 

* Попробуем странную штуку https://github.com/spiwokv/af2cv/tree/main

* Pickle для белка можно сделать в "ColabFold":https://github.com/sokrypton/ColabFold?tab=readme-ov-file с выбором AlphaFold2_mmseqs2 и обязательно поставить галочку Save All. Pickle будет в архиве.

* Вероятно, что при утряске воды система будет взрываться, надо в em.mdp  поменять integrator на steep. И потом:

<pre><code class="bash">

gmx grompp -f em -c p_si -p  -o p_em2 -maxwarn 1

gmx mdrun -deffnm p_em2 -v

gmx grompp -f pr -c p_em2 -p  -o p_pr -maxwarn 1

gmx mdrun -deffnm p_pr -v

</code></pre>

* Скопируйте себе  библиотеку для этой "CV":https://gitlab.sirius-web.org/golovin/protein-design/-/raw/master/pracs/p5/AFDistProb.cpp?inline=false и файл "plumed.dat":https://gitlab.sirius-web.org/golovin/protein-design/-/raw/master/pracs/p5/plumed.dat

* В начале файла надо внести изменения в путь к библиотеке с CV

>  Для варианта  дизайна этого пептида  надо догадаться как поменять номера атомов в двух местах: WHOLEMOLECULES и ATOMS=

* Запускаем 

<pre><code class="bash">

gmx  mdrun -deffnm p_md -ntomp 4 -v  -plumed plumed.dat

</code></pre>

h2.  Анализ метадинамики 

* В ходе расчета у вас получились файлы со структурами : trr и xtc

* Файл HILLS это информация об наложенных потенциалах, это тот самый "песок"

* Файл COLVAR это файл в который сохраняли значения CV и сумарное наложенный bias. 

* Посомотрите на результат в PyMol

<pre>

		gmx trjconv -f p_md -s p_md -o view1.pdb -skip 20 -pbc mol

		или

		trjconv -f p_md.xtc -s p_md.tpr -o view2.pdb -skip 20 -fit rot+trans

</pre>

* Воспользуйтесь plumed sum_hills для востановления поверхности свободной энергии

* Используя matplotlib и подобное постройте график

* Попробуйте придумать способ для оценки сходимости

* Давайте попробуем построить график FES  на основе перевзвешенных данных. 

<pre>

		 plumed driver --plumed rew.dat --noatoms

</pre>

* Постройте график и сравните с sum_hills

* Возможно это не самая лучшая переменная, давайте попробуем превзвесить на радиус инерции (предположительное распрямление белка)

* Расчитаем RG c помошью этих директив в отдельном dat файле 

<pre>

		rg: GYRATION TYPE=RADIUS ATOMS=1-284

		PRINT ARG=rg STRIDE=1 FILE=rg

</pre>

* Расчёт новой переменной будет сделан в новый файл

<pre>

		plumed driver --plumed rg.dat --mf_trr p_md.trr

</pre>

* Теперь надо скомбинировать данные meta.rbias из COLVAR и значения радиуса гирации в один новый файл формата COLVAR , обратите внимание на разницу в шаге. Вот некоторые функции в помощь

<pre>

		def GetPlumedColumns(f):

			with open(f,'r') as dat:

			arr = dat.readline().split()

			return arr[2:]

		def LoadNumpy(f,names):

			a = np.loadtxt(f)

			df = pd.DataFrame()

			for i,n in enumerate(names):

			df[n] = a[:,i]

			return copy.deepcopy(df)

		def combinevars(.......):

			c = LoadNumpy(file, cnames )

			gr = LoadNumpy(file, danames )

			out =  pd.merge(crf, da, on='time')

</pre>

* Расчитайте график свободной энергии в завимисомти от радиуса инерции, надо сделать новый rev.dat

* Попрбуйте сделать анализ  от двух переменных RG и alfascore

* Может есть силы поискать другие переменные которые могут лучше описывать процесс : https://www.plumed.org/doc-v2.9/user-doc/html/_colvar.html

* Для выбранной переменной сделайте анализ мутантного белка и ранее стабилизированного белка, данные лежат выше

* 

h2.  Подсказка 

> Общий пайплан построения FES : Подготовить COLVAR (слить, фильтровать) - Адаптировать rew.dat - Запустить plumed driver --plumed rew.dat --noatoms

* На hub запуск Plumed: 

<pre>

		/home/anur/miniconda3/envs/plumed/bin/plumed driver --plumed rew2d.dat --noatoms

</pre>

* Для работы дома надо подготовить все возможные COLVAR и положить их в свою директорию на hub