福建水产设备联盟

GMX模拟多肽-蛋白相互作用(四)

分子模拟之道 2020-08-12 06:32:44

第四部分 分子动力学模拟数据的分析(二)

结构分析: 构象派生的性质

说明:

  • 当提示选择时, 如果教程没有明确说明如何选择, 或者没有遵循教程自身的逻辑, 请选择蛋白Protein组.

  • 如果需要, xmgrace可以使用滑动平均方法对数据进行平滑(去除噪声): Data -> Transformation.

  • 确认模拟已经收敛到平衡态后, 就可以进行一些真正的分析了, 模拟数据的分析可以分为几种类型. 第一类包括对单个构象进行解释, 在每个时间点根据一些函数获得一个值或多个值. RMSD和回旋半径就是例子. 这样的性质, 可称为构象依赖或瞬时性质. 此外, 也可以在时域对过程进行分析, 例如, 通过对一段时间内的平均化得到(自)相关或涨落. 在本部分会进行一些常见的分析, 其中的每种分析都会得到直接由轨迹(随时间变化的坐标)导出的某个值的时间序列. 问题可以参考程序运行时的屏幕输出或图像.

    1. 氢键

    氢键数目是一类信息丰富的性质, 无论是内部氢键(蛋白-蛋白)还是蛋白和周围溶剂之间的氢键. 氢键存在与否可以通过氢键施体-H-受体之间的距离和施体-H-受体之间角度来推断. hbond命令可计算模拟过程中分子间或组间的氢键数目以及氢键距离或角度的分布. 使用下面的命令, 然后查看得到的输出文件

    echo 1 1 | gmx hbond -f traj_nojump.xtc -s topol.tpr -num hydrogen-bonds-intra-protein.xvg
    echo 1 12 | gmx hbond -f traj_nojump.xtc -s topol.tpr -num hydrogen-bonds-protein-water.xvg

    讨论两种情形下氢键数目的关系, 每种情形下的氢键数目的波动情况.

    特定的氢键可以使用包含待研究原子编号的索引文件来考察. 查看多肽前一半和后一半所涉及的氢键. 你必须看看confout.gro文件来检查残基编号并将多肽大致分为两半. 假定你的多肽含有14个残基, 其编号始于22终于35. 你想将它分为两部分, 22-28和29-35. 下面的第一个命令将在组选择菜单中创建两个新的条目, part_1part_2. 第二个命令是一个使用索引文件运行hbond的通用命令. 你可能想要看看N端和C端之间的氢键, 例如, 可以用来监测β发卡的形成.

    echo “r 22-28\nname 19 part_1\nr 29-35\nname 20 part_2\nq” | gmx make_ndx -f confout.gro -o my_index.ndx
    gmx hbond -f traj_nojump.xtc -s topol.tpr -num hydrogen-bonds-strand.xvg -n my_index.ndx

    注意, 你必须根据自己的多肽替换上面命令中的数字.

    根据氢键分析, 你的多肽在模拟过程中是否形成类似β发卡的构象?

    2. 二级结构

    判定蛋白结构的最常用参数是指定二级结构元素, 如α螺旋, β片层. 能提供这一信息的一个程序是dssp, 但前提是你的电脑已经安装了dssp程序. 此程序并不是GROMACS发行的一部分, 但能够从CMBI, Radboud University获得. 下载后请解压, 然后将环境变量DSSP设定为程序的路径, 如/home/student/dssp.

    GROMACS提供了一个dssp的接口, 可以计算轨迹中每帧的二级结构.

    首先, 需要生成消除跳跃的轨迹

    gmx trjconv -f traj.xtc -o traj_nojump.xtc -pbc nojump -dt 50

    然后运行以下命令

    gmx do_dssp -f traj_nojump.xtc -s topol.tpr -o secondary-structure.xpm -sc secondary-structure.xvg

    secondary-structure.xvg文件包含一个时间序列, 列出了每帧中与每一二级结构类型相关的残基数目. 更多的细节在.xpm文件中, 它使用颜色编码了每个残基在一段时间内的二级结构. .xpm文件可以使用类似Gimp的程序查看, 但可以使用GROMACS工具xpm2ps添加一些有用的元数据, 得到的结果可以使用gview查看, 或转换为pdf后使用xpdf或evine查看.

    gmx xpm2ps -f secondary-structure.xpm -o secondary-structure.eps -size 4000
    ps2pdf secondary-structure.eps
    evince secondary-structure.pdf

    讨论二级结构的变化, 如果有的话.
    比较不同蛋白二级结构的稳定性.

    3. 拉氏图(Ramachandran (phi/psi) plots)

    足球直播免费高清无插件直播骨架的phi和psi扭转角是两个洞察足球直播免费高清无插件直播结构特性的有用参数. phi对psi的绘图称为拉氏图, 该图中某些特定区域反映了蛋白二级结构元素或氨基酸的特性, 而(这些区域之外的)其他区域被认为是禁阻的(不可到达). phi/psi随时间的变化可以体现出结构的转变. 这些角度可以通过rama程序进行计算, 尽管结果有些粗略, 因为程序将所有角度都输出到单一的文件中. 若想研究单个残基, 可通过Linux程序grep将其从图中选出来.

    gmx rama -f traj_nojump.xtc -s topol.tpr -o ramachandran.xvg

    此文件中包含了全部氨基酸的所有phi/psi角. 要抽取某个残基的角度, 例如LEU-24, 可键入

    egrep “@|LEU-24” ramachandran.xvg > rama-LEU-24.xvg

    抽取每个残基(除边界外)的拉氏图, 并用xmgrace对其进行可视化, 描述它们的主要相似性与差别.

    根据拉氏图, 讨论每个残基的构象稳定性.

    动力学和时间平均性质的分析

    4. 合并轨迹

    每组学生模拟了一个从不同构象开始的多肽, 这样做的目的在于增加采样, 或者增加模拟能够覆盖的结构空间. 为了对轨迹和多肽构象性质有一个完整的观点, 我们必须将轨迹拼合到一起.

    由于大多数操作/计算在执行时都需要一个拓扑文件, 为避免使用哪个文件, 从那个轨迹开始之类的问题, 请从这里下载你的多肽文件. 这一文件包含了多肽分子位于原始蛋白-多肽复合物中的直角坐标. 你也可以使用这一结构来计算相对于结合结构的结构相似性. 下载pdb文件后, 你必须将其转换为GROMACS的坐标文件(.gro), 这是教程中第一步骤的重复操作. 别忘了选择合适的力场和水模型, 尽管后一选择对分析并无影响.

    gmx pdb2gmx -f XXXX_peptide.pdb -o XXXX_peptide.gro -ignh

    这样, 使用新生成的这个文件来创建索引文件, 这样我们可以进行更多的选择. 在这一步, 与前面不同, 我们不需要指定任何特定的残基组, 因此, 我们简单地推出程序(键入q)

    echo “q” | gmx make_ndx -f XXXX_peptide.gro -o XXXX_peptide_index.ndx

    现在我们可以对我们新生成的轨迹进行截断, 使用trjconv去除前10 ns. -b选项设定要创建的新xtc文件的开始时间, 单位为ps, 这意味着我们需要新文件从10 ps开始. -dt选项指定我们要保留的时间精度. 当提示选择时, 选择Protein来生成xtc文件.

    gmx trjconv -f traj_helix.xtc -o traj_helix_10-50ns.xtc -n XXXX_peptide_index.ndx -pbc nojump -dt 100 -b 10000

    为什么我们只分析最后40 ns的轨迹?

    下面的命令会输出一个单一的.xtc文件, 其中包含以下列顺序排列的四个轨迹, bound-helix-polypro-extended. 这一顺序很重要因为后面要进行比较, 因此需要特别注意. -settime选项指定在拼合轨迹时使用连续的时间, 这意味着首帧时间为10 ns的第二个轨迹, 会恰好放在第一个轨迹的后面50 ns处. 如果你忘记了这个选项, 所拼合的轨迹具有重复的时间戳(即每一轨迹的原点都是10 ns). 我们要选择选项congtinuec, 指示trjcat从前一轨迹的最后时间戳处开始新的轨迹.

    gmx trjcat -f traj_bound_10-50ns.xtc traj_helix_10-50ns.xtc traj_polypro_10-50ns.xtc traj_extended_10-50ns.xtc -o traj_concat_10-50ns.xtc -cat -settime

    5. 又是RMSD

    我们在前面已经计算过RMSD, 并用其来检查模拟的收敛性, 但它也可以用于更进一步的分析. RMSD是两个结构之间的表征. 如果我们对轨迹文件中每一对结构的组合计算RMSD, 就可以看到是否有属于同一类型或具有相同特征的结构组. 属于同一组的结构其RMSD值较低, 而与其他组结构的RMSD值更高. 利用矩阵来表示RMSD值, 可以用于识别转变状态.

    要建立RMSD矩阵, 可使用rms处理两条轨迹. 如果你要单独考虑组(簇)及其在不同轨迹之间的转变, 可以将所有的四条轨迹合并为一条, 再使用rms生成交叉RMSD矩阵. 所有的RMSD计算时, 选择主链组Mainchain.

    gmx rms -s XXXX_peptide.gro -f traj_concat_10-50ns.xtc -f2 traj_concat_10-50ns.xtc -m rmsd-matrix.xpm -tu ns

    得到的矩阵是灰阶图. 为了显示更清楚, 可以使用彩虹梯度图.

    gmx xpm2ps -f rmsd-matrix.xpm -o rmsd-matrix.eps -rainbow blue
    ps2pdf rmsd-matrix.eps
    evince rmsd-matrix.pdf

    分析时为什么选择主链组?
    你看到多少簇?
    在不同的轨迹中, 你是否对相同构象进行了采样? 你是如何发现的? 什么是轨迹的重叠, 即对相同构象空间进行采样?
    你能发现多少转变?
    从这个分析你能得到什么结论? 这是你预期的结果么? 请证实你的观点.

    6. 聚类分析

    基于结构间的距离RMSD, 可以将结构归并为反映构象可及性范围及其相对权重的一组组团簇. 这可以通过聚类算法来完成, cluster命令实现了一些聚类算法. 这一程序会生成好几个输出文件. 检查该程序的帮助文档, 了解每一个它们的含义, 然后运行程序. 注意, 我们已经计算了RMSD矩阵, 可以将它作为cluster的输入.

    gmx cluster -h
    gmx cluster -s XXXX_peptide.gro -f traj_concat_10-50ns.xtc -dm rmsd-matrix.xpm -dist rmsd-distribution.xvg -o clusters.xpm -sz cluster-sizes.xvg -tr cluster-transitions.xpm -ntr cluster-transitions.xvg -clid cluster-id-over-time.xvg -cl clusters.pdb -cutoff 0.2 -method gromos

    聚类算法使用的RMSD截断值为多大? 这一值代表什么, 如何影响得到的团簇数目?
    共有多少团簇? 最大的两个其尺寸多大?

    更改RMSD的截断值, 以获得适当数目的团簇(10到100之间).
    使用PyMOL打开cluater.pdb, 比较前两个团簇的结构.

    disable all
    splitstates clusters
    delete clusters
    /for i in range(3,100): cmd.delete( “clusters
    %04d” %i )
    dss show cartoon
    util.cbam clusters_0002
    align clusters_0001 and ss h, clusters_0002 and ss h

    前两个团簇之间是否存在可觉察的差别?

    7. 距离RMSD

    采用RMSD比较结构的一个不足之处在于, 它包含了最小二乘叠合, 这会影响结果. 但是, 一个蛋白的结构也可以使用一系列的原子间距来表示. 这可以用于获得一个比较的表征量, 并且不依赖于叠合, 这就是距离RMSD(dRMSD). 可以利用rmsdist命令来计算dRMSD.

    gmx rmsdist -s XXXX_peptide.gro -f traj_concat_10-50ns.xtc -o distance-rmsd-all-atoms.xvg

    dRMSD何时收敛? 收敛值多大?
    所得图像与标准RMSD相比有何区别? 查阅GROMACS手册, 看二者是如何计算的. 试着根据二者的计算方法解释你观察到的区别.

    我们现在要回答构象采样的问题, 在我们的轨迹中, 我们是否对多肽的结合构象进行了采样? 为此, 我们要对多肽的结合构象再次运行rms. 我们可以使用合并的轨迹来进行这一分析, 而不是计算并生成四个不同的RMSD图形. 记住你生成的合并轨迹中构象的顺序, 这很重要. 当叠合计算RMSD值时, 选择主链Mainchain组.

    gmx rms -f traj_concat_10-50ns.xtc -s XXXX_peptide.gro -o rmsd_concat-mainchain-vs-bound.xvg
    xmgrace rmsd_concat-mainchain-vs-bound.xvg
    gmx rmsdist -f traj_bound_10-50ns.xtc -s XXXX_peptide.gro -o dist-rmsd_bound-mainchain-vs-bound.xvg
    xmgrace dist-rmsd_bound-mainchain-vs-*.xvg

    如果在计算时我们包含整个蛋白(即所有原子), 而不是只选择主链Mainchain原子, 结果会有什么变化?
    RMSD和距离RMSD图形有区别么?
    在不同轨迹中, 你是否对多肽的结合构象进行了采样? 对你研究的情况, 关于对蛋白多肽识别的构型采样假说你能得出什么结论?

    8. 自由能形貌图

    本分析的最后一步是计算自由能形貌(FEL, free energy landscape), 它由不同多肽构象开始的轨迹采样所得. 除了看起来很酷以外, 它能显示出 在整个模拟过程中多肽所经历的自由能的谷和丘. 此外, 在对接计算中需要选择有代表性的结构时, FEL会变得更有意义.

    FEL表示一个映射, 分子在模拟过程中所经历的所有可能构象到相应能量的映射. 典型的能量是Gibbs自由能. 正如你所想象的, 使用直角坐标来代表不同的构象是不合适的(你能想象4维形貌么?). 因此, FEL通常使用两个变量来表示, 它们反映了体系的特定性质, 并表征了构象变化. 例如你可以使用绕一根特定键的扭转角, 或分子的回旋半径, 或相对于天然状态的RMSD来作为这两个变量. 第三个变量是自由能, 可以从体系相对前面所选变量的分布(概率分布)来估计. 当使用三维表示时, 形貌图中的谷表示低自由能区域, 代表体系的亚稳定构象, 丘表示连接这些亚稳定状态的能量势垒.

    FEL的计算可以使用GROMACS命令sham, 为了体验更好, 我们提供了一些额外的脚本将FEL载入Mathmatica(建议使用版本9), 并漂亮地显示出来. 此外, Mathmatica还可以用于识别形貌中的谷, 并确定哪些坐标匹配这些特定的状态. 利用这些信息可以追溯sham的输入, 用于提取轨迹中的时间帧, 进而得到代表性的结构.

    我们将计算两个变量用于表示自由能. 为表示多肽的构象变化, 我们计算其回旋半径, 以及相对于平均结构的RMSD. 我们不使用天然状态, 因为那会破坏后面的对接模拟. 并且, 我们可以基于合并的轨迹计算这些变量. 计算FEL的一个前提是充分采样, 或者使用不同的初始速度重复进行长时间的模拟, 或者如在我们所做的, 从不同的构象开始进行模拟.

    我们如何从MD轨迹中抽取代表性结构用于对接模拟? 给一个前面用来分析多肽构象变化的方法, 并说明这种方法可用于此目的.

    因为FEL分析依赖于在多肽构象空间的充分采样, 并基于构象的分布进行Gibbs自由能估计, 我们需要准备新的轨迹文件(.xtc), 其中的帧数是我们到目前为止所用的十倍. 注意我们没有使用-dt选项. 对所有四条原始轨迹重复下面的步骤, 注意要更改名字.

    gmx trjconv -f traj_helix.xtc -o traj_helix_10-50ns-highresolution.xtc -pbc nojump -n XXXX_peptide_index.ndx -b 10000

    使用和前面完全一样的做法合并高分辨率轨迹(当提示时间戳时选择ccontinue)

    gmx trjcat -f traj_bound_10-50ns-highresolution.xtc traj_helix_10-50ns-highresolution.xtc traj_polypro_10-50ns-highresolution.xtc traj_extended_10-50ns-highresolution.xtc -o traj_concat_10-50ns-highresolution.xtc -cat -settime

    现在我们需要生成FEL计算所需的数据. 先计算合并轨迹的回旋半径Rg, 方法如前:

    gmx gyrate -f traj_concat_10-50ns-highresolution.xtc -s XXXX_peptide.gro -o rg-concatenated_traj.xvg

    为计算相对于平均结构的RMSD, 首先必须重新计算整条轨迹的平均结构. 再次使用rmsf命令生成平均结构(-ox选项), 然后使用平均结构做参考, 利用rms计算RMSD, 方法如前. 当提示时, 对叠合和计算都选择蛋白Protein组.

    gmx rmsf -f traj_concat_10-50ns-highresolution.xtc -s XXXX_peptide.gro -ox average-structure-concat.pdb -o rmsf-per-residue.xvg
    gmx rms -f traj_concat_10-50ns-highresolution.xtc -s average-structure-concat.pdb -o rmsd-vs-average-concat.xvg

    sham命令需要一个文件, 其中包含多个列, 每一列代表不同的坐标. 为生成一个正确的输入文件, 我们使用Perl脚本sham.pl, 并将我们刚生成的两个xvg文件作为它的输入. Perl是一种很类似Python的编程语言. 这个脚本将两个文件作为输入, 并假定数据处于文件的列中. xvg文件的第一列通常是时间, 第二列是感兴趣的坐标. 与Python一样, Perl的计数也是从零开始, 因此选择列数时小心. 下载脚本sham.pl并执行下面的命令. 输出文件只是两个xvg文件的简单合并, 第一列代表时间, 第二列和第三列是特定时刻的Rg值和相对于平均结构的RMSD.

    perl sham.pl -i1 rg-concatenated_traj.xvg -i2 rmsd-vs-average-concat.xvg -data1 1 -data2 1 -o gsham_input.xvg

    现在我们有了正确的sham输入文件, 可用于生成FEL. 如果你还记得, Gibbs自由能可以根据分布概率计算出来, 并且依赖于指定的温度. 使用sham的选项-tsham来指定正确的体系温度(如果你不记得这个值, 可检查成品模拟所用mdp文件中的值).

    gmx sham -f gsham_input.xvg -ls free-energy-landscape.xpm

    如果使用sham计算FEL时, 指定非常高的温度, 会出现什么情况?

    sham的输出可以使用xpm2ps输出到ps文件, 再用pdf阅读器查看. 但是, 当选择代表性结构时这种二维等值线图帮助不大, 需要查看数据时也很麻烦. 为此, 我们准备了一个Mathmatica文件, 可用于FEL的2D或3D可视化和检查. 由于Mathmatica不支持xpm文件, 我们准备了一个Python脚本用于将任意的xpm文件转换为3列数据的文本文件, 这样更易操控. 你可以下载脚本xpm2txt.py, 用它转换xpm文件

    python xpm2txt.py -f free-energy-landscape.xpm -o free-energy-landscape.txt

    下载Mahmatica脚本, 然后打开, 遵照其中说明, 更改开头的文件路径. 如果一切正常, 你可以看到类似下面的图像

    查看FEL并找到其最小点(谷)的位置. 选择你认为具有代表性的一些最小点(5). 你可以通过在2D等值线图上右键点击, 选择Get Coordinates选项来获知对应的坐标. 之后, 当你在图像上移动鼠标时, 会显示每个特定点的坐标. 记录这些坐标, 然后返回原始的sham.pl脚本输出文件查找相应的时间戳. 使用sham生成FEL的过程涉及坐标的分格, 因此对原始值有所近似, 不要期待你能找到精确的对应值.

    从FEL中选择你认为能代表多肽状态的五个点, 记下它们的坐标.

    为了在150000行中查找特定的内容(尽管使用gedit的搜索选项很方便), 我们提供了一个搜索时间戳的脚本. 可以使用它来得到5个时间戳. 你需要提供sham.pl的输出, 并给出两个坐标(顺序要正确). 下面是一个示例

    python get_timestamp.py -f gsham_input.xvg -1 0.657 -2 3.1415

    如果你总是得到No timestamp found...这样的错误信息, 试着使用附近稍有区别的点, 或打开脚本增加变量nval的值到100或更大.

    将代表性结构的时间帧与以前根据RMSD矩阵获得的相比, 这些代表性结构是否与结构簇相符?

    一旦你有了5个点的时间戳, 你可以使用它们从合并的轨迹中抽取代表性结构. 下面的例子抽取45 ns时刻的结构:

    gmx trjconv -f traj_concat_10-50ns-highresolution.xtc -s XXXX_peptide.gro -dump 45000 -o representative_45000ps.pdb

    在PyMOL中打开代表性结构以及结合结构, 并比较二者.

    其他分析

    盐桥

    除了氢键之外, 足球直播免费高清无插件直播不同的带电残基之间也常形成盐桥. 它对足球直播免费高清无插件直播的结构起着重要的稳定作用, 尤其是当它们处于憎水环境中时, 例如足球直播免费高清无插件直播核心. 但是盐桥也能在足球直播免费高清无插件直播暴露的表面形成, 这对于介导足球直播免费高清无插件直播的识别过程往往很重要. 残基间的盐桥分析可以用saltbr命令进行. 程序会输出一系列xvg文件, 给出-/-, +/-(最关注的)和+/+残基间的距离. 当需要时, 通过设置-sep选项, 这个程序可以为每对相反的带电残基产生一个输出文件, 这些残基位于轨迹中的某点, 彼此处于一定的截断距离范围内(这里是0.5 nm, 通过-t选项设置). 这将产生许多文件, 所以分析时最好建立一个单独的目录. 执行以下命令:

    mkdir saltbridge
    cd saltbridge
    gmx saltbr -f ../traj_nojump.xtc -s ../topol.tpr -t 0.5 -sep

    为了更清晰, 删除与钠离子和氯离子有关的文件:

    rm -f CL- NA+ #*

    看看以下残基之间的相互作用:

  • GLU-56 (resp. ASP-56) 和 LYSH-60

  • LYSH-60 和 ASP-88

  • 对于这些相互作用你有什么想法? 残基K60和D88高度保守, 仅仅将UbcH8中的残基D56突变为E56都会改变蛋白的相互作用方式. 其可能原因何在?

    溶剂可及表面积(SAS)

    一个可能感兴趣的性质是溶剂可触及的足球直播免费高清无插件直播表面的面积, 通常称为溶剂可及表面(SAS, solvent accessible surface)或溶剂可及表面面积(SASA, solvent accessible surface area). 溶剂可及表面积可还可进一步细分为亲水性SAS和疏水性SAS. 此外, 结合一些经验参数, SAS还可以用于估计溶剂化自由能. 这四个参数都可用sasa命令计算, 它还可以计算每个残基或原子一段时间内的平均SAS. 输键入下面的命令, 要计算SAS的组和输出组都选择蛋白Protein, 然后查看输出文件.

    gmx sasa -f traj_nojump.xtc -s topol.tpr -o solvent-accessible-surface.xvg -oa atomic-sas.xvg -or residue-sas.xvg

    哪个残基是最容易被溶剂触及的?

    原子间距离的分析和NOE

    前面用到的rmsdist命令也可用于进一步的距离分析. 特别地, 为了解结构及其稳定性, 查看原子间的平均距离及其波动可能会有用处. 使用下面的命令对足球直播免费高清无插件直播中每对原子间的平均距离及其波动进行计算获得矩阵, 然后用上面的步骤重新着色, 并显示.png文件的结果(rmsmean.pngrmsdist.png).

    gmx rmsdist -s topol.tpr -f traj.xtc -rms rmsdist.xpm -mean rmsmean.xpm -dt 10

    简单地解释两个图像: rmsmean表示结构, rmsdist表示柔性/稳定性. 回想前面分析得到的信息并查看结构.

    从实验角度看这些距离同样重要. 原子间距离的边界值可以由NMR实验推断——利用核的Overhauser效应(NOE)——这也是NMR结构计算间的主要推动力. 如果足球直播免费高清无插件直播模型正确, 这个预期得可以到的NOE信号可通过MD模拟来计算. 这些信号与距离密切相关, 特别是r-3和r-6权重的距离. 这些信号也可以用rmsdist计算.

    gmx rmsdist -s topol.tpr -f traj.xtc -nmr3 nmr3.xpm -nmr6 nmr6.xpm -noe noe.dat -dt 10

    nmr3.xpmnmr6.xpm重新着色, 查看这些矩阵, 也查看noe.dat文件.

    模拟战, 最小的1/r3和1/r6平均距离为多少?

    弛豫和序参量

    计算向量的弛豫计算及其自相关. 对足球直播免费高清无插件直播, 通常包括碳骨架N-H或侧链C-H向量. 这种自相关给出了向量能保持其方向多长时间的量度, 因而为表征可变性和稳定性提供了指示. 序参量是自相关的长程限制. 如果一个分子能够自由旋转, 序参量将不可避免地减小到零; 但在分子框架内(内部参考框架), 序参量常有一个明确的值, 这个值表明总体稳定性. 通过叠合足球直播免费高清无插件直播, 这个参考框架可以固定下来, 因此序参量就确定了.

    N-H序参量可以用chi命令计算. 这个程序可以写出一个.pdb文件, 把序参量加入了b因子列, 更容易查看. 该程序也计算J耦合参数, 并可以和NMR结果相比较, 或用于指导NMR实验.

    gmx chi -s topol.tpr -f traj.xtc -o order-parameters.xvg -p order-parameters.pdb -jc Jcoupling.xvg

    看看.xvg文件中的序参量, 并用PyMOL查看.pdb文件, 根据b因子的值给残基着色.

    记下起始与终止的残基, 具有最高序参量值的两个区域的平均值.
    序参量与波动(RMSF)相比如何?

    构象主成分分析

    一个常用的, 但常常理解不深的分析方法是对轨迹进行主成分分析(PCA, PrincipalComponents Analysis). 这种方法有时候也称为本征动力学(ED, essential dynamics), 目的在于识别原子的大尺度集约运动, 从而揭示隐藏在原子波动后面的结构信息, 帮助确定哪些运动方式对足球直播免费高清无插件直播的整体动力学贡献最大.

    在含有N个原子的体系中, 存在3N-6个可能的内部运动方式(另外6个自由度用来描述体系的外部整体的平动和转动). 在MD模拟中, 粒子的波动是互相关联的, 因为粒子彼此之间存在相互作用. 关联的程度有大有小, 但那些通过键直接相连或者彼此位置接近的粒子会产生显著的运动一致性. 粒子运动之间的相关性导致了体系总体波动的结构性, 对大分子而言, 这种结构常常与其功能或(生物)物理特性直接有关. 因此, 研究原子波动的结构性可以为了解这些大分子的行为提供有价值的洞见. 然而, 确实需要一定程度的线性代数方法和多元统计方面的知识才能来解释结果并认识到该方法的缺陷. 特别的, PCA的目标是用新的变量来描述原始数据, 这些新变量是原始变量的线性组合. 这也是PCA存在的最主要问题: 它仅仅使用原子运动间的线性关系来解释问题.

    PCA的第一步是构建协方差矩阵, 它表征了每对原子之间原子运动的共线性程度. 协方差矩阵从定义上说是一个对称矩阵. 接着将这个矩阵对角化, 得到一个特征向量矩阵和特征值的对角矩阵. 每个特征向量描述了粒子的集约运动, 其中的向量值表示相应原子参与运动的程度. 通常, 体系中的大多数(>90%)运动是由10个以下的特征向量或主成分来描述的. 与特征向量相应的特征值等于以集约运动中的每个原子描述的波动的总和, 因此是与特征向量相关的总运动的一个度量, 可用于比较不同情况下足球直播免费高清无插件直播的可变性, 但对不同大小的足球直播免费高清无插件直播进行比较时, 就难以得到有意义的解释. 更多信息请参考Leach的9.14.

    协方差分析会生成很多文件, 因此最好在一个新的目录中运行:

    mkdir COVAR
    cd COVAR

    协方差矩阵的构建和对角化可使用covar命令. 键入下述命令进行分析:

    gmx covar -s ../topol.tpr -f ../traj.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm ascii covariances.dat

    对PCA, 我们主要关心足球直播免费高清无插件直播骨架上的原子, 因此提示时选择backbone组. 构建和对角化协方差矩阵可能需要一些时间.

    协方差矩阵的维数多少? 特征值的总和多少?

    现在, 看看covar.xpm文件中的协方差矩阵.

    xview covar.xpm

    矩阵显示了原子间的协方差. 红色表示两个原子运动一致, 而蓝色表示它们彼此向相反的方向运动. 红色的深度表示波动振幅的大小. 对角线上的值对应于前面得到的RMSF图.

    查看除端基外运动最剧烈的两个部分, 它们之间如何相对运动, 它们各自相对于蛋白的其他部分如何运动?

    从协方差矩阵能得到一组组相关或反相关运动的原子, 从而可以将其集约运动重新写入总运动. 我们前面提到过, 特征值保存在eigenvalues.xvg文件中, 通过相应特征值表示出总波动.

    使用文本编辑器查看eigenvalues.xvg文件, 计算头五个特征向量对总运动的百分比以及累积百分比.

    典型地, 最初五个特征值将捕获主要运动, 这相当于>80%的总运动. 如果解释的总波动较低, 就说明没有明确的集约运动.

    为了解特征值的实际意义, 可用anaeig命令作进一步的分析. 为了更近地看看前两个特征值, 键入以下命令

    gmx 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
    gmx 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中的极端结构显示为蓝-绿色而特征值2为黄-红色.

    spectrum count
    dss hide everything
    show cartoon

    用PyMOL的align命令, 能画出表示两种构象差异的小条.

    align ev1_0001 and (n. c,n,ca),ev1_0002 and (n. c,n,ca),object=diff1
    align ev2_0001 and (n. c,n,ca),ev2_0002 and (n. c,n,ca),object=diff2

    对特征值1而言, 极端结构之间的最大区别是什么? 对特征值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 -n 1500 ev1-vs-ev2.dat > last-7.5ns.dat
    tail -n 1000 ev1-vs-ev2.dat > last-5.0ns.dat
    xmgrace ev1-vs-ev2.dat last-7.5ns.dat last-5.0ns.dat

    投影的形状如何? 它们相互依赖么(分布有所重叠)?
    如果只对最后7.5 ns进行分析, 是否得到相同的特征向量(轴)? 使用最后5.0 ns呢?