CPPTRAJ聚类/团簇分析教程

by Daniel R. Rose
原文地址

从分子动力学模拟中确定结构布居(structure populations)的一种方法是聚类/团簇分析。聚类/团簇是指将部分具有相似性的数据点从其他数据中区分出来的一个集合。在分子模拟的背景下，这意味着将相似的构象分组在一起。其间的相似性由一个距离标准来度量——距离越小，结构越相似。 一种常用的距离度量便是基于坐标的RMSD。

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

在此示例中，我们将使用CPPTRAJ进行聚类/团簇分析和组合聚类/团簇分析（作为查明是否收敛的方法）。CPPTRAJ支持多种聚类/团簇算法，距离度量，聚类/团簇度量和输出选项。该实施例将使用多维副本交换动力学（MREMD，24×Temperture，8×aMD）产生的四核苷酸rGACC的轨迹。尽管轨迹最初是使用显式的TIP3P水模型生成的，但为了减小轨迹的大小，这里我们将去除溶剂分子。 该数据和分析由Roe等人发表。 Roe et al., J. Phys. Chem. B, 2014, 118(13), pp 2543-3552.

文件

• rGACC.nowat.parm7: 拓扑文件, rGACC + 3 Na+ ions.
• rGACC.MREMD1.nowat.nc.40: rGACC M-REMD 轨迹 @ 300 K, run1.
• rGACC.MREMD2.nowat.nc.40: rGACC M-REMD 轨迹 @ 300 K, run2.
• 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: 通常情况下，生成成对距离矩阵（每帧到每隔一帧的距离）是聚类/团簇计算中非常耗时耗力的部分。“筛分”是一种通过使用“总/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），表明两个轨迹相对较好地收敛。