- 2017-12-14 07:45:45
在以前整理的GROMACS分子动力学模拟教程: 多肽-蛋白相互作用中, 提到了一种使用gmx sham计算自由能形貌图的简单方法. 这种方法适用于计算那种和两个变量相关的自由能形貌图, 并且容易作图. 有些文章中给出的就是这种图, 比如很多蛋白的自由能形貌图中给出了自由能与回旋半径以及RMSD的关系, 并以彩色填充图的形式给出, 看起来很高端的样子. 实际上做这种图还是比较简单的, 下面我就给个示例.
模拟
我们使用一个28肽的小蛋白, BBA折叠多肽1fme来做示例:
作为示例, 我们只跑真空中的模拟, 100万步, 每步都保存到xtc文件. 这里要注意的是, 计算自由能时要尽量多保存轨迹, 越多越好, 因为构象越多, 得到的自由能形貌图越光滑.
计算相关性质
得到轨迹之后, 我们来计算回旋半径和RMSD(根据自己的需要选择分组):
gmx gyrate -f traj_comp.xtc
gmx rms -f traj_comp.xtc
这样就得到了gyrate.xvg, rmsd.xvg, 它们包含了每一步的总回旋半径Rg和RMSD值.
准备计算自由能的输入文件
要使用gmx sham计算自由能与Rg和RMSD的关系, 我们需要为它准备一个输入文件, 文件名称设为Rg-RMSD.xvg, 其中列出每一步对应的时间, Rg, RMSD值. 准备这种文件的简单方法是使用列编辑功能, 直接切出几列, 放到同一个文件中即可. 但如果每列的数据行数过多的话(在我们的例子中有100万行), 许多文本编辑器的列编辑功能可能会崩溃. 所以, 更稳定的方法是使用bash的paste命令, 它可以将两个文件按行合并在一起. 注意使用这个命令时, 文件的换行符应该是Linux下的换行符.
先将gyrate.xvg和rmsd.xvg文件中前面的注释和设置部分删掉, 只留下数据部分, 然后执行下面的命令
paste gyrate.xvg rmsd.xvg > Rg-RMSD.xvg
这样就得到了合并好的文件. 不过, 文件中还含有我们不需要的数据, X/Y/Z三个方向的Rg值和RMSD对应的时间值, 将这些列直接使用列编辑功能删除就好了.
有了输入文件后, 就可以计算自由能了:
gmx sham -f Rg-RMSD.xvg -ls -b 20 -ngrid 200 200 200 -nlevels 100
这里指定从20 ps开始计算, 是因为从RMSD曲线很容易看出, 20 ps后RMSD才稳定. 如果不忽略前面未平滑部分, 计算出的自由能范围可能过大, 对后面的作图不利. 后面的选项根据你自己的测试来设置就好.
作图
得到的gibbs.xpm文件就包含了自由能与Rg和RMSD的关系, 用它就可以作图了. 如果是作二维图的话, gmx xpm2ps就可以了. 要作三维填色图的话, 麻烦点, 需要使用Origin, MatLab之类的作图软件. 但这些软件在作图的时候无法将自由能形貌图和蛋白一起渲染. 要达到这个目的, 我们可以使用PyMOL.
我以前曾经示例过如何利用PyMOL制作反应势能面示意动画, 将那里的脚本稍加修改就能满足我们的需要. 不过在作图之前, 我们需要将xpm中的数据提取出来. 这可以使用前面提到的那篇教程中提供的一个脚本, 但那个脚本只适用于水平数低于80的情况. 我们这里使用的水平数为100, 所以需要对脚本稍加修改. 此外, 如果得到的图形尖峰太多, 可以对数据进行简单的平滑, 这样图形看起来更漂亮.
然后再载入蛋白, 调整好位置, 将其与形貌图放在一起渲染.
看, 蛋白在形貌图上有投影的.
当然, 你还可以将其他有代表性的构象放在自由能形貌图的不同位置, 这样信息量就更大了. 其他的很多设置都可以试着调整. 这些就不赘述了.
不过这种作图方法的一个缺点就是不好加坐标轴, 需要额外的代码才可以.
