GROMACS是一个功能强大的分子动力学的模拟软件，我们在前面的推文中为大家分享了一些使用GROMACS运行分子动力学的教程，上期内容也为大家分享了分子动力学之后的一些常用的分析方法，本期内容继续带大家看看在GROMACS中，还有哪些常用的分析方法。
1.导出固定的帧
很多时候我们需要导出分子动力学轨迹中某一时刻的结构构象，比如比较动力学前后的构象差异，就需要导出轨迹中的第一帧和最后一帧。在Gromacs中，导出固定帧十分容易，常用命令如下：
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o first.pdb -dump 0
以上命令可输入轨迹中的第一帧，即0 ps时的结构，通过“-dump”可指定时间，单位为ps。
2.间隔固定时间导出帧
除了单帧导出外，可能还需要间隔帧导出，比如制作分子动力学轨迹动画时，直接使用原轨迹制作，轨迹文件较大，容易造成卡顿，就可以间隔一定的时间导出帧，比如每隔100 ps导出一帧到PDB文件，可输入以下命令：
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o traj.pdb -dt 100
以上命令导出的全部帧都保存为一个PDB文件，如果需要将每一帧都输出为独立的PDB文件，可以通过“-sep”命令来指定。
3.计算平均结构
“gmx covar”主要用来计算并对角化（质量加权的）结构协方差矩阵，也可以用来计算平均结构。
gmx covar -s md_0_10.tpr -f md_0_10_center.xtc -b 500 -e 800 -av traj_aveg.pdb
以上命令可以得到分子动力学模拟轨迹中500到800 ps的分子平均结构，结果输出为traj_aveg.pdb。
此外，上期内容为大家介绍了使用“gmx rmsf”进行RMSF分析，计算分子结构的均方根涨落，也可使用该命令获取某一时间段的平均结构，如输入以下命令获得分子动力学模拟轨迹中500到800 ps的分子平均结构：
gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -b 500 -e 800 -ox traj_aveg.pdb
4.叠合两结构并输出RMSD
“gmx confrms”可将两个分子结构进行叠合，并且计算RMSD值，值得一提的是，叠合的两个结构的原子数可以不相同，只要用于叠合的两个索引组一样即可。常用命令为：
gmx confrms -f1 structure_1.gro -f2 structure_2.gro -o fit.pdb -n1 index_1.ndx -n2 index_2.ndx
通过以上命令，叠合的两个结构会写入至“fit.pdb”文件中。
5.分子动力学轨迹聚类
聚类（Clustering） 是按照某个特定标准（如距离）把一个数据集分割成不同的类或簇，使得同一个簇内的数据对象的相似性尽可能大，同时不在同一个簇中的数据对象的差异性也尽可能地大。
在分子动力学中，可以依据RMSD对轨迹进行聚类，把相似的构象聚集在一起，聚类的常用命令为：
gmx cluster -f md_0_10_center.xtc -s md_0_10.tpr -cutoff 0.8 -b 1000 -wcl 10 -method gromos -sz
其中：
-b：指定从多少ps开始分析，一般情况下轨迹前期结构变化较大，可以从稳定后如1000 ps开始分析；
-method：指定算法，默认为linkage，可以选jarvis-patrick、monte-carlo、diagonalization、gromos；
-cutoff：截断值，可以根据实际情况需要反复尝试，cutoff设置的越大，簇里的帧之间的RMSD允许差异就越大，聚类的数目越少。
通过以上命令，每个簇当中最有代表性的结构一起被输出为多帧文件cluster.pdb，可以通过各类可视化软件如VMD、PyMOL等进行观看。
6.氢键分析
氢键在分子动力学中相互作用的分析中十分常用，在Gromacs中可通过“gmx hond”命令对氢键的数量、稳定性等进行分析。在氢键分析中需要指定两个不同的组，它们必须完全相同或者彼此之间无任何重叠。
分析氢键常用的命令为：
gmx hbond -s md_0_10.tpr -f md_0_10_center.xtc -n index.ndx -num hbond_num.xvg -dt 10 -life hbond_life.xvg -ac hbond_ac.xvg
通过以上命令可输出hbond_ac.xvg、hbond_life.xvg、hbond_num.xvg文件。其中：
hbond_ac.xvg：输出了氢键的平均存在周期（forward lifetime），可以作为衡量氢键稳定性的一个指标；
hbond_life.xvg：主要需要关注p（t）；p（t）是t的概率密度函数，就是统计不同的t出现的频率得到的结果，即氢键从0能维持到t时刻的概率；
此外，t*p（t）对时间积分，用统计学语言说这个值就是时间的期望值，也就是平均寿命；
hbond_num.xvg：第二列表示的是氢键数目，而第三列显示的是Pairs within 0.35nm，表示距离在0.35nm不考虑氢键键角判据的氢键数，实际意义不大。
需要注意的是，hbond程序使用几何准则来对氢键加以判定，当氢键供受体距离小于3.5埃且氢-供体-受体所成角度小于30度时，即认为其为一个氢键。而不同软件如VMD、PyMOL在默认情况下对于氢键的判定与Gromacs有所差异，所以如果使用不同的软件进行氢键分析，会获得不同的结果。
相关链接：
1. https：//www.gromacs.org/
2. https：//manual.gromacs.org/current/index.html
转自：“叮当学术”微信公众号
如有侵权，请联系本站删除！
