- 2018-06-24 23:03:19 翻译: 李传舜(江南大学)
在本教程中, 我将以我的一个研究项目为例, 向你介绍使用AMBER进行分子动力学模拟的基本步骤. 我们将要模拟的分子是HIV-1整合酶的核心领域. 但在开始任务之前, 让我快速给你一些关于这种特定分子的背景信息.
整合酶是复制HIV-1病毒所需的三种必需酶之一. 它的主要功能是将病毒DNA插入宿主基因组. 蛋白质全长288个残基, 但由于其在溶液中的低溶解性, 很难通过实验的方法证明整合酶的全长结构. 幸运的是, 蛋白水解分析表明, 这种酶可分为三个独立折叠的结构域: N端结构域, C端结构域和催化核心结构域. 已经报道了每个单独结构域的实验确定的结构, 包括催化核心结构域的几个高分辨率晶体结构. 从这些晶体结构中, 我们发现在保守的活性位点附近存在无序的表面环. 由于蛋白质的表面环通常参与催化, 理解环的构象动力学对于理解蛋白质的催化机制是重要的. 在整合酶的情况下, 诱变研究表明, 环上有几个关键残基可以显着影响整合酶的催化活性. 上图显示了10纳秒分子动力学轨迹的结果. 整合酶的核心结构域用卡通形式表示来描绘, 而DNA底物则用表面形式来表示. 在表面环上, 保守残基Tyr143以蓝色, 绿色和红色突出显示出不同阶段的环的构象变化. 以表面表示描绘了停靠在结合位点中的DNA底物.
下面让我们看看如何用AMBER获得上面的信息.
概述
在本教程中, 我们将首先执行标准的MD模拟. 以下是我们将要模拟HIV-1整合酶循环运动的步骤概述:
- 输入文件准备
- 能量最小化
- 平衡数据的平衡和分析
- 成品模拟
- 分析
- 下载
1. 输入文件准备
AMBER套件由60多个程序组成. 但是, 要进行传统的分子动力学模拟, 只需要知道两个程序: tleap和sandder. 除了这两个, 你也可能会使用另外两个辅助程序carnal和ptraj. 这两个程序是AMBER套件的数据分析主力, 但如果你知道其他可以与AMBER轨迹文件一起工作的分析程序, 你也可以使用这些程序.
像所有其它电脑一样, 运行计算机模拟都是关于输入和输出的. sander是整个套件的主要模拟引擎; 它需要两个输入文件来描述要模拟的分子系统, 一个控制文件指定模拟的条件, 并根据这些信息计算经典的分子动力学轨迹. tleap是一个帮助程序, 它可以获取预定的坐标文件, 例如pdb, 并生成一个 拓扑 文件和一个 重启 文件. 当然, 我在这里仅简单介绍一下. tleap实际上不仅仅是这一点, 而且它不是AMBER中唯一可以采用pdb文件并为sander生成输入文件的程序. 另一个帮助程序叫做xleap. 它们之间的区别在于xleap使用图形界面, 而tleap使用基于文本的界面. 当我第一次使用AMBER时, 更倾向于使用图形化的xleap, 并拒绝”落后”的基于文本的tleap. 而随着时间的推移, 我逐渐认识到基于文本的界面的价值. 如果你打算做很多模拟, 花费一点额外的精力来适应基于文本的界面是值得的.
整合酶核心结构域(pdb id 1qs4)的最新x-射线结构有三种晶体学单体. 在环状区域中有几个残基丢失. 在将pdb文件提供给tleap之前, 我们必须先解决这些问题. 如何在pdb文件中建立丢失残基模型非常依赖于你的研究目的和你使用的特定模型. 如果你刚刚开始建模并不知道如何做这些事情, 可以使用南缅因州大学Gale Rhodes教授在互联网上免费提供的SwissPDB Viewer程序. 在这项研究中, 我使用另一种整合酶晶体结构(pdb id 1bis)模拟了缺失的残基. 除了缺失的残基之外, 在133和185位还有两个点突变. 这些突变主要是为了帮助蛋白质更好地形成晶体. 然而, 后来认为F185K突变会影响附近第二个表面环的动力学. 由于我们对表面环的动力学感兴趣, 所以我们应该将它从Lys改回Phe. E133W突变不在蛋白质的任何关键区域, 并且pdb文件实际上表示pdb序列和相应的序列数据库之间存在某些冲突, 因此我们可以仅保留晶体坐标. 现在, 计算机中的变异残基比试管中的变异要容易得多. 有许多免费的图形分子编辑程序可供您使用. 我使用SwissPDB Viewer来做到这一点, 因为它是免费的, 也因为它有一个内置的能量最小化设施来清理模型化的结构. 为了节省您的时间, 您可以从Download页面下载已清理的结构, 或单击此处.
好吧, 废话说多了. 让我们开始吧.
下载完文件后, 看看里面的内容. 你会注意到没有氢原子. 如果你不想使用我准备好的文件, 但想准备你自己的pdb文件, 现在一定要去掉所有的氢原子. 根据您用来构建初始模型的软件, 最终可能会得到不兼容的氢名称. 将它们剥离并让tleap添加氢原子更容易.
在命令行上, 启动tleap程序. 你应该看到这样的东西:
tleap
-I: Adding /amber/dat/leap/prep to search path.
-I: Adding /amber/dat/leap/lib to search path.
-I: Adding /amber/dat/leap/parm to search path.
-I: Adding /amber/dat/leap/cmd to search path.
Welcome to LEaP!
Sourcing leaprc: /amber/dat/leap/cmd/leaprc
Log file: ./leap.log
Loading parameters: /amber/dat/leap/parm/parmME.dat
Loading library: /amber/dat/leap/lib/all_nucleic94.lib
Loading library: /amber/dat/leap/lib/all_aminoME.lib
Loading library: /amber/dat/leap/lib/all_aminoctME.lib
Loading library: /amber/dat/leap/lib/all_aminontME.lib
Loading library: /amber/dat/leap/lib/ions94.lib
Loading library: /amber/dat/leap/lib/water.lib
>
注: 首先要加载力场文件, 这里应该是自己设置的默认值.
不要担心目录的差异. 这可以通过使用您自己的leaprc文件来指定. 您可以稍后了解如何定义您的leaprc文件. 目前, 只接受默认值.
如果您在启动时遇到问题, 请检查您的环境变量是否已正确设置. 如果您使用的是bash shell, 请确保您的~/.bash_profile中包含这两行
export AMBERHOME=/amber
export PATH=$PATH:/amber/exe
如果您是C-Shell用户, 请将以下内容添加到~/.cshrc文件中.
setenv AMBERHOME "/amber"
setenv PATH "${PATH}:/amber/exe"
其中/amber是您AMBER发行版的位置.
接下来, 我们将把pdb文件加载到tleap中, 并将其分配给一个我们称之为mol的变量
...
...
Loading library: /amber/dat/leap/lib/water.lib
> mol = loadpdb wt1mg.pdb
Loading PDB file: ./wt1mg.pdb
Unknown residue: MG number: 154 type: Terminal/last
..relaxing end constraints to try for a dbase match
-no luck
Added missing heavy atom: .R<CGLN 154>.A<OXT 18>
Creating new UNIT for residue: MG sequence: 155
Created a new atom named: MG within residue: .R<MG 155>
total atoms in file: 1189
Leap added 1192 missing atoms according to residue templates:
1 Heavy
1191 H / lone pairs
The file contained 1 atoms not in residue templates
>
如果一切正常, 你不应该从tleap得到任何错误信息. 通常情况下, 第一次尝试时可能会遇到某些错误的信息. 像上面一样, tleap解释MG是一个未知的残基. 这是因为当我们启动程序时, MG不包含在我们加载的参数库中. 如果您本地的AMBER的安装已经包含MG参数, 那么你不会看到这个问题. 如果没有, 您可以从下载页面下载参数文件的副本, 或者单击此处. 退出tleap, 将下载的文件放入当前工作目录, 再次启动tleap输入命令loadoff MG.off, 然后再次loadpdb. 你的屏幕现在应该是这样的:
Welcome to LEaP!
Sourcing leaprc: /amber/dat/leap/cmd/leaprc
Log file: ./leap.log
Loading parameters: /amber/dat/leap/parm/parmME.dat
Loading library: /amber/dat/leap/lib/all_nucleic94.lib
Loading library: /amber/dat/leap/lib/all_aminoME.lib
Loading library: /amber/dat/leap/lib/all_aminoctME.lib
Loading library: /amber/dat/leap/lib/all_aminontME.lib
Loading library: /amber/dat/leap/lib/ions94.lib
Loading library: /amber/dat/leap/lib/water.lib
> loadoff MG.off
Loading library: ./MG.off
> mol = loadpdb wt1mg.pdb
Loading PDB file: ./wt1mg.pdb
Added missing heavy atom: .R<CGLN 154>.A<OXT 18>
total atoms in file: 1189
Leap added 1192 missing atoms according to residue templates:
1 Heavy
1191 H / lone pairs
>
tleap程序非常聪明足以识别C端缺少OXT原子, 因此它会为我们自动添加. 因此这不是错误消息.
此时, 我们已经加载了pdb文件, 添加了氢原子和丢失的重原子, 并为所有原子分配了参数. 下一步计算强度将会有点大. 我们将添加一个水盒子并添加相反离子来抵消这些电荷.
使用tleap添加水和离子
现在我们将使用tleap为分子和与之电荷相反的离子添加一个水盒子来完成系统模型.
> solvateBox mol WATBOX216 10
Solute vdw bounding box: 49.995 53.684 37.063
Total bounding box for atom centers: 69.995 73.684 57.063
Solvent unit box: 18.774 18.774 18.774
Total vdw box size: 73.439 76.851 60.149 angstroms.
Volume: 339472.593 A^3
Total mass 162173.797 amu, Density 0.793 g/cc
Added 8064 residues.
>
如上所示, 添加水盒子的命令称为solvateBox. 还有其他几种添加水的方法, 但在这儿我们只使用最直接的方式. WATERBOX216是TIP3P水分子的预平衡盒子. 数字10是盒子边缘与蛋白质之间的缓冲距离(埃). 在这里, 您必须对使用的缓冲区的大小进行一些判断. 如果你使用的数字太大, 你最终会得到一个大的水盒子, 并浪费大量不必要的计算时间在无意义的水分子上. 但是, 如果你使用的水盒子太小, 在模拟过程中, 分子可能会发生构象变化, 部分分子可能会跑在盒子外面. 如果你想模拟接近实验条件, 并想要居中系统, 你可以使用这些信息来确定你需要的水盒子大小, 并明确地设置, 否则, 我认为10是一个合理的数字来启动.
接下来我们需要做的是添加相反离子. 在我们输入addions命令之前, 我们需要弄清楚我们的系统是带正电还是负电. 如果它是正电荷的, 我们会想要加上带负电的Cl-来抗衡它, 如果它带负电, 那么我们会添加Na+来对抗它. 要计算我们系统的电荷, 我们可以使用如下命令charge:
> charge mol
Total unperturbed charge: 2.00
Total perturbed charge: 2.00
>
因此我们看到我们的系统是带正电的. 我们将添加Cl-来平衡电荷.
AMBER实际上提供了两种算法来添加离子. addions中实施的第一种方法是简单地在溶质周围绘制网格, 并将离子放置在能量最低的网格点处. 这种方法将忽略水分子定位离子的位置, 如果选定的位置与水分子重叠, 则水被删除并被离子取代. 如果我们使用这种算法, 我们最终会得到Cl-离子, 而不是我们想要的Mg2+. 命令addions2实施的第二种方法除了将溶剂分子与溶质相同外, 与addions几乎完全相同. 我们将使用addions2来确保Cl-离我们的分子有一段距离, 所以它的电荷不会人为地扭曲我们的系统. 我们命令结尾处的数字0意味着我们希望找出合适数量的反离子来中和整个系统.
正如你所看到的, 这个计算在我的Pentium IV个人电脑上花费了990秒, 所以如果你的计算机上一段时间似乎没有发生任何事情, 只需在抵达重置按钮之前锻炼一点耐心.
> addIons2 mol Cl- 0
2 Cl- ions required to neutralize.
Adding 2 counter ions to "mol" using 1A grid
Grid extends from solute vdw + 2.47 to 8.47
Resolution: 1.00 Angstrom.
grid build: 4 sec
Calculating grid charges
charges: 990 sec
Placed Cl- in mol at (-18.42, 1.87, 29.98).
Placed Cl- in mol at (0.58, -17.13, 29.98).
Done adding ions.
>
最后, 我们准备输出结果并将系统保存为输入文件.
> saveAmberParm mol wt1mg.parm7 wt1mg.crd
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 468 improper torsions applied
Building H-Bond parameters.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:
res total affected
CGLN 1
NCYS 1
WAT 8064
)
(no restraints)
>
你现在可以退出tleap.
2. 能量最小化
在创建初始模型后, 最好对模型进行可视化检查以确保一切看起来合理. 您可以使用命令ambpdb将拓扑和重新启动文件转换为pdb格式:
ambpdb -p wt1mg.parm7 < wt1mg.crd > wt1mg_solvated.pdb
| New format PARM file being parsed.
| Version = 1.000 Date = 06/14/03 Time = 14:53:33
然后, 您可以使用您最喜爱的可视化软件检查pdb文件. 我再次使用SwissPDB Viewer作为示例. 当您将pdb文件加载到查看器中时, 您应该看到类似这样的内容.
在SwissPDB查看器中, 在Select菜单下有一个名为aa making clashes的选项. 该命令将突出显示与其邻居进行糟糕接触的氨基酸残基. 要查看冲突, 请从Color菜单中选择by selection. 你现在应该看到像这样的东西:
由蓝色染色的区域是接触不好的残基. 还不错, 只有三个.
在我们开始分子动力学模拟之前, 我们需要删除这些不好的接触. 原因是, 如果我们用这些不好的接触开始分子动力学研究, 那么该区域的能量将会不切实际地高, 并且会导致模拟崩溃或导致轨迹以不切实际的方向前进.
我们通过执行能量最小化来移除这些不好的接触. 即使没有明显的不好的接触, 运行一个短暂的能量最小化来放松结构也是一个好主意.
我们将分两个阶段进行能量最小化. 在第一阶段, 我们只会将水分子最小化, 并保持蛋白质和Mg2+固定. (系统不是添加的氯离子么？) 以下是我用于此计算的控制输入文件:
| bash | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 | Minimization with Cartesian restraints for the solute
&cntrl
imin=1, maxcyc=200,
ntpr=5,
ntr=1,
&end
Group input for restrained atoms
100.0
RES 1 155
END
END
|
将这些控制信息保存在一个名为min.in的文件中, 如下所示将其供给sander:
sander -O -i min.in -p wt1mg.parm7 -c wt1mg.crd -r wt1mg_min_water.rst -o wt1mg_min_water.out -ref wt1mg.rst
注意命令行末尾的-ref选项. 该选项通常不是必需的. 在这种情况下, 我们需要它, 因为我们正在做一个限制性的能量的最小化, 并且这个信息是sander程序能够定位我们在min.in中指定的残基选择所必需的. 文件wt1mg.rst与我们在准备结束时生成的原始重新启动文件是相同的文件. 只需将原始文件wt1mg.crd复制到wt1mg.rst即可. (我们为同一个文件创建另一个副本的原因仅限于文本保存, 如果需要, 可以使用相同的wt1mg.crd.)
该命令要求sander程序将min.in作为控制输入, wt1mg.parm7作为系统的参数文件, wt1mg.crd作为输入坐标. 它将写出wt1mg.rst作为重启文件, wt1mg_min_water.out作为输出文件.
计算完成后, 检查输出以查看它是否如下所示:
计算可能需要几分钟才能完成, 故这是一段很好的休息时间.
能量最小化(续)
之前的最小化计算放松了溶质分子周围的水分. 如果我们想要非常细致, 现在我们可以做相反的事情, 并且在固定水分子的同时最小化溶质, 然后一起放松整个系统. 但是, 请记住, 我们的目标只是消除不好的接触, 因此不需要过分小心. 我们将直接从整体上最小化整个系统:
| bash | |
|---|---|
1 2 3 4 5 | Minimization of the entire molecular system
&cntrl
imin=1, maxcyc=200,
ntpr=5,
&end
|
再次, 将此输入保存在一个文件中, 并将其命名为min_all.in, 然后将其提供给sander程序.
注意: 这提醒我告诉你, 在进行模拟时, 文件数量可以快速增加. 如果你习惯于保留一个单独的笔记本来记录事情, 那很好, 但是如果你和我一样, 很少写下任何东西在纸上, 那么对于你创建的文件有一个很好的命名系统是非常重要的. 拥有一个好的命名系统可以使你的工作或多或少地自我记录, 而不好的命名系统就会使你头疼数周.
sander -O -i min_all.in -p wt1mg.parm7 -c wt1mg_min_water.rst -r wt1mg_min_all.rst -o wt1mg_min_all.out
这将需要5到10分钟的时间在一台快速的Pentium IV机器上, 所以现在是时候休息一下了.
当你回来时, 不要忘记通过从重启文件中生成一个pdb文件来检查你的结果, 就像我们以前一样.
3. 平衡
在MD模拟中, 大分子和周围溶剂的原子在系统达到静止状态之前经历数十或数百皮秒的松弛. 在模拟轨迹的初始非平稳段通常在平衡性质的计算中被丢弃. MD模拟的这个阶段称为平衡阶段.
平衡的方案在很大程度上仍然是个人喜好的问题. 一些方案要求非常复杂的程序, 涉及以多步方式逐渐增加温度, 而其他更好的方法只需使用线性温度梯度并将系统加热至所需温度.
在我们的例子中, 我们将遵循AMBER手册中建议的方案并执行两阶段平衡. 在第一阶段, 我们将从100 K的低温启动系统, 并在10皮秒的模拟时间内逐渐升温至300 K. 我们将在保持体积不变的情况下执行这个平衡阶段. 输入文件如下所示:
| bash | |
|---|---|
1 2 3 4 5 6 7 8 | Heating up the system equilibration stage 1
&cntrl
nstlim=5000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=5000, ntwx=5000,
tempi =100.0, temp0=300.0, ntt=1, tautp=2.0, ig=209858,
ntb=1, ntp=0,
ntc=2, ntf=2,
nrespa=2,
&end
|
让我们仔细的一行一行的看一下.
-
第1行: 我们指定我们要运行5000步(
nstlim = 5000), 时间步长为2 fs(dt = 0.002); 我们正在开始新的运行, 所以没有先前的速度信息(ntx = 1, irest = 0), 我们将每5000步打印一次能量输出(ntwr = 5000)并每5000(ntwx = 5000)保存坐标. -
第2行: 我们希望初始温度为100K(
tempi = 100.0), 参考温度为300K(temp0 = 300.0), 使用Berendsen耦合算法保持恒定的温度. 温度耦合的时间常数为2 ps(tautp = 2.0). 对于初始速度分配, 需要随机数种子(ig = 209858). 如果你使用相同的种子, 你会得到完全相同的轨迹. 如果要运行具有不同初始速度分配的多个轨迹, 则需要使用不同的随机种子编号. (gromacs是根据进程号生成伪随机种子数) -
第3行: 我们将使用具有恒定体积(
ntb = 1)的周期边界, 并且没有压力控制(ntp = 0) -
第4行: 我们将使用SHAKE算法来仅限制涉及氢的键(
ntc = 2)并省略这些键的力评估(ntf = 2). 我们还通过重置坐标来设置更严格的容忍度 -
第5行: 为了加快计算速度, 我们还采取了更大的时间步长(
nrespa = 2)来评估力场中的缓变项. 在此设置中, 时间步长等于nrespa * dt(2 * 2 fs = 4 fs. )
注意: 只有AMBER 7和更高版本才能使用该nrespa选项. 如果您使用的是旧版本的AMBER, 请记住在您的输入中排除nrespa.
好的, 将此输入文件保存为eq_v.in. 现在让我们运行sander.
sander -O -i eq_v.in -p wt1mg.parm7 -c wt1mg_min_all.rst -r wt1mg_eq_v.rst -x wt1mg_eq_v.crd -o wt1mg_eq_v.out
这个计算在双CPU的Pentium III上花费了2个小时, 所以你可以在这里再喝一杯咖啡. 或者, 您可以从下载页面下载我的结果.
平衡(第2部分)
计算完成后, 您的目录中应包含更多文件. 我们来看看这些结果.
wt1mg_eq_v.rst: 我们下一次计算的重启文件wt1mg_eq_v.crd: 保存的轨迹文件wt1mg_eq_v.out: 包含各种能量项的输出文件以及是否有温度项.
回顾我们这部分平衡的目标是将系统温度从100K提高到300K. 为了检查结果, 我们从输出文件中提取温度信息并将其绘制在图上.
在命令行上, 执行以下操作:
grep TEMP wt1mg_eq_v.out | awk'{print $ 6, $ 9}'> temp.dat
这将从输出文件中提取时间和温度信息, 并将其保存在名为temp.dat的新文件中. 然后, 您可以将这些数据绘制在您最喜欢的图形或电子表格程序中. 我使用Excel生成以下图表. 请记住排除最后两个数据项, 因为它们是平均值和均方根偏差, 而不是实际的数据点.
从这张图中, 我们可以看到系统的温度在10 ps模拟时间内逐渐增加. 现在还不够300K, 但对我们来说, 它已经足够接近了.
接下来, 我们将使用压力和温度控制来平衡系统, 来将水的密度调整为实验值.
输入文件如下所示:
| bash | |
|---|---|
1 2 3 4 5 6 7 8 | Constant pressure constant temperature equilibration stage 2
&cntrl
nstlim=5000, dt=0.002, ntx=5, irest=1, ntpr=500, ntwr=5000, ntwx=5000,
temp0=300.0, ntt=1, tautp=2.0,
ntb=2, ntp=1,
ntc=2, ntf=2,
nrespa=1,
&end
|
请注意, 除了高亮显示的项目之外, 新的输入文件与前一个文件几乎相同. 第一个变化是反映我们现在接着前面的模拟的继续的事实, 我们将使用重启文件中的速度信息(ntx = 5, irest = 1). 我们还打开了各向同性位置缩放(ntp = 1)的恒定压力(ntb = 2), 并去除了未使用的参数(tempi和ig). 对于恒压模拟, 我们要关闭respa选项.
让我们继续开心的操作:
sander -O -i eq_pt.in -p wt1mg.parm7 -c wt1mg_eq_v.rst -r wt1mg_eq_pt.rst -x wt1mg_eq_pt.crd -o wt1mg_eq_pt.out
平衡(第3部分)
现在, 我相信你已经开始掌握它了. 再次, 我们已经使系统平衡了10 ps, 并从计算中获得了一些新的输出文件. 此时, 我们的系统应该具有约300K的温度和约1g/ml的密度.
非常好! 我们已经将我们的系统从100 K安全地提升到300 K, 并将水盒子的密度从约0.8调整到约1.0克/毫升. 为了更加肯定, 我们将延长恒温和恒压平衡时间:
| bash | |
|---|---|
1 2 3 4 5 6 7 8 | Constant pressure constant temperature equilibration stage 3
&cntrl
nstlim=50000, dt=0.002, ntx=5, irest=1, ntpr=500, ntwr=5000, ntwx=5000,
temp0=300.0, ntt=1, tautp=2.0,
ntb=2, ntp=1,
ntc=2, ntf=2,
nrespa=1,
&end
|
好吧, 我懂. 我们现在要执行100ps的平衡运行. 从这个角度来看, 除非你有一台非常快的PC或者超级计算机, 否则你可能不想坐在这里等待结果. 但是, 您仍然可以通过设置计算的运动, 并让计算运行一段时间, 以确保您所做的一切都正确. 然后, 我会向您提供结果, 以便你继续下一步.
平衡(第4部分)
按照之前的例子, 我从输出文件中提取温度并将它们组合成一个连续的图. 我们可以看到, 在初始20ps的平衡之后, 温度肯定达到了理想的300K. 对于余下的100ps的平衡, 温度保持相当恒定.
在MD的初始阶段, 如果分子构象不处于稳定, 那么势能会漂移, 因此势能通常用作决定生产数据可以被收集分析的平衡的时间. 用来帮助做出这一决定的另一种方法是相对于初始结构的RMSD. 由于势能值, 就像温度值一样, 可以直接从输出文件中提取, 故我将继续(直接)向您展示结果.
接下来, 我将向您展示如何使用程序ptraj进行RMSD计算.
平衡(第5部分)
数据分析是所有的分子动力学模拟项目中最不标准的部分. 很多时候, 人们不得不编写他们自己的程序, 以达到他们想要的方式分析轨迹. 但是对于轨迹的标准处理, carnal和ptraj在处理这些需求方面非常棒. 在本节中, 我们将介绍ptraj的一些常见用途.
1. 修整AMBER轨迹:
在这一点上, 我们已经收集了大约100 ps的MD轨迹. 您现在可以使用UIUC的VMD程序可视化轨迹. 但是当你玩播放部动画时, 你会注意到水盒会在一段时间后变得分散. 这可以通过使用ptraj重新对轨迹进行重新成像修整.
初始:
最后:
重新绘制轨迹输入:
trajin wt1mg_eq.crd
center :1-154
image center familiar
rms first out wt1mg_eq_rms.out :3-152@CA
trajout wt1mg_eq_nice.crd nobox
保存此输入并将其称为ptraj.in, 然后按如下所示运行ptraj:
ptraj wt1mg.parm7 ptraj.in
注: 最新的使用的是cpptraj程序
基本上, 这将导致ptraj读取sander创建的原始轨迹文件, 并将坐标重新定位到蛋白质质量中心, 然后重新对周期性盒子进行很好的成像并且形成正方形的盒子. 最后, 它将使用第一帧作为参考计算每帧的rmsd, 并且将叠加到C-α原子上的每个末端排除两个残基. 然后将结果叠加的帧写出到一个新的轨迹文件中, 该文件称为wt1mg_eq_nice.crd, 每帧结束时删除框信息.
在屏幕上, 你应该看到这样的东西:
....
....
....
Processing AMBER trajectory file wt1mg_eq.crd
Set 1 ..........
PTRAJ: Successfully read in 10 sets and processed 10 sets.
Dumping accumulated results (if any)
PTRAJ RMS: dumping RMSd vs time data
这里显示的RMSD图显示骨架构象在大约20ps后已经达到静止状态. 我们现在可以开始收集生产数据.
4. 成品模拟
恭喜! 经过这些艰苦的工作, 你终于到了容易的部分.
在此阶段, 您基本上使用与平衡的最后阶段相同的方案, 并且只需继续简单进行, 直到您满意或耗尽计算机时间.
再次, 输入的文件为:
| bash | |
|---|---|
1 2 3 4 5 6 7 8 | Constant pressure constant temperature production run
&cntrl
nstlim=500000, dt=0.002, ntx=5, irest=1, ntpr=500, ntwr=5000, ntwx=5000,
temp0=300.0, ntt=1, tautp=2.0,
ntb=2, ntp=1,
ntc=2, ntf=2,
nrespa=1,
&end
|
请注意, 我已将时间步数增加到500,000步. 你需要在这里做出的一个决定是你想要运行多久? 目前, 大众的工作通常需要纳秒级的模拟. 如果您打算收集几纳秒的生产数据, 那么最好先对自己进行调整并抵制一次性运行整个模拟的冲动. “提交和忘记程序”很容易. 对你的模拟进程保持警惕或许看起来有很多工作要做, 但最终它会让你免于浪费大量时间和不必要的痛苦.
5. 分析
正如我前面提到的, 分子动力学轨迹的分析是一项高度定制的工作. 预测所有可能的分析方法并压缩在一个文件当中实施是不不切实际的(尽管我不得不承认它听起来像一个有吸引力的想法). 在上一节中, 我们已经看到了如何使用ptraj来处理轨迹文件并计算RMSD.
I. 通过回溯RMSD分析确定初始弛豫时间
在本节中, 我们将重新讨论平衡问题. 最初, 我们进行了一次有教育的猜测, 并将系统平衡了100ps. 然后, 我们检查了势能和rmsd, 并确定轨迹稳定, 其数据十可以进行收集的. 在论文 J. Chem. Phys. Vol 109(23),10115, Stella等人指出使用初始结构的均方根偏差来确定初始弛豫时间的一般做法经常高估了真实值并导致浪费有价值的模拟数据. 他们提出了一种更好的方法来确定多少轨迹真的属于初始松弛阶段, 因此提供了一种方法来恢复我们最初可能抛弃的一些数据. 他们的方法包括使用许多后期阶段结构作为参考来计算rmsd值, 然后将它们与使用初始结构作为参考计算出的rmsd值进行比较. 通过向后追溯曲线, 我们可以更好地识别轨迹达到静止状态的点. 我们将把这个分析应用到我们收集的1纳秒轨迹上, 并包括第一个100 ps的恒温恒压的平衡轨迹. 我已经将轨迹文件连接成一个单独的crd文件并在这里为您重新映射. 为了减小文件大小, 我还从轨迹中去除了水和反离子. 如果你打算自己做, 不要忘记使用tleap为剥离的轨迹生成一个新的拓扑文件. 您可以通过3个简单的步骤完成此操作:
1. 采取初始拓扑并重新启动文件并从它们生成一个pdb文件
ambpdb -p wt1mg.parm7 ＜wt1mg.rst＞ wt1mg_water.pdb
2. 编辑wt1mg_water.pdb文件并删除所有水分子和Cl-离子
3. 使用tleap并加载这个没有水和离子的pdb文件并保存新的拓扑文件
tleap
> mol = loadpdb wt1mg_water.pdb
> saveAmberParm mol wt1mg_dry.parm7 wt1mg_dry.rst
> quit
用于清理轨迹的ptraj输入文件是:
trajin wt1mg_eq.crd
trajin wt1mg_md.crd
center :1-154
image center familiar
strip :WAT
strip :Cl-
rms first :3-152@CA
trajout wt1mg_dry.crd
go
您也可以从下载页面下载这些文件.
我写了一个快速而肮脏的shell脚本去做这个工作:
| csh | |
|---|---|
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 | #!/bin/csh
### take snapshots at 100ps generate reference files
cat << EOF > ptraj.in
trajin wt1mg_dry.crd 1 100 10
trajout wt1mg_ref restart
go
EOF
ptraj wt1mg_dry.parm7 ptraj.in
### calculate the rmsd for each reference structure
foreach ref (`ls wt1mg_ref*`)
set ofile=$ref.rms
cat << EOF > rms.ptraj
trajin wt1mg_dry.crd 1 100
reference $ref
rms reference out $ofile :1-154
go
EOF
ptraj wt1mg_dry.parm7 rms.ptraj
end
### clean up the files and merge them into a single file
cp wt1mg_ref.2.rms rms1
set i=10
while ($i <= 100)
cut -b12-20 wt1mg_ref.$i.rms > tmp
paste rms1 tmp > rms
cp rms rms1
@ i=$i + 10
end
rm *.rms rms1 tmp
mv rms relax.rms
|
运行上面的脚本后, 你应该得到一个名为relax.rms的输出文件. 看看下面的图, 并作出自己的结论;-)
II. B因子计算
另一个常见的分析是从MD轨迹计算B因子并将其与晶体学的B因子进行比较. 这对于ptraj来说非常简单
ptraj输入文件是:
trajin wt1mg_dry.crd
rms first out wt1mg_ca.b :1-154@CA
atomicfluct out wt1mg.bfactor @CA byatom bfactor
go
这会给你每个C-α原子的B因子. 我将结果与来自pdb文件1QS4, 链B(这是我们用作我们的初始模型的结构)的实验性B因子作图. 您可能注意到, MD模拟的B因子比晶体的B因子低得多, 除少数关键残基外. 这是可以预料的, 因为与晶体结构相比, 1-ns并不能提供足够的采样.
III. 可视化
最后, 怎么少的了一个MD动画项目?
有几个可以播放AMBER格式轨迹(PyMol, VMD等)的免费程序. 我发现VMD在将AMBER轨迹作为动画方面做得非常好, 并提供了许多用于交互式处理数据的有用工具. 我已经使用VMD在这里为您愉快的生成了一个1ns轨迹的小电影.
你可以从MD模拟中得到什么样的见解？ 就像挖掘黄金一样, 它们都隐藏在那里等着你去发现.
