MMPBSA计算工具使用

g_mmpbsa的使用教程

(2020.01) 新增g_mmpbsa Win版本(调用外部apbs程序)
(2020.02) 新增g_mmpbsa_with_intra_apbs(内置apbs1.4并行,不再需要额外安装apbs) Win版本
(2020.02) 新增g_mmpbsa(内置apbs1.4) Linux版本
(2020.03) 新增g_mmpbsa(内置apbs1.4) Win版本,直接支持到最新GROMACS2019.6
(2020.05) 新增g_mmpbsa 2019.6(内置apbs1.4) Win版本,采用sse2指令集,适合旧电脑
(2020.02) 功能新增:支持单个残基MM能量分解输出为静电能和范德华能
(2020.02) 功能新增:提供两个bash脚本,与官网Python功能一致,其中awk_mmpbsadecomp.bsh支持输出残基静电能和范德华能。具体选项采用-h查看
(2020.04) 功能新增:新增awk_sum_energy2xvg.bsh脚本,所得xvg文件可直接使用xmgrace/qtgrace进行绘图显示
(2020.03) 新增修改版GMXPBSAtool

说明:
  如上面所述,最近更新了一些工具。为什么要更新呢?我简单的讲一下:其一: 由于g_mmpbsa的作者不再更新该工具,因此有些用户喜欢的功能无法加入,最常见的就是残基能量进一步分解为静电能和范德华能。其实这个功能说简单也特别简单,说难也可能对于编程一窍不通的人来说有点儿,这个能量分解直接在代码中添加输出就可以了。因此上面的更新中默认也会产生分解能量的文件。其二: 平台,我一直认为Linux和Windows平台各有特色,因此也致力于把一些某一平台特有的工具进行共有化。如上面说的,原作者并没有给出Windows平台g_mmpbsa,其实也是因为该软件依赖的库比较多,而在Windows平台又没有现成的,比如gromacs.lib,就需要通过VS编译获得,比较麻烦。因此上面给出了两种平台下编译好的二进制g_mmpbsa,注意都采用的AVX2_256指令集(我觉得现在的CPU都应该具备了吧,这也是优势)。其三: 兼容性。什么软件都应该与时俱进,由于本人能力也比较有限,不可能对代码进行大量更改,此次更新也包括了GROMACS5.X版本的最后一个版本[已更新到GROMACS2019.6],我觉得兼容能力(对不同版本的tpr读取)应该是更强的。另外也特别提供了awk版本的后处理工具,支持残基分解能的计算。


废话不多说,因为网上有许多教程写过g_mmpbsa的使用教程,具体见最下方参考资料[1-2]。此处我所要讲的是g_mmpbsa网站自带的一个做氨基酸突变的python脚本。因为官网没有给出具体的使用实例,所以导致很多人可能不知道或者也没有使用过,我也侥幸使用了一番,还挺有意思,因此记下使用过程。

脚本功能介绍

根据脚本中所提及的,该脚本的功能主要有两个:其一,丙氨酸扫描;其二,氨基酸突变

脚本使用方式

python环境配置

py脚本当然是使用python。此脚本使用python2写的,因此你必须安装python2,这里我建议用anaconda环境,安装和环境切换十分的方便。初次使用python mutate_traj_tpr.py -h查看帮助信息会出错,提示from modeller import *出错,不难看出这个脚本需要调用modeller来处理氨基酸突变。解决方法也十分简单,安装modeller,官网可以通过学术机构邮箱获得秘钥,安装方式我就不讲了,网上有很多。安装完成以后。再次使用python mutate_traj_tpr.py -h可以看到下图所示:

一眼就可以看出上半部分是设置氨基酸突变/丙氨酸扫描用的参数,下半部分用过GROMACS软件的人一眼就能看出是能量极小化的mdp参数部分。

辅助软件

因为该脚本是通过modeller将突变完成以后的蛋白直接进行能量最小化以备GROMACS软件模拟用,因此有几点要注意的地方:

  • 蛋白体系中不能含有其他配体或者辅酶等物质,因为GROMACS的pdb2gmx无法直接处理这样的体系。
  • 本文使用的是GROMACS 4.X版本,与脚本中出现的命令保持一致。如果你要使用最新的GROMACS软件,就需要将所有使用4.x版本的命令加上gmx前缀。另外,其中有一个4.x版本的trjcat命令可能需要删除,因为如果你直接换成新版本的gmx trjcat命令,会报错对于-f选项无法使用.gro文件.

具体实例

成功安装和配置好环境以后你就可以使用该脚本做自己的研究了。在分子动力学模拟或者生化领域,做氨基酸突变主要是为了考察该氨基酸的功能,并且对突变以后的蛋白结构和功能造成怎样的影响。本文以1AKI(鸡蛋白溶菌酶)蛋白为例,考察定点突变。
比如,我想考察一号残基位的赖氨酸(LYS)突变成半胱氨酸(CYS)以后对整个蛋白体系的结构稳定性和功能的影响。那么我们可以输入以下命令:

1
python mutate_traj_tpr.py -drWT C:\Users\liuyujie714\Desktop\test -drMT C:\Users\liuyujie714\Desktop\MT -rr "A:1"  -rnm CYS  -ff 6

这里,几个参数需要解释一下:
-drWT指的是你蛋白pdb文件所在的文件夹位置。
-drMT指的是你需要存储输出文件的文件夹位置,即突变以后的氨基酸文件。
-rr指定的是需要突变的氨基酸编号,其中 "A:1"表示A链的1号氨基酸,本例中指的是赖氨酸,如果含有其他链可以相应更改。
-rnm指定的是突变成哪种氨基酸,其中CYS指的是半胱氨酸。
-ff指定GROMACS力场。此处6表示amber99sb-ildn.ff力场。
其余参数保持默认,包括em.mdp中的参数。
好了,等待完成即可,屏幕输出如下图:

这样你会在-drMT指定的文件夹下得到如下文件:

其中,1aki_em.gro为突变后经过能量极小化以后的结构文件,你可以用它进行下一步的动力学模拟。

结构检查

过程虽然很顺利,但是还是应该检查一下突变的位置是否正确,通过可视化可以看出突变前后的变化,如下图:

看起来还是正确的突变。。
先记录到这里,以后有添加再补充。

补充(2020.02)

详见置顶说明,下载不同平台更新的g_mmpbsa:
Windowsg_mmpbsa_win
Linuxg_mmpbsa_linux

注意:对于Windows版本,我自己在Win7Win1064位系统上都测试过没问题;Linux版本仅限于Ubuntu16.04

补充(2020.03)

将GMXPBSAtool更改为适用于gmx版本>5.1,修改输出文件排版(原排版真是辣眼睛)。当然也修复了Windows下使用的一些bug。下载地址如下,包括了1EBZ实例。
gmxpbsatool_modified_by_liuyujie

补充(2020.05)

gmx_mmpbsa工具是李继存老师写的,也比较方便,需要提前安装好gmxapbs,测试对于新版apbs3.0也是没问题的。
值得注意的是如果新手运行后发现PBSA项均为0,那么很可能是因为你没有正确安装apbs或者脚本里面的路径设置不正确。因为这个脚本里面没有检查程序纠错,我就在原脚本中加了几行来确定gmxapbs是否正常工作的异常处理报错,下载链接

参考资料

[1] g_mmpbsa 安装与使用笔记
[2] GROMACS教程:使用GROMACS计算MM-PBSA结合自由能