GROMACS氢键分析工具hbond的使用及扩展

类别:    标签: gmx bash   阅读次数:   版权: (CC) BY-NC-SA

前一篇文章介绍了VMD氢键插件的使用, 这篇文章就再接再厉, 说明下GROMACS氢键分析工具hbonds的使用, 顺便发布一个gmx_hbdat脚本, 使其能够整理氢键数据.

无论准备使用GMX的哪个工具, 首先要做的是查看其文档, 以及GMX手册中的相关说明, 确保能理解工具给出的结果. 对有些工具, 文档和手册中的说明太多简略, 有时候还必须查看相关的源代码. 这是使用GMX工具的必修课, 切记.

GMX的氢键工具hbond选项众多, 细究起来还是比较复杂的, 我们这里只介绍最基本的功能, 也就是分析一组原子内, 或两组原子间的氢键. 这需要指定两个索引组. 如果两个索引组完全相同, 那得到的就是这个索引组内部的氢键; 如果两个索引组完全不同(即两组不会同时包含相同的原子), 那得到的就是两个索引组之间的氢键. 这是目前仅支持的两种情况. 如果两个索引组既不完全相同, 却又部分相同, 那就无法计算, 结果当然也不可靠.

分析结果主要涉及三个文件,

如果查阅MD文章中的氢键分析部分, 可以发现很多文章给出的结果表格都是AMBER格式的, 其中会给出每个氢键所涉及的残基/原子, 每个氢键的占有率. 我们利用GMX的拓扑文件以及后两个文件, 容易整理得到这样的氢键数据. hbond.ndx中给出了原子编号, 根据编号查看gro文件中的信息, 或者拓扑中的数据, 就可以得到相关原子的残基, 原子类型, 名称等数据. 处理下hbmap.xpm文件就可以计算出每个氢键的占有率.

简单示例

继续以上篇文章的体系做示例吧. 这是一个包合物体系, 我们计算两个分子间的氢键. 先使用make_ndx做个索引文件, 其中包含两个分子各自的索引组. 然后运行命令

gmx hbond -f traj.xtc -s topol.tpr -n index.ndx -num -hbn -hbm

提示时, 选两个分子各自的索引组, 我用的是MOLA, MOLB.

运行完, 就得到了三个文件输出文件.

hbnum.xvg
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
@    title "Hydrogen Bonds"
@    xaxis  label "Time (ps)"
@    yaxis  label "Number"
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Hydrogen bonds"       # 同时满足距离标准, 角度标准的氢键数
@ s1 legend "Pairs within 0.35 nm" # 只满足距离标准的氢键数
         0           2           9
      0.02           2           8
      0.04           2           6
      0.06           2           6
      0.08           2           6
       0.1           2           6
      0.12           2           6
      ... ...
hbond.ndx
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
[ MOLA ]                  # MOLA的所有原子
    1     2     3     4     5     6     7     8     9    10    11    12    13    14    15
    ... ...
[ donors_hydrogens_MOLA ] # MOLA的供体原子, 氢原子
    8    9
   12   13
   22   23
   ... ...
[ acceptors_MOLA ]        # MOLA的受体原子
    3     8    12    16    22    24    29    33    37    43    45    50    54    58    64
    ... ...
[ MOLB ]                  # MOLB的所有原子
  148   149   150   151   152   153   154   155   156   157   158   159   160   161   162
    ... ...
[ donors_hydrogens_MOLB ] # MOLB的供体原子, 氢原子
  148  201
[ acceptors_MOLB ]        # MOLB的受体原子
  148   149
[ hbonds_MOLA-MOLB ]      # 整段分析轨迹中MOLA和MOLB之间所有可能的氢键
     64     65    149
    146    147    148
    148    201     43
hbmap.xpm
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
/* XPM */
/* This file can be converted to EPS by the GROMACS program xpm2ps */
/* title:   "Hydrogen Bond Existence Map" */
/* legend:  "Hydrogen Bonds" */
/* x-label: "Time (ps)" */
/* y-label: "Hydrogen Bond Index" */
/* type:    "Discrete" */
static char *gromacs_xpm[] = {
"101 3   2 1",
"   c #FFFFFF " /* "None"    */, # 无氢键, 空白
"o  c #FF0000 " /* "Present" */, # 有氢键, 红色o
/* x-axis:  0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16 1.18 1.2 1.22 1.24 1.26 1.28 1.3 1.32 1.34 1.36 1.38 1.4 1.42 1.44 1.46 1.48 1.5 1.52 1.54 1.56 1.58 */
/* x-axis:  1.6 1.62 1.64 1.66 1.68 1.7 1.72 1.74 1.76 1.78 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 */
/* y-axis:  0 1 2 */
# 三条氢键, 与前面`hbond.ndx`中的 [ hbonds_MOLA-MOLB ] 对应
"ooooooooooooooooo        o ooooo oo oooooooo oo ooooooooooooooooooooooooooooooooooooooooooooooooooooo",
"                 oooooo           ooo  oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo",
"ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo"

使用bash脚本hbdat.bsh获取类似AMBER的氢键数据

bash hbdat.bsh

得到hbdat.dat

hbdat.dat
1
2
3
4
#Donor             Hydrogen           Acceptor            Occupancy%
64#"MOLA"_MOL@OH   65#"MOLA"_MOL@HO   149#"MOLB"_MOL@O_3       100.0
146#"MOLA"_MOL@OH  147#"MOLA"_MOL@HO  148#"MOLB"_MOL@OH         70.3
148#"MOLB"_MOL@OH  201#"MOLB"_MOL@HO  43#"MOLA"_MOL@OH          87.1

其中原子的标识为 总体编号#"分子名称"_残基名称@原子名称. 此数据可以与前一篇文章中的数据对比, 可以看到氢键所涉及的原子完全一致, 但由于VMD和GMX的氢键判断标准无法精确地相互转换, 所以占有率只是基本一致, 无法完全一致.

其他

分析氢键时的另一个常见问题是, 如何分析GMX不支持原子间的氢键, 如涉及卤素原子的氢键. 这可以通过修改拓扑文件, 将卤素原子改为ON, 骗过GMX来达到目的. 因为GMX不会检查实际原子是什么元素, 只看其名称. 修改了原子名称GMX就会以为原子真的变了, 即便名不符实. 但注意不要把拓扑搞混了, 这种修改过的拓扑只能用来计算氢键, 不要用它来做其他计算.

◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: VMD氢键插件及其应用
后一篇: 使用gnuplot评估蛋白二级结构

访问人次(2015年7月 9日起): | 最后更新: 2022-12-17 03:04:57 UTC | 版权所有 © 2008 - 2022 Jerkwin