搜索
Gromacs中文教程（II）：结果分析淮海一粟MD结果分析模拟结束后，就可以分析数据了。分析包括三个阶段。首先，有必要对模拟的质量进行检查，如果检查结果表明模拟良好，那么就可利用该模拟对所研究的问题作出回答；最后，不同的模拟结果可以合并。注：文件名应该反映文件内容，这根据你的模拟系统不同而不同。这里我们假定使用默认文件名，那么就会产生以下文件：topol.tpr:模拟开始时包含完整系统描述的输入文件；confout.gro:结构稳健，包含最后一步的坐标和速度；traj.trr:全部精确轨迹，包括位置、速度和随时间变化的力traj.xtc:压缩的轨迹文件，只包含低精度(0.001nm)的坐标信息；ener.edr:随时间变化的有关能量数据md.log:模拟过程的日志附注：许多分析工具都能生成.xvg文件。这些文件能用xmgr或xmgrace查看，也可用python脚本程序xvg2ascii.py在终端显示出来。Eachgroupwritesonereport.Forgeneralquestionsasingleanswershouldbegiveninthereport.Questionsspecifictoeachsimulationshouldbegivenintable,indicatedwith(T).Answerstoquestionsfromonesectionshouldbecombinedinasingletableifpossible.先检查一下结果在开始分析前，要确定模拟是否正确地完成。有很多原因会导致模拟中断，尤其是与力场和系统平衡不充分引起的问题。为了检查模拟是否正确完成，运行gmxcheck程序：gmxcheck-ftraj.xtc看看是否执行了10纳秒的模拟。==Q==Howmanyframesareinthetrajectoryfileandwhatisthetimeresolution?(T)另一个重要的信息源是日志文件。在文件md.log的末尾有模拟过程的统计数据；包括内存和CPU的使用情况和模拟时间。看看日志文件的末尾，使用'less'命令时，你可以用'G'(shift-g)命令，跳到文件末尾。==Q==Howlongdidthesimulationruninrealtime(hours),whatwasthesimulationspeed(ns/day)andhowmanyyearswouldthesimulationtaketoreachasecond?(T)==Q==Whichcontributiontothepotentialenergyaccountsformostofthecalculations?Don'tbeafraidtousetheGromacsonlinemanual,tosearchintheGromacsmailinglistoreventouseGoogletogetinformationabouttheterminologyusedbyGromacs.结果可视化好玩的环节开始了。虽然很多分析都能归结为从轨迹文件中提取图像，MD当然首先要关注系统的移动。来看看轨迹文件。首先用gromacs提供的查看器ngmx来看看。虽然该软件的完善程度和视觉效果不及其他查看器，但它能够在拓扑文件信息的基础上画出键。其他查看器可能隐含远程键，这可能导致这些键被认为太长而不画出，或者会在非常接近的原子之间画出键。这是对模拟结果分析的一个常见错误源。使用ngmx载入拓扑和轨迹文件：ngmx-stopol.tpr-ftraj.xtc看看程序菜单，试试不同的选项。播放动画。观看过程通过右边的选项控制。右击或左击选择选项来改变查看。==Q==Whathappensiftheproteindiffusesovertheboundaryofthebox?为了视觉美观，我们将从轨迹文件中提取1000帧(-dt10)并忽略水分子（当软件询问时，选择Protein）。而且，我们还将忽略边界的跨越，作出连续轨迹（-pbcnojump）。为了做这些工作，我们使用瑞士军刀般的gromacs工具trjconv，该工具有1001个选项组合。我们用它写出一个多模型的pdb文件，从而能在Pymol中观看。trjconv-stopol.tpr-ftraj.xtc-oprotein.pdb-pbcnojump-dt10在Pymol中提取轨迹文件：pymolprotein.pdb当所有的帧装入完毕，播放动画。Mplay动画播放时，其他控制键仍在运行，可以用鼠标旋转、放大或缩小图像，也可以改变分子外观。SpectrumShowcell如果没错的话，你现在能看到蛋白质扩散、翻转跳跃。但我们对内部运动更感兴趣，而不是总体行为。在Pymol中，你能使用命令intra_fit将其他所有帧与第一帧对齐。随后，你可以用定向工具设定蛋白质中心：intra_fitproteinand(namec,n,ca)Orient现在，所有的帧应该都对齐了，你可以看到蛋白质的哪一部分移动得更厉害。这些差异将在以后定量分析。当然，在cartoon模式下，蛋白质看起来更舒服，试试这个命令：showcartoon因为.pdb文件里面没有二级结构信息，你可能会看到碳骨架是个粗管状，而看不到正确的二级结构。Pymol可以自动计算蛋白质的二级结构，但只计算一帧，并将其映射到其它的帧。例如，下述命令可以计算第一帧的二级结构：dss通过设定状态，用于计算的帧可以改变：dssstate=1000最后，同时查看所有帧并检查蛋白质的柔性和刚性部分。setall_states=1请随便练习Pymol的使用。试试放大柔性或刚性区域，并检查侧链构象。使用'ray'和'png'制作一份图像，即使浪费点CPU时间也不要紧。但记住，如果图像太复杂，可能会导致pymol的插件ray-tracer崩溃，这种情况下，你可以直接用'png'得到屏幕上的图像。Thefollowingpart,upto'qualityassurance',isoptionalanditmaybebesttofirstfinishtheothersections.如果你有足够机时做此教程，或者你已经做完了前面的功课，那就来做个不错的电影吧。你可能注意到，这些轨迹噪音非常多，那是热噪音，是正常的；但是对于制作好的电影会有影响。我们可以屏蔽这些高频运动，只保留低速和平滑的总体移动。为了达到这个目的，使用g_filter程序：g_filter-stopol.tpr-ftraj.xtc-olfiltered.pdb-fit-nf5现在，在Pymol中倒入轨迹文件。计算二级结构(dss)并显示(showcartoon)，隐藏碳骨架上的侧链(hidelines,not(namec,n,o))并给蛋白质上色，然后orient分子设好视角。现在开始制作电影了：viewport640,480setray_trace_frames,1mpngframe_.png现在退出Pymol(quit)并显示文件路径(ls)。你能看到，文件的数量多了好些，包括1000个图像。以每秒30帧的速度，将会产生30秒的电影。下载mpeg_encode程序和参数文件movie.param，用它来产生单帧图像的电影（你可能需要编辑参数文件来改变文件名）：mpeg_encodemovie.param质量确认进行了最初的轨迹图像查看后，该对模拟的质量进行彻底检查了。这个质量检查(QA)包括热动力参数（如温度、压力、势能和动能）的收敛情况。更一般地说，QA试图评价（模拟）是否达到平衡状态。结构上的收敛也需要检查，这个用起始结构和平均结构的均方差(RMSD)来表示。随后，还要检查邻近的周期性图像之间没有互相作用，因为这将导致非物理学效应。最后，QA检查包括原子的均方差，这个可以与晶体学数据b-factors进行比较。能量收敛我们首先从能量文件中提取一些热动力学数据。研究以下性质：温度、压力、势能、动能、单位盒子体积、密度和盒子尺寸。这些大部分性质已在系统准备步骤中检查过了。用工具g_energy进行能量分析。该程序读出能量文件，也就是模拟过程中产生的扩展名为.edr的文件。g_energy程序将会问需要提取什么参数并将产生一幅图像。输入下面的命令：g_energy-fener.edr-otemperature.xvg此命令将产生一列能量及其参量，这些参数存贮在.edr能量文件中。本教程的能量文件可能含有68个参量，每个都可以提取并画出图像。最开始九个对应于力场中的不同能量。还应注意，从第47个开始同时列出了蛋白和非蛋白的参量，及两者之间的相互作用。为了提取温度，输入140(Gromacsversion4.0.5)，回车。用xmgrace程序看图，看看在指定温度附近(300K)的温度如何波动。从波动也可以计算体系热容，热容在g_energy输出文件的结尾。xmgrace-nxytemperature.xvg==Q==Whatistheaveragetemperatureandwhatistheheatcapacityofthesystem?(T)通过调用能量参量名，可以自动运行能量文件。使用'echo'和管道命令(|)，实现从一个程序到另一个程序的数据传输，g_energy可以自动回应。为了提取多个参量，每个参量以\n划分。拷贝并粘贴，或者输入以下命令行提取其他参量。不幸的是，能量参量必须用数字指定。echo150|g_energy-fener.edr-opressure.xvgecho-e11\n12\n130|g_energy-fener.edr-oenergy.xvgecho200|g_energy-fener.edr-ovolume.xvgecho210|g_energy-fener.edr-odensity.xvgecho-e17\n18\n190|g_energy-fener.edr-obox.xvg逐个查看这些文件，看看数值的收敛情况。如果数值没有收敛，就表明模拟尚未达到平衡，需要延时才能进行进一步分析。而且，平衡附近的数据不能用于分析。这里，为了简便起见，我们不管他了，直接用这些数据分析。==Q==Whatarethetermsplottedinthefilesenergy.xvgandbox.xvg==Q==Estimatetheplateauvaluesforthepressure,thevolumeandthedensity.(T)一些参量比其他的收敛慢。特别地，温度很容易收敛而系统各部分的相互作用收敛可能较慢。看看蛋白质和溶剂之间的相互关系：echo-e51\n530|g_energy-fener.edr-ocoulomb-inter.xvgecho-e52\n540|g_energy-fener.edr-ovanderwaals-inter.xvg周期性图像建的最小距离对于QA，一个最重要的需要检查的事项是，周期性图像之间不应该有相互作用；因为周期性图像是有独立身份的，这些相互作用在物理学上不应该发生。设想双极性蛋白质图像存在直接的相互作用，那么同一蛋白质在盒子边界的两个末端就会产生作用力，这将影响蛋白质的本身行为并使模拟结果失效。为了确认这些相互作用没有发生，我们用g_mindist命令计算周期性图像每次的最小距离：g_mindist-ftraj.xtc-stopol.tpr-odminimal-periodic-distance.xvg-pi==Q==Whatwastheminimaldistancebetweenperiodicimagesandatwhattimedidthatoccur?(T)==Q==Whathappensiftheminimaldistancebecomesshorterthanthecut-offdistanceusedforelectrostaticinteractions?I