发布时间 : 星期日 文章Gromacs教程II-MD结果分析 - 图文更新完毕开始阅读469eaee7dd3383c4ba4cd243
g_rmsdist -s topol.tpr -f traj.xtc -rms rmsdist.xpm -mean rmsmean.xpm -dt 10
==Q== Briefly explain the images: rmsmean in terms of structure and rmsdist in terms of flexibility/stability. Recall the information from earlier analysis and viewing the structure
从实验角度讲，这些距离同样重要。原子间距离的边界可以由NMR实验推断——利用核的Overhauser效应 (NOE)——这也是NMR结构计算间的主要推动力。如果蛋白质模型正确，这个应该可以预期得到的NOE信号，可以通过MD模拟来计算。这些信号与距离密切相关，特别是r-3 和 r-6 权重的距离。这些信号也可以用g_rmsdist来计算。
g_rmsdist -s topol.tpr -f traj.xtc -nmr3 nmr3.xpm -nmr6 nmr6.xpm -noe noe.dat -dt 10
给nmr3.cpm 和 nmr6.xpm重新上色，看看这些矩阵 ，也看看文件noe.dat。 ==Q== What are the smallest 1/r^3 and 1/r^6 averaged distances in the simulation? ( T )
松弛和序参量
向量的松弛的计算与其自身相关。对蛋白质，通常包括碳骨架N-H 或侧链C-H向量。这种自相关给出了一个向量能保持其方向多长时间的量度，因而为表征可变性和稳定性提供了指示。序参量是自相关的长程限制。如果一个分子能够自由旋转，序参量将不可避免地减小到零；但在分子框架内（内部参考框架），序参量常有一个明确的值，这个值表明总体稳定性。通过固定（fitting mistake for fixing?）蛋白质，这个参考框架可以固定下来，因此序参量就确定了。 N-H序参量可以用g_chi计算。这个程序可以写出一个.pdb文件，该文件把序参量加入了b-factor场，使观看更容易。该程序也计算J-couplings，这个参数可以和NMR 的结果相比较，或者用于指导NMR实验。
g_chi -s topol.tpr -f traj.xtc -o order-parameters.xvg -p order-parameters.pdb -jc Jcoupling.xvg
看看.xvg文件中的序参量，并用Pymol 查看.pdb文件，根据在b-factor场中的值，给残基上色。
==Q== Write down the start and end residues, and the average value for the two regions having highest order parameters. ( T )
==Q== How do the order parameters compare to the fluctuations (RMSF)?
构象主成分分析
See Leach Chapter 9.14 for more information on the following section.
一个常用的，但常常对之理解不深的方法是轨迹的主成分分析法。这种方法——有时候也指“本征动力学（ED）”——目的在于鉴别大量原子的整体运动并因此揭示隐藏在原子波动后面的结构信息。在MD中，粒子波动确实互相关联，因为粒子间的彼此作用。关联的程度有大有小，但需注意的是，那些通过键直接相连或者彼此位置接近的粒子会产生谐动。粒子之间的运动相关性产生系统总体波动中的结构，对于大分子来说，这种结构常常与其功能和生物物理学特性直接有关。因此，原子波动结构的研究可以对了解这些大分子的行为提供有价值的信息。然而，它确实需要相当程度的线性代数和多变量统计方面的知识，来解释结果并鉴别该方法的缺陷。特别是，PCA的目标是用新的变量来描述原始数据，这些新变量与原始变量线性关联。这也是PCA最主要的问题：它仅仅依据原子运动间的线性关系来解释问题。
PCA的第一步是构建协方差矩阵，它将捕捉每对原子间的原子运动相关度。协方差矩阵从定义上说，是个对称矩阵。这个矩阵随后将用于分析，产生一个特征向量矩阵和特征向量的对角线矩阵。每个特征向量描述了粒子的集体运动，这里向量的值表示相应原子参与了多少运动。给出的相关的特征值等于波动的总和，这些波动以集体运动中每个原子来描述，因此是与特征向量有关的全部运动的一个度量值。通常，系统中的大多数 (>90%) 运动是用少于10个特征向量或主成分来描述的。
协方差分析产生相当多的文件，因此最好在一个新的子文件夹中运行： mkdir COVAR cd COVAR
协方差矩阵的构建和对角线化可用g_covar程序。输入下述命令进行分析。注意，这可能要花费一些时间。
g_covar -s ../topol.tpr -f ../traj.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm
==Q== What are the dimensions of the covariance matrix and what is the sum of the eigenvalues? ( T )
特征值之和是系统总运动性的一个度量。它可用于比较不同情况下蛋白质的可变性，但是对不同大小的蛋白质进行比较时，就难以得到有意义的解释。
现在，看看协方差矩阵自身，它在covar.xpm文件里。 xview covar.xpm
矩阵显示了原子间的协方差。红色表示两个原子一起移动，而蓝色表示它们彼此向相反的方向运动。对角线化相应于早前得到的RMSF图，红色的深度表示波动的振幅大小。
==Q== Look at the two most moving parts, excluding the termini. How do they move with respect to each other and to the rest of the protein? 从协方差矩阵能看到一组组相关或反相关运动的原子。这允许往原子组的集体运动中重新写入总运动。我们提到过，特征值储存在eigenvalues.xvg文件中，通过相应特征值表示出总波动。
==Q== Look at the file eigenvalues.xvg with a text editor. Calculate the percentage and cumulative percentage of the motion explained for the first five eigenvectors. ( T )
典型地，最初五个特征值将捕获主要运动，这相当于 >80% 的总运动。如果解释的总波动较低，就提示没有明确的集体运动。
为了了解特征值的实际意义，可用g_anaeig工具作进一步的分析。为了更近地看看前两个特征值，输入以下命令：
g_anaeig -s ../topol.tpr -f ../traj.xtc -v eigenvectors.trr -eig eigenvalues.xvg -proj proj-ev1.xvg -extr ev1.pdb -rmsf rmsf-ev1.xvg -first 1 -last 1
g_anaeig -s ../topol.tpr -f ../traj.xtc -v eigenvectors.trr -eig eigenvalues.xvg -proj proj-ev2.xvg -extr ev2.pdb -rmsf rmsf-ev2.xvg -first 2 -last 2
特征值对应于运动方向。选项-extr沿着选定的特征值从轨迹中提取极端结构。把这些结构导入PyMOL看看： pymol ev?.pdb
把pdf文件中的模型分开，删除原始结构。 split_states ev1 split_states ev2 delete ev1 or ev2
给模型上色。特征值1中的极端结构显示为蓝-绿色而特征值二为黄-红色。 spectrum count dss
hide everything show cartoon
用PyMOLs的'align'命令，能画出表示两种构象差异的小条。 align ev1_0001 and c,n,ca),object=diff1 align ev2_0001 and c,n,ca),object=diff2
(n.
c,n,ca),ev1_0002
and
(n. (n. c,n,ca),ev2_0002 and (n. ==Q== What is the largest difference between the extreme structures for eigenvector 1? And for eigenvector 2?
为了理解特征值的意义，想象一下欧洲城市间卖货人的移动所发生的事情。这种移动可用地球坐标体系解释，这种情况下我们需要对每个位置采用三个坐标。虽然这是对的，但如果你只想解释卖货人的移动，就不是最佳的。因为理论上，任何坐标体系都一样好，我们可以定义一个新的（坐标体系）来解释卖货人的移动。实际上，因为地球表面可以用二维空间代表，我们只需要两个坐标而无须三个。 直观看来，有人会拿南北极和东西轴说事，但也可以从移动中推断出轴。比方说卖货人在伦敦-雅典轴上走的最多。这个轴作为第一个特征值；第二个特征值与第一个正交。用这种方式，卖货人在欧洲的每个位置就可以用在这两个特征值上的投影唯一地确定下来。对它们的投影作图，就可以显示出旅行路线。沿着第一个坐标的极端投影对应于Reykjavik和莫斯科，即使它们实际上不在这个轴上。
蛋白质构象也和上面所述的一模一样。你看到的极端投影并不一定对应于物理结构，但是它们允许看到沿轴的运动和总的运动程度。为了解蛋白质沿构象空间的移动，我们可以画出特征值2对特征值1的投影图。为此，从两个.xvg文件中提取投影数据并且合并到文件ev1-vs-ev2.dat中。注意'>'表示输出由屏幕重定向到一个文件中，所以你看不到任何屏幕输出。
grep -v \proj-ev1.xvg | awk '{print $2}' > proj-ev11.dat
grep -v \proj-ev2.xvg | awk '{print $2}' > proj-ev12.dat
paste proj-ev11.dat proj-ev12.dat > ev1-vs-ev2.dat 为了解释模拟过程，我们也提取了最后7.5 ns（最后1500点）和5.0 ns（最后1000点）的投影。然后把它们导入xmgrace。
tail -1500 ev1-vs-ev2.dat > last-7.5ns.dat tail -1000 ev1-vs-ev2.dat > last-5.0ns.dat
xmgrace ev1-vs-ev2.dat last-7.5ns.dat last-5.0ns.dat ==Q== What is the shape of the projections? Are these mutually independent (oval distribution)?
==Q== Would the same eigenvectors (axes) be obtained if the analysis were performed on the last 7.5 ns? And on the last 5.0 ns?
结语
Write a concluding paragraph, comparing the results obtained for the different proteins. Also reflect on the overall stability and the probability that the structure deposited in the PDB properly reflects the solution structure. In other words, does the structure stay close to the starting structure or does it drift away and how much? Can you suggest a mechanism that would explain the differences in UbcH6 and UbcH8 interaction profiles and why a single conserved mutation can restore the rich interaction profile of UbcH6 when starting from UbcH8? What are the limits of the simulations to explain such a complex biological process?