• 2019-07-31 10:14:58

缘起

很多年后, 当我再次更新这个工具, 就会记起, 很多年前, 我开始研究GROMACS中MMPBSA实现方法的情境. 那时的我就已经发现已有了两种解决方式: GMXPBSA脚本与g_mmpbsa程序. 我粗略地测试一下, 还算可用, 但未能尽美, 可改进之处也颇多. 因为我不经常使用, 便也就无所谓了, 暴虎冯河, 用过即忘, 倒也合适.

时间流逝, 形势在变, 心境也在变. 当我觉得真正需要测试一下MMPBSA方法时, 又回想起这两种解决方式, 重新捡拾起先前的操作, 才发现颇不顺手. GMXPBSA对GROMACS版本有要求, 脚本也颇得意大利人的风骨, 一如Italic的中文, 让我厌烦. 出起错来却又如意大利面, 纠纠结结, 纠缠不清, 让人无可措手, 只想烈马快刀, 一斩了之. g_mmpbsa呢, 印度人的作品, 传至东土已有些时日, 受者众多, 却仍有水土不服之虞. 安装麻烦, 二进制包Windows下也无法使用. 外加只能使用特定版本的GROMACS和APBS, 而无法兼容最新版本. 我本想好好参阅代码真经, 希冀修饬一二, 传与他人, 也算得我功德一件. 可总不得五蕴皆空, 六根清净. 既已如此, 又何须执著? 倒不如我也来写一段, 看看传得传不得.

简介

gmx_mmpbsa用于计算GROMACS轨迹的MMPBSA结合能, 并进行能量分解. 计算时, MM部分由脚本自行完成, PBSA部分借助APBS程序完成. 因此, 在使用gmx_mmpbsa前, 必需安装好GROMACS和APBS程序. 如果你已经需要使用这个脚本了, 那说明你已经得到了GROMACS的轨迹, GROMACS自然是已经安装好了的. 而APBS的安装也很简单, 官方网站提供了二进制版本, Linux和Windows下的都有, 拖过来就可以用.

gmx_mmpbsa计算MM和PBSA所需的原子参数(电荷, 半径, LJ参数等)来自GROMACS的.tpr文件, 因此, 需要调用gmx dump程序处理.tpr文件, 以便抽取所需的参数. 采用这种方法, 可以避免版本不兼容的问题, 因此也就可以支持任意版本的GROMACS. 对于APBS也是一样, 只要APBS的输入文件格式没有变化, 那gmx_mmpbsa生成的APBS输入文件就可以用于任意版本的APBS. 这样我们就可以安装GROMACS和APBS的最新版本而不用担心兼容性问题.

gmx_mmpbsa还使用了gmx trjconv程序, 用于处理轨迹的周期性叠合问题. 当然, 这一步你可以自己完成, 而无需借助脚本, 特别是在需要对周期性进行特殊处理的情况下.

在Windows下使用时, 还需要一个bash环境, 你可以使用msys2, Git Bash, 或cygwin.

使用流程

1. 获得输入文件

1. 生成tpr文件, 预平衡

2. 运行成品模拟, 得到xtc轨迹文件

3. 生成ndx文件, 其中需要定义三个组: 复合物(com), 蛋白(pro), 配体(lig). 名字虽然是pro, lig, 但其实可以代表任意分子, 比如两个有机小分子, 而不一定就是蛋白和配体.

4. 处理轨迹. 如果存在二聚体, 团簇等情况, 确保组成的原子间没有分离. 其他情况脚本可自动处理.

2. 设定gmx_mmpbsa计算参数

程序路径

• apbs='c:/apbs1.5/bin/apbs.exe' # APBS: APBS可执行文件的完整路径. Linux系统下如: apbs='/home/users/APBS/bin/apbs'.

基本参数

• ff=AMBER: 力场类型, AMBER, OPLS, CHARMM

• trj=1EBZ.xtc: 轨迹文件

• tpr=1EBZ.tpr: tpr文件

• ndx=index.ndx: 索引文件

• com=System: 复合物索引组

• pro=Protein: 蛋白索引组

• lig=BEC: 配体索引组

将trj, tpr, ndx等变量指向相应文件, 将com, pro, lig指向索引文件中相应的索引组.

其他参数主要是APBS计算所用的极性参数, 非极性参数, 网格参数, 一般无需更改太多.

3. 运行脚本

1. 脚本首先处理轨迹: 1. 完整化; 2. 居中叠合. 之后将构型输出到pdb文件.

2. 脚本其次抽取tpr中的原子信息, 存放在qrv文件中. 主要是复合物中每个原子的电荷, 半径, LJ参数以及残基信息.

3. 脚本根据pdb文件中原子的坐标获取APBS网格参数, 并将每帧构型输出到APBS所需的pqr文件, 同时生成APBS输入文件*.apbs. 然后调用APBS计算每帧构型对应的apbs文件, 并计算极性PB, 非极性SA部分的贡献, 再计算MM贡献, 同时进行残基分解, 输出结果.

由于脚本在计算时是分步进行的, 因此你可以先将所有apbs文件都产生出来然后并行计算它们, 最后再统计结果. 当然这只有在计算构型较多的情况下才值得尝试.

文件说明

用户反馈

与GMXPBSA相比

• GMXPBSA只能给出总的结果, 而gmx_mmpbsa可给出MM, PB, SA总的结果和残基分解结果;

• GMXPBSA使用时若体系中存在自定义拓扑或残基, 则需要提供相应的拓扑或力场文件;

• GMXPBSA原有脚本生成并使用的Mm.mdp文件较老, 需要手动修改print_files.dat脚本;

• GMXPBSA出错时需要进入子目录查看错误原因, 重新运行时又要删除子目录, 否则不能运行, 因而调试较为繁琐;

• GMXPBSA需先后运行三个脚本, 而gmx_mmpbsa只需要运行一次;

• GMXPBSA可实现甘氨酸扫描(CAS), gmx_mmpbsa尚未实现.

与g_mmpbsa相比

• g_mmpbsa需要安装, 而gmx_mmpbsa无需安装, 简单修改后可直接使用;

• g_mmpbsa需手动分别运行计算MM, PB, SA的能量, 而gmx_mmpbsa只需要运行一次.

• g_mmpbsa计算速度较慢, 而gmx_mmpbsa脚本计算速度较快, 因可以使用最新版本的APBS程序, 支持并行.

• g_mmpbsa可实现任意残基突变, gmx_mmpbsa尚未实现.

• g_mmpbsa支持其他非极性模型, 如SAV, WCA, gmx_mmpbsa尚未实现.

有待改进

• 结果的误差统计

• GB方法

• 其他非极性模型

• 调用其他PB程序, 如AMBER的PBSA

• 熵的贡献

• 形变能的贡献

• 甘氨酸扫描, 任意残基突变

• APBS mg-para方法处理超大蛋白

测试: g_mmpbsa的1EBZ示例

MM能量

计算方法: 直接对势累加, 考虑溶质(蛋白)的介电常数pdie, 不使用周期性, 不使用截断, 不对构型进行优化.

注意: GMXPBSA未考虑pdie, 且使用周期性, 使用截断, 对构型进行了优化, 但所得结果差距不大.

MM能量分解

计算方法: 直接累加

PB(极性)能量, Polar solvation enrgy

计算方法: 水相能量减去真空中的能量, PB=Esol-Evac

与g_mmpbsa比较, PB能量有差距, 原因在于

• 所用网格不同, g_mmpbsa对每个分子构型使用基于其自身的网格, 而gmx_mmpbsa使用基于所有构型的网格, 所有APBS计算使用相同的网格.

• 有些原子的半径不同, g_mmpbsa根据蛋白原子名称推断半径的方式有些怪异, 与论文所给数据不一致. 对配体, 由于原子名称可以任意取, 所以没有办法直接根据名称推断半径. 下面是g_mmpbsa所给蛋白中原子的半径, 可以看到有很多不一致之处.

    CA      1.7
    CB      1.7
    CD      1.7
    CD1     1.7
    CD1     1.77
    CD2     1.7
    CD2     1.77
    CE      1.7
    CE1     1.77
    CE2     1.77
    CE3     1.77
    CG      1.7
    CG      1.77
    CG1     1.7
    CG2     1.7
    CH2     1.77
    CZ      1.7
    CZ      1.77
    CZ2     1.77
    CZ3     1.77
    H       1.2
    H1      1.2
    H2      1.2
    HA      1.2
    HA1     1.2
    HA2     1.2
    HB      1.2
    HB1     1.2
    HB2     1.2
    HB3     1.2
    HD1     1.
    HD1     1.2
    HD11    1.2
    HD12    1.2
    HD13    1.2
    HD2     1.
    HD2     1.2
    HD21    1.2
    HD22    1.2
    HD23    1.2
    HD3     1.2
    HE      1.2
    HE1     1.
    HE1     1.2
    HE2     1.
    HE2     1.2
    HE21    1.2
    HE22    1.2
    HE3     1.
    HE3     1.2
    HG      1.2
    HG1     1.2
    HG11    1.2
    HG12    1.2
    HG13    1.2
    HG2     1.2
    HG21    1.2
    HG22    1.2
    HG23    1.2
    HH      1.2
    HH11    1.2
    HH12    1.2
    HH2     1.
    HH21    1.2
    HH22    1.2
    HZ      1.
    HZ1     1.2
    HZ2     1.
    HZ2     1.2
    HZ3     1.
    HZ3     1.2

• 换用与g_mmpbsa相同的网格与半径后, 所得PB作用能差距很小.

• GMXPBSA所用网格更大, 且使用了来自LJ参数的半径, 并对构型进行了优化, 所得结果也更大.

PB能量分解

计算方法: 残基处于复合物中时的能量减去处于其所属分子(蛋白/配体)中的能量, G(res@com)-G(res@pro)

SA(非极性)能量, APolar solvation energy

计算方法: 每个原子的溶剂可及表面积乘以表面张力, 加上常量

SA能量分解

计算方法: 直接累加

能量分解图示