Gromacs——教程学习(1)

1. 分子动力学模拟（Molecular Dynamics）全流程

所有的xvg格式文件，都可以使用大神编写的python DuIvyTools脚本可视化，很方便，只要你的电脑配置了python或者anaconda或者miniconda

pip install DuIvyTools

dit xvg_show -f xxx.xvg

gmx mdrun -deffnm md_0_1 -nb gpu
#-nb gpu调用GPU进行计算，可加速

B站“杜艾维”老师的课程

2. Gromacs蛋白质动力学模拟入门简明步骤-v3

#gmx energy模块来分析能量最小化的结果
gmx energy -f em.edr -o potential.xvg
#提示时, 输入10 0来选择势能Potential(10), 并用零(0)来结束输入

#图形化查看能量最小化的势能变化
xmgrace potential.xvg

NVT温度平衡

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
#平衡结果分析
gmx energy -f nvt.edr -o temperature.xvg
#在提示符下键入“16 0”，选择温度代号16，选择0退出
#图形化查看温度平衡情况
xmgrace temperature.xvg

NPT压力平衡

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
gmx energy -f npt.edr -o pressure.xvg
#根据提示输入18 0 获得压力的数据结果
xmgrace pressure.xvg

密度分析

gmx energy -f npt.edr -o density.xvg
#根据提示 输入 24 0 获得结果
#图形化查看
xmgrace density.xvg

正式的动力学模拟

#执行grompp生成模拟的tpr文件：
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
#执行模拟
gmx mdrun -deffnm md_0_1

分子动力学模拟结果分析

#选择1 (Select 1 ("Protein") as the group to be centered and 0 ("System") for output)
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
#Choose 4 ("Backbone") for both the least-squares fit and the group for RMSD calculation
gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
#计算相对于模拟之前晶体的结构差异，可以使用下面的命令, 也是选择4 ("Backbone")：
gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns
#可以把两个RMSD拟合一下查看：
xmgrace rmsd.xvg rmsd_xtal.xvg
#计算回旋半径：
#选择group 1 (Protein)，执行完成后新生成了gyrate.xvg
gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg
#图形化查看
xmgrace gyrate.xvg

3. 分子动力学模拟-gromacs的基本使用

#gromacs推荐使用conda安装:
conda install -c bioconda gromacs

(我没有用conda装过，用的是集群统一安装的，不过一般是源码安装吧)

对于盒子的设置，本例使用的是正方体，但是正十二面体的菱形可能是一个更好的选择（能够使用更少的水分子作为填充）

之前的溶剂和溶质的优化是分割的，所以我们还要使得溶剂和溶质同时优化到我们想要的状态（温度，压力）
温度平衡，NVT系综
压力平衡，NPT系综
MD模拟
结果分析

4. 分子动力学模拟（Molecular Dynamics）全流程【二】

2、模拟结果分析：
参考诸多教程和文献，关于分子动力学模拟的结果分析与展示，大多数是观察RMSD稳定状态、计算RMSF、计算回旋半径；并基于RMSD和回旋半径，绘制自由能景观形貌图；

2.1 蛋白定位至中心
在计算RMSD以及回旋半径过程，大分子（蛋白）以及小分子会在盒子中不断变换位置和构象，甚至会出现跨盒子的现象；所以在整个计算开始之前，我们需要将整个体系定位到盒子中心；
使用gmx trjconv命令完成操作：

#选择蛋白质，1；1
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center

选择蛋白质，1
gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o center.pdb -dt 1000

运行后得到md_0_1_noPBC.xtc文件，这是居中后的体系文件，用于后续分析，这里我们可以可视化一下刚刚生成的center.pdb文件（pymol），这里记录蛋白运动1000帧的动画，时长00:10

2.2 去除平衡转动
上面的动画我们发现整个蛋白复合物运动幅度都比较大，所以下一步需要去除体系中的平衡转动；

#选择蛋白骨架，4，选择蛋白，1
gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o fit.xtc -fit rot+trans

#选择蛋白质，1
gmx trjconv -s md_0_1.tpr -f fit.xtc -o fit.pdb -dt 1000

去除平衡转动后，我们再继续可视化一下fit.pdb，并使用fit.xtc进行后续的计算时长00:11

2.3 计算RMSD、RMSF、回旋半径
使用fit.xtc作为输入文件，进行下一步的计算：

gmx rms -s md_0_1.tpr -f fit.xtc -o rmsd.xvg -tu ns

选择蛋白骨架，4，随后选择蛋白，1

dit xvg_show -f rmsd.xvg

从RMSD的结果来看，整个体系在7ns-10ns基本保持稳定，如果计算步长更长，时间更久，可能之后一直趋于稳定（个人观点，请批评指正！）；
我们可以使用模拟前的tpr文件同样做一次计算，对比结果：

gmx rms -s npt.tpr -f fit.xtc -o rmsd_xtal.xvg -tu ns
选择蛋白骨架，4，随后选择蛋白，1
dit xvg_compare -f rmsd.xvg rmsd_xtal.xvg -c 1 1

接下来计算RMSF，使用gmx rmsf命令完成：
计算RMSF

gmx rmsf -s md_0_1.tpr -f fit.xtc -o rmsf.xvg
dit xvg_show -f rmsf.xvg

计算回旋半径，使用gmx gyrate进行计算：

gmx gyrate -s md_0_1.tpr -f fit.xtc -o gyrate.xvg
选择蛋白，1
dit xvg_show -f gyrate.xvg

2.4 绘制自由能景观形貌图
因为这是基于RMSD和回旋半径绘制的，所以在此之前，我们需要将上述得到的结果（rmsd.xvg、gyrate.xvg进行合并），先观察数据，直接用文本编辑器打开
之后，使用gmx sham命令生成图谱：

gmx sham -tsham 310 -nlevels 100 -f output.xvg -ls gibbs.xpm -g gibbs.log -lsh enthalpy.xpm -lss entropy.xpm
dit xpm_show -f gibbs.xpm
python xpm2png.py -f gibbs.xpm -o gibbs.png -ip yes
#观察谱图

2.5 通过主成分绘制自由能形貌图
同样的使用主成分分析，绘制自由能形貌图同样需要去除体系的平衡转动：
协方差计算（这里不太懂）：

gmx covar -s md_0_1.tpr -f fit.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covapic.xpm
python xpm2png.py -f covapic.xpm -o covapic.png -ip yes

接下来进行主成分PCA分析，绘制形貌图并可视化：

##生产主成分数据
gmx anaeig -s md_0_1.tpr -f fit.xtc -v eigenvectors.trr -first 1 -last 1 -proj pc1.xvg
gmx anaeig -s md_0_1.tpr -f fit.xtc -v eigenvectors.trr -first 2 -last 2 -proj pc2.xvg

注意计算完PCA,同样需要对两个结果进行组合，重命名为pc_output.xvg
计算自由能形貌图

gmx sham -tsham 310 -nlevels 100 -f pc_output.xvg -ls pca12_gibbs.xpm -g pca12_gibbs.log -lsh pca12_enthalpy.xpm -lss pca12_entropy.xpm
python xpm2png.py -f pca12_gibbs.xpm -o pca12_gibbs.png -ip yes

从形貌图我们可以看到整个体系能量最小的构象（即深蓝色），通过下面两个文件（bindex.ndx、pca12_gibbs.log）寻找最低能量帧数……
确认帧数之后，可以通过gmx trjconv导出pdb构象

gmx trjconv -s md_0_1.tpr -f fit.xtc -b 30 -e 30 -o test.pdb

当然部分文章也会进行氢键分析，代码如下：

#氢键分析
gmx hbond -s md_0_1.tpr -f fit.xtc -n bindex.ndx -num hbond_num.xvg -dt 100 -life hbond_life.xvg -ac hbond_ac.xvg
gmx hbond -f fit.xtc -s md_0_1.tpr -num hb_num.xvg -ac hbond_ac.xvg -dist hbdist.xvg -hx hbhelix.xvg