NAMD入门教程(二)

3. 分析方法

在实际工作中，分子动力学模拟体系的设计和模拟计算的进行可能并不困难。关键之处在于如何分析动力学模拟得出的结果，说明一定的问题。因此我们专门设置一章讲解动力学模拟的分析方法。

本章中我们将讲解四个例子：残基RMSD值，平均能量，mb分布和温度分布。每一个例子都很典型，代表了一定的分析思想。但是由于某些例子需要较长的计算时间，读者没有必要自己进行动力学模拟。在NAMD教程文件中已经给出了动力学模拟得到的结果。

3.1 每个氨基酸残基的RMSD值

上一章中，我们已经计算了整个蛋白质的RMSD值变化。我们初步提到了，RMSD 值反映的是偏离平均位置的程度，是原子运动幅度的大小的衡量。这里我们将计算每个氨基酸残基的RMSD值，按照RMSD值对蛋白着色，并在最后进行问题讨论。本节的目的是使读者明白进行RMSD值分析的基本方法和意义所在。

3.1.1 RMSD值的计算与蛋白着色显示

我们将要使用一个VMD提供的脚本计算每个残基的平均RMSD值，并把这个值储存在VMD提供的用户存储区中（User Field），然后对蛋白进行着色。

1、首先，打开VMD，选择File→New Molecule菜单项，在弹出的文件浏览中找到common目录下的文件usq_wb.psf，点击按钮Load。这样我们就载入了蛋白质的结构信息。（注意图形窗口中是空白的，并没有显示出蛋白结构，因为psf文件不存储各个原子的坐标）

2、此时Molecule File Browser窗口上端一栏应该显示Load files for: 0: ubq_wb.psf （图）说明如果再次载入新文件，所载入的文件将与usq_wb.psf结合起来一同处理显示。我们将载入原子运动轨迹文件，和psf结构文件配合即可显示出蛋白质原子的运动轨迹。因此点击Browse，找到1-2-sphere文件夹下的ubq_wb_eq.dcd，点击Load，这时关闭Molecule File Browser窗口。

图载入ubq_wb.psf后的Molecule File Browser

3、打开VMD TKConsole，使用cd命令改变当前目录到namd-tutorial/2-1-rmsd/

4、在这一目录下有我们准备好的脚本residue_rmsd.tcl（可能读者也注意到了，VMD 中脚本文件的后缀名均为“.tcl”）。在TKConsole中输入：

source residue_rmsd.tcl

这一命令并不进行RMSD值的求算，而是载入了该脚本以供用户调用其中的命令。求算RMSD值所用的命令是：rmsd_residue_over_time。其格式为：

rmsd_residue_over_time 目标分子所选残基

5、这里我们将选择蛋白质的所有氨基酸残基进行计算。选择方法是输入：

set sel_resid [[atomselect top “protein and alpha”] get resid]

这一命令的作用是：定义一个变量sel_resid，储存氨基酸残基的序号列表，方便下一步调用。在TKConsole中可以看到列出的氨基酸残基序号列表（图），由1至76依次编号。

6、对蛋白的所有残基（1号至76号）进行RMSD值计算。输入命令：

rmsd_residue_over_time top $sel_resid

输入后VMD调用脚本计算每个残基的RMSD值。计算公式和2.6节中介绍的公式相同。此时应当能够在窗口中看到程序处理过程（图）：

待计算完成，输出结果应该如下图：

与此同时，程序会将每个原子的RMSD值存储在该原子的VMD用户存储区中。我们将使用VMD对每个氨基酸残基按照RMSD值着色，以清楚地显示各个残基在动力学模拟时的运动剧烈程度。

7、选择Graphics→Representations菜单项。

8、在Atom Selection window一项中输入protein后回车，然后在Drawing Method 下拉列表中选择Tube，在Coloring Method下拉列表中选择User。然后单击Trajectory 选项卡，在Color Scale Data Range中输入0.40 和1.00，单击Set。

以上设置方法如图：

现在在图形窗口中，应该可以看到按照残基RMSD 值着色的蛋白骨架（图）。

在这幅彩图中，蓝色的氨基酸残基是运动幅度最大的，红色残基是幅度最小的，这和我们通常的着色方式恰好相反。我们要进行一下更改：

9、选择Graphics→Colors菜单项。选择Color Scale 选项卡，在Method下拉列表中选择BWR。现在红色是活动剧烈的残基，蓝色是运动幅度较小的残基，白色介于二者之间（图）。转动分子仔细观察，氨基酸残基运动的幅度和二级结构有什么关系？（如果不方便观察二级结构，可以选择Gr aphics→Representation，将Drawing Method改变为Cartoon。

3.1.2 氨基酸RMSD值分布图

我们下面将使用Excel处理数据，作出RMSD值分布图

10、关闭VMD，打开Excel，选择“文件”→“打开”，找到2-1-rmsd目录，选择文件residue_rems.dat打开，注意“文件格式”下拉列表要选择“所有文件”，否则将看不到我们所要的文件。

11、在弹出的“文件导入向导”中，首先单击“下一步”，然后单击“完成”。可以看到A列是氨基酸残基的序号（1至76），B列是每个残基对应的RMSD值。

12、选择“插入”→“图表”，在“图表类型”中选择“XY 散点图”，在“子图表类型”中选择左下角的“折线散点图”，点击“完成”。

得到的图表（图）就是泛素的氨基酸残基RMSD值变化分布。

当然，必须记住我们不是为了作图而作图。我们作图的目的是分析RMSD 值的变化，发现问题和规律。仔细观察你得到的RMSD 值图，可以看到一些有意思的问题。

以下是我们提出的两个问题，供读者思考：

1、我们将蛋白质的二级结构大体分成三类：α螺旋，β片层和无规则卷曲（loop ）。使用前面讲过的DeepView 可以得到泛素二级结构的氨基酸残基分布：

1-7 片层；11-17 片层；23-33 螺旋；40-45 片层；48-50 片层；57-59 螺旋；65-72 片层，其余部分认为是无规则卷曲（loop ）。

在Excel 中计算一下哪一种二级结构的平均RMSD 值较高？哪一种二级结构的平均RMSD 值较低？注：二级结构的平均RMSD 值 = 该种二级结构的总RMSD 值 / 该二级结构氨基酸残基数。思考一下这种规律背后的原因是什么？

参考答案： β片层：0.537； α螺旋：0.562； loop ：0.744

2、从我们作出的图中可以看到，C 末端最后4个氨基酸的RMSD 值明显较高。 在VMD 中观察C 末端loop

的构象，有什么特色？

前面提到，泛素的功能是通过标记特定蛋白质，使得被标记的蛋白特异性降解（参见。你认为泛素最有可能是通过哪一端和被降解蛋白相连？C 末端还是N 末端？如果已知泛素和标记蛋白的半胱氨酸相连，推测一下泛素末端的哪一个氨基酸负责连接？

参考答案：C 末端；Arg

思考过上面两个问题后，可能读者已经初步理解了RMSD 值分析的意义。我们从X 射线衍射得到的蛋白结构仅仅是一个静态的“死”模型，这一结构并不能直接告诉我们各个氨基酸残基的运动特性。但是进行动力学模拟时，我们手中的蛋白质模型就由“死”变“活”，运动起来了。而通过RMSD 值分析，我们可以初步了解各个氨基酸残基的空间位阻大小，运动的方式，以及不同的二级结构的运动特点等等信息。这些信息往往同蛋白质的功能相联系（比如上面第二个问题）。总而言之，RMSD 值分析可以为我们阐明结构和功能的关系提供思路和突破口。

3.2 能量分析

上一节中我们所分析的RMSD值反映的是体系中不同原子的运动情况。在这一节中，我们将对体系能量进行分析。下面我们将计算整个体系中不同形式的能量（如动能，键能，静电势能等）在动力学模拟时的平均值。

我们将计算上一章完成的立方水体分子动力学模拟过程中的能量平均值。计算时我们又一次需要使用一个脚本文件：namdstats.tcl。这个脚本文件在进行能量平均值的计算时是十分有用的。

1、打开VMD，打开TKConsole，使用cd命令改换当前目录至2-3-energies/

2、在TKConsole中输入：

source namdstats.tcl

data_avg ../1-3-box/ubq_wb_eq.log 101 last

第一行命令将载入脚本文件namdstats.tcl, 第二行命令将提取日志文件中记录能量的那一部分并进行计算。101 last表示计算动力学模拟第101步至最后一步过程中的平均能量值。

3、结果会直接输出在TKConsole中（图）。

图计算能量平均值的输出结果

上图中所使用的缩写：BOND键能，ANGLE键角能DIHED二面角能ELECT静电能VDW范德华力能KINETIC动能TEMP温度。在这里仅需要读者明白这些平均能量计算的方法，熟悉namdstats.tcl脚本的使用即可。

本节到此就算完成了。但是关于脚本文件namdstats.tcl我们还需要多了解一些，因为这一脚本在进行动力学模拟结果分析时是非常常用的。

脚本namdstats.tcl定义了两个命令：data_avg和data_time。

?data_avg 的功能是计算特定时间段内整个体系各种能量的平均值。调用方法是：data_avg <日志文件> <起始> <终止>

其中，<日志文件>表示用于计算的”.log”文件的位置。<起始> 和<终止>分别代表计算。第一个timestep 可以用first代替，最后一个timestep可以用last代替。因此我们在例子中输入的101 last，表示的是从第101步计算到最后一个timestep。如果省略<起始>和<终止>，仅仅输入data_avg <日志文件>，VMD将默认对于动力学模拟的整个过程进行计算。

?data_time 的功能是输出特定时间段内整个体系某个特定能量随时间的连续变化。调用方法是：

data_time <能量> <日志文件> <起始> <终止>

NAMD将会从日志文件中提取<能量>这一参数设定的能量值从起始timestep到终止timestep的连续变化过程，并输出到<能量>.dat这一文件中。<起始>和<终止>的设定方法同上。

在动力学模拟的多种分析方法中，能量分析是最基础最经典的，通过能量分析可以直接或间接说明许多问题。但是进行分析时需要过硬的理论功底，因此本书仅作简单介绍，其他具体的分析实例请参照相关的文献专著。

上两节我们分别对体系中原子的运动和体系的能量进行了简单分析。我们的分析有助于阐明泛素的结构和功能的关系，以及泛素自身的某些性质，因此我们得出的结果具有明确的生物学意义。

在下面两节中，我们将抛开我们所研究体系的生物学意义，对系统的统计力学性质进行分析。读者可能会问：我们做的是生物学研究，为什么不分析体系的生物学性质，却要考察体系的纯物理学性质？这里有必要作一下说明。

我们知道，物理学是公理化的科学，经过几百年的发展，理论体系相对较为完备，纯粹理论推导出的结果能够很好地与实验符合。但生命科学依然是描述性的科学。得出的许多基本规律，比如中心法则，遗传密码，其本质是经验性的，没有理论推导可言。而在分子动力学模拟领域，虽然力学规律是严格的微分方程，基本的相互作用力场参数（如CHARMM，Amber等）却是半经验式的，依靠大量实验得出的，仅在统计意义上成立，而不能保证动力学模拟的结果一定正确。

除了力场参数引入的误差，我们在进行模拟时还必须人为地做出一些假设。读者可以再仔细阅读一下2.3.1一节，配置文件中我们的设定：我们假定范德华力和静电力只在某一范围内起作用，假定计算时某些原子的相互作用被忽略等等，都使得体系越来越不精确。而在真正的研究工作开展时，我们可能需要在体系中放入离子以模拟真实的细胞质基质，需要使用格点计算静电力场分布……这些人为引入的假设和近似能否符合生命体系中的真实过程？从基础理论的推导我们得不出答案。

那么，如何验证我们得出的结果的正确性？一个方法当然是通过实际实验检验，但实验可能需要大量时间，并且要受条件的限制。因此，另一个方法是：使用物理学基本规律进行证伪。前面说到了，物理学规律是建立在严密的理论基础之上的，其结论是经历过实验的考验的。因此我们对所得结果进行分析，如果符合物理学基本规律，并不能证明结果正确——因为还不能说明是否符合生命体系中的真实过程；但是如果不符合，那么一定能说明结果是错误的。因此，物理学基本规律可以作为一个只能证伪的有力判据。

具体到本例，我们将检验分子动力学模拟得出的结果是否符合统计物理中的两个结论：mb分布和温度分布。如果不符合，说明我们的结果不正确，原因可能是体系有问题，力场参数不适合，边界条件设置不当等等。

同时，完成下面两节后，读者应该掌握一些重要的数据处理思想和方法。

3.3 麦克斯韦-波尔兹曼（Maxwell-Boltzmann ）能量分布

下面我们将通过对分子动力学模拟所得结果的分析，验证动力学模拟的结果是否符合Maxwell-Boltzmann能量分布。本节内容较多，也较为综合，并将详细讲述Excel回归分析，拟合数据的基本方法。

3.3.1 计算每个原子的动能

计算动能分布自然首先要得到各个原子的动能，而动能的计算需要得到某一瞬间各个原子的速率。在2.6节我们提到，NAMD 动力学模拟的输出文件包括“.vel”文件，该文件记录原子的瞬时速率。我们这里使用的就是立方水体中泛素分子动力学模拟时原子的速率文件。

虽然我们在2.5节已经进行了立方水体的动力学模拟，但是在模拟时我们设定了Rigid Bonds （见2.3.1）,步长是2fs ，这样的近似设定使得体系不能精确反映各个原子的运动速率。因此我们取消Rigid Bonds ，步长1fs 重新进行动力学模拟，并将结果输出文件提供给读者。

1、首先，我们仍旧要载入泛素的结构信息。打开VMD ，选择File→New Mol ecule 菜单项，找到common 目录下的文件ubq_wb.psf ，载入该文件。不要关闭Molecule File Browser 窗口，窗口第一项应该显示Load file for: 0:ubq_wb.psf 。

2、下面载入泛素中各原子的分子速率信息。再次单击Browser 按钮，找到2-2-maxwell 目录下，载入ubq_wb_eq_1fs.restart.vel 文件。在Determine file type ：下拉列表中选择NAMD Binary Coordinates ，然后再点击Load 。现在可以关闭Molecule File Browser 。

图 Molecule File Browser 载入速率文件时的设置

图载入速率文件后图形窗口的显示

可以在图形窗口中看到一个海胆一样浑身是刺的怪物！这是因为VMD把各个原子的速度当作原子的坐标处理（瞬时速度是以x,y,z三个方向的分速度记录的，记录格式和座标一样）。不过这并不影响我们下一步的操作。我们将使用VMD得到只存储各个原子质量和瞬时速度的文件，便于我们计算动能。

3、在VMD TKConsole窗口中，用cd命令改变目录到namd-tutorial/2-2-maxwell。

4、然后输入：

set all [atomselect top all]

这样我们建立了一个变量all ，储存系统中的所有原子。（系统返回：atomselect0，表示这一选择的序号是0）

5、下面在当前目录下新建一个文件energy.dat : set fil [open energy.dat w]

6、下面计算每个原子的动能：2

12

k mv ε=

,并把计算结果写入文件energy.dat 。 foreach m [$all get mass] v [$all get {x y z}] {

puts $fil [expr 0.5 * $m * [vecdot $v $v] ] }

(注意：连续输入既可，不需要每一行回车) 7、关闭文件： close $fil 然后退出VMD

这样我们就获得了文件energy.dat ，存储所有的原子及其相应的动能。

3.3.2 绘制实验数据柱状图

下面我们将首先使用Excel 处理原始数据，做出柱状图，作为我们下一步拟合的准备工作。

1、打开Excel ，选择“文件”→“打开”，找到2-2-maxwell 目录下的energy.dat 文件并打开它，注意“文件类型”设置为“所有文件”

2、此后会弹出文本导入向导，直接点击“完成”，这样在A 列显示的就是各个原子的动能值。

下面我们将把这些动能值处理成柱状图，并和Maxwell-Boltzmann 分布拟合，最后可以得出系统的平均温度。

1、我们首先要做的是：统计所有原子的动能值的分布频数。打个比喻，我们现在得到了一个班全体同学身高的列表，身高的分布频数就是统计有多少人的身高在[150cm, 155cm], [155cm,165cm],[165cm, 170cm]...区间之内。显然，为了将动能分布等分成区间，我们首先要知道这些动能值中的最大值和最小值。

在B1中输入 =max （A1：A7051）回车。max 是Excel 提供的求指定范围内最大值的函数。我们指定的范围是A1到A7051，因为向下拖动滚动条可知一共有7051个原子的动能值。同样地，在B2中输入 =min （A1：A7051）回车，可得到这一范围内的最小值。结果如图

图利用max和min函数得到最大最小值

因为最大值是6.5左右，最小值在0附近，我们将0到7以0.1为间距等分成70个区间，可以包含所有的动能值。

将B1和B2单元格中的公式删除，然后在B1中输入0，B2中输入=B1+0.1 回车，然后再次选中B2，单击B2右下角的小黑点后，不要松开鼠标左键，一直向下拖动到B71后松开左键。可以看到单元格已经被自动填充上了我们想要的区间（图）

图得到分布区间

2、下面我们将使用专门的分析工具统计动能值在这70个区间中的分布频数

单击“工具”→“加载宏”，在弹出的对话框中选择“分析工具库”，点击“确定”。再次打开“工具”菜单，应该出现“数据分析”这一项（图）。

图“加载宏”对话框

图“工具”中新增的“数据分析”菜单项

图“直方图”对话框

单击打开“数据分析”。选择“直方图”，在弹出的“直方图”对话框中填入：

输入数据：$A$1:$A$7051

接受区域：$B$1:$B$71

单击“确定”后会自动新建一个新的工作表，A列为“接受”，B列为“频率”。我们可以保

存一下以防工作丢失。不要单击，选择“文件”→“另存为”，“保存类型”选择第一项“Microsoft

Office Excel 工作簿”，重新命名为rmsd，单击“保存”。

3、保存之后，我们开始对数据进行作图。我们希望获得的是标准化的频率分布图。

所谓标准化（normalization）就是使得频率分布柱状图的面积之和为一。为了达到这一点，

需要先对数据进行如下处理：

在C2单元格中输入公式=B2*0.1 回车，然后按照上面用过的方法单击并拖动C2单元格的右下角一直到第72行，将C2到C72填充该公式。

4、在D2单元格输入公式=sum（C2：C72）回车。sum是Excel提供的求和函数。这样在D2中求得了C2至C72之中所有数据的和。

5、在E2单元格中输入公式=B2/D$2回车。同样地，单击并拖动E2的右下角一直到第72行。这样数据处理过程就进行完毕了。完成后的工作表顶端应当如图：

图完成后的工作表顶端

6、下面开始作图。选择“插入”→“图表”，在弹出的对话框中选择图表类型：“柱形图”，子图表类型：左上角第一项“簇状柱形图”，单击“下一步”，选择“系列”选项卡（图）

图系列选项卡

7、然后不断单击按钮“删除”，删除掉所有的已有系列。再按“添加”新增一个系列，单击“值”这一项最右边的按钮（图）。这个按钮的作用是允许用户开始选定特定的区域。因此我们以后称这个按钮为“区域选择按钮”

图单击区域选择按钮

8、此时图表向导对话框消失，并出现一浮动条（图）。此时鼠标变成十字指针，用鼠标选定E2到E72单元格，单击区域按钮回到图表向导（图）。

图：图表源数据选择

图：再次单击按钮回到图表向导

图完成以上过程后的图表向导对话框

9、然后以同样的方法设置“分类X轴标志”的区域：单击最右侧的区域选择按钮，用鼠标选定A2到A72单元格，再次单击按钮回到图表向导对话框。完成后对话框应该如图：

10、下面可以点击按钮“完成”，得到最初的柱状图（图）。

3.3.3 绘制mb分布理论曲线

上一小节，我们使用Excel绘制出了实验数据的柱状图。下面我们将根据mb分布方程绘制理论曲线，以便直观的看到拟合前后发生的变化。

1、在F2单元格输入以下公式：

= (2 / SQRT(pi() * G$2^3)) * SQRT(A2) * EXP(-A2/G$2)

回车后，单元格内显示错误，这是因为我们还需要在G2单元格中输入公式的一项参数。输入0.5（注意：这个值是随意的，仅仅是作为一个初值。读者可以输入0.3，0.7等等），作为初始参数。然后选中F2单元格，单击按住单元格右下角的黑点向下拖动，一直到F72单元格松开鼠标，将这些单元格内填充上同样的公式。这样我们就在F列中得到了理论计算求得的Maxwell-Boltzmann分布。

在拟合之前，我们先直观地看一下理论求得的Maxwell分布和我们的实验值（E列）的偏差有多大。

3、在刚刚生成的柱状图上先单击左键，再单击右键，这时弹出的右键菜单第一项应当是“数据系列格式”（图）。如果不是，说明刚才单击左键的位置不对。注意一定要左键单

击有柱形图存在的区域，不要单击空白区域。

4、右键菜单中选择“数据系列格式”，找到最后一个选项

卡“选项”，更改“分类间距”为0（图）。确定后，可以看到柱形

图中一个个条形之间已经没有间距（图）。

图：更改数据系列格

式所弹出的右键菜单