Gromacs教程II-MD结果分析.docx

1、Gromacs 中文教程（II）：结果分析淮海一粟MD结果分析模拟结束后，就可以分析数据了。分析包括三个阶段。首先，有必要对模拟的质量进行检查，如果检查结果表明模拟良好，那么就可利用该模拟对所研究的问题作出回答；最后，不同的模拟结果可以合并。注：文件名应该反映文件内容，这根据你的模拟系统不同而不同。这里我们假定使用默认文件名，那么就会产生以下文件： topol.tpr: 模拟开始时包含完整系统描述的输入文件； confout.gro: 结构稳健，包含最后一步的坐标和速度； traj.trr: 全部精确轨迹，包括位置、速度和随时间变化的力 traj.xtc: 压缩的轨迹文件，只包含低精度(0.0

2、01 nm)的坐标信息； ener.edr: 随时间变化的有关能量数据 md.log: 模拟过程的日志附注：许多分析工具都能生成.xvg文件。这些文件能用xmgr或xmgrace查看，也可用python脚本程序xvg2ascii.py在终端显示出来。Each group writes one report. For general questions a single answer should be given in the report. Questions specific to each simulation should be given in table, indicated wi

3、th (T). Answers to questions from one section should be combined in a single table if possible. 先检查一下结果在开始分析前，要确定模拟是否正确地完成。有很多原因会导致模拟中断，尤其是与力场和系统平衡不充分引起的问题。为了检查模拟是否正确完成，运行gmxcheck程序：gmxcheck -f traj.xtc看看是否执行了10纳秒的模拟。=Q= How many frames are in the trajectory file and what is the time resolution? (T)

4、 另一个重要的信息源是日志文件。在文件md.log的末尾有模拟过程的统计数据；包括内存和CPU的使用情况和模拟时间。看看日志文件的末尾，使用less命令时，你可以用G (shift-g) 命令，跳到文件末尾。=Q= How long did the simulation run in real time (hours), what was the simulation speed (ns/day) and how many years would the simulation take to reach a second? (T) =Q= Which contribution to the p

5、otential energy accounts for most of the calculations? Dont be afraid to use the Gromacs online manual, to search in the Gromacs mailing list or even to use Google to get information about the terminology used by Gromacs. 结果可视化好玩的环节开始了。虽然很多分析都能归结为从轨迹文件中提取图像，MD当然首先要关注系统的移动。来看看轨迹文件。 首先用gromacs提供的查看器ng

6、mx来看看。虽然该软件的完善程度和视觉效果不及其他查看器，但它能够在拓扑文件信息的基础上画出键。其他查看器可能隐含远程键，这可能导致这些键被认为太长而不画出，或者会在非常接近的原子之间画出键。这是对模拟结果分析的一个常见错误源。使用ngmx载入拓扑和轨迹文件：ngmx -s topol.tpr -f traj.xtc看看程序菜单，试试不同的选项。播放动画。观看过程通过右边的选项控制。右击或左击选择选项来改变查看。=Q= What happens if the protein diffuses over the boundary of the box? 为了视觉美观，我们将从轨迹文件中提取100

7、0帧 (-dt 10)并忽略水分子（当软件询问时，选择 Protein）。而且，我们还将忽略边界的跨越，作出连续轨迹（-pbc nojump）。为了做这些工作，我们使用瑞士军刀般的gromacs工具trjconv，该工具有1001个选项组合。我们用它写出一个多模型的pdb文件，从而能在Pymol中观看。trjconv -s topol.tpr -f traj.xtc -o protein.pdb -pbc nojump -dt 10在Pymol中提取轨迹文件：pymol protein.pdb当所有的帧装入完毕，播放动画。Mplay动画播放时，其他控制键仍在运行，可以用鼠标旋转、放大或缩小图像

8、，也可以改变分子外观。SpectrumShow cell如果没错的话，你现在能看到蛋白质扩散、翻转跳跃。但我们对内部运动更感兴趣，而不是总体行为。在Pymol中，你能使用命令intra_fit将其他所有帧与第一帧对齐。随后，你可以用定向工具设定蛋白质中心：intra_fit protein and (name c,n,ca)Orient现在，所有的帧应该都对齐了，你可以看到蛋白质的哪一部分移动得更厉害。这些差异将在以后定量分析。当然，在cartoon模式下，蛋白质看起来更舒服，试试这个命令：show cartoon因为.pdb文件里面没有二级结构信息，你可能会看到碳骨架是个粗管状，而看不到正确

9、的二级结构。Pymol可以自动计算蛋白质的二级结构，但只计算一帧，并将其映射到其它的帧。例如，下述命令可以计算第一帧的二级结构：dss通过设定状态，用于计算的帧可以改变：dss state=1000最后，同时查看所有帧并检查蛋白质的柔性和刚性部分。set all_states=1请随便练习Pymol的使用。试试放大柔性或刚性区域，并检查侧链构象。使用ray和 png制作一份图像，即使浪费点CPU时间也不要紧。但记住，如果图像太复杂，可能会导致pymol的插件ray-tracer崩溃，这种情况下，你可以直接用png 得到屏幕上的图像。The following part, up to quali

10、ty assurance, is optional and it may be best to first finish the other sections. 如果你有足够机时做此教程，或者你已经做完了前面的功课，那就来做个不错的电影吧。你可能注意到，这些轨迹噪音非常多，那是热噪音，是正常的；但是对于制作好的电影会有影响。我们可以屏蔽这些高频运动，只保留低速和平滑的总体移动。为了达到这个目的，使用g_filter程序：g_filter -s topol.tpr -f traj.xtc -ol filtered.pdb -fit -nf 5现在，在Pymol中倒入轨迹文件。计算二级结构 (ds

11、s)并显示(show cartoon)，隐藏碳骨架上的侧链 (hide lines, not (name c,n,o) 并给蛋白质上色，然后orient分子设好视角。现在开始制作电影了： viewport 640,480set ray_trace_frames,1mpng frame_.png现在退出Pymol (quit) 并显示文件路径 (ls)。你能看到，文件的数量多了好些，包括1000个图像。以每秒30帧的速度，将会产生30秒的电影。下载mpeg_encode程序和参数文件movie.param，用它来产生单帧图像的电影（你可能需要编辑参数文件来改变文件名）：mpeg_encode m

12、ovie.param质量确认进行了最初的轨迹图像查看后，该对模拟的质量进行彻底检查了。这个质量检查(QA)包括热动力参数（如温度、压力、势能和动能）的收敛情况。更一般地说， QA试图评价（模拟）是否达到平衡状态。 结构上的收敛也需要检查，这个用起始结构和平均结构的均方差 (RMSD)来表示。随后，还要检查邻近的周期性图像之间没有互相作用，因为这将导致非物理学效应。最后，QA检查包括原子的均方差，这个可以与晶体学数据b-factors进行比较。 能量收敛我们首先从能量文件中提取一些热动力学数据。研究以下性质：温度、压力、势能、动能、单位盒子体积、密度和盒子尺寸。这些大部分性质已在系统准备步骤中检

13、查过了。用工具g_energy进行能量分析。该程序读出能量文件，也就是模拟过程中产生的扩展名为.edr的文件。g_energy程序将会问需要提取什么参数并将产生一幅图像。输入下面的命令：g_energy -f ener.edr -o temperature.xvg此命令将产生一列能量及其参量，这些参数存贮在.edr能量文件中。本教程的能量文件可能含有68个参量，每个都可以提取并画出图像。 最开始九个对应于力场中的不同能量。还应注意，从第47个开始同时列出了蛋白和非蛋白的参量，及两者之间的相互作用。为了提取温度，输入14 0 (Gromacs version 4.0.5)，回车。用xmgrace

14、程序看图，看看在指定温度附近(300 K)的温度如何波动。从波动也可以计算体系热容，热容在g_energy输出文件的结尾。 xmgrace -nxy temperature.xvg=Q= What is the average temperature and what is the heat capacity of the system? (T) 通过调用能量参量名，可以自动运行能量文件。使用echo和管道命令 (|)，实现从一个程序到另一个程序的数据传输，g_energy可以自动回应。为了提取多个参量，每个参量以n划分。拷贝并粘贴，或者输入以下命令行提取其他参量。不幸的是，能量参量必须用数字

15、指定。echo 15 0 | g_energy -f ener.edr -o pressure.xvg echo -e 11n12n13 0 | g_energy -f ener.edr -o energy.xvg echo 20 0 | g_energy -f ener.edr -o volume.xvg echo 21 0 | g_energy -f ener.edr -o density.xvg echo -e 17n18n19 0 | g_energy -f ener.edr -o box.xvg 逐个查看这些文件，看看数值的收敛情况。如果数值没有收敛，就表明模拟尚未达到平衡，需要延时

16、才能进行进一步分析。而且，平衡附近的数据不能用于分析。这里，为了简便起见，我们不管他了，直接用这些数据分析。=Q= What are the terms plotted in the files energy.xvg and box.xvg =Q= Estimate the plateau values for the pressure, the volume and the density. (T) 一些参量比其他的收敛慢。特别地，温度很容易收敛而系统各部分的相互作用收敛可能较慢。看看蛋白质和溶剂之间的相互关系：echo -e 51n53 0 | g_energy -f ener.edr -

17、o coulomb-inter.xvg echo -e 52n54 0 | g_energy -f ener.edr -o vanderwaals-inter.xvg 周期性图像建的最小距离对于QA，一个最重要的需要检查的事项是，周期性图像之间不应该有相互作用；因为周期性图像是有独立身份的，这些相互作用在物理学上不应该发生。设想双极性蛋白质图像存在直接的相互作用，那么同一蛋白质在盒子边界的两个末端就会产生作用力，这将影响蛋白质的本身行为并使模拟结果失效。 为了确认这些相互作用没有发生，我们用g_mindist命令计算周期性图像每次的最小距离： g_mindist -f traj.xtc -s

18、topol.tpr -od minimal-periodic-distance.xvg -pi =Q= What was the minimal distance between periodic images and at what time did that occur? (T) =Q= What happens if the minimal distance becomes shorter than the cut-off distance used for electrostatic interactions? Is it the case in your simulations? 小

19、距离的事件发生是偶然性的还是持续性的，也会有影响。如果是持续性的，就可能影响蛋白质动力学；但是如果只是偶尔发生，就一般不会有影响。=Q= Run now g_mindist on the C-alpha group, does it change the results? What does is mean for your system? 不只是直接相互作用需要关注，那些由水分子介导的间接效应，也会产生问题。例如，蛋白质可能影响最初四层水分子的排列，这相当于1 nm的距离；理想的情况是，最小距离不应该小于2nm。 波动的均方根除了能量本身的检查外，还应该检查模拟过程中，松弛后的蛋白质向平衡状

20、态收敛的情况。通常，松弛只是考虑了与参考结构（例如，晶体结构）的Euclidean距离。这个距离以 均方根偏差 (RMSD)表示。然而，我们推荐再检查一下与平均结构的松弛情况，也就是说，与平均结构的RMSD。个中原因，将在下一段说明。但为了计算与平结构的RMSD，需要首先得到平均结构。这个结构可在计算波动均方根(RMSF)时附加得到。RMSF的捕获每个原子相对平均位置的波动。这能深入观察蛋白质柔性部位的性质，相应于晶体结构测定中的b-factors (温度因子)。通常，我们希望RMSF与b-factors相近，这能表明模拟结果与晶体结构相适应。RMSF (和平均结构)用g_rmsf计算。-oq

21、选项运行计算b-factors，并把它们添加到参考结构中。我们最关心每个残基的波动，可用选项-res设定。g_rmsf -f traj.xtc -s topol.tpr -o rmsf-per-residue.xvg -ox average.pdb -oq bfactors.pdb -res 用xmgrace看看RMSF，找出柔性和刚性区域。=Q= Indicate the start and end residue for the most flexible regions and the maximum amplitudes. (T) =Q= Compare the results fro

22、m the different proteins. Are there differences? If yes, which is the most flexible and which least? 为了获得这些结果关联性的印象，这里有个人类prion蛋白质1qlz 的RMSF，图中会导致CJD疾病的突变残基被标示出来。现在来看看两个pdb文件，把它们读入Pymol。然后根据b-factors给结构文件bfactors.pdb上色，检查柔性区域。平均结构是非物理结构。看看其中的一些侧链，并注意平均结构对其构想的影响。pymol average.pdb bfactors.pdb spectru

23、m b, selection=bfactors 颜色分布在b-factor范围内，蓝色代表最低（最稳定），红色代表最高（最易波动）。你可以减少最大值来调整颜色范围，如设为500: Q = 500; cmd.alter(all, q = b Q and Q or b); spectrum q 以下图像显示的是人类野生型UbcH8蛋白根据模拟计算得到的b-factors上的颜色。蓝色是低，红色是高。白色区域表示目前已知的能够反转蛋白质相互作用特征的残基。在图像右侧，你可以看到那些在Helix 2前后的loops有较高的b-factors。第二个图像是helix 2前loop的放大显示。 RMSD的

24、收敛当心！因为你的蛋白质可能会跳出盒子外， 我们必须重新制作轨迹来把中间周期性图像中的粒子拽回来。为此，运行以下命令： trjconv -f traj.xtc -o traj_nojump.xtc -pbc nojump 由于RMSF计算也给出平均结构，我们现在计算均方根偏差（RMSD）。RMSD通常被用作指示到平衡状态的收敛情况。前面说过，RMSD仅仅是个距离单位，值越小越好。RMSD用g_rms程序计算。首先计算所有原子的RMSD，使用起始结构作为参考：g_rms -f traj_nojump.xtc -s topol.tpr -o rmsd-all-atom-vs-start.xvg =

25、Q= If observed, at what time and value does the RMSD reach a plateau? (T) 再次计算RMSD，但只计算骨架原子：g_rms -f traj_nojump.xtc -s topol.tpr -o rmsd-backbone-vs-start.xvg 这次，RMSD达到了一个低值，这是由于排除了侧链原子（侧链原子往往更容易发生运动）。两个RMSD最后都应到达一个平台值。这意味着蛋白质结构达到了与参考结构有一定距离并或多或少保持那个距离。然而，随着距离的增加，可能构象的数目也在增加。这意味着虽然RMSD达到了一个平台值，结构仍然

26、可能在向平衡状态收敛。出于这个原因，建议同时检查一下向平均结构的收敛情况。echo 1 | trjconv -f traj_nojump.xtc -s topol.tpr -o protein.xtc g_rms -f protein.xtc -s average.pdb -o rmsd-all-atom-vs-average.xvg g_rms -f protein.xtc -s average.pdb -o rmsd-backbone-vs-average.xvg 比较与上次图像的差别。注意RMSD值的平衡点。=Q= Briefly discuss two differences betw

27、een the graphs against the starting structure and against the average structure. Which is a better measure for convergence? 回旋半径的收敛QA的最后一步，我们将计算回旋半径。这个值表示每次分子形状的变化。将回旋半径与试验得出的回旋半径相比较。可以用g_gyrate计算回旋半径。该程序也会给出个性因子（individual components），相应于惯性矩阵的本征值。这意味着，第一个individual component对应于原子的最长轴，最后一个individual

28、 component相应于最小轴。这三个轴能说明分子的形状。输入以下命令： g_gyrate -f traj.xtc -s topol.tpr -o radius-of-gyration.xvg 看看回旋半径和individual components，注意这些值如何达到平衡值。=Q= At what time and value does the radius of gyration converge? (T) 这里，我们完成了分析的第一部分，包括图像检查和质量检查。现在该深入一点了，挖掘一下蛋白质内部发生的情况。第二部分的分析包括结构特性，这些特性可用蛋白质形状，例如氢键数量、溶剂可接近表

29、面或特定原子-原子之间的距离等。Lets go。结构分析：由形状得出的特性如果提示需要选择，选择 Protein；如果没有特别提示或者没有 if no selection is specifically stated or does not follow logically from the text. To get rid off the noise, please use the Running Average method in Data-Transformation to smooth your graphs with xmgrace 如果确认了模拟已经收敛，就可以进行真正的分析了。对

30、模拟数据的分析，可以分为几种类型。第一种包括对单个图形的解释，通过一些函数逐个点进行计算得到一个值或者一些值；例如RMSD和回旋半径的计算。这些特性可称为构型，依赖或瞬间特性。另外一种是可以通过一定时间范围内的平均化来分析的特征，比如相互关系或波动。本部分将进行一些通常的分析，每个都能得到直接来自于轨迹（不同时间下的坐标）的一个时间序列值。这些问题可以参考程序运行时的屏幕输出或者图形。溶剂可及表面面积一个可能感兴趣的特性是蛋白质表面可被溶剂到达的面积，一般指溶液科技表面 (SAS)或者溶剂可及表面面积 (SASA)。还可细分为亲水性SAS和疏水性SAS。除此之外，SAS可以和一些经验参数一起使

31、用，得到一个溶剂化自由能的估计。 所有这四个参数都可用g_sas程序完成。本程序也允许计算每个残基或原子一段时间内的平均SAS。输入下面的命令，设置需要计算SAS的基团和输出基团，查看输出文件。 g_sas -f traj_nojump.xtc -s topol.tpr -o solvent-accessible-surface.xvg -oa atomic-sas.xvg -or residue-sas.xvg=Q= Focusing on the loop1 (residues 53-62), which residues are the most accessible to the so

32、lvent? 氢键另一可能能提供很多信息的性质是氢键，内部氢键（蛋白质-蛋白质）或蛋白质与其周围的溶剂之间的氢键都是这样。氢键的存在与否可以通过氢键受体和供体对的距离和键角来推断。为了计算氢键，使用如下命令（并用Xmgrace查看输出文件）：echo 1 1 | g_hbond -f traj_nojump.xtc -s topol.tpr -num hydrogen-bonds-intra-protein.xvgecho 1 12 | g_hbond -f traj_nojump.xtc -s topol.tpr -num hydrogen-bonds-protein-water.xvg=Q= Discuss the r