CADD课程学习（13）-- 研究蛋白小分子动态相互作用-II(水中的溶菌酶 GROMACS)

案例一：模拟体系：水盒子中的蛋白质（溶菌酶，lysozyme）及离子

本案例为一个蛋白质（溶菌酶，lysozyme）加上离子在水盒子里的模拟过程。
官方教程链接：http://www.mdtutorials.com/gmx/lysozyme/index.html
非官方中文翻译链接：http://jerkwin.github.io/GMX/GMXtut-1/
下载地址：http://public-ehs.oss-cn-hangzhou.aliyuncs.com/packages/Lysozyme.targz本教程适用于2018.x以上GROMACS版本

整体步骤

本案例需要的外部文件及跑完分子模拟的结果文件：
链接：https://pan.baidu.com/s/1jgoJR_kBSH4tdrTpBw4D0Q
提取码：asd8

第一步：准备拓扑文件

PDB 网站下载溶菌酶（PDBID：1AKI）

GROMACS基础：
获得GROMACS某一模块的帮助信息。可使用：
gmx help (module)
或
gmx (module) -h
使用时只要将其中的（module）替换为要查询的命令的实际名称即可.相关信息会输出到终端，其中包括可用的算法，选项，需要的文件格式/已知的缺陷和限制，等等，对GROMACS的新用户来说，查看常用命令的帮助信息是了解每个命令功能的有效方式。

处理pdb文件
查到是是否有断链，去掉结晶水分子，保存新的PDB（1AKI_clean.pdb）
方法一：在pymol中只保存蛋白部分
方法二：文本编辑器（如Notepad中，不要用Word，易出错），直接删掉PDB文件中与结晶水相关的行（残基为“HOH”的行）即可
方法三：grep -v HOH 1aki.pdb>1AKI_clean.pdb
注意：注意意检查pdb文件中MISSING注释下面列出的项，这些项代表了在晶体结构文件中缺失的那些原子鼓残基：在模拟中，未端区域的缺失可能并不会引起问题，但缺失原子的不完整内部序列或任何氨基酸残基都将导致pdb2gmx程序运行失败，必须使用其他软件对这些缺失的原子/残基进行建模并补充完整（如薛定谔或MOE等准备蛋白）。

使用：pdb2gmx模块
可产生：

分子拓扑文件
位置限制文件
后处理结构文件
拓扑文件（默认为topol.top）包含了定义模拟中质用分子的质有必需信品，包适非续参数（原子类型和电荷）和键合参数（键长，键角和二面角）

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

力场包含了将要写人到拓扑文件中的信息。这个选择非常重要！
你必须仔细了解每个力场，并决定哪个力场最适用于你的体系。
在本教程中，我们选用全原子OPLS力场，因此在命令提示行中输入15，然后回车

上一条命令的其他参数：
pdb2gmx程序还接受很多其他选项，下面列出常用的几个：

-ignh:忽略PDB文件中的氢原子，对NMR结构非常有用，否则，如果存在氢原子，它们的名称和顺序必须与GROMACS力场中的完全一样，对氢原子存在各种不同的约定所以处理起来有时让人很头疼！如果你需要保留初始氢原子的坐标但需要重命名，可以试试Linux的sed命令
-ter：交互式地为N未端和C末端分配电荷态
-inter：交互式地为Glu(谷氨酸)，Asp(天冬氨酸)，Lys(赖氨酸)，Arg(精氨酸)和His指定电荷态，选择涉及二硫键的Cys(胱氨酸)
现在你已经得到了三个新的文件1AKI_processed.gro，topol.top和posre.itp 1AKI_processed.gro 是GROMACS格式的结构文件，包含了力场中定义的所有原子即，已经将氢原子加到蛋白质中的氨基酸上了）.topol.top文件是体系的拓扑文件（稍后会解释）posre.itp文件包含了用于限制重原子位置的信息（后面解释）。

第二步：检查拓扑文件

用文本编辑器打开topol.top文件，即：vi topol.top
注释：此行调用了OPLS-AA力场的参数，它位于文件的开头，这意味着下面的所有参数都来自OPLS-AA力场。

注释：“Protein_A”定义了分子名称，这是因为这个蛋白质在PDB文件中被标定为A链，对键合近邻的排除数为3.

**注释：**离子的参数

**注释：**最后是体系级别的定义。[system]指令给出了体系的名称，在横拟中此名称将被写入到输出文件中。[molecules]指令列出了体系中的所有分子

[molecule]指令的几个关键注意点：

列出分子的顺序必须与坐标（本例中为.gro）文件中的分子顺序完全一致
对每一初待，列出的名称必须与[molcculetype]中的名称一致，而不是残基名称或其他名称
任何时候如果不能满定这些明确的要求，运行grompp程序（后面介绍）时都会产生致命错误如mismatched names，molecules not being found或者其他各种错误。

**注释：**以下定义了蛋白质中的[atoms]，信息按列始出

以上信息的解释如下：

nr:原子序号
type：原子类梨
resnr: 氨基锻残基序号
residue：氨基酸残基名
注意这里的残基在原来的PDB文件中为"LYS"，使用.rtp中的"LYS"项意味着这是质子化的残基（中性pH时的占主导）
atom：原子名称
cgnr: 电荷组序号
电得组定义了整数电荷单元，可加速计算。
charge.电荷
gtot: 为对分子总电荷的持续累加
mass.原子质量
typeB，chargeB，massB：用于自由能微扰
后面几节包括[bonds]，[pairs]，Iangles]和[dihedralsl的定义。

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

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

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

使用editconf来定义盒子：
执行代码gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic
注开以上命令将蛋白质置于盒子的中心(-c)，并且它到盒子边缘的距离至少为1.0nm（-d 1.0），盒子类型是立方体(-bt cubic).

使用solvate模块添加溶剂：
执行代码`gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top``
注释：以上使用的蛋白质构型文件（-cp）来自前面editconf步骤中的输出文件，而溶剂的构型文件（-cs）来自标准安装的GROMACS.
我们使用的spc216.gro是通用的已平衡的三位点溶剂模型。你也可以使用spc216.gro作为SPC，SPC/E或TIP3P水模型的溶剂构型，因为它们都是三位点的水模型.输出文件的名称为1AK1_solv.gro，并且我们为solvate模块指定了拓扑文件的名称（topol.top），这样它就能修改拓扑文件.注意topol.top中[molecules]的变化：

第四步添加离子

以上获得的溶液体系带有8个正电荷：topol文件[atoms]的最后一部分：gtot 8.需要添加离子使总电荷为0.
使用：genion工具
为了用gromoo产生tpr文件，我们还需要一个扩展名为mdp(molecauar dynamics parameter）的输入文件，grompp会将坐标和拓扑信息与mdp文件中设定的参数组合起来生成tpr文件

在genion命令中，我们以结构/状态文件（-s）作为输入，以.gro文件作为输出（-o）。
以拓扑文件（-p）来反映去除的水分子和增加的高子，并且定义阳离子和阴离子的名称（分别使用-pname 和-nname），告诉genion只需要添加必要的离子来中和蛋白质所带的净电荷，添加的阴离子数目为8（-nn 8）。对于genion，除了简单地中和体系所带的净电荷以外，你也可以同时指定-neutral和-conc选项来添加指定浓度的离子，关于如何使用这些选项，

使用下面的命令来产生.tpr文件
执行代码gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr
以上代码得到一个原子级别的.tpr文件，将此文件用于genion，
执行代码gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral
注释：出现提示时，选择第13组"SOL"用于嵌入离子。

第五步：能量最小化

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

调用mdrun来进行能量最小化：
执行代码gmx mdrun -v -deffnm em
注释：使用-v选项可快连看到运行结果：它使mdrun输出更多信息，这样就会在屏幕上输出每步运行的情兄-deffnm选项定义了输入文件和输出文件的名称可获得如下文件：

em.log ASCII文本日志，记录能力最小化过程
em.edr 二进制能量文件
em.trr 全精度的二进制轨迹文件
em.gro 能量最小化后的结构

有两个里要的指标来决定能量最小化是否成功：

第一个是势能（在能量最小化过程的最后输出即使你未使用-v选项。Epot应当是负值，根据体系大小和水分子的多少，大约在10510^5105-10610^6106的数量级（对水中的单个蛋白质而言）。
第二个重要的指标是力的最大值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
注释；提示输入10 0来选择势能Potential（10），并用（0）来结束输入

第六步：

第五步的能呈最小化可保证我们的初始结构在几何构型和溶剂分子取向等方面都合理

为了开始真正的动力学模拟，我们必须对蛋白质周围的溶剂和离子进行平衡.

如果我们在这时就尝试进行非限制的动力学模拟体系可能会崩溃.

原因在于我们基本上只是优化了溶剂分子自身，而没有考虑溶质.我们需要将体系置于设定的模拟温度下以确定溶质（蛋白质）的合理取向.

达到正确的温度（基于动能）之后，我们要对体系施加压力直到它达到合适的密度。

因此需要使用posre.itp文件，其目的在于对蛋白质中的重原子（非氢原子）施加位置限制（position restraining）力.这些原子不会移动，除非增加非常大的能量，位置限制的用途在于，我们可以平衡蛋白质同用的溶剂今子，而不引起蛋白质结构的变化.

平衡往往分两个阶段进行：
第一个阶段在NVT系综（粒子数，体积和温度都是恒定的）下进行.这个系综也被称为等温等容系综或正则系综.这个过程的需要的时间与体系的构成有关但在NVT系综中，体系的温度应达到预期值并慕本保持不变，如果温度仍然没有稳定，那就需要更多的时间.通常情况下，50ps到100ps就足够了，因此在本例中我们进行100ps的NVT平衡.根据机器的不同，运行可能需要一段时间.需要的.mdp参数配置文件

genvel-yes产生初始速度使用不同的随机数种子（gen_seed）会得到不同的初始速度，因此从一个相同的初始结构开始可进行多个（不同的）模拟。
tcoupl -V -rescale 速度重缩放控温器改进了Berendsen弱耦合方法，后者不能始出正确动能系综。
pcoupl - no 不使用压力耦合.

使用grompp和mdrun，像在能呈最小化过程中做的一样

执行代码gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
执行代码gmx mdrun -deffnm nvt

使用energy模块分析温度变化情况：

执行代码gmx energy -f nvt.edr -o temperature.xvg
注释：提示时输入16 0来选择体系温度并返出，得到的结果应该和下面的差不多：

上图显示：体系的温度很快就达到了目标温度（300K），并在平衡过程中后面的时间内保持稳定对于这个体系，更短的平衡时间（50ps）也足够了！

第七步：NPT平衡

第二个阶段在NPT系综**（粒子数，压力和温度都是恒定的）下进行.
前一步的NVT平衡稳定了体系的温度.在采集数据之前，我们还需要稳定体系的压力（因此还包括密度）.
压力平衡是在NPT系综下进行的，其中粒子数，压力和温度都保持不变这个系综也被称为等温等压系综，最接近实验条件**。
100 ps NPT平衡的.mdp参数配置文件：
该文件与NVT平衡时所用的参数文件没有太大不同，注意添加的压力耦合部分，其中使用了Parrinello-Rahman控压器其他几项改动如下：

continuation =yes：我们将从NVT平衡阶段开始继续进行模拟

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

注释：使用-t选项以包括NVT平衡过程中的产生的检查点文件，这个文件包含了继续模拟所需要的所有状态变量，为使用NVT过程中得到的速度我们必须包含这个文件.坐标文件(-c)是NVT模拟的最终输出文件。
执行代码gmx mdrun -deffnm npt

使用energy模块分析温度变化情况：
执行代码gmx energy -f npt.edr -o pressure.xvg
注释：提示时输入18 0来选择体系压力并退出，结果应与下图类侧：

上图显示：在100ps的平衡过程中压力值涨落很大。这并不意外。
圈中的红线为数据的移动平均值.在整个平衡过程中，压力的平均值为1.05 bar.

使用energy模块分析密度变化情况：
执行代码gmx energy -f npt.edr -o density.xvg

注释：提示时输入24 0来选择体系压力并退出，结果应与下图类似

上图显示：跟压力一样，红线是密度的移动平均值，100ps过程中密度的平均值为998.3kg m-3，比较接近实验值1000kgm-3与SPC/E水模型的值1008kg m-3.

SPC/E水模型的参数给出的密度值接近水的实验值，在整个过程中密度值都很稳定，意味着体系的压力和密度下都平衡得很好。

第八步：成品MD

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

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

我们要进行一个1ns的MD模拟，所用的参数配置文件mdp

使用grompp和mdrun，类似前面优化过程中做的一样

执行代码gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
**注释：**使用-t选项以包括NVT平衡过程中的产生的检查点文件.这个文件包含了继续模拟所需要的所有状态变量为使用NVT过程中得到的速度我们必须包含这个文件，坐标文件(-c)是NVT模拟的最终输出文件

执行代码gmx mdrun -deffnm md_0_1

使用GPU运行GROMACS
执行代码 gmx mdrun -deffnm md_0_1 -nb gpu

第九步：分析

现在已经完成了对蛋白质的模拟，我们应该来分析一下我们的体系.
使用trjconv工具来处理体系中的任何周期性，从而获得“修正”后的轨迹：

执行代码gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
再次输出1 0
注释：Select 1（“Protein”）as the group to be centered and 0（“System”）for output.
执行代码gmx rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns
再次输出4 4
注释：Choose 4（“Backbone"）for both the least-squares fit and the group for RMSD calculation.
执行代码gmx rms -s em.tpr -f md_0_1_noPBC.xtc -o rmsd_xtal.xvg -tu ns
再次输出1 1
注释：计算相对于晶体结构的RMSD.

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

上面两个图都显示出RMSD大约是0.1 nm（1A），这表示蛋白质的结构非常稳定.两图之间的微小差异意味着，当t=0ns时的蛋白质的结构与晶体结构稍有不同.这是预期结果，因为它已经进行了能量最小化，而且如我们前面讨论的，位置限制并不是100%完美的.

将初始构型与模拟后的构型进行比较，这样可以更直观地看出二者的区别：

蛋白质的回旋半径Rg可衡量其密实度，如果蛋白质的折叠很稳定，其Rg将保持一个相对稳定的值，如果蛋白质去折叠，它的Rg将随时间变化.我们来分析一下模拟的溶菌酶的回旋半径。

执行代码gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg
再次输出1
注释：Chcose group 1（Proteinmrtor analysis.

上图显示：Rg值基本不变，这预示着在温度为300K时，1ns的时间内蛋白质很稳定，处于紧密（折叠）的形式，这一结果并非意外，但说明了GROMACS具有先进的分析功能