CADD课程学习（13）-- 研究蛋白小分子动态相互作用-III(蛋白配体复合物 GROMACS)_cadd教程

CADD课程学习（13）-- 研究蛋白小分子动态相互作用-III(蛋白配体复合物 GROMACS)

官方教程：http://www.mdtutorials.com/gmx/complex/09 analysis.html
参考中文教程：http:/jerkwin.github.io
本教程适用于2018x以上GROMACS版本

案例二：蛋白配体复合物

全部资料，包括本案例使用的文件，及案例结果文件：
链接：https://pan.baidu.com/s/1KdITmBWChKl3bNSd34MTWQ
提取码：9whx

第一步：准备蛋白拓扑文件

下载3HTB文件

查到是是否有断链，去掉结晶水分子，保存新的PDB（3HTB_clean.pdb）

• 方法一：在pymol中只保存蛋白和配体部分（2-propylphenol，2-丙基苯酚）
• 方法二：文本编组器（如Notepad中，不要用word，易出错），直接删掉PDB文件中与结晶水相关的行（残基为"HOH"等的行）即可
• 方法三：
  grep -v HOH 3HTB.pdb>3HTB1.pdb
grep -v PO4 3HTB1.pdb>3HTB2.pdb
grep -v BME 3HTB2.pdb>3HTB_clean.pdb
grep -v JZ4 3HTB_clean.pdb>3HTB_protein.pdb

注意：注意意检查pdb文件中MISSING注释下面列出的项，这些项代表了在晶体结构文件中缺失的那些原子或残基。在模拟中，末端区域的缺失可能并不会引起问题，但缺失原子的不完整内部序列或任何氨基酸残基都将导致pdb2gmx程序运行失败，必须使用其他软件对这些缺失的原子/残基进行建模并补充完整（如薛定谔蛋白准备）。

下载CHARMM力场： http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs

根据已安装python版本下载相应的脚本，此教程下载了：
cgenff_charmm2gmx_py3_nx1.py

解压缩包，获取力场文件：
执行代码：tar -zxvf charmm36-mar2019.ff.tgz

执行代码：gmx pdb2gmx -f 3HTB_protein.pdb -o 3HTB_processed.gro
将处理结构，输出一些相关信息后，提示你选择一个力场：

在本教程中，我们选用CHARMM36全原子力场，因此在命令提示行中输入1，然后巴牛提示水分子模型，选择1, Tip3p模型

第二步：准备配体拓扑文件

使用grep命令获取配体信息：
执行脚本grep JZ4 3HTB_clean.pdb > jz4.pdb

用Avogadro软件给配体加氢 ：https://sourceforge.net/projects/avogadro/
下载安装后，use Avogadro open jz4.pdb，from the"Build"menu，choose"Add Hydrogens."
Avogadro will build all of the H atoms onto the JZ4 ligand.Save a.mol2 file（File->Save As.and choose Sybyl Mol2 from the drop-down menu）named"jz4.mol2."

用文本编辑器打开jz4.mol2，并进行如下修改：

如下：

下载sort_mol2bonds.pl工具，以保证构建正确拓扑文件进行坐标匹配时，化学键的排序是正确的。
执行代码perl sort_mol2_bonds.pl jz4.mol2 jz4_fix.mol2

利用CGenFF产生配体Z4的拓扑文件：上传jz4fix.mol2，生成jz4.str文件
将z4_fix.mol2提交到CGenFF服务器（需要注册账号）https://cgenff.umaryland.edu/

利用之前下载好的cgenff_charmm2gmx_py3_nx1.py脚本处理配体文件获得配体拓扑文件
执行代码：python cgenff_charmm2gmx_py2.py JZ4 jz4_fix.mol2 jz4.str charmm36-mar2019.ff
注意以上脚本需要NetworkX package，version 1.11，和Networkx2.x不兼容，会出错。出现下面的提示，表示转换成功，获得jz4.itp文件和jz4.prm文件。

**注释：**该配体引入了新的结合参数，这些参数不是现有力场的一部分，并且将这些参数写入到名为"JZ4.prm"的文件中，该文件的格式为GROMACS.itp。

构建复合物（Build the Complex）：
生成配体gro文件：
执行代码：gmx editconf -f jz4_ini.pdb -o jz4.gro

复制3HTB_processed.gro为complex.gro文件
执行代码：cp 3HTB_processed.gro complex.gro

构建复合物（Build the Complex）：
将配体jz4.gro文件中的信息复制到complex.gro中的末尾部分：


最后结果

注意：由于我们在complex.gro文件中添加了22个以上的原子，因此请增加complex.gro的第二行以反映此更改。
现在，坐标文件中应该有2636个原子。

准备复合物拓扑文件：
在系统拓扑中加入JZ4配体的参数信息。topol.top文件开头和末尾部分，只需在topol.top中插入两行，是显示#include ”jz4.itp”。
文件开头部分：

准备复合物拓扑文件：
在系统拓扑中加入Z4配体的参数信息。Topol.top文件开头和末尾部分，只需在topol.top中插入一行，显示#include “jz4.itp”。
文件结尾

第三步：定义单位盒子并填充溶剂

定义一个模拟用的盒子并添加溶剂要分两步完成

1. 使用editconf模块定义盒子的尺寸
2. 使用solvate模块（以前的版本中称为genbox）向盒子中填充水

使用editconf来定义盒子：
执行代码： gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

执行代码： gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

使用vmd查看solv.gro文件：

第四步：添加文件

以上获得的溶液体系带有6个正电荷：topol文件[atoms]的最后一部分：gtot 6，需要添加离子使总电荷为0.
使用： genion工具

注释：genion的功能是读取拓扑信息，然后将体系中的一些水分子替换为指定的离子。

genion需要的输入文件称为运行输入文件，扩展名为.tpr这个文件可使用GROMACS的grompp（GROMacs Pre-Processor）
模块产生而且后面我们运行模拟时也会用它.

grompp的功能是处理坐标文件和拓扑（它描述了分子）以产生原子级别的输入文件(.tpr)如tpr文件包含了体系中所有原子的所有参数

为了用grompp产生t如文件，我们还需要一个扩展名为.mdp(molecular dynamics parameter）的输入文件、grompp会将坐标和拓扑信息与mdp文件中设定的参数组合起来生成tpr文件

使用下面的命令来产生**.tpr**文件
执行代码：gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
以上代码得到一个原子级别的.tpr文件，将此文件用于genion.

执行代码gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

注意：charmm-2022的版本离子格式改变了， “Na+”会变成“SOD”，“K+”会变成“POT”，“Cl-”会变成“CLA”
注释出现提示时，选择第15组”SOL“用于嵌入离子。

第五步：能量最小化

使用grompp处理这个参数文件，以便得到二进制的输入文件：
执行代码 ``gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

调用mdrun来进行能量最小化：
执行代码 gmx mdrun -v -deffnm em

**注释：**使用选项可快速看到运行结果它使mdrun输出更多信息，这样就会在屏幕上输出每步运行的情况，-deffnm选项定义了输入文件和输出文件的名称

有两个重要的指标未决定能县最小化是否成功；

• 第一个是势能（在能量是小化过程的最后输出，即使你未使用-v选项. Epot应当是负值，根据体系大小和水分子的多少，大约在10^5-106的数量级（对水中的单个蛋白质而）。
• 第二个重要的指标是力的最大值Fmax. 我们在minim.mdp中设置的目标是emtol=1000.0，这表示Fmax的目标值不能大于1000 KJ mol-1 nm-1.能量最小化完成后，有可能得到一个合理的Epot，但Fmax>emtol.如果是这样，用于模拟时你的体系可能不够稳走， 思考一下为什么会这样，可能需要更改一下能量最小化的参数设置（integrator，emstep等），再试试重新进行能量最小化过程

使用energy模块来分析任何一个edr文件：
执行代码：gmx energy -f em.edr -o potential.xvg
注释：费示输入11 0来选择势能Potential（11），并用零（0）来结束输入

执行代码: xmgrace potential.xvg

第五步：NPT平衡

对配体施加约束
首先，为JZ4创建一个仅包含其非氢原子的索引组（在index_jz4.ndx文件中为group 3):

执行代码： gmx make_ndx -f jz4.gro -o index_jz4.ndx
> 0 & ! a H*
> q
执行代码：gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000
> 3

1.对配体施加约束
将该信息添加到top文件中：
topol.top文件结尾部分：（这里我们只做限制配体）

或者

1. 温度耦合基团处理
   Do not couble every single specie in your system separately.
   执行代码： gmx make_ndx -f em.gro -o index.ndx
   合井"Protein"and JZ4"groups：
   1 | 13
q
   然后可设置tc-grps=Protein JZ4 Water and ions
   

第六步：NVT平衡

使用grompp和mdrun，像在能量最小化过程中做的一样
执行代码： gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
执行代码： gmx mdrun -deffnm nvt

第七步：NPT平衡

使用grompp和mdrun，像在能量最小化过程中做的一样
执行代码：gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr -maxwarn 100
**注释：**使用-选项以包括NVT平衡过程中的产生的检查点文件，这个文件包含了继续模拟所需要的所有状态变量。为使用NVT过程中得到的速度我们必须包含这个文件，坐标文件（-c）是NVT模拟的最终输出文件。
执行代码： gmx mdrun -deffnm npt

第八步：MD成品

随着两个平衡阶段的完成，体系已经在需要的温度和压强下平衡好了

我们现在可以放开位置限制并进行成品MD以收集数据了.这个过程跟前面的类似.运行grompp时，我们还要用到检查点文件（在这种情况下，其中包含了压力耦合信息）.

我们要进行10ns的MD模拟，所用的参数文件md.mdp如下图：

使用grompp和mdrun类似前面优化过程中做的一样，运行10ns的成品MD
执行代码： gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr
注释： 使用-t选项以包括NVT平衡过程中的产生的检查点文件。这个文件包含了继续模拟所需要的所有状态变量，为使用NVT过程中得到的速度我们必须包含这个文件，坐标文件（-c）是NVT模拟的最终输出文件。
执行代码： gmx mdrun -deffnm md_0_10

第九步：分析

提取能量数据观察动力学能量是否达到平衡，代码如下：
温度：Temperature
gmx energy -f md_0_10.edr -o md_0_10_Temperature.xvg
15 0
压力：Pressure
gmx energy -f md_0_10.edr -o md_0_10_Pressure.xvg
16 0
密度：Density
gmx energy -f md_0_10.edr -o md_0_10_Density.xvg
22 0
势能：Potential
gmx energy -f md_0_10.edr -o md_0_10_Potential.xvg
11 0
总能量：Total-Energy
gmx energy -f md_0_10.edr -o md_0_10_TotalEnergy.xvg
13 0

xmgrace可观察动力学是否达到能量平衡，代码如下：
xmgrace md_0_10 _Potential.xvg
xmgrace md_0_10_ TotalEnergy.xvg
xmgrace md_0_10__Temperature.xvg
xmgrace md_0_10_ Pressure.xvg
xmgrace md_0_10_ Density.xvg

现在已经完成了对蛋白质的模拟，我们应该来分析一下我们的体系。

1. 修正轨迹-1
   使用trjconv工具来处理体系中的任何周期性，从而获得“修正”后的轨迹：
   执行代码：gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact
   注释：Select 1(“Protein”）as the group to be centered and 0（“System"）for output.
   修正轨迹-2
   使用triconv工具来处理体系中的任何周期性，从而获得“修正”后的轨迹：
   执行代码： gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0
   **注释：**提取轨迹的第一帧（t = 0 ns）
   选 0

修正轨迹-3
使用tricomv工具来处理体系中的任何周期性，从而获得“修正”后的轨迹：
为了使可视化更加流畅，执行旋转和平移拟合可能会有所帮助
执行代码： gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans
注释：Choose “Backbone” to perform leastsquares fitting to the protein backbone，and “System” for output.

2. 回旋半径分析蛋白的密实度（Radius of gyration）
   执行代码：gmx gyrate -s md_0_10.tpr -n index.ndx -f md_0_10.xtc -o md_0_10_Gyrate-Protein.xvg
   Protein

3. 均方根偏差分析蛋白构象变化（RMSD）
   执行代码：gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg
   Protein_JZ4
   **注意：**The calculated RMSD should be about 0.1 nm (1A).indicating only a very small change in the ligand’s position.

4. 均方根涨落分析蛋白原子构象变化（RMSF）
   执行代码：gmx rmsf -s md_0_10.tpr -f md_0_10.xtc -o rmsf.xvg -n index.ndx
   Protein

5. 蛋白配体构象聚类（To cluster the protein least squaies fit and RMSD/output）
   执行代码：gmx cluster -s md_0_10.tpr -n index.ndx -f md_0_10.xtc -cl 3HTB_Clusters.pdb -dist 3HTB_Cluster_RMSD-Dist.xvg -b 800 -e 1000
   Protein_JZ4
   Protein_JZ4
   EOI

6. 蛋白配体相互作用氢键分析（Hbond）Number of interactions between two groups
   执行代码: gmx hbond -s md_0_10.tpr -n index.ndx -f md_0_10.xtc -num md_0_10_HB-Num_Protein_JZ4.xvg
   Protein 13
   EOI

7. 蛋白配体相互作用能 Protlin-Ligand Interaction Energy
   为了量化JZ4和T4溶菌酶之间相互作用的强度，计算这两个物种之间的非键相互作用能可能是有用的。GROMACS能够分解任意数量的已定义基团之间的短距离非结合能。重要的是要注意，该度量值不是自由能或结合能。
   
   执行代码： gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr
   执行代码： gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpugmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu
   执行代码： gmx energy -f ie.edr -o interaction_energy_Coul-SR.xvg
51 0
   执行代码： gmx energy -f ie.edr -o interaction_energy_LJ-SR.xvg
52 0
   The terms we are interested in are Coul-SR:Protein-Jl4 and uronirroelmet.

生成动态轨迹文件

在任何用周期性边界条件进行的模拟中，分子可能出现“断裂”或可能在盒子中来回“跳跃”。为了重新定位蛋白质并重新包裹单元细胞内的分子，以恢复所需的菱形十二面体形状，调用 trjconv：

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact
1

选择“Protein”作为中心，选择“System”输出。注意将配合物(蛋白质-配体，蛋白质-蛋白质)定位于中心对于包含许多跨越周期边界的长时间模拟可能是困难的。在这些情况下(特别是在蛋白质-蛋白质复合物中)，可能需要创建一个定制的索引组（对应于复合物中一个蛋白质的活性位点或一个单体的界面残基），用于定心。

为了更平滑的可视化，进行适当的旋转和平移拟合可能是有益的（请注意，同步PBC重新包装和拟合坐标在数学上是不兼容的。如果您希望执行拟合（这对可视化很有用，但对大多数分析例程不是必需的），请按此处所示，分别执行这些坐标操作）。按如下方式执行trjconv：

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans

1
2

选择“Backbone”进行蛋白质主链的最小二乘拟合，选择“System”进行输出。

由于原始文件太大，我们间隔抽取运动轨迹

gmx trjconv -s md_0_10.tpr -f md_0_10_fit.xtc  -o whole.pdb  -dt 500 -b 1000 -e 20000
1

选择“Protein_JZ4”进行输出。
注意：没有“Protein_JZ4”的话，需要做索引,

gmx make_ndx -f em.gro -o index.ndx
1

gmx trjconv -s md_0_10.tpr -f md_0_10_fit.xtc  -o whole.pdb  -dt 500 -b 1000 -e 20000 -n index.ndx
1

whole.pdb 导入Pymol, 做平滑处理， 输入“intra_fit complex”， “smooth”