4 Trajectory Analysis：4.4 Introduction to Principal Component Analysis with CPPTRAJ

##PCA浅谈

    N个原子的柔性大体系运动轨迹需要3N维笛卡尔坐标描述，然而这样高维度的数据很难直观分析，尤其是随着计算能力越来越高，模拟的体系也越来越大、模拟的时间越来越长、内容越来越复杂，这就需要发展更多的分析方法去挖掘出海量信息中感兴趣的信息，而使数据“噪音”被过滤，比如从复杂轨迹中了解蛋白的折叠/去折叠主要路径、酶重要域的大尺度运动行为、体系构象分布倾向、配体/残基突变/质子化等效应对体系构象的影响等等。一种解决途径是将坐标维度约化得尽可能低，同时尽可能使原始高维空间下的信息最大程度保留，这样就方便人直观考察体系的运动过程。统计学上的Principal component analysis(PCA)正适合解决这个问题，一般是将变量选为3N个笛卡尔坐标（也有其它选择方式，诸如二面角PCA等），要处理的数据集就是各帧结构。首先由MD轨迹构建协方差矩阵，然后计算其本征值和本征向量。这3N个本征向量也就是由原始笛卡尔坐标组合而成的新坐标，它们彼此间从协方差的意义上说是不相关的（也就是把轨迹转换到这些新坐标上后，协方差均为0）。本征向量对应的本征值越大说明体系的运动行为越多地体现在这个本征向量坐标上，若最大的几个本征值的和除以协方差矩阵的迹等于0.85，说明在这几个维度上就可以描述体系85%的运动信息。这样有高本征值的本征向量也称主成分(Principal component, PC)。

    实际上一般不将大分子全部N个原子都考虑，比如研究蛋白骨架行为就考虑alpha-C即可，研究某一个区域的运动就选择这个区域的原子即可，不仅计算会省时、省内存，而且如果把过多无关原子也考虑进去，若它们在模拟中运动也比较明显，则感兴趣的原子的运动反倒被掩盖了。另外，做PCA前必须将整个轨迹向指定帧的结构做align以消除平动和转动的影响，前文所谓“运动”皆是指体系的内运动，否则分子内运动行为会被整体运动掩盖住。（实际上即便做了align，整体运动的成分也并没彻底消除，因而出现了基于内坐标的PCA方法，即dihedral PCA）。限于本文篇幅和目的不再对PCA进行更具体讨论，可以参考Principal Components Analysis: A Review of its Application on Molecular Dynamics Data。下载地址：/usr/uploads/file/20150610/20150610200236_48141.pdf

Principal component analysis with CPPTRAJ

PCA的输入是从3D位置坐标的时间序列计算出的坐标协方差矩阵，因此PC将代表系统经历的某些运动模式，其中第一个PC代表主导运动。

需要注意的是，尽管PCA对于考察系统的动力学是有用的，但是系统真实的运动往往是多个PC的集合。因此一个单独的PC可能代表了某种运动，但是系统如何经历这种运动是不清楚的。

步骤一：计算坐标协方差矩阵

如上所述，PCA的输入将是一个坐标协方差矩阵。这个矩阵的条目是每个原子的X、Y和Z分量之间的协方差，因此最终的矩阵大小为[3 * # selected atoms] X [3 * # selected atoms]。这意味着为了正确地填充这个矩阵，我们将需要至少与我们的行/列(即3 * #选择原子)一样多的输入帧来计算坐标协方差矩阵。

CPPTRAJ脚本如下：

#由于轨迹可能具有不同数量的原子，我们需要为每个轨迹加载拓扑文件(prmtop)信息。我们用名称处理程序[name]功能指定每个拓扑。

parm cpu/cpu.prmtop [cpu]

parm gpu/gpu.prmtop [gpu]

#现在可以用' [cpu] '或者‘[gpu]’来代替完整的拓扑文件名。在加载了拓扑之后，我们现在希望用它们各自的拓扑加载轨迹。在这个例子中，每个轨迹都是10001帧长，所以我们的合并数据集将是20002帧，前10001帧对应于cpu帧，从10002到20002对应于GPU帧。

trajin cpu/cpu.nc parm [cpu]

trajin gpu/gpu.nc parm [gpu]

#现在，如果我们只是从原始轨迹数据中计算坐标协方差矩阵，我们将不仅捕获内部运动，而且还捕获系统的全局旋转和平移。由于在这个例子中，我们只对内部动力学感兴趣，我们需要去除旋转/平移运动，我们将通过对参考结构进行坐标RMS最佳拟合来实现，在我们的例子中，这将是平均坐标。为了生成平均坐标，我们首先通过对第一帧进行RMS拟合，使用除氢原子外的所有DNA原子（残基1-36），将各帧放在一个共同的参考中。

rms first :1-36&!@H=

#然后，我们在加载的整个框架上创建一个平均结构，并将坐标保存为 "cpu-gpu-average"，该坐标随后可作为参考结构使用。请注意，如果我们愿意，我们也可以将平均坐标写入CPPTRAJ支持的任何格式的文件中。

average crdset cpu-gpu-average

#CPPTRAJ有一个 "数据集 "的概念，可以有多种格式。这里我们创建了一个坐标数据集，它将保存所有的输入帧。这使我们以后可以对整个数据集采取行动，而不必从磁盘上重新读取轨迹。我们把加载的帧坐标数据集称为cpu-gpu-trajectories。

createcrd cpu-gpu-trajectories

#上面的命令将生成我们想要作为参考的平均结构。现在要运行上述程序，在不退出程序的情况下，我们输入运行命令。

run

#现在我们已经生成了平均坐标，'cpu-gpu-average'，同时也保存了输入轨迹的帧。现在我们可以将保存的轨迹帧与平均坐标进行RMS拟合，以消除全局旋转/平移运动。这是用crdaction命令完成的。

crdaction cpu-gpu-trajectories rms ref cpu-gpu-average :1-36&!@H=

#现在我们使用矩阵命令来生成坐标协方差矩阵，我们将其命名为 "cpu-gpu-covar"。

crdaction cpu-gpu-trajectories matrix covar \

  name cpu-gpu-covar :1-36&!@H=

#请注意，反斜杠'\'字符可以用来延续一个行；这对于使复杂的输入更加可读很有用。

步骤二：计算主成分和坐标投影

#现在我们有了矩阵，我们可以通过对角线化来获得PC；这将给我们提供特征向量（即PC）和特征值（即每个PC的 "权重"）。首先，我们将获得前三个特征向量。

runanalysis diagmatrix cpu-gpu-covar out cpu-gpu-evecs.dat \

    vecs 3 name myEvecs \

    nmwiz nmwizvecs 3 nmwizfile dna.nmd nmwizmask :1-36&!@H=

#runanalysis命令告诉cpptraj立即运行'diagmatrix'，而不是将其加入分析队列。nmwiz和相关的关键字产生的输出可以用VMD的'nmwiz'插件来实现主成分数据的可视化。

#一旦这个命令完成，文件'cpu-gpu-evecs.dat'和数据集'myEvecs'将包含特征向量（PC）和特征值（这些数据被统称为 "特征模式数据"）。将这些数据写到文件中通常是有用的，因为它们可以在以后被读回来进行进一步分析。我们现在可以沿PC投射轨迹坐标，看每一帧的坐标沿每个主成分 "匹配 "的程度。我们可以对来自原始单个轨迹的帧进行这一操作，这将使我们能够比较每个单个轨迹的运动的匹配程度。

#请注意，用于投影的帧与用于生成坐标协方差矩阵的帧是一样的，这一点很关键

crdaction cpu-gpu-trajectories projection CPU modes myEvecs \

  beg 1 end 3 :1-36&!@H= crdframes 1,10001

crdaction cpu-gpu-trajectories projection GPU modes myEvecs \

  beg 1 end 3 :1-36&!@H= crdframes 10002,last

#在这种情况下，来自cpu轨迹的投影被命名为 "CPU"，来自gpu轨迹的投影被命名为 "GPU"。一旦这些数据被生成，我们就可以用hist分析命令对每个轨迹的三个计算出的投影制作归一化的直方图，然后用运行命令来实际操作

hist CPU:1 bins 100 out cpu-gpu-hist.agr norm name CPU-1

hist CPU:2 bins 100 out cpu-gpu-hist.agr norm name CPU-2

hist CPU:3 bins 100 out cpu-gpu-hist.agr norm name CPU-3

hist GPU:1 bins 100 out cpu-gpu-hist.agr norm name GPU-1

hist GPU:2 bins 100 out cpu-gpu-hist.agr norm name GPU-2

hist GPU:3 bins 100 out cpu-gpu-hist.agr norm name GPU-3

run

数据集索引（如':1'）指的是主成分，因此'CPU:1'是CPU等的第一个主成分。

图略

步骤三：主成分的可视化

#现在，这一阶段的分析已经完成，我们可以发出clear all命令，清除所有存储的数据，这样我们就可以用一块 "干净的石板 "做进一步分析。

clear all

#我们的下一步是将特征模式的波动可视化。要做到这一点，我们使用readdata命令读入生成的带有特征向量的文件。

readdata cpu-gpu-evecs.dat name Evecs

#加载一个拓扑结构，并对其进行修改，使其与坐标协方差矩阵（以及随后的特征向量）的计算方式相匹配。

parm cpu/cpu.prmtop

parmstrip !(:1-36&!@H=)

parmwrite out cpu-gpu-modes.prmtop 

#创建一个沿第一个PC运动的NetCDF伪轨迹文件。注意！需要通过观察PC投影的直方图来选择最小pcmin和最大值pcmax，如果设置的过大，则运动方向延伸过大，轨迹会很夸张。

runanalysis modes name Evecs trajout cpu-gpu-mode1.nc \

  pcmin -100 pcmax 100 tmode 1 trajoutmask :1-36&!@H= trajoutfmt netcdf

#现在你可以在Chimera / VMD中打开文件，观看第一种运动模式的电影

实践例子：

parm ../../../6z5r_wm5_final_U2H.prmtop

trajin ../../../hold_npt4/hold_npt4.dcd 1 last 1

trajin ../../../prod/prod.dcd 1 4900 1                #导入拓扑和轨迹

rms first :925-1005@CA                         

average crdset sec1-average

average sec1_average.pdb                     #获得平均结构

createcrd sec1-trajectories                    #获得轨迹集合（的代名词）

run

crdaction sec1-trajectories rms ref sec1-average :925-1005@CA     #将轨迹集合与平均结构进行叠合，消除旋转和平移

crdaction sec1-trajectories matrix out sec1_matrix.dat covar \

  name sec1-covar :925-1005@CA                       #生成坐标协方差矩阵

runanalysis diagmatrix sec1-covar out sec1-evecs.dat \

  vecs 3 name myEvecs \

  nmwiz nmwizvecs 3 nmwizfile sec1.nmd nmwizmask :925-1005@CA     #矩阵对角化并获得前3个PC（特征向量）和每个PC对应的权重（特征值）

crdaction sec1-trajectories projection sec1 modes myEvecs \

  beg 1 end 3 :925-1005@CA crdframes 1,5000             #获得轨迹投影并命名为sec1

hist sec1:1 bins 100 out sec1-hist.agr norm name sec1-1

hist sec1:2 bins 100 out sec1-hist.agr norm name sec1-2

hist sec1:3 bins 100 out sec1-hist.agr norm name sec1-3      #根据投影制作归一化的直方图，横跨X轴的范围越大说明越占主导（即PC1最大），这个横跨x轴的范围在后面做特征模式的运动伪动画的时候会使用

run

clear all

readdata sec1-evecs.dat name Evecs      #使用readdata命令读取带有特征向量的文件

parm ../../../6z5r_wm5_final_U2H.prmtop

parmstrip !(:925-1005@CA)

parmwrite out sec1-modes.prmtop        #修改拓扑结构，使其与坐标协方差矩阵相匹配

runanalysis modes name Evecs trajout sec1-mode1.nc \

  pcmin -6 pcmax 4 tmode 1 trajoutmask :1-81@CA trajoutfmt netcdf      #生成沿PC1运动的NetCDF伪轨迹文件，需要通过观察PC投影的直方图来选择最小值pcmin和最大值pcmax。

#####################################################

# Now you can open the files:                       #

# cpu-gpu-modes.prmtop                                         #               

# cpu-gpu-modes.nc                                         #

# in Chimera / VMD and watch the movie                    #

# which shows the first mode of motion                       #

#####################################################

#获取某一个PC对原子运动信息的体现

#cpptraj不输出矩阵的所有本征值，所以无法使用 PC1/all_PC 来获取某一个PC对整体运动的贡献。

#这里使用python计算，输出文件第一行是矩阵的迹（所有本征值的和），后面列表是每一个本征值大小。

#!coding=utf-8
import numpy as np

f = open('../sec1_matrix.dat')  
lines = f.readlines()  
list_all = []
for line in lines:
    list = map(float,(line.strip('\n').split()))
    list_all.append(list)
print ("list_all over")
#print (list_all)
A = np.matrix(list_all)
print ("matrix_all over")
#print (A)
b = np.linalg.eig(A)
print ("eig_all over")
#print (type(b[0]))
ff = open("eig.dat",'w')
ff.write(str(np.sum(b[0]))+'\n')
ff.write(str(b[0]))
ff.close()
 

参考：

http://sobereva.com/73

https://amberhub.chpc.utah.edu/introduction-to-principal-component-analysis/