amber教程4.6：对体系氢键分析

一、简介

利用amber的cpptraj模块对显式溶剂模型进行H键作用分析，对得到的溶质-溶质间H键存在与否的变化进行gnuplot绘图，对体系整体H键的变化,体系的RMSD变化进行xmgrace绘图，对溶质内H键寿命变化进行计算统计。

二、下载教程提供的准备文件

1、top文件：

wget https://amber.utah.edu/AMBER-workshop/London-2015/Hbond/trpzip2.ff10.tip3p.parm7

2、轨迹文件：

wget https://amber.utah.edu/AMBER-workshop/London-2015/Hbond/trpzip2.unfold.nc

3.cpptraj程序的input文件：

wget https://amber.utah.edu/AMBER-workshop/London-2015/Hbond/cpptraj.hbond.in

三、执行cpptraj模块

导入top文件，赋予标签：

cpptraj
parm trpzip2.ff10.tip3p.parm7 [top]

导入轨迹文件：

trajin trpzip2.unfold.nc parm [top]

计算轨迹中溶质和溶质、溶质和溶剂(WAT)之间的H键还有溶质、溶剂和溶质之间形成的桥接作用：

#这里U特指溶质，V特指溶剂
hbond All out All.hbvtime.dat solventdonor :WAT solventacceptor :WAT@O \
  avgout All.UU.avg.dat solvout All.UV.avg.dat bridgeout All.bridge.avg.dat
run

 ALL：是生成数据集的名字，类似的有下面的Backbone

solventdonor：指定溶剂H供体（也可以是离子），这里限制为残基WAT

solventacceptor：指定溶剂H受体（同样可以是离子），这里限制为残基原子O（......）

这里没有限制溶质的H供受体哦

avgout：将溶质-溶质H键平均值写入到后接着的文件

solvout：将溶质-溶剂H键的平均值写入到后接的文件

bridgeout： 将溶质-溶剂桥接作用写入到后接的文件

生成的ALL.UU.avg.dat文件内容如下：

Acceptor、DonorH和Donor是参与组成H键的原子，Frames是这个H键存在的帧数，2099就是出现了2099次，Frac是这个次数占总次数的百分比。AvgDist是根据这些帧计算出的H键平均距离，AvgAng则是三个原子组成的平均角度。

 生成的ALL.UV.avg.dat文件：

Count也是出现的帧数总值，为什么会高出实际总帧数（......）

生成的All.hbvtime.dat文件内容：

All[ID] 组成桥接的三个原子序号x(x+x+)，若有多组桥接作用。组间由逗号分隔

这个文件中，行数是帧数，表示一帧中存在多少H键和溶剂桥接作用

生成的ALL.bridge.avg.dat文件：

计算轨迹中溶质蛋白主链上的H键：

hbond Backbone :1-12@C,O,N,H avgout BB.avg.dat series uuseries bbhbond.gnu
run

Backbone是一个数据集，series 以每帧有哪些残基组成H键出现的格式保存文件，uuseries 将溶质-溶质之间的H键写入到后接着的文件（series uvseries [filename]则是溶质-溶剂之间的H键）

用gnuplot和生成的gnu文件进行绘图：

gnuplot bbhbond.gnu
在win下的版本：
load "bbhbond.gnu"

图中， 部分残基组成的H键在2000帧后存现，同时部分残基形成的H键消失，和下图RMSD观察到的蛋白构象变化有关。

利用数据集的内容生成新文件：

# Create file containing only number of solute-solute and solute-solvent
# hydrogen bonds, and number of solute-solvent-solute bridges.
create nhbvtime.agr All[UU] Backbone[UU] All[UV] All[Bridge]
run

  这个文件的数据是接成同一列的：

所以这里实际的列数>帧数（3290*4），还要加上注释

xmgrace绘图：

xmgrace nhbvtime.agr

计算主链的RMSD：

# Calculate RMSD to first frame.
rms BBrmsd :1-12@C,CA,N out BBrmsd.agr
run

xmgrace绘图：

xmgrace BBrmsd.agr

 结合轨迹H键数目的变化图，在2000Frames后，H键数目减少，溶质的RMSD发生明显的改变。

H键寿命分析：

# Perform lifetime analysis on backbone hydrogen bonds
lifetime Backbone[solutehb] out backbone.lifetime.gnu
runanalysis

寿命lifetime的定义：一个物质存在的时长，就是寿命。这里研究的对象是H键存在的时长，若一个H键在连续的10帧中都存在，这个H键的寿命时长就是10

我们设置一个截断距离，但一帧内H键时长大于这个值，就判断其存在：

cut <cutoff>
#默认的cutoff是0.5(单位？？）

标准化：对数据计算均值，所有的数值都进行减掉均值后在除以标准差，得到一组均值为0，方差为1的标准数据，这个过程称为标准化。

[rawcurve]
#添加这个选项后，不会对数据进行标准化处理，在out [filename] 得到的文件不会有crv.[filename]

out 会输出多个文件，一个backbone.lifetime.dat(这里简称为[file])，一个crv.[file],这个文件是对前者进行标准化后得到的数据。

若是添加参数window（指定计算的帧数范围），还会多出文件max.[file]和avg.[file]，前者是这个帧数范围内的最大寿命，后者是这个帧数范围内的平均寿命。

得到的*.dat文件：

这里的lifetime_00027是寿命期的个数，lifetime_00027[max]是存在的最长寿命期的长度，lifetime_00027[avg]是平均长度，lifetime_00027[frames]是存在的帧数。

e.g.第一个是11号残基TRP和1号残基SER形成的H键（lifetime_00027[name]列），这个H键分别出现了两次（lifetime_00027列），最长的一次寿命长度是1帧（lifetime_00027[max]列），平均帧时长是1.0000（lifetime_00027[avg]列）

四、利用gromacs和gunplot绘制H键热图：

利用gromcas完成MD，用gmx hbond计算轨迹中的h键数据：

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

这里的traj.xtc是要分析的轨迹文件，也可以是gro。在回车这个命令后，会出现选择对应的组所以，分别选择要分析的两组相互作用的对象。例如1（回车）2（回车）就是分析1和2之间H键作用。

得到的一个xpm文件，将文件用xmp2all.bsh转化为xyz格式文件,格式内容大概如下：

1 0 0
2 0 1
3 0 0
4 0 0

1 0 0
2 0 1
3 0 1
4 0 0

1 0 1
2 0 0
3 0 0
4 0 0

第1列表示x，是时间。第2列是y，表示，第3列是z轴，只有1和0两个值，表示H键存在或不存在。

这里每个区块表示1个H键随时间的变化，如第1列，表示一个H键在第2ns出现，后第3ns消失，依次类推。

转化的命令如下，在win下可以借用git-bash软件完成：

bash xpm2all.bsh hbond.xpm

利用下面的脚本hbond.gun，用命令load “hbond.gun”，要修改文件名，完成画图：

# 设置输出文件格式为PNG
set terminal png size 800,600 enhanced font "Arial,12"

# 设置输出文件名
set output "heatmap.png"

# 设置标题和轴标签
set title "Heatmap of XYZ Data"
set xlabel "X-axis (0 to 400)"
set ylabel "Y-axis (1 to 5)"

# 设置颜色映射
set palette defined (0 "white", 1 "blue")

# 设置颜色范围
set cbrange [0:1]

# 设置不显示色标
unset colorbox

# 设置绘图区域
set xrange [0:400]
set yrange [1:5]

# 设置绘图样式为热图
set style line 1 linecolor rgb "black" linewidth 0.5
set style fill solid 1.0 border -1

# 读取数据并绘制热图
plot "data.xyz" using 1:2:3 with image

这里data.xyz要修改为用xpm2all.bsh转化后的xyz格式文件 

参考链接：

hydrogen-bond-analysis-within-a-protein

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