巴拉巴拉分子模拟里用到的PCA

简单介绍主成分分析

就是那些吧，计算协方差矩阵，把有关联的N个变量根据方差投射成无关联的N个变量，后者就是主成分，变量值换算为主成分的投影值。经过特征分析，得到特征值和特征向量，绘制碎石图，取涵盖原数据80%以上的成分（数目小于N）进行分析，达到降维的目的。

主成分（principal component）

不是很懂，主成分应该是根据协方差矩阵得到的。

协方差矩阵（covariance matrix）

矩阵x的行数=轨迹的帧数；列数由3N个原子的笛卡尔坐标组成。

一般进行PCA的变量是体系的Cα原子或是所有的非H原子。

通常，用轨迹中蛋白质的Cα的笛卡尔坐标进行PCA分析（简称cPCA）在一些特殊情况下，如柔性蛋白在模拟过程中发生折叠，cPCA会产生虚假结果。

①计算第i个原子在轨迹里的平均坐标<xi>

②计算第i个原子和第j个原子的协方差：

(3)式里的n是指样本数

③P*P个cij组成一个矩阵，就是协方差矩阵，协方差矩阵属于对角线矩阵。

cij的意义：

若是原子xi在轨迹里沿着一个方向（e.g.单纯沿着x轴）发生激烈的变化，xi-<xi>将是一个正值，同样，若是xj在轨迹里也沿着相同的方向闹腾，xj-<xj>也会是一个正数，那cij就是一个胖胖的正数（也叫positive）

一直纠结：笛卡尔坐标有xyz三维，算多维变量吧，计算的协方差是一个数值吗？

突然理解，是一个数值，因为向量相乘是标量（数量积）！｡:.ﾟヽ(｡◕‿◕｡)ﾉﾟ.:｡+ﾟ

#比如：
"样本1的第i个ca的坐标是(1,2,3)，它的平均坐标是(1,1,1)"
"同样样本1里的第j个ca的坐标是(2,3,4),它的平均坐标是(2,2,2)"
"那么,它俩的协方差："
[xi1-<xi>]*[xj1-<xj>]=(0,1,2)*(0,1,2)=1+4=5
"然后，是n个样本计算n个数量积，相加除以n得到其协方差"

高中三个数学老师得殴死我...

特征分析（eigenanalysis）

看文献里有一段筛选指数，大意就是具体问题具体分析，PCA选择所有原子不能很好捕获局部运动，这个问题再xxx文献里展示过...阿巴阿巴

找到协方差矩阵里的特征向量和特征值，就算完成PCA的第一步分析

"M*λ=v*λ"
#这里M代表矩阵，v代表特征向量，λ代表特征值

 就是一个矩阵换算，具体的例子如下（虽然它不是一个协方差矩阵）：

这里，3就是一个PC的特征值λ，(1,1,1)就是一个特征向量，其维度=变量数目P

一般，λ越大，代表对应的PC涵盖的信息越多。

绘制碎石图（creating a scree plot）

绘制碎石图可以清楚的查看各个PC涵盖的数据信息量。

①在碎石图里，x轴表示维度（PC的index）

②y轴表示特征值占比（PC1的y值是λ1\sum(λ)）

y值的总和是1，这样能更容易确定数据方差占相关本征向量的百分比。

 文献里介绍的两类碎石图情况，对于图A，方差下降较平缓，这样数据信息近乎平摊给各个PC，降维效果很差。对于图B，方差下降先陡峭，后平稳，前3个PC的曲线下S大于80%，说明这3个PC涵盖数据的大部分信息，我们取前2个或前3个PC进行分析，合理达到数据降维的目的。

自由能面图（Free Energy Landscape）

挑取PC后，下一步就是数据可视化。

我们要算每帧构象在挑选主成分下的投影：

 注意：Zi是一个数值！由于Xi是一个笛卡尔坐标，投影后，都会变成（0,0,1）这种形式（投影到一条轴上），累积后，Zi是一个（0,0,n）的值。

下图很形象说明投影的效果，（不过原子坐标是3D）

若是我们轨迹里有k帧构象，就能获得（投影得）k个PC下的投影值。

要绘制的2D平面图里，x、y轴一般取PC1、PC2（最大的两个特征向量）

计算轨迹中各帧的PC1投影值和PC2投影值，作为一个构象的（x,y）值。若是A帧中的构象a和B帧中的构象b结构很接近，那这俩的（x,y）值就会挨得很近。

若是轨迹中某个构象出现帧数很多，其坐标值就会挨得很近，若是散点图，就表现为密集的点。

图中，轨迹中的构象分摊成三堆，说明这个大分子模拟过程中大体向着三个方向变形。

PS:文献里进行PCA的是一个体系里的两个残基。

上面的图，只是散点图，而自由能井图描述大分子各种构象的自由能大小。

概念是想通的，在能量盆地中，多个极小值被聚集到被更高的能量屏障隔开的区域中，在散点图里，是堆积在一个PC1和PC2围成的圈内的高密度点，这些点代表着低能量构象。

要多加一步的，就是表现出构象（图中的点）的自由能。

通过Boltzmann关系估算构象x的自由能：

"G(x)=-kT*ln(P(x))+常数项"
#关心的是构象之间的相对自由能，而不是绝对自由能;
#故常数项的具体数据并不重要，可以取任意值。
#k常量：1.380649 × 10^-23 J/K
#T取模拟的温度，一般是310K
#KT单位的转换：J,若要转换为kj/mol，要除以4186再乘上阿伏伽德罗常数(6.023*10^23)

要求的，就是这里的P(x) 啦

在cpptraj里，有个hist命令，可对不止一个PC下的投影值进行统计，和归一化（计算分布概率）。

输出的就是组段下的概率值。

 

 上图，就是输出的gnu文件里，PC1&PC2都分成25段，统计概率值。

相应的，需要对第3列的数据代入公式换算，变成构象的相对自由能。

常数取换算后的负的最大值，效果是把最大值变成0

将cpptraj得到的gnu文件里的概率用脚本转化为相对自由能，绘制填色图：

左图是用概率绘制的填色图，右图是修饰后的填色图。

是的，这是同一份数据做出来的图，Σ(っ °Д °;)っ。

红色的区域表示概率为0的区域，也可以理解为难跨越的能垒。左图是蓝色的区域，概率值越高，表示构象落在该区域的次数越多。对应右图的蓝色区域，就是相对自由能越低。

参考文献：

浅谈PCA与g_covar+g_anaeig+ddtdp+sigmaplot做自由能面图的方法

Principal Components Analysis: A Review of its Application on Molecular Dynamics Data