gromacs PCA 做自由能景观图

在文献中看到有人做自由能景观图，记录下自己的操作过程
参考http://www.kangsgo.com/49.html
http://sobereva.com/73
https://www.zhihu.com/question/19895566

跑完动力学后蛋白比对

gmx trjconv -f md.xtc -s  md.tpr -fit rot+trans -o mdfit.xtc

选择backbone 以及 system用于输出

然后生成协方差

gmx covar -s md.gro -f mdfit.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covapic.xpm

选择C-alpha

这里可以用xmgrace 查看生成的

投影本征向量

gmx anaeig -f mdfit.xtc -s md.tpr -v eigenvectors.trr -first 1 -proj pc1.xvg
gmx anaeig -f mdfit.xtc -s md.tpr -v eigenvectors.trr -first 2 -last 2 -proj pc2.xvg

然后利用perl脚本合并

perl sham.pl -i1 pc1.xvg -i2 pc2.xvg -data 1 -data2 1 -o gsham_input.xvg

这里perl 脚本只是一个实现pc1.xvg pc2.xvg 合并的功能，也可以使用awk 命令快速实现

#!/bin/bash

# 提取PC2.xvg的第二列
awk '{print $2}' pc2.xvg > pc2_values.txt

# 合并PC1.xvg和提取的PC2列
paste pc1.xvg pc2_values.txt > gsham_input.xvg

再使用sham工具计算Gibbs自由能：

gmx sham -f gsham_input.xvg -ls FES.xpm

用脚本xpm2txt.py转换为文本文件用origin软件做图

python2.7 xpm2txt.py -f FES.xpm -o free-energy-landscape.txt

用origin 软件作图
打开free-energy-landscape.txt，x y z 分别对应pc1 pc2 能量项,

然后点击Plot —> 3D XYZ —> 3D Scatter就好了(这步没看出来有什么用）

数据列表转换为矩阵表(Matrix)
Worksheet —> Convert to Matrix —> Direct —> Open Dialog

在弹出的对话框中的Data Format 选项中选择X across column，X values in选项中选择 First Data Row，勾选Y Valuesin First Column，保持Even Spacing Relative Tolerance (%)，
这里需要注意一下，有时候Origin程序会在转换矩阵表的时候报错，一般是“ Your X data are not evenly spaced, it has fluctuation of 11.11%.Please try changing the tolerance setting.If the fluctuations are very large, you should transform yourdata to XYZ column format by Worksheet: Convert to XYZ menu and then use the matrix gridding tool. ”
遇到这种情况，一般的解决方式是增大Even Spacing Relative Tolerance (%) 值就可以了，但是往往这时候数据点的间隔会有所增大，不过一般来说不会有太大的影响。

按照教程自己的报错，重新选择了XYZ gridding,成功转换成 矩阵

Origin菜单栏上的Plot —> 3D Surface —>Color Map Surface
这样做出来粗糙的一个3D自由能景观图。

接下来是修改参数
双击绘制好的平面，在弹出的对话框Plot Details 中，选择左侧对话框中Graph - Layer【注意，是在Layer激活的状态下，如下图所示】，然后在右侧的对话框中选择Size/Speed，在下面的Worksheet data, Skip Points if needed下的数据框中填写数字：数值越大，图中的平面约光滑。自己修改为300 200 200

然后点击level, 出现set level, 修改increment, 这里自己修改为1

然后再点击Surface, 修改 Enable Grids–None,同时在surface/projections取消勾选contour line

然后再修改色带legend, 双击legend,弹出color scale control对话框

点击Increment,这里自己修改为20，这样色带标签就只有三个数字了，画面更简洁。

修改level 中的 Increment 时值小点，一般为0.04-0.08，total 值越多，曲线越光滑。
在色带label 中的 increment 设置为20-40（label 将会是 level 中的 0.04乘以这里的值，这样色带label 稀疏，值可以清楚的显示出来。

同时也可以做一个2D map 图，这里2Dmap 参数设置需要同3D相同，这样两者结果一致

接着我们要从自由能景观图中获取点的坐标，首先获取最低势能构象的X Y Z坐标，获取后再origin的sheet中找到，然后在gsham_input.xvg文件中找到相对应的XY坐标，gsham_input.xvg有三列值，第一列为时间（ps）,第二列为X坐标，第三列为Z坐标，这样匹配到时间戳后(假设这里第45000ps结构为最低势能），从VMD中提取相对应的复合物坐标，对该结构进行详细的分析。
提取的脚本为
gmx trjconv -f md_fit.xtc -s md.tpr -dump 45000 -o representative_45000ps.pdb

看到有人问xpm2txt.py的文件，这里贴出来

#!/usr/bin/env python

import sys

"""
Utility tool to convert xpm files generated by GROMACS to a 3-column text file.
"""

USAGE = "USAGE: xpm2txt.py -f <input xpm file> -o <output txt file> [-s]\n"
USAGE+= "Options:\n"
USAGE+= "\t-s\t(int)\tSorts the output by a given column"
USAGE+= "\n" # always keep this line

# Parse arguments
read_input, read_output, sort = False, False, False
xpm_file, out_file, column_sort = None, None, None
for arg in sys.argv[1:]:
    if read_input:
        read_input = False
        xpm_file = arg
    elif read_output:
        read_output = False
        out_file = arg
    elif sort:
        sort = False
        column_sort = int(arg)
    if arg[0