AMBER 2015伦敦研讨会 - 组合聚类分析示例

类别:    标签: amber   阅读次数:   版权: (CC) BY-NC-SA

使用CPPTRAJ进行聚类/团簇分析

从分子动力学模拟中确定结构布居数目(structure populations)的一种方法是聚类/团簇分析. 聚类/团簇化是一种对数据进行分区的方法, 也就是将部分具有相似性的数据点从其他数据点中区分出来形成一个群集/团簇, 这样群集内数据点之间的相似性要比它们与群集外数据点之间的相似性更大. 在分子模拟的背景下, 这意味着将相似的构象分组在一起. 其间的相似性由一个距离标准来度量——距离越小, 结构越相似. 一种常用的距离度量便是基于坐标的RMSD.

值得注意的是, 并没有一种所谓”正确”的方法来进行聚类/团簇分析. 有许多不同的算法和距离度量可用, 并且不同的方法组合对于某些系统可能有更好的效果. 聚类/团簇分析往往需要不断地试错才能获得满意的结果. 本教程只是聚类/团簇分析的一个示例. 如需更深入地讨论MD轨迹的聚类/团簇分析, 可以阅读Shao等人2007年关于这一主题的论文.

在此示例中, 我们将使用CPPTRAJ进行聚类/团簇分析以及组合的聚类/团簇分析(作为确定收敛性的一种手段). CPPTRAJ支持多种聚类/团簇算法, 距离度量, 聚类/团簇度量和输出选项. 该示例将使用四核苷酸rGACC的轨迹, 它是由多维副本交换动力学(MREMD, 24个Temperture, 8个aMD)产生的. 尽管最初的轨迹中包含了显式的TIP3P水分子, 但我们将去除溶剂分子以减小轨迹的大小. 该数据和分析首先(部分地)由Roe等人发表在J. Phys. Chem. B, 2014, 118(13), pp 2543-3552

文件

请注意, 虽然提供了输入文件, 我们仍鼓励用户使用命令行输入命令以便更好地熟悉CPPTRAJ工作流程和命令选项.

使用CPPTRAJ进行聚类/团簇分析

第一步: 加载拓扑与轨迹文件

我们将首先对单个轨迹进行聚类/团簇分析. 在终端中输入cpptraj, 并输入以下命令加载拓扑与轨迹文件:

Parm rGACC.nowat.parm7
Trajin rGACC.MREMD1.nowat.nc.40

由于这些轨迹中包含离子, 而我们只对rGACC本身感兴趣, 因此我们可以使用strip命令删除离子(确保它们不会出现在任何输出结构中). 我们还将使用outprefix关键字, 以便将去除离子后相应的拓扑文件命名为noions.rGACC.nowat.parm7. 请注意, 这是可选步骤, 并不是后续聚类/团簇分析所必需的.

strip :Na+ outprefix noions

第二步: Clustering命令

接下来我们说明一下Clustering命令, 该命令下有许多选项, 主要分为四大类: 聚类/团簇算法, 距离度量, 输出以及坐标输出.

cluster C0 \
        dbscan minpoints 25 epsilon 0.9 sievetoframe \
        rms :1@N2,O6,C1',P,:2@H2,N6,C1',P,:3@O2,H5,C1',P,:4@O2,H5,C1',P \
        sieve 10 random \
        out cnumvtime.dat \
        sil Sil \
        summary summary.dat \
        info info.dat \
        cpopvtime cpopvtime.agr normframe \
        repout rep repfmt pdb \
        singlerepout singlerep.nc singlerepfmt netcdf \
        avgout Avg avgfmt restart

我们将简要介绍一下各选项的一些细节:

聚类/团簇选项

距离度量选项

输出选项

坐标输出选项

当读入输入文件后, 输入run命令开始对轨迹进行处理与分析. 即便对相当先进的CPU来说, 这个过程也将运行大约一分钟, 所以请保持耐心. 当聚类/团簇分析结束后可以输入quit来退出CPPTRAJ. 现在你便得到了许多包含聚类结果的输出文件. 检查summary.dat文件便可发现共有16个聚类/团簇. 有关输出选项的更详细介绍参见Amber手册.

通常我们最感兴趣的输出是cpopvtime.arg(xmgrace格式)文件中聚类/团簇布居数目随时间的变化. 可以看到, 在轨迹的开始阶段聚类/团簇布居随时间快速变化, 随后, 聚类/团簇布居逐渐趋于稳定, 最终达到约为70000帧. 这表明聚类/团簇布居达到了平衡, 其结果可能适用于热力学等方面的分析. 不过, 确定这一点的更好方法是将当下结果与独立运行的其他结果进行比较.

使用CPPTRAJ进行组合聚类/团簇分析

“收敛”是分子模拟中的一个重要概念, 它是一种衡量构象势能面采样有效性的方法. 从单个点(即构象)开始的模拟在一段时间后, 动力学上可能会陷于某个局部能量极小点, 从单次模拟中, 我们很难确定这种情况是否发生. 然而, 如果以不同的初始条件和/或不同的初始构象(理想情况下两者都不同)对体系进行其他模拟, 最终应该会对相同的构象进行采样, 并得到相同的布局分布. 一旦发生这种情况, 我们便认为模拟已经”收敛”. 当两个或多个模拟已经收敛且结构布居不再发生显著变化时, 我们便可以从中计算具有一定确定性的热力学量.

聚类/团簇分析可以被用来比较不同运行模拟结果的结构布居, 并作为确定是否收敛的判据. 但比较不同轨迹之间的聚类/团簇分析结果是项较有挑战性的工作. 不同轨迹中可能存在完全不同的布居, 有些特殊构象也可能只存在于某一轨迹中. 为了便于比较不同轨迹之间的聚类/团簇, CPPTARJ程序提供了”组合聚类/团簇分析”功能, 其中聚类/团簇分析可以在两个(或者更多)组合轨迹上进行, 随后基于原始轨迹对结果进行划分.

第一步: 载入拓扑与轨迹文件

由于我们在聚类/团簇分析之前需要将两个轨迹组合在一起, 因此我们需要知道每个轨迹中分别有多少帧, 以便使cluster命令明确何处去划分聚类/团簇分析的结果. 这可以通过一下命令实现:

cpptraj -p rGACC.nowat.parm7 -y rGACC.MREMD1.nowat.nc.40 -tl

这里-p选项载入拓扑文件, -y选项载入轨迹文件, -tl选项用以输出帧数. 此处的输出显示为:

Frames: 83843

所以我们应该在组合轨迹的第83843帧处将其划分开.

在命令行中再次输入cpptraj启动CPPTRAJ. 显示交互信息后, 按照之前同样的方式载入拓扑与轨迹文件(并去除离子):

parm rGACC.nowat.parm7
trajin rGACC.MREMD1.nowat.nc.40
trajin rGACC.MREMD2.nowat.nc.40
strip :Na+ outprefix noions

请注意, 这是可选步骤, 并不是后续聚类/团簇分析中所必需的.

第二步: 组合聚类/团簇分析中的Clustering命令

现在开始输入clustering命令选项, 大部分与之前的相同, 除了我们加入了进行组合分析的关键词.

cluster C1 \
        dbscan minpoints 25 epsilon 0.9 sievetoframe \
        rms :1@N2,O6,C1',P,:2@H2,N6,C1',P,:3@O2,H5,C1',P,:4@O2,H5,C1',P \
        sieve 20 \
        out combined.cnumvtime.dat \
        info combined.info.dat \
        summarysplit split.dat splitframe 83843

这里有两个新的关键词:

组合聚类/团簇分析的主要输出信息在split.dat文件中, 文件首两行说明了组合轨迹是如何被分割的以及每部分的总帧数.

# 1st <= 83843 < 2nd
# 1st= 83843  2nd= 108887

接下来几列分别为:

第一和第二轨迹之间的聚类/团簇布居非常一致(布居分数的最大绝对差值约为0.015), 表明两个轨迹收敛得相对较好.

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 三维空间随机旋转
后一篇: GROMACS输入文件模板

访问人次(2015年7月 9日起): | 最后更新: 2018-11-16 03:57:15 UTC | 版权所有 © 2008 - 2018 Jerkwin