仓库源文，站点原文

layout: post title: GROMACS分析教程：氢键分析 categories:

科 tags:
gmx
vmd
氢键

2017-01-09 16:56:34

卢天给过一个VMD的tcl脚本, 用于计算不同z位置水能形成的平均氢键数. 但轨迹大了以后, VMD的tcl脚本分析起来有点吃力, 且得到的结果与GROMACS的默认氢键标准存在差距. 这里我以此为例, 展示下如何组合GROMACS自带的工具进行复杂一点的分析, 并且利用gmx trjconv的外挂脚本功能让分析自动化. 示例所用GROMACS版本为5.1.4.

运行模拟

示例文件用的是GROMACS自带的spc216.gro, 只不过沿z方向将盒子扩大了一些, 留出一定的真空层. 运行NVT模拟.

gmx grompp
gmx mdrun

预处理轨迹

得到轨迹后, 先对轨迹进行预处理

gmx trjconv -f traj_comp.xtc -pbc whole -o traj_pbc.xtc

选择0, 对整个体系进行PBC校正, 保证分子完整. 如果需要, 根据情况进行其他特殊处理.

分析单帧构型

在对整个轨迹进行分析前, 先选择一帧构型测试命令和脚本的正确性, 这样可以节约时间.

gmx trjconv -f traj_pbc.xtc -o traj.gro -dump 0 -sep

选择0, 输出整个体系, 得到traj0.gro.

选择需要的原子组

我们的最终目的是分析不同高度层内水分子所成的平均氢键数目, 要使用的命令是gmx hbond. 使用这个命令分析氢键时, 需要指定两个组, 它们可以相同也可以不同. 如果相同, 分析的是本组分子之间形成的氢键; 如果不同, 分析的则是两组分子之间形成的氢键.

对于我们的目的而言, 我们需要两个组, 组1 Wsel是某高度层内的水分子, 组2 Wexc是排除Wsel后的其余水分子. 这样我们就可以计算组1之间的氢键数目N(Wsel-Wsel), 组1和组2之间的氢键数目N(Wsel-Wesc), 而Wsel中的水分子所成的平均氢键数目 Nhb = ( 2*N(Wsel-Wsel)+N(Wsel-Wexc) )/N(water). N(Wsel-Wsel)之前的因子2是因为gmx hbond报告同组分子所成的氢键数目时会排除重复, 而我们不需要这样做.

选择不同高度层内的分子, 选择方式有几种策略, 或者根据氧原子的位置, 或者根据分子的质心位置, 几何中心位置, 或者根据任意原子位置等, 具体的语法看参考gmx select的语法及用法.

选中一定高度层内的水分子, 最简单的策略是根据氧原子的位置, 先选中符合条件的氧原子, 然后扩展到与氧原子所属分子相同的原子. 下面是的语句选择z坐标位于3到3.2之间的水分子以及其余水分子:

"Wsel" same mol as (name OW and (3<z and z<3.2))
"Wexc" not same mol as (name OW and (3<z and z<3.2))

如果利用质心条件的话, 可以使用下面的语句

resname SOL and (3<res_com z and res_com z<3.2)
或 resname SOL and (res_com z 3 to 3.2)

如果要同时选择Wsel和Wexc的话, 最好先定义一个变量, 这样写起来简洁, 执行起来效率也更高

Wsel = resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel

最终, 我们可以使用下面的命令获得gmx hbond所需要的分组及其原子数目

gmx select -f traj0.gro -s -select 'Wsel = resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel' -os traj0.xvg -on traj0.ndx

traj0.xvg最后一行的第二列是Wsel中的原子数目, 其1/3就是我们所需要的N(water); traj0.ndx中保存了我们需要的两个分组.

对原子组进行分析

获得了原子组的索引文件, 就可以用它进行氢键数目分析了.

gmx hbond -f traj0.gro -n traj0.ndx -num traj0.xvg

运行上面的命令时会提示选择两个分组, 手动选择不利于自动化, 我们可以利用echo命令和管道来替代手动选择

<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">echo 0 0 | gmx hbond -f traj0.gro -n traj0.ndx -num traj0_00.xvg echo 0 1 | gmx hbond -f traj0.gro -n traj0.ndx -num traj0_01.xvg </pre></div>

traj0_00.xvg最后一行的第二列是我们需要的N(Wsel-Wsel), traj0_01.xvg最后一行的第二列则是我们需要的N(Wsel-Wexc).

整理报告分析结果

有了这两个数据, 以及前面的水分子数目, 我们就可以计算出每个水分子所成的平均氢键数目了. 当然还是需要用脚本来自动计算, 这样后面才能实现自动化.

<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">Nsel=$(tail -n 1 traj0.xvg) N00=$( tail -n 1 traj0_00.xvg) N01=$( tail -n 1 traj0_01.xvg) echo $Nsel $N00 $N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}' </pre></div>

上面的脚本先获取每个文件的最后一行, 然后使用管道将其传awk, awk计算出所需要的值并保存到HB.xvg中.

通用单帧构型分析脚本

将前面的命令写到一个bash脚本中, 就可以自动执行上面的分析了

<table class="highlighttable"><th colspan="2" style="text-align:left">hb.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">i=0 file=traj$i gmx select -f $file.gro -s -select 'Wsel = resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel' -os $file.xvg -on $file.ndx echo 0 0 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg echo 0 1 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_01.xvg Nsel=$(tail -n 1 $file.xvg) N00=$( tail -n 1 ${file}_00.xvg) N01=$( tail -n 1 ${file}_01.xvg) echo $Nsel $N00 $N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}' rm -rf $file.gro $file.xvg $file.ndx ${file}_00.xvg ${file}_01.xvg </pre></div> </td></tr></table>

我们在最前面定义了帧号, 这样脚本就很容易用于其他帧了, 只要改变帧号对应的变量i就可以了. 脚本的最后我们删除了用到的中间文件, 这样既可以避免重复运行GROMACS工具时因备份文件太多而导致的错误, 也可以让我们的目录更清爽.

分析整条轨迹

既然已经可以对一帧构型进行分析, 并完成了一个通用的分析脚本, 那么将脚本用于多个构型就比较简单了. 最直接的方式就是先使用gmx trjconv输出所有的构型, 然后循环处理每帧构型. 比如我们有帧号为0到100的多个构型文件, 那么在bash中可以使用

<table class="highlighttable"><th colspan="2" style="text-align:left">hb.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">for i in {0..100}; do file=traj$i gmx select -f $file.gro <略> -os $file.xvg -on $file.ndx echo 0 0 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg echo 0 1 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_01.xvg Nsel=$(tail -n 1 $file.xvg) N00=$( tail -n 1 ${file}_00.xvg) N01=$( tail -n 1 ${file}_01.xvg) echo $Nsel $N00 $N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}' rm -rf $file.gro $file.xvg $file.ndx ${file}_00.xvg ${file}_01.xvg done </pre></div> </td></tr></table>

这种方式的缺点在于需要先输出一大堆构型文件, 当处理的帧数很多时, 就凌乱了. 好在 gmx trjconv支持一个外挂脚本的选项, -exec "命令", 可以在输出每一帧后对此帧构型执行指定的命令, 脚本或程序, 且以帧号作为命令行参数. 举例来说, 如果选项为-exec "cmd", 那么对输出每帧构型后对其执行的命令就是cmd 帧号. 利用这一功能, 我们就不再需要自己写循环了, 只要将我们上面的脚本改为以帧号为输入参数就可以了

这样我们直接使用

gmx trjconv -f traj_pbc.xtc -o traj.gro -sep -exec "bash hb.bsh"

就可以自动对整条轨迹进行分析, 得到所需要的数据了. 对我们的测试轨迹, 得到的平均数目为3.5左右, 符合预期.

总结说明

gmx trjconv的-sep选项只支持gro或pdb格式, 不支持二进制格式构型文件.
这种组合GROMACS已有工具, 以及简单数据处理小脚本的分析方式不需要自己编写代码处理轨迹, PBC, 分析等问题, 只需要一些代码处理输出输入文件, 流程具有通用性, 但难度远小于自己从头写代码处理轨迹.
这种处理方法在分析时每次只处理一帧构型, 对于需要多帧构型才能计算的物理量就不适用了. 在那种情况下, 只能先输出所有的信息, 再统一处理了.
这种方式每次只处理一帧构型, 效率不是太高, 但对机器的要求低, 不像VMD那样要先载入整条轨迹, 需要很大内存.
如果你对于效率很在意, 那么可以采用并行的方式进行处理: 先输出所有的构型文件, 使用脚本并行处理左右构型文件, 最后再将这些文件的结果合并. 至于如何并行处理, 这里就不再细说了, 可参考我以前的两篇博文Bash脚本实现批量作业并行化, GNU Parallel.
如果你想试验一下这种处理方式, 那可以试着完成我前面提出的一个问题, 蛋白口袋或纳米通道内水分子的个数统计. 当然, 这些口袋或通道都是柔性的, 否则的话, 就没有必要这么处理了.

附: 类似功能的VMD tcl脚本

使用VMD的tcl脚本进行分析时, 思路也是类似的. 但值得注意的是, VMD的默认氢键标准与GROMACS不同, 这是因为目前存在多种判断氢键的标准. 简言之, VMD的 3.5埃-40度标准大致和GROMACS的标准一致. 此外, VMD还没有考虑PBC的问题. 详情可见前一篇博文GROMACS的默认氢键标准.

<table class="highlighttable"><th colspan="2" style="text-align:left">HB.tcl</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">#mol delete all #mol new conf.gro traj_pbc.xtc set sel "same resid as resname SOL and z 30 to 32 " set Zset "z 30 to 32 and x 4 to 14 and y 4 to 14" set Wsel [atomselect top "same resid as name OW and $Zset"] set Wexc [atomselect top "same resid as name OW and not($Zset)"] set Nfrm 100 set Rcut 3.5 set Angle 40. puts "#Frm #Wsel #Wexc #sel-sel #sel-exc #exc-sel #HB/Wat" for {set i 0} {$i<=$Nfrm} {incr i} { $Wsel frame $i $Wsel update $Wexc frame $i $Wexc update set Nsel [expr [$Wsel num]/3.] set Nexc [expr [$Wexc num]/3.] set Nss [llength [lindex [measure hbonds $Rcut $Angle $Wsel] 0]] set Nse [llength [lindex [measure hbonds $Rcut $Angle $Wsel $Wexc] 0]] set Nes [llength [lindex [measure hbonds $Rcut $Angle $Wexc $Wsel] 0]] puts "$i $Nsel $Nexc $Nss $Nse $Nes [expr (2.*$Nss+$Nse+$Nes)/$Nsel]" } </pre></div> </td></tr></table>