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

• 作者: Daniel R. Roe, 原文地址
• 2018-10-14 17:15:49 翻译: 刘恒江; 校对: 李继存

使用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

文件

• rGACC.nowat.parm7: 拓扑文件, rGACC和3个Na+离子.
• rGACC.MREMD1.nowat.nc.40: rGACC M-REMD 轨迹, 300 K, 运行1.
• rGACC.MREMD2.nowat.nc.40: rGACC M-REMD 轨迹, 300 K, 运行2.
• cpptraj.cluster.in: 用于单轨迹聚类/团簇分析的CPPTRAJ输入文件
• cpptraj.combined.cluster.in: 用于组合聚类/团簇分析的CPPTRAJ输入文件

请注意, 虽然提供了输入文件, 我们仍鼓励用户使用命令行输入命令以便更好地熟悉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

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

• C0: 聚类/团簇输出数据集合的名称.

聚类/团簇选项

• dbscan: 使用DBSCAN(基于密度)的团簇算法
• minpoints: 形成聚类/团簇的最小点数
• epsilon: 形成聚类/团簇的截断距离
• sievetoframe: 通过与所有聚类/团簇帧(而不仅仅是质心)进行比较来恢复经过筛选的帧.

距离度量选项

• rms: 使用原子的RMSD作为距离度量
• sieve 10: 通常情况下, 生成成对距离矩阵(每帧到其他帧的距离)是聚类/团簇计算中非常耗时耗力的部分. “筛分”是一种通过使用total/10帧进行初始聚类的方法, 可以减少此步骤的开销. 然后将筛分后的帧添加到那些簇中, 作为附加步骤. 请注意, 在大多数情况下, 您还需要使用random关键字来使用随机(而不是有序)筛选; 为保证可重复性, 本教程中并没有使用这个关键词.

输出选项

• out <file>: 将聚类/团簇数量随时间的变化写入file文件. 注意, 由于DBSCAN算法具有”噪声”概念, 因此任何噪声帧都将分配给簇-1(无团簇).
• summary <file>: 将全部聚类/团簇计算摘要写入file
• info <file>: 将详细的聚类/团簇结果(包括DBI, pSF等)写入file
• cpopvtime <file> normframe: 将聚类/团簇布居数目随时间的变化写入file, 布局数目根据总帧数进行归一化.

坐标输出选项

• repout <prefix> repfmt pdb: 将聚类/团簇代表以pdb的格式写入prefix.cX.fmt, 其中X为聚类/团簇编号, fmt为默认格式拓展名
• singlerepout <file> singlerepfmt netcdf: 以NetCDF格式将所有聚类/团簇代表写入file
• avgout <prefix> avgfmt restart: 以Amber重启文件格式将每个聚类/团簇中所有帧的平均值写入prefix.cX.fmt

当读入输入文件后, 输入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

这里有两个新的关键词:

• summarysplit <file>: 将分割的聚类/团簇信息写入文件file
• splitframe 83843: 将聚类/团簇信息于83843帧处分割(即第一段轨迹结束的地方) 当读入输入文件后, 输入run命令开始对轨迹进行处理与分析. 保持耐心, 这将比之前的分析花费更多的时间. 当聚类/团簇分析结束后可以输入quit来退出CPPTRAJ.

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

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

接下来几列分别为:

• #Cluster: 聚类/团簇数
• Total: 聚类/团簇中的总帧数
• Frac: 聚类/团簇的帧数占总帧数的比例
• C#: Grace程序可兼容的聚类/团簇的颜色代码, 从1开始, 14及以上的聚类/团簇都被指定为颜色15
• Color: Grace程序兼容的颜色名称(使用Grace默认方案)
• NumIn1st NumIn2nd: 分别落入第一轨迹和第二轨迹的聚类/团簇的帧数
• Frac1 Frac2: 分别落入第一轨迹和第二轨迹的聚类/团簇的比例
• First1 First2: 最先在第一轨迹或第二轨迹中形成聚类/团簇的帧

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