Gromacs教程IIMD结果分析.docx

Gromacs教程IIMD结果分析.docx

《Gromacs教程IIMD结果分析.docx》由会员分享，可在线阅读，更多相关《Gromacs教程IIMD结果分析.docx（23页珍藏版）》请在冰点文库上搜索。

Gromacs教程IIMD结果分析

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

动画播放时，其他控制键仍在运行，可以用鼠标旋转、放大或缩小图像，也可以改变分子外观。

Spectrum

Showcell

如果没错的话，你现在能看到蛋白质扩散、翻转跳跃。

但我们对内部运动更感兴趣，而不是总体行为。

在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,480

setray_trace_frames,1

mpngframe_.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.xvg

echo-e"11\n12\n130"|g_energy-fener.edr-oenergy.xvg

echo200|g_energy-fener.edr-ovolume.xvg

echo210|g_energy-fener.edr-odensity.xvg

echo-e"17\n18\n190"|g_energy-fener.edr-obox.xvg

逐个查看这些文件，看看数值的收敛情况。

如果数值没有收敛，就表明模拟尚未达到平衡，需要延时才能进行进一步分析。

而且，平衡附近的数据不能用于分析。

这里，为了简便起见，我们不管他了，直接用这些数据分析。

==Q==Whatarethetermsplottedinthefilesenergy.xvgandbox.xvg

==Q==Estimatetheplateauvaluesforthepressure,thevolumeandthedensity.（ T ）

一些参量比其他的收敛慢。

特别地，温度很容易收敛而系统各部分的相互作用收敛可能较慢。

看看蛋白质和溶剂之间的相互关系：

echo-e"51\n530"|g_energy-fener.edr-ocoulomb-inter.xvg

echo-e"52\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?

Isitthecaseinyoursimulations?

小距离的事件发生是偶然性的还是持续性的，也会有影响。

如果是持续性的，就可能影响蛋白质动力学；但是如果只是偶尔发生，就一般不会有影响。

==Q==Runnowg_mindistontheC-alphagroup,doesitchangetheresults?

Whatdoesismeanforyoursystem?

不只是直接相互作用需要关注，那些由水分子介导的间接效应，也会产生问题。

例如，蛋白质可能影响最初四层水分子的排列，这相当于1nm的距离；理想的情况是，最小距离不应该小于2nm。

波动的均方根

除了能量本身的检查外，还应该检查模拟过程中，松弛后的蛋白质向平衡状态收敛的情况。

通常，松弛只是考虑了与参考结构（例如，晶体结构）的Euclidean距离。

这个距离以均方根偏差（RMSD）表示。

然而，我们推荐再检查一下与平均结构的松弛情况，也就是说，与平均结构的RMSD。

个中原因，将在下一段说明。

但为了计算与平结构的RMSD，需要首先得到平均结构。

这个结构可在计算波动均方根（RMSF）时附加得到。

RMSF的捕获每个原子相对平均位置的波动。

这能深入观察蛋白质柔性部位的性质，相应于晶体结构测定中的b-factors（温度因子）。

通常，我们希望RMSF与b-factors相近，这能表明模拟结果与晶体结构相适应。

RMSF（和平均结构）用g_rmsf计算。

-oq选项运行计算b-factors，并把它们添加到参考结构中。

我们最关心每个残基的波动，可用选项-res设定。

g_rmsf-ftraj.xtc-stopol.tpr-ormsf-per-residue.xvg-oxaverage.pdb-oqbfactors.pdb-res

用xmgrace看看RMSF，找出柔性和刚性区域。

==Q==Indicatethestartandendresidueforthemostflexibleregionsandthemaximumamplitudes.（ T ）

==Q==Comparetheresultsfromthedifferentproteins.Aretheredifferences?

Ifyes,whichisthemostflexibleandwhichleast?

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

现在来看看两个pdb文件，把它们读入Pymol。

然后根据b-factors给结构文件bfactors.pdb上色，检查柔性区域。

平均结构是非物理结构。

看看其中的一些侧链，并注意平均结构对其构想的影响。

pymolaverage.pdbbfactors.pdb

spectrumb,selection=bfactors

颜色分布在b-factor范围内，蓝色代表最低（最稳定），红色代表最高（最易波动）。

你可以减少最大值来调整颜色范围，如设为500:

Q=500;cmd.alter（"all","q=b>QandQorb"）;spectrumq

以下图像显示的是人类野生型UbcH8蛋白根据模拟计算得到的b-factors上的颜色。

蓝色是低，红色是高。

白色区域表示目前已知的能够反转蛋白质相互作用特征的残基。

在图像右侧，你可以看到那些在Helix2前后的loops有较高的b-factors。

第二个图像是helix2前loop的放大显示。

RMSD的收敛

当心！

因为你的蛋白质可能会跳出盒子外，我们必须重新制作轨迹来把中间周期性图像中的粒子拽回来。

为此，运行以下命令：

trjconv-ftraj.xtc-otraj_nojump.xtc-pbcnojump

由于RMSF计算也给出平均结构，我们现在计算均方根偏差（RMSD）。

RMSD通常被用作指示到平衡状态的收敛情况。

前面说过，RMSD仅仅是个距离单位，值越小越好。

RMSD用g_rms程序计算。

首先计算所有原子的RMSD，使用起始结构作为参考：

g_rms-ftraj_nojump.xtc-stopol.tpr-ormsd-all-atom-vs-start.xvg

==Q==Ifobserved,atwhattimeandvaluedoestheRMSDreachaplateau?

（ T ）

再次计算RMSD，但只计算骨架原子：

g_rms-ftraj_nojump.xtc-stopol.tpr-ormsd-backbone-vs-start.xvg

这次，RMSD达到了一个低值，这是由于排除了侧链原子（侧链原子往往更容易发生运动）。

两个RMSD最后都应到达一个平台值。

这意味着蛋白质结构达到了与参考结构有一定距离并或多或少保持那个距离。

然而，随着距离的增加，可能构象的数目也在增加。

这意味着虽然RMSD达到了一个平台值，结构仍然可能在向平衡状态收敛。

出于这个原因，建议同时检查一下向平均结构的收敛情况。

echo1|trjconv-ftraj_nojump.xtc-stopol.tpr-oprotein.xtc

g_rms-fprotein.xtc-saverage.pdb-ormsd-all-atom-vs-average.xvg

g_rms-fprotein.xtc-saverage.pdb-ormsd-backbone-vs-average.xvg

比较与上次图像的差别。

注意RMSD值的平衡点。

==Q==Brieflydiscusstwodifferencesbetweenthegraphsagainstthestartingstructureandagainsttheaveragestructure.Whichisabettermeasureforconvergence?

回旋半径的收敛

QA的最后一步，我们将计算回旋半径。

这个值表示每次分子形状的变化。

将回旋半径与试验得出的回旋半径相比较。

可以用g_gyrate计算回旋半径。

该程序也会给出个性因子（individualcomponents），相应于惯性矩阵的本征值。

这意味着，第一个individualcomponent对应于原子的最长轴，最后一个individualcomponent相应于最小轴。

这三个轴能说明分子的形状。

输入以下命令：

g_gyrate-ftraj.xtc-stopol.tpr-oradius-of-gyration.xvg

看看回旋半径和individualcomponents，注意这些值如何达到平衡值。

==Q==Atwhattimeandvaluedoestheradiusofgyrationconverge?

（ T ）

这里，我们完成了分析的第一部分，包括图像检查和质量检查。

现在该深入一点了，挖掘一下蛋白质内部发生的情况。

第二部分的分析包括结构特性，这些特性可用蛋白质形状，例如氢键数量、溶剂可接近表面或特定原子-原子之间的距离等。

Let’sgo。

结构分析：

由形状得出的特性

如果提示需要选择，选择"Protein"；如果没有特别提示或者没有ifnoselectionisspecificallystatedordoesnotfollowlogicallyfromthetext.

Togetridoffthenoise,pleaseusethe'RunningAverage'methodin'Data->Transformation'tosmoothyourgraphswithxmgrace

如果确认了模拟已经收敛，就可以进行真正的分析了。

对模拟数据的分析，可以分为几种类型。

第一种包括对单个图形的解释，通过一些函数逐个点进行计算得到一个值或者一些值；例如RMSD和回旋半径的计算。

这些特性可称为构型，依赖或瞬间特性。

另外一种是可以通过一定时间范围内的平均化来分析的特征，比如相互关系或波动。

本部分将进行一些通常的分析，每个都能得到直接来自于轨迹（不同时间下的坐标）的一个时间序列值。

这些问题可以参考程序运行时的屏幕输出或者图形。

溶剂可及表面面积

一个可能感兴趣的特性是蛋白质表面可被溶剂到达的面积，一般指溶液科技表面（SAS）或者溶剂可及表面面积（SASA）。

还可细分为亲水性SAS和疏水性SAS。

除此之外，SAS可以和一些经验参数一起使用，得到一个溶剂化自由能的估计。

所有这四个参数都可用g_sas程序完成。

本程序也允许计算每个残基或原子一段时间内的平均SAS。

输入下面的命令，设置需要计算SAS的基团和输出基团，查看输出文件。

g_sas-ftraj_nojump.xtc-stopol.tpr-osolvent-accessible-surface.xvg-oaatomic-sas.xvg-orresidue-sas.xvg

==Q==Focusingontheloop1（residues53-62）,whichresiduesarethemostaccessibletothesolvent?

氢键

另一可能能提供很多信息的性质是氢键，内部氢键（蛋白质-蛋白质）或蛋白质与其周围的溶剂之间的氢键都是这样。

氢键的存在与否可以通过氢键受体和供体对的距离和键角来推断。

为了计算氢键，使用如下命令（并用Xmgrace查看输出文件）：

echo11|g_hbond-ftraj_nojump.xtc-stopol.tpr-numhydrogen-bonds-intra-protein.xvg

echo112|g_hbond-ftraj_nojump.xtc-stopol.tpr-numhydrogen-bonds-protein-water.xvg

==Q==Disc