• Ross Walker 原始文档

• 2018-05-16 19:56:42 翻译: 张少波; 校对: 李继存

初始线型结构

色氨酸肽笼的NMR结构(1L2Y)

本教程旨在作为一个案例研究, 向您展示如何重现以下文章中讨论的工作:

Simmerling, C., Strockbine, B., Roitberg, A.E., J. Am. Chem. Soc., 2002, 124, 11258-11259 DOI

建议您在开始本教程之前阅读这篇文章. 对于氨基酸序列和其他相关信息, 可参考上文中引用的文献.

本教程由本介绍页面和七个小节组成:

1. 建立结构

2. 创建prmtop和inpcrd文件

3. 能量最小化结构

4. 加热系统

5. 成品模拟

6. 分析结果

7. 测试NMR结构的稳定性

警告: 本教程中的一些计算可能需要运行很长时间. 我使用了一台1.3GHz共享内存的SGI Altix机器中的16个CPU来运行计算, 总共运行了27小时. 尽管我强烈建议您运行这些模拟以熟悉SANDER, 但我仍然提供相关的输出文件, 这样即使您没有足够的计算资源, 仍然可以按照教程进行操作.

警告2: 按教程运行这些模拟并不能保证将得到与我完全相同的结果. 机器架构的差异会导致计算中的舍入误差, 致使不同机器上的不同模拟探索的相空间区域不同, 但是平均性质应该是具有可比性的. 正如您在进行本教程时将要看到的那样, 没有办法能保证您可以重现论文中描述的结果. 事实上, 即使我们设定了与论文几乎完全相同的模拟, 我们也不能得到他们发表的结果. 这可能是因为我们的模拟时间不够长, 或者陷入了局域能量最低点. 无论如何, 本教程展现了一些有趣的结果.

背景

论文描述了肽折叠模拟, 使用全原子经典模拟和略微修改的AMBER FF99力场. “Trpcage”是由华盛顿大学的Andersen课题组优化的含20个残基的氨基酸序列. 它是目前能显示两种不同折叠状态性质的最小的蛋白质, 且在室温下稳定. 这种蛋白质的尺寸小, 使其成为计算折叠模拟的理想模型. 初始折叠模拟是在实验测定结构工作前完成的, 因此预测并没有参考实验. 当实验结构测定完成后, 发现与预测结构的RMSD在1.4埃范围内, 对于从没有限制的直链线性序列开始的模拟来说这是非常好的结果.

在本教程中, 我们将尝试重现论文获得的结果. 所使用的模拟流程非常相似, 尽管由于时间限制, 我们只能运行一个约50 ns的模拟. 这应足以重现折叠结果, 但尽管如果要得到 新 结果, 我们需要延长模拟并进行其他类似的模拟以验证我们的结果以及预测结构的稳定性. 事先声明, 由于模拟的长度不同, 机器和/或处理器数量的不同, 模拟结果会有差异. 这是分子动力学的原理决定的, 执行顺序的小变化以及浮点计算中的舍入意味着不同机器采样的轨迹将随时间发散. 这不是误差或漏洞, 也不代表一个模拟比另一个更 正确, 只是简单地说明两个模拟正在探索相空间的不同区域. 如果我们对结果进行平均, 而且我们已经运行了足够长的模拟, 那么两台机器应该给出相同的结果, 只不过它们是通过稍微不同的路径达到的. 因此, 我们很难在本教程中精确地重现论文的结果, 但是我们可以尝试重新创造重要的结果, 也就是使用AMBER预测20个氨基酸多肽的折叠结构.

所以, 记住上面的话, 让我们开始吧.

步骤1: 建立初始结构.

在之前的所有教程中, 我们或者有一个可用的晶体结构, 或者体系足够小, 可以手动建立结构. 在本教程中, 结构太复杂无法手动建立, 而且也没有可用的pdb. 因此我们需要一种方法来构建线性肽链. 所幸LEaP包含了一个叫做sequence命令的工具.

论文给出的Trpcage的氨基酸序列为:

NLYIQWLKDGGPSSGRPPPS

它是以单个字母符号表示的. 在将其输入LEaP中之前, 我们需要把它转换成sequence命令能识别的标准3字母标记形式. 转换表如下:

所以序列可以重写为:

ASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO SER

这还没有完成, 因为LEaP不会自动终止链的两端, 所以我们必须在前面加上N来指定N端, 加上C来指定C端. 因此在LEaP中使用的序列是:

NASN LEU TYR ILE GLN TRP LEU LYS ASP GLY GLY PRO SER SER GLY ARG PRO PRO PRO CSER

这会得到一个合理的线性链结构, 我们将使用它开始折叠模拟. 所以让我们开始吧. 启动xleap, 我们将使用FF99力场:

    $AMBERHOME/bin/xleap -s -f $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99

(在使用XLEAP时谨记关闭数字锁定)

然后使用sequence命令来构建我们的结构(查找关于该命令用法请输入: help sequence). 注意: 下面的命令已被分成三行, 以便于阅读. 在xleap中, 你应该把它作为一行输入.

>TC5b = sequence { NASN LEU TYR ILE GLN TRP LEU LYS					ASP GLY GLY PRO SER SER GLY ARG					PRO PRO PRO CSER }

该结构会创建为TC5b, 我们可以使用edit命令来查看它.

>edit TC5b

可以看出我们得到了一个很好的线性初始结构. 由于存在许多空间重合, 所以在开始分子动力学之前, 我们需要对这个结构进行段时间的能量最小化. 让我们将这个结构保存为库文件以便将来使用. 我们将其保存为pdb文件, 以便在VMD中查看.

>saveoff TC5b TC5b_linear.lib>savepdb TC5b TC5b_linear.pdb

下载文件 TC5b_linear.lib TC5b_linear.pdb

步骤2: 创建prmtop和inpcrd文件.

现在我们已经有了结构, 可以创建prmtop和inpcrd文件了. 在开始之前, 我们应该确保使用的参数和模拟设置与论文中的相同. 论文第3段写道:

我们只使用Trpcage TC5b2氨基酸序列(N20LYIQWLKDGGPSSGRPPPS39)开始模拟, 扩展的初始构象由AMBER 6.0.4的LEaP模块构建, 所有的分子动力学(MD)模拟都完全没有使用限制, 并使用SANDER模块进行正则系综模拟, 我们修改该模块以提高用于所有计算的Linux/Intel PC群集的性能. 模拟采用ff99力场^5, 但重新拟合了[phi/psi]二面角参数^6(见补充材料)以改善与丙氨酸四肽构象的从头算相对能量^7的一致性. 参数并没有拟合trpcage的数据. 使用广义Born模型^8考虑溶剂效应, 与AMBER一样^9.

所以他们首先在xleap中建立了一个线性结构, 像我们刚刚完成的那样. 然后他们在恒温(正则系综)隐式广义Born溶剂中进行不受限制的MD模拟. AMBER支持许多不同的广义Born模型, 在本教程写作时最先进的是由A. Onnfriev, D. Bashford和D. A. Case修改的模型, 使用了II型半径(IGB = 5), 如在AMBER手册的GB部分引用的论文中所述. 在论文中, 他们没有指定使用的具体的GB模型, 但很可能是因为他们当时只能使用IGB = 1. 为了尽量保证我们的模拟设定与原文相近, 我们也使用IGB = 1. LEaP使用的默认半径与IGB = 1兼容, 因此我们不需要指定要使用的半径.

他们还指出, 他们使用了FF99力场, 像我们一样, 但修改了phi/psi二面角参数. 这些参数是对FF99中过度支持α-螺旋结构的phi/psi参数的校正. 为了重现论文中的结果, 我们在构建prmtop文件时必须包含这些更改. 不幸的是, 论文的补充信息中并未明确指出哪些参数已更改, 他们只是提供了修改后的parm99.dat文件. 我相信这是由于AMBER 6没有内置FF99力场, 它是后期加入的. 因此当时有着包含不同错误的几个力场版本存在. 通过在论文中重复整个parm99.dat文件, 作者确保人们看到他们使用的是官方的FF99力场. 不幸的是, ACS以PDF格式提供这些数据, 而不是文本文件, 这使得提取信息格外困难. 幸运的是, 这些修改已经包含在AMBER8版本$AMBERHOME/dat/leap/parm/目录中的frcmod.mod_phipsi.1中. 在Amber9中, 这些修正已被删除, 因为它会导致甘氨酸参数化的问题. 它已经基本上被深入的重新参数化的FF99所取代, 被称为FF99SB. 如果你正在考虑自己的模拟, 建议使用FF99SB代替mod_phipsi.1修改的版本. 然而, 我们的目标是重现JACS论文中的模拟, 其作者使用了mod_phipsi.1, 所以我们也会使用这个版本. 对于Amber 9用户, 我在下面列出了该文件的链接.

frcmod.mod_phipsi.1

from Simmerling, Strockbine, Roitberg, JACS 124:11258, 2002. Modifies parm99.
MASS

BOND

ANGL

DIHEDRAL
N -CT-C -N 1 0.700 180.000 -1.
N -CT-C -N 1 1.100 180.000 2.
C -N -CT-C 1 1.000 0.000 1.

NONB

如上所示, 只有修改了3个二面角参数. 所以我们可以简单地把这个文件加载到xleap中, 它会覆盖相关的参数. 因此, 如果您关闭了xleap, 请重新打开它并将结构加载回来:

$AMBERHOME/bin/xleap -s -f $AMBERHOME/dat/leap/cmd/oldff/leaprc.ff99>loadoff TC5b_linear.lib

现在加载更改后的二面角参数:

>loadamberparams frcmod.mod_phipsi.1

我们现在可以保存prmtop和inpcrd文件:

>saveamberparm TC5b TC5b.prmtop TC5b.inpcrd

这是得到的文件: TC5b.prmtop, TC5b.inpcrd

步骤3: 能量最小化结构

在开始模拟之前, 我们需要对初始结构进行短时间的能量最小化. 这不会导致我们的系统折叠, 因为最小化会保持在起始结构附近的局部最小值. 它可以去除结构的不合理之处, 修正氢原子位置等, 这样开始模拟时体系可以稳定.

这是输入文件:

1
2
3
4
5
6

Stage 1 - minimisation of TC5b
 &cntrl  imin=1, maxcyc=1000, ncyc=500,  cut=999., rgbmax=999.,igb=1, ntb=0,  ntpr=100
 /

一共进行1000步的能量最小化；前500步用最陡下降法(ncyc=500), 后500步用共轭梯度法(maxcyc-ncyc), 这足以清理初始结构. 注意所用截断值很大(cut = 999), 原因在于我们用的是非周期性模拟(ntb = 0), 所以无法使用PME方法得到无限的静电作用. 使用PME时, 推荐的截断值为8埃, 因为这只能截短极短距离的范德华相互作用. 但是, 如果没有PME, VDW和静电相互作用都会在截断处被截断, 所以截断值应该尽可能地大. 不幸的是, 模拟所花费的时间与截断值的平方成正比-关于模拟时间与截断距离大小的关系, 请参见教程1第5.1.2节. 所幸这个系统足够小, 我们可以设置无截断. 999埃的截断值足够大, 以至于在模拟过程中不会截断任何非键相互作用. 同样, 将rgbmax设置为999埃, 这个参数控制了计算有效Born半径时将要考虑的原子对之间的最大距离. 与截断值一样, 这个值越大越好, 但该值越大模拟时间越长. 由于这个系统只有20个残基的大小, 可以把所有的原子对包括在有效的Born半径计算中, 因此设定的rgbmax值远远大于系统的最大范围.

进行最小化:

$AMBERHOME/bin/sander -O -i min1.in -o min1.out -p TC5b.prmtop						 <span style="color:#666">-c</span> TC5b.inpcrd -r min1.rst

输入文件: TC5b.prmtop, TC5b.inpcrd, min1.in

输出文件: min1.out, min1.rst

在16个1.3GHz SGI Altix的处理器上, 这需要运行大约3.5秒.

使用Amber9或以上的用户注意: 随着AMBER 9的发布, 作为AMBER版本一部分的高性能MD代码PMEMD现在能够运行广义Born模拟以及PME模拟. 因此, 如果您使用的是AMBER 9或更高版本并安装了PMEMD(在AMBER 11之前您必须单独构建它), 则可以使用PMEMD代替SANDER来运行这些模拟. 在大型并行系统和集群上, 这应该比SANDER运行得更快, 同时产生相同的结果.

为了比较最终的最小化结构和输入结构, 请执行以下操作:

$AMBERHOME/bin/ambpdb -p TC5b.prmtop < min1.rst > min1.pdb

用看图软件例如VMD打开这两个文件(min1.pdb和TC5b_linear.pdb)

原始结构以蓝色显示, 而最小化结构以黄色显示. 如图所示, 最小化并没有太多地改变骨架结构, 但它相当大地移动了色氨酸和酪氨酸残基以消除空间冲突. 在分子动力学模拟开始进行的时候, 这些高能热点将会引发问题. 如果你不相信的话, 尝试用我们创建的inpcrd文件运行MD并观察会发生什么.

步骤4: 加热体系

下一阶段是对这个体系进行MD. 我们首先分7个阶段将体系缓慢加热超过50 ps. 分阶段对体系在每个温度下均衡加热可以降低系统崩溃的几率. 另一种方法是使用权重限制来控制加热. 更多信息请参阅AMBER手册. 通常我们会在300 K(室温)下运行MD模拟. 但是, 我们应该参考他们在文章中提出的建议:

在300 K下进行了100 ns的MD模拟, 但是所有的模拟在这个时间尺度上都受到动力学圈闭, 对初始条件有强相关性, 并且不能收敛成类似的构象组合. 因此, 我们将温度升高到325 K.

文中发现有必要将系统加热至325 K以避免动力学陷入局部最小值. 因此我们也应该这样做. 在高温运行时, 需要仔细考虑时间步长, 因为系统中的振荡会更加明显. 这意味着如果以2 fs的时间步长在600 K下运行, 那这个步长对于在300 K下的动态系统也是足够的, 时间步长过长会导致不稳定. 所幸325 K非常接近300 K, 只要对设计氢的键进行限制, 我们就可以采用2 fs的时间步长. 但是, 如果稍后将系统加热到400 K, 可能需要将时间步长减少到1.5 fs左右.

由于我们的初始结构是手动构建的而不是实验晶体结构, 因此在初始阶段跟实验结构相比更不稳定. 为了使初始系统的能量释放得到控制, 在加热阶段采用很小的时间步长0.5 fs, 然后在平衡阶段改为2 fs. 这可能是矫枉过正, 但安全第一. 我们也将坐标(mdcrd)写入频率设置为50(ntwx = 50). 这是一个非常高的写入频率, 并且可能大大降低模拟的效率, 但在加热的最初阶段, 系统很不稳定, 尤其当所用模型为手动建造的结构时. 把坐标写入步长设置减小意味着当出现问题时, 我们能够知道发生了什么并确定问题来自哪里. 当进行正式平衡MD时, 更合适步长为每500到1000步记录一次. 该系统只有304个原子, 所以每个坐标集都非常小, 因此可以在成品阶段每隔500步输入一次而不会影响计算性能. 但如果这是一个在大型群集上并行运行的大型显式溶剂模拟, 那么除非出于某种原因需要更频繁的采样, 否则应该每隔2到5 ps左右写一次坐标. 对于2 fs的时间步长相当于ntwx设置在1000到2500之间.

因此, 加热系统的参数如下:

Step 1 - 10,000 steps, 0.5fs time step (5ps), initial temp = 0.0K, target temp = 50.0K, temperature coupling constant = 1.0ps
Step 2 - 10,000 steps, 0.5fs time step (5ps), target temp = 100.0K, temperature coupling constant = 1.0ps
Step 3 - 10,000 steps, 0.5fs time step (5ps), target temp = 150.0K, temperature coupling constant = 1.0ps
Step 4 - 10,000 steps, 0.5fs time step (5ps), target temp = 200.0K, temperature coupling constant = 1.0ps
Step 5 - 10,000 steps, 0.5fs time step (5ps), target temp = 250.0K, temperature coupling constant = 1.0ps
Step 6 - 10,000 steps, 0.5fs time step (5ps), target temp = 300.0K, temperature coupling constant = 1.0ps
Step 7 - 40,000 steps, 0.5fs time step (20ps), target temp = 325.0K, temperature coupling constant = 1.0ps

这是第一个输入文件:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

Stage 1 heating of TC5b 0 to 50K
 &cntrl  imin=0, irest=0, ntx=1,  nstlim=10000, dt=0.0005,  ntc=2, ntf=2,  ntt=1, tautp=1.0,  tempi=0.0, temp0=50.0,  ntpr=50, ntwx=50,  ntb=0, igb=1,  cut=999.,rgbmax=999.
 /

其他6个输入文件非常相似, 只是相关参数有所改变

在此处下载: (heat2.in, heat3.in, heat4.in, heat5.in, heat6.in, heat7.in)

警告: 在本教程的示例中, 我们不会更改用于随机数生成器的初始随机值. 这由命名列表变量ig控制. 这主要是针对教程设置中结果的可重复性问题. 但是, 在运行成品模拟时, 特别是当用ntt = 2或3(Anderson或Langevin热浴)时, 必须在每次重新启动MD时将初始随机数的默认值更改. 如果使用AMBER 10(bugfix.26或更高版本)或AMBER 11或更高版本, 则可以通过在cntrl命名列表中设置ig = -1自动执行此操作. 否则, 在每次重新开始计算时需要为ig指定一个正的随机数. 有关不执行该操作的缺点的更多细节, 请参阅以下文章:

Cerutti DS, Duke, B., et al., “A Vulnerability in Popular Molecular Dynamics Packages Concerning Langevin and Andersen Dynamics”, JCTC, 2008, 4, 1669-1680, DOI

以下是用来运行模拟的PBS脚本, 可以根据自己的系统进行修改:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

#PBS -l ncpus=16#PBS -l walltime=500:00:00#PBS -l cput=2000:00:00#PBS -j oesetenv AMBERHOME /usr/people/rcw/amber9cd ~rcw/initial_heating
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat1.in -p TC5b.prmtop -c min1.rst -r heat1.rst -o heat1.out -x heat1.mdcrd
gzip -9 heat1.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat2.in -p TC5b.prmtop -c heat1.rst -r heat2.rst -o heat2.out -x heat2.mdcrd
gzip -9 heat2.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat3.in -p TC5b.prmtop -c heat2.rst -r heat3.rst -o heat3.out -x heat3.mdcrd
gzip -9 heat3.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat4.in -p TC5b.prmtop -c heat3.rst -r heat4.rst -o heat4.out -x heat4.mdcrd
gzip -9 heat4.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat5.in -p TC5b.prmtop -c heat4.rst -r heat5.rst -o heat5.out -x heat5.mdcrd
gzip -9 heat5.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat6.in -p TC5b.prmtop -c heat5.rst -r heat6.rst -o heat6.out -x heat6.mdcrd
gzip -9 heat6.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i heat7.in -p TC5b.prmtop -c heat6.rst -r heat7.rst -o heat7.out -x heat7.mdcrd
gzip -9 heat7.mdcrdecho "DONE"

这7个模拟在16个1.3GHz SGI Altix处理器上运行大约需要7分钟. 以下是输出文件, 可以单独下载或打包成tar.gz文件下载:

把mdcrd文件加载到VMD中(先解压)可以观察到加热过程中的变化. 可以看到系统开始折叠. 我们并不担心在这个阶段抽样的结构, 只需要确认结构看起来没问题, 没有不正常的事情发生.

以下是heat7.mdcrd的最终结构:

我们开始看到一些α螺旋形成, 但在开始看到稳定的折叠结构之前, 还有很多步骤.

步骤5: 成品模拟

本教程模拟部分的最后一步是在325 K下进行长时间的成品模拟. 原文在325 K下运行50 ns的模拟, 但只给出前20 ns的结果, 表明其余30 ns与模拟的前一部分没有显着变化. 尽管作者表示系统在5 ns到20 ns之间折叠:

两次独立的模拟在5-20 ns后收敛到基本相同的结构族.

50 ns的模拟分为10个阶段, 每个阶段5 ns. 分10个独立的阶段运行的原因是, 如果我们遇到系统崩溃, 则可以更容易地恢复模拟. 同时这也使得输出文件和mdcrd文件大小合适. 10个阶段的模拟使用相同的输入文件, 如下所示:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

Stage 2 equilibration 1 0-5ns
 &cntrl  imin=0, irest=1, ntx=5,  nstlim=2500000, dt=0.002,  ntc=2, ntf=2,  ntt=1, tautp=0.5,  tempi=325.0, temp0=325.0,  ntpr=500, ntwx=500,  ntb=0, igb=1,  cut=999.,rgbmax=999.
 /

每个阶段将包括2,500,000步(nstlim), 时间步长2 fs(dt), 共5 ns. 注意整段模拟都设定了SHAKE(ntc = 2, ntf = 2), 控温使用Berendsen热浴(ntt = 1), 系统加热好并稳定后, 热浴参数改为0.5 ps(tautp = 0.5)以保持系统接近325 K. 和文中一样, 目标温度设为325 K, 每500步记录一次输出文件和mdcrd文件. 过高的写入频率会产生过大的文件, 每5 ns写入一次产生的mdcrd文件大小为35 MB, 每500步记录一次对于本研究而言刚好合适.

以下是10个模拟的PBS脚本, 可根据自己的研究体系进行修改:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

#PBS -l ncpus=16#PBS -l walltime=500:00:00#PBS -l cput=8000:00:00#PBS -j oesetenv AMBERHOME /usr/people/rcw/amber9cd ~rcw/production
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c heat7.rst -r equil1.rst -o equil1.out -x equil1.mdcrd
gzip -9 equil1.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil1.rst -r equil2.rst -o equil2.out -x equil2.mdcrd
gzip -9 equil2.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil2.rst -r equil3.rst -o equil3.out -x equil3.mdcrd
gzip -9 equil3.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil3.rst -r equil4.rst -o equil4.out -x equil4.mdcrd
gzip -9 equil4.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil4.rst -r equil5.rst -o equil5.out -x equil5.mdcrd
gzip -9 equil5.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil5.rst -r equil6.rst -o equil6.out -x equil6.mdcrd
gzip -9 equil6.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil6.rst -r equil7.rst -o equil7.out -x equil7.mdcrd
gzip -9 equil7.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil7.rst -r equil8.rst -o equil8.out -x equil8.mdcrd
gzip -9 equil8.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil8.rst -r equil9.rst -o equil9.out -x equil9.mdcrd
gzip -9 equil9.mdcrd
mpirun -np 16 $AMBERHOME/bin/sander -O -i equil.in -p TC5b.prmtop -c equil9.rst -r equil10.rst -o equil10.out -x equil10.mdcrd
gzip -9 equil10.mdcrdecho "DONE"

在16个1.3GHz SGI Altix的CPU上, 10个阶段的模拟总共需运行约27小时. 以下是输出文件, 可单独或者打包下载(注意每个mdcrd压缩文件大小超过13Mb)

把mdcrd文件加载到VMD中(先解压, 导入文件需要一小段时间, 也许以后VMD的开发者会支持读取压缩文件), 可以看到完整的轨迹. 采用ribbon样式的效果最好, 以下是以residue type上色的结构效果图:

步骤6: 结果分析

首先根据输出文件对温度, 总能量, 动能和势能进行作图. 检查这些性质以确认系统在模拟过程中没有问题. 我们可以检查温度是否平稳升高到325 K并稳定平衡. 同样的, 动能和势能也应该有同样的趋势. 温度和能量图中的任何突变都说明模拟出现了问题: 有问题的初始结构, 过长的步长, 不合适的参数等.

为了从输出文件中提取数据, 可以用以下的perl脚本自动进行(process_mdout.perl). 用以下的命令行可以一次性地处理输出文件:

mkdir analysiscd analysis../process_mdout.perl ../heat1.out ../heat2.out ../heat3.out ../heat4.out ../heat5.out    ../heat6.out ../heat7.out ../equil1.out ../equil2.out ../equil3.out ../equil4.out ../equil5.out ../equil6.out ../equil7.out ../equil8.out ../equil9.out ../equil10.out

以下是通过process_mdout.perl获得的所有文件压缩包 analysis.tar.gz(2.3 Mb)

我们可以用一些作图软件, 如教程1中的xmgr来对温度和能量进行作图.

首先是温度:

接近升温阶段的温度

如我们所望, 温度在325 K左右平衡. 加热阶段温度上升均匀, 没有突变, 看起来没有问题.

再看能量图. 势能变化图最重要因为需要找出能量最低点的结构, 才能重现文章中RMSD图. 在能量最低点的结构将被标记为 折叠 状态.

总能量(黑), 势能(红), 动能(绿)

由上温度变化是稳定的, 而温度作为动能的一部分, 可以料想到, 动能也是稳定的. 总能量和势能开始有所下降, 最后稳定下来说明结构已经折叠成比初始的线性结构更稳定的状态. 单从势能图上能更明显地看到能量的下降趋势. 也可以像文中那样对10 ps以后的平均势能作图, 现象更加明显. 注意, 在xmgr中添加running average可以通过Data->Transformations->Running averages操作.

势能 (红: 10 ps滑动平均值)

这是论文中的Figure 1(A). 注意我们的初始能量比文中要高, 这可能是由于广义Born项的差异, 因为文章中用的是Amber6, 我们用的是Amber8. Amber7和8与Amber6相比在Born半径上有细微的差别. 也可能是我们的结构折叠状态与文中不同. 后续的分析将进一步解答我们的疑问.

能量最初急剧增加, 这是由于对系统进行加热, 但随着系统折叠能量又开始减少. 总体而言在50 ns内大约减少了40 kcal/mol, 这与文中第4段中的说法是一致的.

下一步是定位能量最低的结构, 我们称之为折叠状态, 并将用其来重复文中Figure 1 B. 可以通过查看process_mdout.perl生成的summary.EPTOT文件来找到能量最低点. 但首先我们需要去除开始50 ps的数据, 因为这是加热阶段. 以下是删除了该部分的文件(summary_EPTOT_stripped.dat.gz). 下面的awk脚本可以很快地找到文件中的最低值. 它会在每次找到新的最小值时输出, 最后输出的值就是最低的能量.

>gunzip summary_EPTOT_stripped.dat.gz>cat summary_EPTOT_stripped.dat | awk '{if($2<min) {min=$2;print $1"   "min}}'

51.000   -453.0616
52.000   -458.9413
53.000   -466.4404
57.000   -471.1558
61.000   -478.0281
71.000   -479.0216
72.000   -482.4780
74.000   -488.0607
122.000   -494.2667
188.000   -506.0259
224.000   -511.8178
510.000   -519.0695
1278.000   -521.5923
1292.000   -526.4179
3903.000   -528.1458
14896.000   -534.8027
16909.000   -537.3546
17745.000   -539.2624
17759.000   -539.3218
19120.000   -539.6930
19719.000   -539.9273
21316.000   -545.5322
26834.000   -545.9479
42505.000   -547.4913

注意: 如果你模拟时, 机器的精度变化不同, 计算的空间区域稍微不同, 那么你的结果可能会有所不同.

因此, 最低能量出现在模拟的42.505 ns, 为-547.4913 kcal/mol, 通过grep可以找到其在输出文件中的位置:

>grep 42505.000 *.out

equil9.out: NSTEP =  1227500   TIME(PS) =   42505.000  TEMP(K) =   367.02  PRESS =     0.0

这出现在第9步平衡中的第1227500步. 由于平衡过程中每500步(ntwx=500)在mdcrd文件中写入一次, 所以这表示的是equil9.mdcrd中的第2455幀(1227500/500). 用ptraj来找出这个结构.

这是ptraj脚本(extract_frame9_2455.trajin):

trajin equil9.mdcrd.gz 2455 2455trajout lowest_energy_struct.pdb pdb$AMBERHOME/bin/ptraj TC5b.prmtop ＜ extract_frame9_2455.trajin ＞ extract_frame9_2455.out

输入文件: extract_frame9_2455.trajin, equil9.mdcrd.gz, TC5b.prmtop

输出文件: extract_frame9_2455.out, lowest_energy_struct.pdb.2455

检查输出文件确保无误.

在VMD中查看该结构的pdb文件.

vmd lowest_energy_struct.pdb.2455

结构看起来很好, 系统很好地与色氨酸折叠成笼状结构, 与论文中Figure 2所示的笼状结构非常相似. 在这一点上, 我们可以将这种结构与晶体结构进行比较, 但当这篇论文发表时, 实验结构尚未确定, 因此我们应该将这一步骤留到最后. 首先, 我们将尝试生成figure 1B, 其显示轨迹结构与最低能量结构的骨架之间的RMSd. 以下是ptraj脚本:

trajin equil1.mdcrd.gztrajin equil2.mdcrd.gztrajin equil3.mdcrd.gztrajin equil4.mdcrd.gztrajin equil5.mdcrd.gztrajin equil6.mdcrd.gztrajin equil7.mdcrd.gztrajin equil8.mdcrd.gztrajin equil9.mdcrd.gztrajin equil10.mdcrd.gzreference lowest_energy_struct.pdb.2455rms reference out rmsd_to_lowest_energy_struct.dat @N,CA,C time 1.0$AMBERHOME/bin/ptraj TC5b.prmtop ＜ measure_rmsd.trajin ＞measure_rmsd.out

输入文件: measure_rmsd.trajin, TC5b.prmtop, lowest_energy_struct.pdb.2455, 上表中的所有mdcrd文件

输出文件: measure_rmsd.out, rmsd_to_lowest_energy_struct.dat

检查ptraj中的输出是否报错. 唯一的报错显示每个轨迹的50001幀都损坏了. 这不是真正的报错因为每个mdcrd文件只有5000幀, 所以我们可以放心地忽略它.

这是rmsd图:

我们的rmsd图比论文中的图波动更加明显, 尽管在16到38ns间也出现了平台. 然而, 并没有如论文第5段中所说的50 ns后变得稳定. 这可能说明了很多问题. 有可能是我们的轨迹没有折叠到原始结构, 而是随机结构, 该结构虽然很稳定, 但在38 ns之后仍会展开. 确实, 如果你在vmd中看equil8.mdcrd, equil9.mdcrd和equil10.mdcrd文件, 发现其非α螺旋部分仍在反复展开和折叠, 系统还在寻找更稳定的结构. 与论文中结果不同的原因还有可能是我们采用了隐式溶剂模型, 或者其他参数的细微差别. 数据结果表明我们的模拟的时间不够长, 系统仍在折叠与展开. 我们采用的初始结构与论文中的结构有细微的差别, 所以他们的结构快速折叠, 而我们的结构困在高的能量低点.

也有可能是我们使用了稍微不同的Born半径(文中使用AMBER 6 Born半径)可能延长了折叠过程. 有很多方法可以进一步延长模拟, 使其正确折叠. 我们可以简单地运行一个更长的轨迹, 等待系统脱离被困的高能盆地, 这可能需要很长时间. 或者可以尝试模拟退火方法, 即将系统加热后再让它冷却. 也可以尝试副本交换模拟, 即有多个模拟副本, 对它们之间的属性互换. 这避免了陷入高能陷阱的问题. 从作者文章给出的结构, 可以确认它比我们的最低能量结构的能量更低. 因此, 很不幸我们的模拟不够长, 系统无法完全折叠. 如果你愿意, 可以尝试延长模拟, 或用上文我提到的其他方法. 这里有一个教训值得学习, 模拟永远不会 太 长. 仅仅因为系统在50 ns内没有折叠并不意味着它在下一个50 ns内不会折叠, 所以当你正在考虑你结果有什么意义时要注意这一点.

尽管我们的系统似乎没有达到完全稳定的结构, 但它至少可以折叠成一个稳定时间超过20 ns的结构. 由于没有时间继续尝试不同的模拟或延长模拟, 我们将用现有的数据继续研究, 对50 ns的轨迹数据进行分析.

首先, 把我们的最低能量结构与现有的NMR结构模型(1L2Y.pdb)进行比较. 像所有核磁共振结构一样, 该pdb文件包含所有可能的结构, 这里共有38个. 理想情况下, 我们应该对这些结构进行平均, 并将平均结构与我们的理论结构进行比较. 但是, 为了提高速度, 我将跳过这一步, 只将文件中的第一个结构与能量最低结构进行比较. 用VMD叠合两个结构并计算RMSD. 如下图所示, 这是合理的, 因为NMR结构3到19号残基是非常保守的, 38个结构只有链的末端之间差别很大.

因此我们只对高度保守的结构进行计算.

注意: 我们使用的是VMD v1.8.3 - 如果你使用不同的版本, 你的菜单可能会略有不同. 请注意, 需要v1.8.2或更高版本才能进行RMSD计算.

以下是含有第一个NMR结构的pdb文件(nmr_struc_1.pdb).

• 第一步: 载入理论结构

• 第二步: 加载第一个NMR 结构 (确认选中New Molecule)

此时结构不会正确地重合, 这是因为坐标系是不同的.

为了比较两个结构, 我们需要进行RMSD拟合并合并他们到一起. 所幸VMD内置的工具让这操作起来很方便. 你可以通过点击Extensions->Analysis->RMSD Calculator找到:

现在可以进行叠合. 我们只对骨架进行叠合, 所以要确保选中Backbone only. 通过在输入框中输入residue 2 to 18选择残基3到19(记住VMD里编号是从0开始的). 然后点击Align对结构进行叠合, 然后点击 RMSD计算两个结构之间的rmsd.

计算可得骨架的 RMSD为3.72埃. 这与论文中的1.4埃相去甚远, 几乎说明我们的结构根本没有折叠到论文中相同的状态. 以下是VMD中显示的叠合状态. (我改变了图像显示模式使其更容易分辨, NMR结构用蓝色表示):

确实, 比较结果看起来很不理想. 色氨酸看起来在非常不同的位置, 与实验结构的位置几乎相差了90°. 但这还没完, 我们对3到19号残基(也就是VMD里的2 to 18号)进行拟合. 如果重复这个过程并且仅拟合第3到9个α螺旋所涉及的残基, 会发生什么情况?

再次打开RMSD计算工具, 这次输入residue 2 to 8, 然后点击align, 再点击RMSD:

我们得到α螺旋骨架的RMSD仅为0.7 埃. 确实这次叠合得非常好. 问题出现在9号残基, 它的骨架扭转方向错了, 从这开始结构变得不同. 有一点需要注意, 仔细观察两个结构链的末端:

理论结构

NMR结构

很明显蛋白链的末端氢键作用有所不同. 我们的结构有2个, 可能3个氢键把链的末端紧密结合. 在NMR结构中氢键作用大不相同. 所以这可能是两个结构不同的一个原因, 可能是这些氢键使得我们的的结构困在能量陷阱中. 让我们来看下如何用AMBER的ptraj通过轨迹文件来看氢键. 这也可以阐明为什么我们的初始结构在轨迹早期形成并保持稳定15 ns, 却没有真正的二级结构(无α螺旋). 通过轨迹你可以看到这个结构, 在大约15到16 ns它突然断开, 使得系统该位置上α螺旋重新折叠. 同样让我们来学习如何用MMTSB工具包对轨迹进行簇分析. 这将使我们能够随着时间变化将结构叠合成不同的类, 以便分析他们的维持时间.

6.1) 用Ptraj进一步分析轨迹

6.1.1) 转换MDCRD文件为Binpos格式

首先把压缩文件格式的mdcrd文件转为binpos格式. 这是轨迹文件的二进制版本. 用mdcrd文件当然可以很容易地进进行后续的步骤, 但用binpos会明显更快. 同时我们把十个mdcrd 文件串成一个单独的binpos文件. 以下是所需的ptraj 脚本:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

mdcrd_to_binpos.ptraj

trajin equil1.mdcrd.gz
trajin equil2.mdcrd.gz
trajin equil3.mdcrd.gz
trajin equil4.mdcrd.gz
trajin equil5.mdcrd.gz
trajin equil6.mdcrd.gz
trajin equil7.mdcrd.gz
trajin equil8.mdcrd.gz
trajin equil9.mdcrd.gz
trajin equil10.mdcrd.gz

trajout equil_1-10.binpos binpos

>$AMBERHOME/bin/ptraj TC5b.prmtop ＜mdcrd_to_binpos.ptraj ＞mdcrd_to_binpos.out

如果成功的话你的工作目录会有一个 174mb 大小的equil_1-10.binpos文件.

6.1.2) 计算平均结构

接下来我们计算整个轨迹的平均结构. (除去加热阶段)

1
2
3

trajin equil_1-10.binpos

average equil_average.pdb pdb

>$AMBERHOME/bin/ptraj TC5b.prmtop ＜average_crd.ptraj

这会生成文件equil_average.pdb

太可怕了… 这说明了什么?

请记住平均结构就是整个模拟轨迹的的坐标平均. 系统在旋转平动, 所以不用指望从这个平均值中获取有用信息. 这不是一个有物理意义的结构. 在做坐标平均时请记住这点.

但是, 如果我们更好地计算平均值, 我们可以获得更有意义的结果. 首先, 我们可以对轨迹中的每一帧相对第一帧的骨架进行质量加权的RMS拟合, 这将消除结构的旋转和平移. 其次, 我们可以看看轨迹的特定部分, 其平均结构可能有意义. 尝试查看28到35 ns之间的部分, 这是轨迹文件的28000到35000帧.

1
2
3

trajin equil_1-10.binpos 28000 35000
rms first mass @C,CA,N
average equil_average_28-35.pdb pdb

>$AMBERHOME/bin/ptraj TC5b.prmtop ＜average_28-35_crd.ptraj

这会生成文件equil_average_28-35.pdb

这比原来的结构更好, 更直观. 尝试在VMD中自己打开并探索该结构. 首先加载TC5b.prmtop文件, 然后再加载pdb文件, 这样成键更合理. 注意一些残基非常小, 特别是一些氢键长度很短. 这表明在这段轨迹期间这些残基是非常活跃的, 发生很多转动. 如果旋转某个键周围的东西, 则平均坐标是该旋转中心的一个点. 但是, 请注意色氨酸看起来几乎完美. 在这段轨迹它肯定不怎么移动. 同样地, 它的骨架形成得很好表明结构的折叠部分在28和35ns之间很稳定, 与RMSD图结果一致.

6.1.3) 氢键作用分析

接下来研究在轨迹当中如何用ptraj来观察氢键的出现.

以下是用到的ptraj 脚本:

trajin equil_1-10.binpos#-- Donors from standard amino acidsdonor mask :GLN@OE1donor mask :GLN@NE2donor mask :ASN@OD1donor mask :ASN@ND2donor mask :TYR@OHdonor mask :ASP@OD1donor mask :ASP@OD2donor mask :GLU@OE1donor mask :GLU@OE2donor mask :SER@OGdonor mask :THR@OG1donor mask :HIS@ND1donor mask :HIE@ND1donor mask :HID@NE2#-- Acceptors from standard amino acidsacceptor mask  :ASN@ND2 :ASN@HD21acceptor mask  :ASN@ND2 :ASN@HD22acceptor mask  :TYR@OH  :TYR@HHacceptor mask  :GLN@NE2 :GLN@HE21acceptor mask  :GLN@NE2 :GLN@HE22acceptor mask  :TRP@NE1 :TRP@HE1acceptor mask  :LYS@NZ  :LYS@HZ1acceptor mask  :LYS@NZ  :LYS@HZ2acceptor mask  :LYS@NZ  :LYS@HZ3acceptor mask  :SER@OG  :SER@HGacceptor mask  :THR@OG1 :THR@HG1acceptor mask  :ARG@NH2 :ARG@HH21acceptor mask  :ARG@NH2 :ARG@HH22acceptor mask  :ARG@NH1 :ARG@HH11acceptor mask  :ARG@NH1 :ARG@HH12acceptor mask  :ARG@NE  :ARG@HEacceptor mask  :HIS@NE2 :HIS@HE2acceptor mask  :HIE@NE2 :HIE@HE2acceptor mask  :HID@ND1 :HID@HD1acceptor mask  :HIP@ND1,NE2 :HIP@HE2,HD1#-- Backbone donors and acceptors for this particular molecule#   N-H for prolines do not exist so are not in the mask#donor mask @Oacceptor mask  :2-11,13-16,20@N :2-20@H#Terminal residues have different atom namesdonor mask @OXTacceptor mask :1@N :1@H1acceptor mask :1@N :1@H2acceptor mask :1@N :1@H3#-- series hbt is just a placeholder to ensure we get the full analysis. If you don't#have the word series you don't get a full analysis.hbond print .05 series hbt

这是个相当复杂的 ptraj 脚本, 所以需要解释一下. 第一部分列出了标准氨基酸侧链氢键的供体和受体. 如果你想找出所有的氢键的话这部分在任何ptraj氢键脚本中都必须有该部分. 注意, 只要在这里定义的残基氢键才会被找到. 因此如果你的系统中有一个辅酶, 你需要在这个列表中手动添加它的供体和受体.

下一部分(#-- Backbone donors...) 列出了能作为受体和供体的骨架原子. 这里列出的是除了前面列出的特定氨基酸的供体和受体. 骨架原子很棘手, 原因有很多. 首先, 脯氨酸不具有质子化的主链氮:

脯氨酸残基和相邻的骨架

所以必须从潜在受体中排除. 其次, N端和C端具有不同的原子命名规则, 并且在C端还具有额外的原子, 更多的质子和OXT原子. 每个特定系统都需要考虑这些问题. 如以下命令

acceptor mask  :2-11,13-16,20@N :2-20@H

需要根据你的系统中的残基来调整.

最后一行实际上是用来告诉ptraj我们需要什么:

hbond print .05 series hbt

hbond表示 “做氢键分析”.

print .05表示”只输出持续超过5%模拟时间的氢键”.

series hbt“让ptraj 输出与时间相关的信息. 换句话说, 我们可以得到轨迹特定氢键存在的一个很好的图形. hbt 只是这个系列的名字, 在这里只是一个临时的代称”

警告: 氢键分析可能会非常耗费内存. 例如, 对不含水的20个氨基酸的肽链进行这种分析, 需要622MB内存. 在具有1 GB RAM的P4 3.0GHz计算机上运行大约需要1分钟. 如果你的内存小于1 GB, 可能需要更长的时间

>$AMBERHOME/bin/ptraj TC5b.prmtop ＜analyse_hbond.ptraj ＞analyse_hbond.out

这会生成文件analyse_hbond.out

在查看实际结果之前, 请始终检查以确保您的选项真正起作用, 并给出你期望的原子数. 要特别注意任何被忽略的选项, 例如:

PTRAJ: donor mask :HIS@ND1
Mask [:HIS@ND1] represents 0 atoms !!!NO ATOMS DETECTED!!!
WARNING in ptraj, donor: No atoms selected (:HIS@ND1), ignoring...

在这种情况下警告是无害的, 因为我们的TRPcage结构中没有组氨酸残基.

如果确信事情按照预期运行, 我们可以看看结果. 以下是输出文件的最后的部分, 在我们的案例中, 其内容如下:

  Dumping schematic of time series after each h-bond, key follows:
   |          .       -       o      x      *      @    |
	  0-5%   5-20%  20-40%  40-60% 60-80% 80-95% 95-100% occupancy

		DONOR         ACCEPTORH      ACCEPTOR
  atom# :res@atom   atom# :res@atom atom# :res@atom %occupied  distance       angle              lifetime maxocc
|   197  :12@O   |   228   :16@H     227   :16@N   |  25.83  2.871 ( 0.08)  22.38 (12.37)      2.2 (  1.9)     30 |   oooo-..|
|   197  :12@O   |   221   :15@H     220   :15@N   |  25.68  2.889 ( 0.08)  28.88 (11.51)      1.5 (  0.9)      9 |--.------.|
|    56   :3@O   |   118    :7@H     117    :7@N   |  25.11  2.887 ( 0.08)  21.18 (11.42)      1.7 (  1.2)     14 |  .---o---|
|   167   :9@OD2 |   248   :16@HH22  246   :16@NH2 |  23.84  2.818 ( 0.09)  25.29 (10.08)      3.9 (  3.8)     35 |   ..xo.o-|
|   167   :9@OD2 |   245   :16@HH12  243   :16@NH1 |  21.40  2.825 ( 0.09)  27.16 ( 9.58)      3.1 (  2.9)     32 |   ..x-.--|
|   116   :6@O   |   171   :10@H     170   :10@N   |  21.21  2.880 ( 0.08)  27.43 (13.53)      1.7 (  1.4)     28 |  .-------|
|   250  :16@O   |   104    :6@HE1   103    :6@NE1 |  19.92  2.850 ( 0.09)  22.71 (11.50)      2.9 (  2.4)     19 |   -o-x-  |
|    92   :5@O   |   159    :9@H     158    :9@N   |  17.78  2.887 ( 0.08)  24.55 (12.88)      1.5 (  0.9)     10 |  .-----..|
|   166   :9@OD1 |   248   :16@HH22  246   :16@NH2 |  16.96  2.813 ( 0.10)  25.17 (10.07)      4.1 (  4.1)     52 |   .o .o-.|
|    35   :2@O   |    77    :5@H      76    :5@N   |  16.38  2.897 ( 0.07)  31.11 (12.86)      1.4 (  0.9)     21 |  .--..--.|
|   166   :9@OD1 |   245   :16@HH12  243   :16@NH1 |  14.69  2.824 ( 0.10)  27.38 ( 9.72)      3.0 (  2.8)     28 |   .- .o..|
|    35   :2@O   |    94    :6@H      93    :6@N   |  14.46  2.888 ( 0.08)  19.03 (11.03)      1.6 (  1.1)     11 |  .-----  |
|   116   :6@O   |   178   :11@H     177   :11@N   |  14.42  2.885 ( 0.08)  28.88 (14.13)      1.5 (  0.9)     11 |   .-..--.|
|   226  :15@O   |    94    :6@H      93    :6@N   |  10.16  2.871 ( 0.08)  27.90 (13.91)      2.0 (  1.5)     14 |o-.       |
|   183  :11@O   |   241   :16@HE    240   :16@NE  |   8.95  2.876 ( 0.08)  27.59 (12.03)      1.7 (  1.5)     15 |   ..-....|
|   116   :6@O   |   247   :16@HH21  246   :16@NH2 |   8.70  2.853 ( 0.09)  42.07 (12.09)      1.6 (  1.1)     14 |   ..-....|
|    75   :4@O   |   137    :8@H     136    :8@N   |   7.83  2.901 ( 0.07)  21.14 (11.68)      1.2 (  0.6)      8 |  ........|
|   183  :11@O   |   247   :16@HH21  246   :16@NH2 |   7.52  2.860 ( 0.09)  32.51 (11.97)      1.7 (  1.6)     25 |   ... ...|
|   176  :10@O   |   104    :6@HE1   103    :6@NE1 |   7.43  2.868 ( 0.08)  24.14 (11.70)      1.6 (  1.0)     14 |--.       |
|    11   :1@OD1 |    18    :2@H      17    :2@N   |   7.02  2.856 ( 0.09)  44.01 (10.18)      1.3 (  0.7)      9 |....... ..|
|    56   :3@O   |    94    :6@H      93    :6@N   |   6.21  2.900 ( 0.07)  36.78 (13.51)      1.3 (  0.6)      8 |  .. .  .-|
|   135   :7@O   |   171   :10@H     170   :10@N   |   5.71  2.906 ( 0.07)  32.21 (12.11)      1.2 (  0.5)      6 |...       |
|   264  :17@O   |   241   :16@HE    240   :16@NE  |   5.61  2.861 ( 0.09)  26.87 (11.39)      2.2 (  1.8)     25 |  -.    . |
|   264  :17@O   |   104    :6@HE1   103    :6@NE1 |   5.50  2.854 ( 0.09)  29.23 (13.35)      2.2 (  2.2)     25 |    ...  .|
 ---------------- --------------------------------- -------------------------------------

这张表很容易看懂. 一个有趣的部分是在表格的最右侧显示轨迹文件每10次(5纳秒)的占用率. 您马上可以发现在轨迹的前15 ns(30%)内几乎没有氢键. 在第一部分确实存在的那些氢键通常不出现在轨迹的其余部分. 例如仅在轨迹的前15 ns期间存在的两个氢键对, 即原子对226-94(GLY15:O-TRP6:H)和176-104(GLY10:O-TRP6:HE). 在VMD中可将其可视化. 启动VMD并加载TC5b.prmtop和equil_1-10.binpos. 如果关闭分子显示(双击VMD Main中的文件名D), 并且最小化VMD, 则加载速度会更快. 然后转到Graphical Representations并选择以下内容:

转到轨迹中的第一个结构. 然后进入OpenGL显示模式, 按2切换到bond label模式. 然后点击两对H-O原子, 以便标记它们的距离.

现在播放整个轨迹, 看看会发生什么. 你可能需要将步长更改10使得轨迹以合理的速度播放. 在整个轨迹中逐步移动. 这是轨迹的MPEG动图(TRPcage_h_bonds.mpg 14.5 MB). 注意这两个氢键如何在模拟的前15 ns迅速形成并保持. 这两个键阻碍了α螺旋的形成. 最终它们断开, 然后α螺旋非常快就形成. 在此之后, 这两个氢键再没有形成. α螺旋是形成的第一个主要二级结构, 这与论文中描述的一致, 除了他们比我们更早看到螺旋形成以外. 在我们的模拟中, 上面讨论的两个氢键的形成足以使我们困入高能陷阱长达很多纳秒, 使得整个折叠延迟.

让我们看看在轨迹的最后3 ns期间的一些其他氢键. 文章中提到了两个不寻常的氢键:

在Trp25 indole NH 1和残基i + 10(Arg35)的骨架羰基之间, 以及Gly30 HN和残基i-5(Trp25)的羰基之间存在两个不寻常的分子内氢键. 折叠后, 两者在模拟中都存在很久(分别为92%和75%)

查看氢键分析, 看到我们都有这些键:

|   250  :16@O   |   104    :6@HE1   103    :6@NE1 |  19.92  2.850 ( 0.09)  22.71 (11.50)      2.9 (  2.4)     19 |   -o-x-  |
|   116   :6@O   |   178   :11@H     177   :11@N   |  14.42  2.885 ( 0.08)  28.88 (14.13)      1.5 (  0.9)     11 |   .-..--.|

然而, 它们分别仅短暂存在于整个轨迹的约20%和14%. 第一个在最后阶段消失, 第二个显然不是所有的时间都存在. 这表明这些氢键不如论文所说的结构折叠的那么稳定, 分别为92%和72%. 这里的百分比差异可能是由于我们在ptraj分析中使用了比论文中更严格的标准, 但更大可能是由于我们的结构没有正确折叠.

请注意, 我们还发现另一个氢键与论文中提到的第二个氢键一样, 是通过同一个供体形成, . 然而, 我们的受体是相邻的甘氨酸残基(10号残基). 这在轨迹的最后30 ns期间或多或少一致地出现, 可能表明色氨酸羰基位于两个甘氨酸残基之间, 而不是主要固定在第二个. 如有所需, 可以用VMD确认. 这也说明可能是由于折叠错误导致我们的色氨酸取向不同.

|   116   :6@O   |   171   :10@H     170   :10@N   |  21.21  2.880 ( 0.08)  27.43 (13.53)      1.7 (  1.4)     28 |  .-------|

6.1.4) 用 Ptraj进行二面角分析

如果将我们的理论折叠结构与NMR结构进行比较, 一个明显的差异是色氨酸的位置, 特别是在色氨酸环相对于蛋白质骨架的位置上. 虽然我们自认为建立了与论文中相同的初始结构, 但很可能实际上并没有. 尽管最终的折叠结构应独立于初始结构, 但折叠发生的时间尺度应该是相关的. 可惜在本教程中我们没有时间尝试从不同的起始结构运行几个模拟来测试它. 然而, 我们可以做的是看看涉及色氨酸的各种骨架相关的二面角. 根据这些信息, 我们可能可以对折叠途径是否可能与色氨酸的初始位置有关做出一些预测.

以函数时间用ptra观察以下4个角度:

以下是用到的ptraj 脚本:

trajin equil_1-10.binposdihedral trpphi :5@C :6@N :6@CA :6@C out trp_phi.datdihedral trppsi :6@N :6@CA :6@C :7@N out trp_psi.datdihedral trpchi1 :6@N :6@CA :6@CB :6@CG out trp_chi1.datdihedral trpchi2 :6@CA :6@CB :6@CG :6@CD2 out trp_chi2.dat

>$AMBERHOME/bin/ptraj TC5b.prmtop ＜measure_trp_angles.ptraj

生成四个输出文件: trp_phi.dat, trp_psi.dat, trp_chi1.dat, trp_chi2.dat

浏览文件并手动调整角度从-180到+180的所有部分, 使绘图整齐, 我不会详细讨论如何操作, 但在xmgrace中用一些代码就可以了:

dihedral cleanup was done using transformations: evaluate expression:
y = (y < -100) ? y=y+360 : y

以下是对四个角度作图:

这个图包含很多信息. 注意在整个模拟过程中, phi角度大致保持不变. 其次, 请注意psi在大约160度开始, 然后在约12.5 ns后急剧下降到约-50度. 这对应于α螺旋开始形成的点. 因此, 如果使用phi和psi角度相反的结构开始模拟, 可能会得到非常不同的变化. 也许这就是文章中的做法? 如果你有意且计算机资源允许, 可以尝试最后获得什么结果.

chi1和chi2的角度信息量也很大. 以下还是RMSD图, 将它与上面的chi1和chi2图进行比较:

注意chi1和chi2的角度变化似乎与RMSD的变化有关. 这表明色氨酸结构的变化对于TRPcage的每个结构样本的稳定性至关重要. 通过计算RMSD与chi1和chi2角之间的相互关系, 可以进一步探索这种相关性. 如果感兴趣可以一试. Numerical recipes一书中包含用于计算相关函数的代码示例. .

6.2) 用MMTSB Toolset进行簇分析

本教程的最后一种分析为簇分析. 在此我们将使用一种算法, 通过轨迹找到结构相同的簇. 由此建立描述每个簇的质心, 然后对轨迹中的每个结构与相应的簇间的RMSD进行计算.

为了完成该部分的教程, 需要用到MMTSB toolset. 这不是AMBER的一部分, 但可以从http://mmtsb.scripps.edu/software/mmtsbtoolset.html免费下载. 如果你在AMBER workshop上运行这个教程, 那MMTSB toolset应该已经在你的机器上安装好了. 你可以通过检查环境变量 MMTSBDIR是否已经设置好. 如:

>echo $MMTSBDIR

如果没有设置, 那么你可能需要安装MMTSB工具集. 请参阅网站上的使用说明. 首先将工具集解压到/usr/local/mmtsb之类的位置, 将MMTSBDIR设置到这个目录, 然后在你的路径中添加$MMTSBDIR/bin和$MMTSBDIR/perl. 完成后可以继续下面的工作.

警告: 不要尝试通过远程文件系统操作, 否则要等待很长时间!!! 如果你当前的工作目录位于NFS共享上, 我强烈建议你将所需的全部文件复制到本地暂存目录, 然后在完成时将结果复制回来. 另外, 避免在包含pdb文件的目录中执行ls命令, 否则会有没完没了的列表滚动浏览.

• 第一步 – 准备目录

簇分析过程会创建大量文件, 即轨迹中每一帧的pdb文件. 因此, 我们不想让这么多文件 “污染” 我们的工作目录, 所以我们要创建一些子目录来进行这个过程:

>mkdir clustering>cd clustering>mkdir PDBfit

• 第二步 – 用ptraj 从轨迹中提取RMSD’d PDB’s

trajin ../equil_1-10.binpos#-- put all the pdb frames in a subdirectorytrajout PDBfit/trp.pdb pdb#orient all frames best fit to backbone of helix in NMR structurerms first mass :2-20@CA,C,N

>$AMBERHOME/bin/ptraj ../TC5b.prmtop ＜extract_pdb.ptraj

这会在PDBfit目录产生50,000个 pdb 文件, 每一个文件代表已经被RMSD拟合到轨迹中NMR第一个结构的2到20号残基单个快照.

• 第三步 – 调整PDB文件编号

这一步严格上并不是必要的, 但它可以使得后续更方便地找到特定的pdb文件.

Ptraj不会在pdb文件名前导0, 这表示对其列出的目录进行排序时是这样的:

1
10
100
...

而不是顺序排序. 用以下c shell脚本用前导零重命名文件:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

#! /bin/cshset ff = *.pdb.1set fnam = $ff:rset numfil = `ls -1 *.pdb.???? | wc -l`if( $numfil != 0)then
  foreach fnam (*.pdb.????)
	 set fr=$fnam:r	 set fnum=$fnam:e
	 mv $fnam $fr.0$fnum
	 echo $fnam $fr.0$fnum
  end
endifset numfil = `ls -1 *.pdb.??? | wc -l`if( $numfil != 0)then
  foreach fnam (*.pdb.???)
	 set fr=$fnam:r	 set fnum=$fnam:e
	 mv $fnam $fr.00$fnum
	 echo $fnam $fr.00$fnum
  end
endifset numfil = `ls -1 *.pdb.?? | wc -l`if( $numfil != 0)then
  foreach fnam (*.pdb.??)
	 set fr=$fnam:r	 set fnum=$fnam:e
	 mv $fnam $fr.000$fnum
	 echo $fnam $fr.000$fnum
  end
endifset numfil = `ls -1 *.pdb.? | wc -l`if( $numfil != 0)then
  foreach fnam (*.pdb.?)
	 set fr=$fnam:r	 set fnum=$fnam:e
	 mv $fnam $fr.0000$fnum
	 echo $fnam $fr.0000$fnum
  end
endif

chmod +x fix_numering_pdb.csh./fix_numbering_pdb.csh

• 第四步 – 进行簇化

归簇将用到mmtsb工具kclust, 该工具需要一个包含结构列表的文件进行操作, 我把这个文件称为clustfils. 我们将轨迹中的所有结构写入该文件. 如果只想对部分轨迹进行簇分析, 当然也可以在文件中仅写入该部分结构. 让我们创建这个文件:

>cd PDBfit>ls -1 . > ../clustfils  # ls -1 means single column)>vi ../clustfils

因为kclust只能处理<50000个结构, 所以删除第一行. 拟合结构2到50,000. 记住这样做表明pdb文件实际上从二开始计数. 这意味着质心文件中的ID从一结束. 你稍后会明白我的意思. 以下是clustfils文件.

现在可以运行归簇:

>kclust -mode rmsd -centroid -cdist -heavy -lsqfit \ <span style="color:#666">-radius</span> 6 -maxerr 1 -iterate \  ../clustfils > ../Centroid_6>cd ..

归簇命令运行大约3到5分钟. 但是, 它非常占用内存, 如果没有足够的可用内存, 可能会在花费相当长的时间. 事实上, 如果内存不足, 交换空间不足, 分析将无法正常进行. 大概需要略少于1 GB内存. 在P4/3GHz/512 Gb RAM机器上, 几乎不断地翻页, 整个运行需要10分钟. 如果你的机器是分页的(像xosview这样的显示器), 那么在作业运行时, 一定不要在机器上做任何东西, 甚至移动鼠标. 这样可以避免机器在运行和分页群集作业的同时, 交换应用程序的导入和导出.

到目前为止, 这个分析的主要输出是文件Centroid_6(Centroid_6.gz 478kB, 下载后别忘了解压), 包括:

• 质心的 pdb 格式

• 每个质心的 members

• 每个 member 对其质心的rmsd距离

在文本编辑器中快速浏览Centroid_6文件, 查看文本的格式. 请注意ID和pdb编号的数差值:

97 trp.pdb.00098 5.221610

ID 97实际上是编号为98的PDB文件. 留心这一点.

• 第五步 – 处理簇

质心本身并不具有物理意义, 因为它们实际上是基于簇的member的数学结构. 然而, 对每个质心具有最低RMSD的结构是有意义的并且更容易理解. 通过这种方式, 可以查看最接近代表所找到的每个簇的质心的结构.

以下是用来处理Centroid_6文件的一个awk脚本:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

BEGIN{b0=2;}

{centind=index($1,"#Centroid");
   c=$2;   getline;centind=index($0,"#Centroid");
   FIL0 = sprintf("centroid%2.2d.member.dat",c)   while(centind != 1){	 print $1,$3 > FIL0 ;	 getline;centind=index($0,"#Centroid");
   }

   numcent=0;    print $2,$1,NR-b0;
   c=$2;   getline; endrec = index("End",$0);   while( endrec != 1 ){
	 FIL = sprintf("centroid%2.2d.pdb",c);	 print > FIL;	 getline; endrec = index($0,"#End");
   }
   b0=NR+2;
 }

>awk -f extract_centroids.awk Centroid_6 | tee Centroid_6.stats

注意: tee命令可以将输出文件传到Centroid_6.stats并同时输出到屏幕.

此命令为每个质心生成一个pdb文件, 每个质心有一个构成文件和一个包含members数量的统计文件. 所有文件的压缩包可以在此下载: centroid_data.tar.gz (326 kB)

此时可以看到一共有6个簇.

1 #Centroid 11732
2 #Centroid 6179
3 #Centroid 26772
4 #Centroid 2735
5 #Centroid 98
6 #Centroid 2483

其中两个 (1和3)非常密集, 三个较为稀松 (2,4 & 6) , 一个不太明显(5). 每个质心生成单独的pdb文件(centroid ??. pdb)和member文件(centroid ??. member.dat). 在VMD中查看pdb文件, 但记住这些只是平均结构, 不是有物理意义的构型, 它们只是簇的质心.

找到离簇中心(质心)最近的结构, 可以看到与质心在物理上更相关的结构. 在第二列上按照数字排序member文件(离质心的rms距离).

>sort -n -k2 centroid03.member | head -5

21340 1.331332
18240 1.346612
18275 1.352564
34161 1.362107
17425 1.370057

ID 21340的结构对质心3的RMSD最低. 复制该结构作为该质心的bestmember. 谨记我们的ID和PDB文件在数序上相差1, 因为我们从数据中删除了第一个pdb文件. 因此ID 21340结构是指pdb文件trp.pdb.21341:

>cp PDBfit/trp.pdb.21341 centroid03_bestmember.pdb

现在根据你的兴趣对这个 best member的 pdb进行分析, 如氢键分析.

对其他簇重复以上分析. 以下是包含6个best members的压缩文件: best_members.tar.gz

• 第六步–对距质心的RMS距离作图

可用用 xmgrace 或 xmgr打开每个质心的member文件. 需要进行一些处理使得图更清晰明了. 通过以下顺序导入每个member数据集:

(密集程度顺序):

centroid03.member.dat
centroid01.member.dat
centroid02.member.dat
centroid04.member.dat
centroid06.member.dat
centroid05.member.dat

然后在 xmgrace 如下操作 (xmgr 菜单可能稍有不同):

1. 从xmgracePlot菜单中打开Set Appearance窗口选择所有数据集

2. 设置Line properties Type 为 None

3. 设置Symbol properties 为 Type circle 并调整大小为 16

4. 点击Apply (在底部)

5. 选择每个数据集, 并更改Symbol colors如下(在选中每个color后点击Apply):

   G0.S0 black (leave it alone, it is already black) G0.S1 red G0.S2 blue G0.S3 green G0.S4 yellow G0.S5 orange

得到如下所示的图:

簇分析中出现了一些重要的现象. 请注意, 6号簇只出现在模拟非常初始的阶段, 从现在开始可以忽略它. 请注意从一个簇到另一个簇的过程, 其中蓝色是唯一的转换为黑色的结构, 除了中间短暂变为绿色. 最初的亚稳态簇是红色的, 并通过绿色, 黄色和蓝色过度为主要的簇.

通过改变归簇参数以及改变完成归簇的结构范围可以做很多其他的簇分析. 一旦选中轨迹中感兴趣的部分, 可以通过创建包含在该部分的pdb文件列表(clustfils)进行簇分析.

其次, 可以改变归簇半径, 较小的半径会产生更多(更小)的簇, 通常难以管理且各不相同. 然而, 不通过试验和错误, 无法获知合适的半径值是多少.

第三, 可以研究以不同方式叠合的pdb文件的簇. 我们的是对第一个结构的2至20号残基进行叠合. 你可以尝试对结构的其他部分进行叠合, 或其他结构如最低能量或是第一个归簇分析中其中一个最相近的member.

可以对以上生成的簇进行的另一个分析是通过ptraj找出关于该簇的其他信息. 最简单的分析是针对一些结构的rmsd, 例如簇的质心, 但我们不会得到任何新信息, 因为无论如何归簇都会生成该信息.

最好的方法是创建仅包含每个簇member的轨迹, 并对每个轨迹(如rms和hbond)进行分析. 我把这一步留给了更富有热情的科研工作者.

步骤 7: 测试NMR结构的稳定性

在我们的模拟中没有把晶体结构折叠起来, 但是, 在我们的轨迹结束时, 我们发现了一些主要的结构变化. 因此, 我们很可能只是没有给系统足够长的时间来折叠到正确的结构.

因此, 出于兴致看看如果我们从NMR结构开始模拟会发生什么. 我们将使用NMR pdb文件中的第一个结构开始模拟. 你现在应该能够自己流畅运行AMBER进行模拟. 首先, 我们从NMR结构中去除氢, 将其加载到XLeap中并保存成一个prmtop和inpcrd文件. prmtop文件应该与原始的TC5b.prmtop文件相同. 应该只有inpcrd文件不同. 文件如下: nmr_struct.inpcrd.

像运行我们的线性结构一样运行完全一样的模拟. 以下是包含所有输出文件的压缩包: nmr_struct_output.tar.gz(124 mb)

如果对比我们原先模拟的相同的温度和能量图, 我们可以看到温度保持很好, 能量图也相似.

尝试进行第一个结构轨迹的骨架RMS拟合. 这将告诉我们这与NMR结构是否有很大偏差. 我们还将计算最低的能量结构的RMS.

trajin equil1.mdcrd.gztrajin equil2.mdcrd.gztrajin equil3.mdcrd.gztrajin equil4.mdcrd.gztrajin equil5.mdcrd.gztrajin equil6.mdcrd.gztrajin equil7.mdcrd.gztrajin equil8.mdcrd.gztrajin equil9.mdcrd.gztrajin equil10.mdcrd.gzreference nmr_struct.inpcrdrms reference out nmr_sim_rmsd_to_nmr_struct.dat @N,CA,C

下图表示的是对第一个NMR结构(黑)和原始折叠结构的最低能量结构的骨架(红)的 RMSD

如果我们对phi, psi, chi1和chi2角度作图, 得到:

![http://ambermd.org/tutorials/basic/tutorial3/files/angles_new.jpg]

现在, 对上图和我们从折叠模拟中得到的表进行比较:

折叠模拟

NMR 结构模拟

为了得出结论, 可以做更多的分析. 基于你在本教程中掌握的知识, 我将这些内容留给读者自己尝试. 尝试看看氢键作用, 这是如何比较的. 同样, 尝试对NMR结构轨迹运行簇分析, 然后查看最接近各个质心的结构.

其他轨迹分析

在本教程中我们看到的系统未能在50 ns内折叠到天然结构, 这表明它在325 K时可能仍动力学被困. 这与他们在论文中陈述的内容相反, 但它们只在325 K运行了两次模拟, 所以没有足够的数据来确定这一点. 如果你感兴趣, 可以重复本教程中的一些模拟. 一些有趣的尝试包括:

1) 模拟退火

2005 PSC AMBER 研讨会的参会者(http://www.psc.edu/training/workshops/2005/Amber/)已经成功用上这种方法. 你可以把系统加热到 375 K, 该温度下运行几个ns, 然后慢慢冷却系统到 325 K. 这可以避免325 K下动力学被困的问题因为初始运行温度在 375 K 使得体系反复折叠和展开, 去除了对初始结构的关系.

2) Langevin 动力学

你可以尝试的另一个选择是使用Langevin控温(ntt = 3)代替Berendsen控温. 这种热浴是通过模拟随机碰撞而起作用, 因为溶剂中的分子可能会感受到, 而不是像Berendsen控温那样简单地调节速度. 这种方法可以更有效地平衡温度, 并且可以使相空间探测地更快. 看看这种方法是否能产生折叠结构会很有趣.

3) 副本交换 MD

本教程未涵盖此内容, 但可能会后续添加. 这种方法中, 系统的许多副本同时运行, 并且每隔一段时间交换温度. 通过这种方式可以快速寻找势能面.

还有许多其他的事情可以尝试, 比如改变控温的松弛时间(tautp)等.

也可以尝试分析以下轨迹. 在联系论文的作者后, 我收到了以下两个文件: carlos.prmtop, carlos.inpcrd. 这些是作者在论文中使用的prmtop文件, 也就是标记为折叠状态的结构. 我用这些文件进行了一些模拟并提供如下, 你有空可以对他们进行分析.

1. 用作者文中的prmtop文件, 设置tautp=0.5ps对线性初始结构进行50 ns模拟: 50ns_linear_carlos_prmtop_tautp0.5.tar.gz (137 mb)

2. 用作者文中的 prmtop文件, 设置tautp=5.0ps 对线性初始结构进行50 ns模拟: 50ns_linear_carlos_prmtop_tautp5.0.tar.gz (137 mb)

3. 用作者获得的最低能量结构开始50 ns模拟. 使用他们的prmtop文件, tautp = 0.5ps: 50ns_linear_carlos_prmtop_carlos_inpcrd_tautp0.5.tar.gz (137 mb)

步骤8: 结束语

从本教程中你应该学到很多重要的技巧. 你也应该懂得了蛋白质折叠模拟是多么棘手的事情.

希望在学习完本教程后, 你能更好地理解如何运行更高级的MD模拟, 更重要的是如何使用AMBER的工具和MMTSB工具集进行一些高级的分析.

◆本文地址: https://jerkwin.github.io/2018/05/16/AMBER基础教程B3-案例研究-稳定蛋白折叠色氨酸肽笼的全原子结构预测和折叠模拟/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆