第三部分 分子动力学模拟数据的分析(一)

模拟结束后就可以进行数据分析了. 这是一个重要的过程, 包括三个阶段. 首先, 有必要进行一些标准的检查, 以便对模拟质量进行评估. 如果评估结果表明模拟良好, 就可以对每个模拟进行分析, 回答预设的研究问题了. 最后, 来自不同模拟的结果可以综合起来.

注: 文件名应当能反映文件的内容, 根据模拟体系的不同而不同. 这里我们假定使用默认的文件名, 这意味着会有以下文件:

• topol.tpr: 运行输入文件, 包含模拟开始时体系的完整描述

• confout.gro: 结构文件, 包含最后一步的坐标和速度

• traj.trr: 全精度轨迹, 包括随时间变化的位置, 速度和力

• traj.xtc: 压缩轨迹, 轻量, 只包含低精度(0.001 nm)的坐标信息

• ener.edr: 随时间变化的有关能量数据

• md.log: 模拟过程的日志, 包含模拟过程中的信息

另外, 许多分析工具都能生成.xvg格式的文件. 这些文件能用xmgr或xmgrace程序查看, 也可用Python脚本xvg2ascii.py在终端以字符形式显示.

1. 结果的简略核查

在进行其他分析之前, 必须确认模拟正常结束. 有许多原因可能导致模拟中断, 特别是涉及力场或体系平衡不充分的问题, 要检查模拟是否正常结束, 可以使用gmx check程序

gmx check -f traj.xtc

确认模拟运行了50纳秒.

轨迹文件中有多少帧, 时间分辨率为多少?

模拟信息的另一个重要来源是是日志文件. 在md.log的结束部分给出了模拟过程的统计数据, 包括CPU和内存资源的使用情况以及模拟时间. 查看日志文件的结束部分. 如果你使用less, 可以按住G(shift-g)来略过前面部分, 跳到文件结束. (提示: 可以使用sed命令只抽取日志文件的最后部分: sed -ne '/NODE/,$p' md.log)

模拟实际运行了多少时间(小时), 模拟速度为多少(ns/day)? 要模拟1 s需要多少年?
势能的哪部分贡献消耗量大部分计算时间?

2. 结果可视化

现在是有趣的部分, 尽管多数分析都归结为从轨迹文件中提取图像, MD当然最重要是的体系的运动. 我们要看看轨迹.

首先使用GROMACS提供的查看器ngmx来看看. 虽然该程序的完善程度和视觉效果不及其他查看软件, 但它能够根据拓扑文件的信息绘制成键. 其他查看软件可能隐含长程键, 导致这些键被认为太长而不画出, 或者会在非常接近的原子之间画出键. 这是对模拟结果分析的一个常见错误. 使用ngmx载入拓扑和轨迹文件:

gmx ngmx -s topol.tpr -f traj.xtc

看看程序菜单, 试试不同的选项. 播放动画. 观看过程. 通过右边的选项控制. 右击或左击选择选项来改变查看.

为了可视化轨迹, 我们将从轨迹中抽取1000帧(-dt 50), 并去除水分子(当提示选择时, 选择蛋白Protein). 此外, 我们还要消除跨过盒子边界的跳跃形成连续轨迹(-pbc nojump). 要处理这些事情, 我们使用瑞士般的GROMACS工具trjconv, 它有1001个可能的组合选项. 我们使用它输出一个多模型的pdb文件, 用于PyMOL可视化.

当提示选择trjconv的输出组时, 选择蛋白Protein(组1)从而忽略溶剂与离子.(提示: 如前面所述, 可以使用echo来传递选项: echo 1 | trjconv ...)

gmx trjconv -s topol.tpr -f traj.xtc -o protein.pdb -pbc nojump -dt 50

在PyMOL中载入轨迹

pymol protein.pdb

载入所有帧后, 播放动画

mplay 

当播放动画时, 所有其他控制仍可以工作. 你可以使用鼠标旋转或缩放体系, 也可以改变分子的表示模式.

spectrum show cell

如果多肽扩散超过盒子边界, 会发生什么?

如果一切正常的话, 现在你可以看到多肽的扩散, 翻转和扭动. 但我们对内部运动更感兴趣, 而不是整体行为. 在PyMOL中, 你可以使用命令intra_fit将轨迹中的所有其他帧与第一帧对齐, 随后, 你可以使用orient命令设定聚焦点在多肽上

intra_fit protein
orient 

现在所有帧都已经叠合好了, 你可以看到多肽的一些部分比另一些部分运动得更剧烈, 这种运动的差异在后面会进行定量化分析.

当然, 使用卡通模式表示多肽显示效果会更好:

show cartoon

上面的命令会将碳骨架表示为粗管状, 而无法显示正确的二级结构元素, 因为.pdb文件中没有二级结构信息. PyMOL能够计算蛋白的二级结构, 但只计算一帧, 然后将结果应用到所有的帧. 例如, 下述命令可以计算第一帧的二级结构:

dss 

通过指定状态编号, 可以改变要计算的帧

dss state=1000

最后, 让我们同时查看所有帧, 检查多肽的柔性和刚性区域

set all_states=1

请随便摆弄PyMOL. 试着放大柔性或刚性区域, 并检查侧链的构象. 请使用ray和png命令制作图像, 即使浪费点(CPU)时间也不要紧(提示: 将场景输出为POV-Ray格式, 得到的图像可能更酷). 但记住, 如果图像的场景太过复杂, 可能会导致PyMOL的内置光线追踪器崩溃, 这种情况下, 你可以直接使用png保存屏幕上的图像.

如果有足够的机时, 你可以考虑制作一段不错的动画. 你可能已经注意到了, 轨迹的噪声很大. 这基本上是热噪声, 因此只是蛋白正常行为的一部分, 但这些噪声对制作好的动画会有影响. 我们可以滤除这些高频的运动, 只保留更慢更平滑的整体运动. 为此, 可以使用filter命令:

gmx filter -s topol.tpr -f traj.xtc -ol filtered.pdb -fit -nf 5

现在在PyMOL中载入滤过后的轨迹, 设置二级结构(dss), 显示二级结构(show cartoon), 隐藏碳骨架上的侧链(hide lines, not (name c,n,o)), 以你喜欢的颜色进行着色, 然后使用orient命令设置好视角. 现在可以开始制作动画了:

viewport 640,480
set raytrace_frames,1
mpng frame.png

退出PyMOL(quit), 检查目录下的文件(ls), 你会发现多了好些文件, 包括250张图片. 以每秒30帧的速度, 这些图片可以制作大约8秒的动画. 下载mpeg_encode程序和参数文件movie.param, 用它来生成单帧图像的动画(你可能需要编辑参数文件来改变文件名):

mpeg_encode movie.param

3. 质量保证

对轨迹进行了最初的可视化检查后, 该对模拟质量进行一些更彻底的检查了. 质量保证(QA)包括对一些热力学参数收敛程度的测试, 如温度, 压力, 势能和动能等. 更常用的含义, QA试图评价模拟是否达到平衡. 通过起始结构和平均结构的均方根偏差(RSMD)也可以对结构的收敛性进行检查. 接下来必须检查相邻周期性映像之间没有相互作用, 因为这种相互作用会导致一些非物理效应. 最后一步QA测试是计算原子的均方根波动, 并与晶体学数据的b因子进行比较.

3.1 能量项收敛

我们首先从能量文件中提取一些热力学数据, 包括温度, 压力, 势能, 动能, 晶胞体积, 密度以及盒子大小等. 这些量中的大多数已经在前面的准备步骤阶段检查过了. 能量分析使用energy命令进行, 该程序读入能量文件, 也就是模拟过程中生成的扩展名为.edr的文件. energy命令会提示需要从能量文件中抽取哪些项, 并为其生成一个图形. 使用下面的命令

gmx energy -f ener.edr -o temperature.xvg

执行命令后将列出一系列存储在.edr文件中的能量以及相关项. 本教程的能量文件中可能含有59项, 每一项都可以抽取并作图. 最前面的9项对应于力场中的不同能量项, 44以上的项列出对蛋白Protein和非蛋白Non-Protein组划分后的结果, 包括二者之间的相互作用. 键入要提取的性质的序号并回车就可以得到相应的.xvg数据文件. 要抽取温度, 键入temperature 0或12 0并回车(也可不键入0, 但需要回车两次). 使用xmgrace查看下图形, 观察温度如何围绕设定值(310 K)上下波动.

也可以从热力学性质的涨落计算体系的热容. 为此, 除体系温度外, 还必须从.edr能量文件中抽取焓(NPT系综)或总能量Etot(NVT系综)的值. 进一步, 我们必须使用-nmol选项明确指定体系中的分子数目(你可以查看拓扑文件的最后部分获知体系中分子的总数). 这样energy就可以自动计算热容并在输出部分的最后报告这个值. 更多细节请参考手册D.29.

xmgrace -nxy temperature.xvg

体系的平均温度和热容多少?

通过名称引用能量项可实现能量文件的自动处理. 使用echo和管道(|)可以将一个程序输出重定向为另一个程序的输入, 这样energy的选择可以自动完成. 要抽取多个项, 各项之间必须以\n分割,
复制粘贴或键入以下命令行来抽取其他项.

echo 13 0 | gmmx energy -f ener.edr -o pressure.xvg
echo “9\n10\n11 0” | gmx energy -f ener.edr -o energy.xvg
echo 18 0 | gmx energy -f ener.edr -o volume.xvg
echo 19 0 | gmx energy -f ener.edr -o density.xvg
echo “15\n16\n17 0” | gmx energy -f ener.edr -o box.xvg

逐个查看这些文件, 观察对应数值的收敛情况. 如果有的数值没有收敛, 这意味着模拟还没有达到热力学平衡状态, 必须延长模拟时间才能进行进一步的分析. 此外, 达到平衡前的过程是不能用于分析的. 这里, 为简单起见, 我们忽略这些考虑, 直接使用模拟的结果进行分析.

energy.xvg和box.xvg文件中给出的是什么性质?
估计压力, 体积和密度的稳定值.

一些性质收敛慢, 达到平衡所需的时间长于另一些性质. 特别的, 温度很容易收敛而体系各部分间的相互作用可能收敛较慢. 这样会导致温度已经收敛到平衡值, 而体系不同部分之间的相互作用仍需要更长时间进行平衡. 查看多肽与溶剂之间的相互作用能:

echo “48\n50 0” | gmx energy -f ener.edr -o coulomb-inter.xvg
echo “49\n51 0” | gmx energy -f ener.edr -o vanderwaals-inter.xvg

coulomb-inter.xvg和vanderwaals-inter.xvg文件中给出的是什么量?

3.2 周期性映像间的最小距离

质量保证中最重要的检查事项之一就是确保周期性映像之间没有直接的相互作用. 由于周期性映像是全同的, 其间的相互作用是物理上不应该发生的自相互作用, 会导致模拟结果不合理. 设想具有偶极矩的蛋白会有直接的相互作用. 那么同一蛋白处于盒子边界的两个末端之间就会产生吸引, 这将影响蛋白质的自身的行为并导致模拟结果不合理. 我们通过计算每一时刻周期性映像之间的最小距离来证实这种情况不会发生. 这通过mindist来实现:

gmx mindist -f traj.xtc -s topol.tpr -od minimal-periodic-distance.xvg -pi

周期性映像之间的最小距离多大, 何时出现?
在你的模拟过程中, 用于长程非键作用的截断距离是多少? 根据这一距离, 两个周期性映像之间允许的最近距离多大?
当最小距离小于截断距离时, 哪一非键能量项受影响最大, 为什么?
如果最小距离小于截断距离, 会发生什么情况? 你的模拟中是否出现了这种情况?
选择C-alpha组重新运行mindist, 结果有变化么? 对你的体系而言, 这意味着什么?

要注意小距离事件是不时出现还是持续出现. 如果持续出现, 很可能影响蛋白的动力学; 如果只是偶尔出现, 那基本没有影响.

不仅直接相互作用需要担心, 非直接效应, 即以水为媒介的间接相互作用也可能引起问题. 例如, 蛋白可以导致水在离其最近的四个溶剂化层中出现有序性, 这大约对应于1 nm的距离. 理想情况下, 最小距离不应该小于2 nm.

3.3 分析回旋半径

回旋半径取决于某(些)原子质量与分子重心的关系, 可用于表征蛋白质结构的密实度. 作为QA的一部分, 我们将计算回旋半径, 它给出了每一时刻分子形状的信息, 可以与实测得的水动力学半径相比. 可以使用gyrate计算回旋半径, 这一程序将会给出回旋半径的各个分量, 对应于惯性矩阵的本征值. 这意味着第一个单独的分量对应于分子的最长轴, 最后一个对应于最短轴. 实际上, 三个轴的值给出了分子形状的整体描述,

gmx gyrate -f traj.xtc -s topol.tpr -o radius-of-gyration.xvg

查看回旋半径及其分量, 注意这些值如何达到平衡值.

回旋半径收敛了么? 如果是, 什么时候收敛, 收敛值为多少?
不同多肽模拟给出的回旋半径其涨落行为是否类似, 为什么?
你能将回旋半径的涨落与周期性映像的最小距离联系起来么?

3.4 均方根波动

除了能量等性质, 也能够通过结构的变化和弛豫来考察模拟趋向平衡的收敛性. 通常, 这种弛豫仅仅使用结构到参考结构(如晶体结构)的欧几里德距离来衡量. 这一距离被称为均方根偏差(RMSD). 然而, 我们也建议再考察一下到平均结构的弛豫, 即相对于平均结构的RMSD, 个中原因将在下节说明. 但是要计算相对于平均结构的RMSD, 需要首先获得平均结构. 平均结构可以在计算均方根波动(RMSF)时顺便获得.

RMSF计算每个原子相对于其平均位置的涨落, 表征了结构的变化对时间的平均, 给出了蛋白各个区域柔性的表征, 对应于晶体学中的b因子(温度因子). 通常, 我们预期RMSF和温度因子类似, 这可以用于考察模拟结果是否与晶体结构符合. RMSF(和平均结构)使用rmsf命令计算. -oq选项可以计算b因子, 并将其添加到参考结构中. 我们最关心的是每个残基的涨落, 这可使用选项-res设定.

gmx rmsf -f traj.xtc -s topol.tpr -o rmsf-per-residue.xvg -ox average.pdb -oq bfactors-residue.pdb -res

使用xmgrace查看RMSF的图形, 区分柔性和刚性区域.

确定最柔性区域的起始和终止残基编号, 给出其最大振幅.
比较不同多肽构象的结果. 是否有区别? 如果是, 哪个构象柔性最大, 哪个构象柔性最小?

为了对这些结果获得关联性的印象, 这里有个人类朊蛋白质1qlz的RMSF, 图中标示出了会导致CJD疾病的突变残基.

将两个pdb文件载入PyMOL, 根据b因子对结构bfactors.pdb进行着色, 并确定柔性区域. 平均结构是非物理结构. 查看下其中的一些侧链, 注意观察平均对构象的影响.

pymol average.pdb bfactors-residue.pdb
spectrum b, selection=bfactors-residue

下图是根据计算的b因子对1KLU多肽的参考结构进行着色后的结果(提示: 下图的设置如下: as sticks; show spheres; set stick_radius, 0.12; set sphere_scale, 0.25)

颜色分布在b因子值的范围内, 其中蓝色表示最小值(最稳定), 红色表示最高值(波动最大). 你可以通过截断最大值来调整颜色范围, 例如将其设置为350:

Q = 350; cmd.alter(“all”, “q = b > Q and Q or b”); spectrum q, selection=bfactors-residue

如果你很好奇, 可以再次计算b因子, 但这次使用每个原子的值:

g_rmsf -f traj.xtc -s topol.tpr -o rmsf-per-atom.xvg -oq bfactors-atom.pdb

在当前的PyMOL中载入新生成的bfactors-atom.pdb文件, 这样可以直接将其与残基平均的b因子相比较. 如果需要, 不要忘了重新调整颜色范围.

load bfactors-atom.pdb
spectrum b, selection=bfactors-atom

比较和对照两个b因子的结构.

以下图像显示的是根据模拟计算得到人类野生型UbcH8蛋白的b因子着色图. 蓝色对应低值, 红色对应高值.白色区域表示目前已知的能够反转蛋白质相互作用特征的残基. 在图像右侧, 你可以看到那些在螺旋2的前后环区有较高的b因子.第二个图像是螺旋2的前环的放大显示.

3.5 RMSD的收敛

当心! 由于你的多肽可能会跳出盒子外, 我们必须处理轨迹, 将粒子重新置于中心的周期性映像中. 为此, 可使用下面的命令:

gmx trjconv -f traj.xtc -o traj_nojump.xtc -pbc nojump -dt 50

由于计算RMSF时也得到了平均结构, 我们现在可以计算均方根偏差(RMSD)). RMSD通常用于表征结构到平衡态的收敛情况. 如前面所讲, RMSD是结构变化对原子总数的平均, 基本上是一个距离表征, 低的值最有意义. RMSD可以用rms命令计算, 此命令可以实现对不同时刻不同分组的原子进行结构平均. 首先计算所有多肽原子的RMSD, 使用初始结构作为参考结构

gmx rms -f traj_nojump.xtc -s topol.tpr -o rmsd-all-atom-vs-start.xvg

如果观察到了, 在什么时间RMSD达到平台期, 平衡值多少?

再次计算, 但只考虑骨架原子:

gmx rms -f traj_nojump.xtc -s topol.tpr -o rmsd-backbone-vs-start.xvg

这次的RMSD值更低, 这是由于计算时排除了通常更柔性的侧链原子, 在两种情形下. 两个RMSD都应该增大到一个平台值, 这意味着相对于参考结构, 多肽的结构达到了一定的距离, 然后或多或少保持在那个距离. 然而, 随着距离的增加, 可能的构象数目也在增加, 这意味着尽管RMSD达到了一个平台值, 但结构可能仍在趋近于其平衡态. 为此, 我们建议同时也检查趋向于平均结构的收敛情况.

echo 1 | gmx trjconv -f traj_nojump.xtc -s topol.tpr -o traj_protein_nojump.xtc
gmx rms -f traj_protein_nojump.xtc -s average.pdb -o rmsd-all-atom-vs-average.xvg
gmx rms -f traj_protein_nojump.xtc -s average.pdb -o rmsd-backbone-vs-average.xvg

与前面得到的图像进行比较. 注意在哪一点RMSD值开始趋平.

简要讨论下两个(相对于起始结构与相对于平均结构)图像的区别, 哪一个更能表征收敛情况?

到从为止, 我们完成了分析的第一部分, 包括可视化考察和质量确保检查. 现在该深入一点, 发掘一下蛋白质内部的情况. 分析的第二部分包括可根据多肽构象计算的结构性质, 例如氢键数量, 溶剂可及表面或特定的原子-原子间距离等.

待续......