教程B5 模拟绿色荧光蛋白

原文:Jason Swails,Dave Case, Tai-Sung Lee
译: 湖南大学-康文渊 建议移步查看许楠(浙江大学)修订的版本

绿色荧光蛋白的结构，配体分子使用球棍模型展示

介绍

这个教程主要是介绍怎样用Amber进行分子动力学模拟一个带有非标准氨基酸残基的体系。需要注意的是非标准的残基是聚合物的一部分，这不同于其他的有机小分子参数力场构建的例子，所以相对更加复杂，主要涉及电荷的计算和准备恰当的文件用于肽链中的非标准残基。

我们假设你已经运行过仅含标准氨基酸残基的基本模拟体系（同时建议熟悉antechamber Sustiva tutorial ）,所以该教程主要关注非标准氨基酸残基的参数构建。不过本教程同时包含完整的GFP（绿色荧光蛋白）分子动力学流程。

教程的大纲如下:
1.使用AMBER准备PDB文件
2.计算非标准氨基酸残基的电荷和原子类型
3.使用LEaP准备残基文库和力场参数
4.创建用于模拟的拓扑和坐标文件
5.使用生成的文件进行最小化，加热，平衡和正式模拟

让我们开始do it！

第1部分：准备PDB文件

我们开始于1EMA.pdbPDB文件，其包含GFP和CRO生色团。因此我们首先从RCSB网站下载PDB文件（或者点此下载副本）.你可以看到CRO残基位于PDB文件中的896行，标记残基编号为66（残基65和67并入残基66生成生色团，PDB文件头部有描述）。

第一步是在准备使用tleap之前给PDB文件“按摩”。为了使用tleap,Amber包括一个名为pdb4amber的脚本，可以修饰PDB文件。--help标签可以展示所有可用的设置。我们将使用--dry标签来移除晶体水，--reduce标签增加氢原子到它们最理想的位置。

1
pdb4amber -i 1EMA.pdb -o gfp.pdb --dry --reduce

需要提醒的是你要了解pdb4amber到底“做”了什么。现在输出的PDB文件应该每个残基都包含了对应的氢键（该体系中没有二硫键，如果体系中包含二硫键，pdb4amber会在正确的位置为你加上二硫键）。运行pdb4amber之后，你需要修改gfp.pdb文件，将MSE残基（硒代甲硫氨酸，selenomethionine）改为标准的MET（我们需要模拟的是野生型变体，而不是为了结晶而刻意突变的版本）。你可以使用sed或者你喜爱的编辑软件将MSE全局替换为MET。此外，你需要转换硒原子变为硫原子（名称为SD）以及将元素符号SE转为S。建议您自己运行命令，当然也可以下载我们的副本gfp.pdb文件进行比较。

第2部分:计算CRO分解电荷和原子类型

我们通过1EMA.pdb的头部文件了解到CRO是非标准氨基酸残基（实际上是三个氨基酸残基环化形成）。因此，Amber的标准氨基酸残基库不包含该残基的原子类型和电荷。单这些是分子模拟所需要的。 这一部分教程主要关注于如何制作CRO电荷和测定其原子类型。
我们将使用antechamber使用bcc电荷方案生成电荷，如果你运行过介绍中提到过的Sustiva教程方案你应该熟知（操作的方法）。需要注意的是该选择不是唯一的，并且对于所有应用也可能不是最佳方案。另外一些工具，例如R.E.D.Tools提供了更加严谨和可重复的计算电荷，但是对于这个教程，antechamber工具已足够。如果你希望使用R.E.D.，他们拥有一定数量的教程，你可以用来生成电荷和原子类型，之后可以直接查看第4部分教程.

CRO模板我们可以使用Protein Data Bank中的非常有用的”components.cif”(在此查看)，对于每个包含在PDB文件中的化合物，其都包含观察到的和假想的结构。使用这个解决了坐标问题，并且确保了原子和残基名称会和您的PDB文件相匹配。如果你在PDB中搜索“CRO”，选择配体ID的搜索结果，你将会进入一个站点，其可以下载CRO component CIF文件，CRO.cif.(可以点此下载copy)

antechamber项目将会阅读component CIF（ccif）文件，并且生成电荷和原子类型。你应该查询”Antechamber and GAFF”章节的Amber参考手册的描述和例子。当然你也可以在命令行中输入antechamber -help获取关于设置的详细信息。在这里我们告知他原子类型为Amber类型(-at amber)，因为这是一个修饰的氨基酸残基，分子和Amber原子类型参数是相似的（对于更普遍的有机小分子，通常使用gaff原子类型会更好）我们使用如下的命令来运行antechamber（通常需要几分钟的时间）:

1
antechamber -fi ccif -i CRO.cif -bk CRO -fo ac -o cro.ac -c bcc -at amber

在生成的许多其他文件中，你可以看到cro.ac文件，其有点像PDB文件，但是包含键，分解的原子电荷和原子类型参数。需要注意的是antechamber是为gaff原子类型所使用的，所以在制作Amber原子类型的时候会出现一些错误，在这个例子中，我们需要修复起始的N原子，其需要拥有酰胺氮(N)同样的类型。你可以手动简单的改变原子类型NT变为N。(或者当作练习，使用sed,awk,perl等工具修改)，你可以下载cro.ac和您的文件做比较。

现在我们有antechamber分子文件，其包含分解原子电荷和原子类型。下一部分将会讨论怎样创建残基库文件和为LEaP做准备。

提醒:你通常需要使用antechamber来创建mol2文件给LEaP阅读。但是我们随后的流程需要先将结果文件给prepgen进行阅读，所以我们需要使用prepgen能够阅读的ac格式文件。

第3部分:准备给LEaP所使用的残基库和力场参数

我们从PDB下载的“CRO”组件是一个完整的分子，包含了所有的氢原子，这是antechamber（或任何量子力学程序）计算电荷所需要的。 但是我们需要在开头和结尾处去除原子，以便制备一个“氨基酸”样残基，该氨基酸残基可以在其N-端和C-端与其他氨基酸连接。

这部分可以使用prepgen项目和”mc”(主链)文件进行准备，该文件定义哪些原子需要被移除，哪些原子是主链的一部分(例如:骨架原子)。通常需要自己创建该文件。下面的内容是该教程的例子:

1
2
3
4
5
6
7
8
9
10
11
12
13
HEAD_NAME N1
TAIL_NAME C3 
MAIN_CHAIN CA1
MAIN_CHAIN C1
MAIN_CHAIN N3
MAIN_CHAIN CA3
OMIT_NAME H2
OMIT_NAME HN11
OMIT_NAME OXT
OMIT_NAME HXT
PRE_HEAD_TYPE C
POST_TAIL_TYPE N
CHARGE 0.0

HEAD_NAME和TAIL_NAME行的原子将会联系前后氨基酸残基。MAIN_CHAIN行列出连接头和尾原子之间链中的原子,OMIT_NAME行列出了CRO残基中应该从最终结构中去除的原子，因为它们不存在于完整蛋白质中。PRE_HEAD_TYPE和POST_TAIL_TYPE行让prepgen知道周围蛋白质中哪些原子类型将用于共价连接。(需要记住的是prepgen在用于蛋白和多肽以外还可以用于聚合物的制作)CHARGE行设定氨基酸残基的总电荷，prepgen将确保“删除的”原子电荷在剩余的原子之间重新分配，以便总电荷是正确的。(在这个例子中为0)

之后，你可以对prepgen使用-help标签来查看该项目所有可用的设置。随后运行prepgen将会除去N端和C端不需要的原子，假设你的主链文件命名为cro.mc:

1
prepgen -i cro.ac -o cro.prepin -m cro.mc -rn CRO

运行该命令之后你将会拥有cro.prepin文件，包含定义的CRO残基（可以从此下载相互比较）

此刻，我们拥有包含我们修饰的氨基酸残基的原子电荷的残基库。接下来我们需要检查其共价参数(bonds,angles和diherals)，查看他们是否可用。parmchk2项目将会描绘出哪些参数是需要的，并且在标准文件中寻找。如果没有找到，将会进行训练过的猜测，并且放入新参数到文件中，我们命名为”frcmod.cro”。

我们使用如下命令运行parmchk2(同样的可以使用-help标签查看可用的设置):

1
2
parmchk2 -i cro.prepin -f prepi -o frcmod.cro -a Y \
         -p $AMBERHOME/dat/leap/parm/parm10.dat

这里我们指定parm10.dat文件，因为它是我们打算使用的ff14SB力场的主要参数数据库，（这个信息可以在您打算使用的力场的leaprc文件中找到）。 如果我们删除这个标签，那么parmchk2将匹配gaff数据库中的参数，这并不是我们想要的。我们还要求打印所有的参数（即使是在parm10.dat中完美匹配的参数），原因很快就会明了。

完成这一步之后，查看frcmod.cro（你可以比较我们的frcmod.cro副本）。 您应该马上看到了含有标记为ATTN的参数的行的问题需要修改。 这表示parmchk2在parm10.dat数据库中找不到适当的相似参数。 还有许多其他的参数被选择为具有高惩罚（这表明parmchk2的结果是不合适的）。

一个简单的解决办法是通过设置parm10.dat,使用parmchk2搜寻gaff.dat来“填充”不存在或者高惩罚的参数。所以我们需要删除frcmod.cro中标为ATTN:need revision的参数，并且告诉parmchk2在gaff.dat(默认设置)中搜寻。来自parm10.dat内的参数我们想要保留，而其他参数使用gaff参数填充到frcmod1.cro文件中（这也是为什么我们需要使用-a Y标签打印所有参数）。完成这些设置的命令如下:

1
2
grep -v "ATTN" frcmod.cro > frcmod1.cro # Strip out ATTN lines
parmchk2 -i cro.prepin -f prepi -o frcmod2.cro

至此，我们拥有两个frcmod文件，一个参数来自于Amber参数数据库(frcmod1.cro,可以再此下载副本进行比较),另一个参数为gaff原子类型(frcmod2.cro,可以再此下载副本进行比较),在这两个参数之中拥有所有我们需要的内容。如果你想，你可以从frcmod2.cro种提取7个缺失参数（和附加的高惩罚参数），将其添加到frcmod1.cro中。然而，如果你按照教程顺序进行操作，这些参数会在下一部分加载，不需要进行这一步操作。然我们继续下面的步骤来创建prmtop和inpcrd文件。

与往常一样，这里生成的参数 - 尤其是由gaff提供的参数，应被视为一个起点，应该根据可用的实验或高水平QM数据进行验证。 为了本教程的目的，我们将简单地继续我们在这里生成的内容。

第4部分:创建用于模拟的拓扑和坐标文件

我们现在有了我们需要的所有文件来创建用于sander或pmemd的拓扑和坐标文件！我们只需要将这些文件加载到LEaP来创建这些文件。 如果你使用R.E.D. 或其他一些工具来做电荷推导和原子类型，欢迎回到教程。

对于这个例子，我们使用ff14SB力场和igb=8隐性水模型，描述可以查看Amber参考手册。将默认PBRadii设置为“mbondi3”集是指定力场的一部分。

我们之后加载之前创建的cro.prepin文件。在这种情况下，我们需要首先加载frcmod2.cro文件，随后是frcmod1.cro文件，确保gaff参数会被我们更加想用的ff14SB参数所覆盖。然后我们加入我们准备的1EMA PDB文件（为gfp.pdb），并保存这个结果。脚本展示如下，将其保存为tleap.in:

1
2
3
4
5
6
7
8
source leaprc.protein.ff14SB
set default PBRadii mbondi3
loadAmberPrep cro.prepin
loadAmberParams frcmod2.cro
loadAmberParams frcmod1.cro
x = loadPDB gfp.pdb
saveAmberParm x gfp.parm7 gfp.rst7
quit

我们可以通过如下命令运行tleap:

1
tleap -f tleap.in

这一步骤以后，你应该有完整的拓扑和坐标文件，你已经准备开始模拟啦！你可以下载我们创建的副本和你的进行比较.进入部分5！

第5部分：模拟；最小化，加热，平衡和正式（模拟）

由于本教程的目的是通过对修饰聚合物“link”进行参数化，本部分仅作简要介绍。 在你自己的项目中，你当然可以自由选择显式溶剂模型和比我们在这里使用的更小心的最小化，加热和平衡程序（也许利用位置限制来防止结构扭曲）。

最小化

我们使用的最小化输入文件如下。我们将该文件命名为min.in:

1
2
3
4
 simple generalized Born minimization script
&cntrl
  imin=1, ntb=0, maxcyc=100, ntpr=10, cut=1000., igb=8, 
/

我们可以使用sander模块运行最小化:

1
sander -O -i min.in -p gfp.parm7 -c gfp.rst7 -o min1.out -r min1.rst7

如同往常一样，我们建议可视化查看结果结构(mini1.rst7)来确保观察没有糟糕的事情发生，并且输出文件确保所有的看起来都正常（例如:结构完整，在最小化期间总能量最大梯度稳步下降）。我们创建的输出文件可以作为压缩包的一部分进行下载，压缩包包含本节最后的计算过程中生成的大部分文件。

加热

下面显示了我们用于加热的输入文件。 我们将这个文件命名为heat.in，并且将在从10K到300K的200ps的过程中线性地改变目标温度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Implicit solvent initial heating mdin
 &cntrl
   imin=0, irest=0, ntx=1,
   ntpr=1000, ntwx=1000, nstlim=100000,
   dt=0.002, ntt=3, tempi=10,
   temp0=300, gamma_ln=1.0, ig=-1,
   ntp=0, ntc=2, ntf=2, cut=1000,
   ntb=0, igb=8, ioutfm=1, nmropt=1,
 /
 &wt
   TYPE='TEMP0', ISTEP1=1, ISTEP2=100000,
   VALUE1=10.0, VALUE2=300.0,
 /
 &wt TYPE='END' /

需要注意的是，我们加热结构时，我们需要使用上一步（最小化）生成的结构来加热。所用使用sander来进行加热操作的命令看起来如下:

1
2
sander -O -i heat.in -p gfp.parm7 -c min1.rst7 -o heat.mdout \
       -x heat.nc -r heat.rst7

值得注意的是使用sander可能需要很长世间，我们使用pmemd.cuda运行我们的模拟，在我们的GTX680上运行能够更加快速的完成。与往常一样，使用您最喜爱的可视化工具检查所得结构和轨迹，以确保没有发生任何明显的不良事件。和以前一样，我们生成的文件将在本节末尾的压缩包中提供。

正式模拟

！！！这一段翻译的不是很好
我们成功的加热了我们的结构！现在我们已经准备开始正式的模拟。请注意，许多人称之为“平衡”的过程实际上只是在分析过程中忽略的正式模拟的一部分，因为要么就是稳定约束的结构，要么就是让系统向平衡位置迁移。为了本教程的目的，我们不区分这两个阶段，我们使用相同的输入对待两者。

我们将使用的输入文件如下，命名为md.in:

1
2
3
4
5
6
7
8
9
Implicit solvent molecular dynamics
 &cntrl
   imin=0, irest=1, ntx=5,
   ntpr=1000, ntwx=1000, nstlim=500000,
   dt=0.002, ntt=3, tempi=300,
   temp0=300, gamma_ln=1.0, ig=-1,
   ntp=0, ntc=2, ntf=2, cut=1000,
   ntb=0, igb=8, ioutfm=1,
 /

你可以使用随后的命令来运行模拟(再一次，我们将使用pmemd.cuda):

1
2
sander -O -i md.in -p gfp.parm7 -c heat.rst7 -o md1.mdout \
       -x md1.nc -r md1.rst7

所有的都完成了！当然，你仍然需要分析你的模拟，以便测试通过首先运行计算得到的结果。现在你已经知道怎样进行一个新的参数化，修饰包含寡聚链的单体单元（例如修饰核苷酸或者氨基酸残基）。这个侧链可以应用你想应用到的任何寡聚单元。

正如我们的承诺，你可以下载我们生成的版本。值得注意的是你的输出绝大部分时候和我们的都会不同（甚至整体参数都会不同，因为1ns的轨迹实在太短）