4.9.2. 主成分分析(PCA) MDAnalysis.analysis.pca

作者:

约翰·德特尔夫斯

年:

2016

版权所有:

GNU公共许可证v3

在 0.16.0 版本加入.

本模块包含线性降维方法主成分分析(PCA)。主成分分析将模拟分成3N个方向的递减方差，N是原子数。这些方向称为主成分。仅通过查看第一主成分的几个投影来降低要分析的维度。要了解如何运行主成分分析，请参阅 PCA教程 。

通过求解协方差矩阵的特征值问题来解决主成分分析问题 \(3N \times 3N\) 矩阵，其中元素 \((i, j)\) 是坐标之间的协方差 \(i\) 和 \(j\) 。主分量是该矩阵的特征向量。

对于每个特征向量，其特征值是该特征向量所解释的方差。存储在 PCA.results.cumulated_variance ，直到索引的每个特征向量的比率 \(i\) 以快速找出需要多少主成分来解释由这些主成分所反映的方差 \(i\) 特征向量。对于大多数数据， PCA.results.cumulated_variance 对于某些人来说大约等于1 \(n\) 这比组件的总数要少得多。这些是主成分分析给出的感兴趣的成分。

从这里，我们可以将轨迹投影到这些主成分上，并尝试从我们的高维数据中检索一些结构。

有关该模块的基本介绍，请参阅 PCA教程 说明如何执行主成分分析。

4.9.2.1. PCA教程

该示例使用作为MDAnalysis测试套件的一部分提供的文件(在变量中 PSF 和 DCD )。本教程介绍如何使用PCA类。

首先加载所有模块和测试数据：：

import MDAnalysis as mda
import MDAnalysis.analysis.pca as pca
from MDAnalysis.tests.datafiles import PSF, DCD

给定一个包含轨迹数据的宇宙，可以使用类来执行主成分分析 PCA 和检索主成分。：：

u = mda.Universe(PSF, DCD)
PSF_pca = pca.PCA(u, select='backbone')
PSF_pca.run()

检查组件以确定要保留的主要组件。这个选择是随意的，但当95%的差异可以用成分来解释时，我就不说了。按分量计算的累积方差方便地存储在一维数组属性中 PCA.results.cumulated_variance 。第i个索引处的值 PCA.results.cumulated_variance 是从0到I的方差之和。：：

n_pcs = np.where(PSF_pca.results.cumulated_variance > 0.95)[0][0]
atomgroup = u.select_atoms('backbone')
pca_space = PSF_pca.transform(atomgroup, n_components=n_pcs)

从这里开始，检查 pca_space 并将从数据中得出的结论留给用户。

4.9.2.2. 类和函数

class MDAnalysis.analysis.pca.PCA(universe, select='all', align=False, mean=None, n_components=None, **kwargs)[源代码]

MD弹道的主成分分析。

在对宇宙或原子群进行初始化和调用方法后，主成分将对原子坐标数据进行减方差排序，以供分析。例如：

pca = PCA(universe, select='backbone').run()
pca_space = pca.transform(universe.select_atoms('backbone'), 3)

生成原子组主干的主分量，然后根据这些方差的方向转换这些原子组坐标。请参阅 PCA教程 以获取更详细的说明。

参数:

  • universe (Universe) -- 宇宙

  • select (string, optional) -- 用于从原子组中选择原子子集的有效选择语句。

  • align (boolean, optional) -- 如果为True，则轨迹将与参考结构对齐。

  • mean (array_like, optional) -- 用作协方差矩阵平均值的可选参考位置。

  • n_components (int, optional) -- 要保存的主成分数量，默认为保存所有主成分

  • verbose (bool (optional)) -- 如果设置为，则显示计算的详细进度 True 。

results.p_components

特征空间的主成分，表示数据中最大方差的方向。列矢量p_Components [:, i] 是与方差对应的特征向量 [i] 。

在 2.0.0 版本加入.

类型:

数组，(n_原子*3，n_组件)

p_components

的别名 results.p_components 。

自 2.0.0 版本弃用: 将在MDAnalysis 3.0.0中删除。请使用 results.p_components 取而代之的是。

类型:

数组，(n_原子*3，n_组件)

results.variance

由协方差矩阵的每个特征向量解释的原始方差。

在 2.0.0 版本加入.

类型:

数组(n_Components，)

variance

的别名 results.variance 。

自 2.0.0 版本弃用: 将在MDAnalysis 3.0.0中删除。请使用 results.variance 取而代之的是。

类型:

数组(n_Components，)

results.cumulated_variance

由所选分量及其前面的分量之和解释的差异百分比。如果没有选择分量的子集，则存储所有分量，并且累积方差将收敛到1。

在 2.0.0 版本加入.

类型:

数组，(n_组件，)

cumulated_variance

的别名 results.cumulated_variance 。

自 2.0.0 版本弃用: 将在MDAnalysis 3.0.0中删除。请使用 results.cumulated_variance 取而代之的是。

类型:

数组，(n_组件，)

备注

通过提供预先计算的平均位置，可以加快计算速度。

在 0.19.0 版本发生变更: 开始帧在执行选择和计算平均位置时使用。以前总是使用第0帧。

在 1.0.0 版本发生变更: n_components 现在限制正确的轴 p_components 。 cumulated_variance 现在准确地表示每个主成分的贡献，并且在以下情况下不会更改 n_components 是被给予的。如果 n_components 不是一个也不是小于 p_components ， cumulated_variance 将不会等于1。 align=True 现在可以正确对齐轨迹并计算正确的均值和协方差矩阵。

在 2.0.0 版本发生变更: mean_atoms 删除，因为这没有可靠地包含平均头寸。 mean INPUT现在接受坐标数组，而不是原子组。 p_components ， variance 和 cumulated_variance 现在存储在一个 MDAnalysis.analysis.base.Results 实例。

cumulative_overlap(other, i=0, n_components=None)[源代码]

计算子空间中向量的累积重叠。

这是不对称的。累积重叠测量在此实例中所选向量在 other 子空间。

请举出 [Yang2008] 如果您使用此功能。

参数:

  • other (PCA) -- 另一节PCA课程。这肯定已经运行过了。

  • i (int, optional) -- 待分析的特征向量指标。

  • n_components (int, optional) -- 中的组件数量 other 以计算累计重叠。 None 计算所有这些参数。

返回:

在此实例中选择的向量与 other 子空间。0表示它们是相互正交的，1表示它们是相同的。

返回类型:

float

参见

cumulative_overlap()

在 1.0.0 版本加入.

project_single_frame(components=None, group=None, anchor=None)[源代码]

计算将结构投影到选定PC上的函数

将逆PCA变换应用于PCA原子组。(可选)计算每个残基一个位移向量，以将变换外推到不在PCA原子组中的原子。

参数:

  • components (int, array, optional) -- 要投影到的组件。默认设置 None 映射到所有组件。

  • group (AtomGroup, optional) -- 包含要投射的原子的原子组。该预测适用于中的整个残基 group 。PCA类中的原子不受此参数的影响。默认设置 None 不会将投影外推到非PCA原子。

  • anchor (string, optional) -- 用于选择其位移向量应用于残数中的非PCA原子的PCA原子的字符串。这个 anchor 选择将应用于 group 得到的原子选择必须在每个残基中恰好有一个PCA原子 group 。默认设置 None 不会将投影外推到非PCA原子。

返回:

得到的函数f(Ts)将a作为输入 Timestep Ts，并返回具有投影结构的ts。警告：：转换函数使用 :class:``Timestep` 作为输入，因为这是 :ref:``transformations 。然而，逆-PCA变换被应用于用于PCA的宇宙的原子。它是 *expected* 那就是 `ts 来自同一个宇宙，但目前还没有被选中。

返回类型:

function

备注

当为原子组运行PCA类时，将缓存主要组件。然后可以将轨迹投影到这些主分量中的一个或多个上。由于主成分是按解释方差递减的顺序排列的，所以前几个成分捕捉到了基本的分子运动。

如果N是PCA基团中的原子数，则每个组分的长度为3N。PCA评分 \(w\_i\) ，沿组件 \(u\_i\) ，是为一组坐标计算的 \((r(t))\) 由相同的原子组成。然后使用PCA分数来改变结构， \((r(t))\) 一步一个脚印，回到原来的空间。

\[\begin{split}W_{i}(T)=({\extbf r}(T)-\bar{{\extbf r}})\cdot {\extbf u}_i\\ {\extbf r‘}(T)=(w_{i}(T)\cot{\extbf u}_i^T)+ \bar{{\extbf r}}\end{split}\]

对于每个残基，通过将主成分分析原子的位移向量应用于残基中的所有原子，可以将投影扩展到不是主成分分析一部分的原子。这可能有助于保持残基中PCA原子和其他非PCA原子之间的键距离。

如果总共有r个残基和n个非PCA原子，则位移向量的大小为3r。这需要被广播到3n大小。使用外推技巧来整形阵列，因为对每一帧的每个残差进行检查可能是昂贵的。对锚点的位移矢量进行花式索引，计算出非主元分析原子的位移矢量。 index_extrapolate 保存哪些原子属于哪些锚。如果第一个锚的残基中有两个非PCA原子，第二个锚的残基中有三个非PCA原子， index_extrapolate 是 [0、0、1、1、1]

示例

在使用此函数之前运行PCA类。对于主干PCA，运行：：

pca = PCA(universe, select='backbone').run()

获得将主干轨迹投影到第一主分量上的变换函数：

project = pca.project_single_frame(components=0)

要投影到前两个组件，请运行：：

project = pca.project_single_frame(components=[0,1])

或者，变换可以应用于PCA原子，并根据CA原子在每个残基中的平移外推到其他原子：

all = u.select_atoms('all')
project = pca.project_single_frame(components=0,
                                   group=all, anchor='name CA')

最后，将转换函数应用于TimeStep：：

project(u.trajectory.ts)

或将投影应用于宇宙：

u.trajectory.add_transformations(project)

在 2.2.0 版本加入.

rmsip(other, n_components=None)[源代码]

计算子空间之间的均方根内积。

仅当两个实例的零部件数量相同时，这才是对称的。RMSIP有效地测量此实例的向量与 other 。

请举出 [Amadei1999] 和 [Leo-Macias2004] 如果您使用此功能。

参数:

  • other (PCA) -- 另一节PCA课程。这肯定已经运行过了。

  • n_components (int or tuple of ints, optional) -- 要为内积计算的组件数。 None 计算所有这些参数。

返回:

所选子空间的均方根内积。0表示它们是相互正交的，1表示它们是相同的。

返回类型:

float

示例

您可以比较同一轨迹的不同间隔之间的RMSIP。例如，要比较前三个主成分中的相似性：

>>> first_interval = pca.PCA(u, select="backbone").run(start=0, stop=25)
>>> second_interval = pca.PCA(u, select="backbone").run(start=25, stop=50)
>>> last_interval = pca.PCA(u, select="backbone").run(start=75)
>>> first_interval.rmsip(second_interval, n_components=3)
0.38147609331128324
>>> first_interval.rmsip(last_interval, n_components=3)
0.17478244043688052

参见

rmsip()

在 1.0.0 版本加入.

run(start=None, stop=None, step=None, frames=None, verbose=None, *, progressbar_kwargs={})

执行计算

参数:

  • start (int, optional) -- 分析的起始框架

  • stop (int, optional) -- 停止分析框架

  • step (int, optional) -- 要在每个分析的帧之间跳过的帧数

  • frames (array_like, optional) -- 用于切片轨迹的整数或布尔数组； frames 只能使用 instead 的 start ， stop ，以及 step 。设置 both frames 并且至少有一项 start ， stop ， step 设置为非默认值将引发 ValueError 。。。版本已添加：：2.2.0

  • verbose (bool, optional) -- 启用详细信息

  • progressbar_kwargs (dict, optional) -- 带有有关进度条位置等的自定义参数的ProgressBar关键字；请参见 MDAnalysis.lib.log.ProgressBar 查看完整的列表。

在 2.2.0 版本发生变更: 添加了分析任意帧的功能，方法是在 frames 关键字参数。

在 2.5.0 版本发生变更: 增列 progressbar_kwargs 参数，允许修改tqdm进度条的描述、位置等

transform(atomgroup, n_components=None, start=None, stop=None, step=None)[源代码]

在轨迹上应用降维

参数:

  • atomgroup (AtomGroup or Universe) -- 包含要进行PCA变换的原子的原子群或宇宙。

  • n_components (int, optional) -- 要投影到的组件的数量。默认设置 None 映射到所有组件。

  • start (int, optional) -- 用于PCA变换的开始帧。默认设置 None 变为0，即第一个帧索引。

  • stop (int, optional) -- 停止PCA变换的帧索引。默认设置 None 成为轨迹中的总帧数。迭代停止 在此之前 该帧编号，这意味着轨迹将被读取直到结束。

  • step (int, optional) -- 包括每个 step PCA变换中的帧。如果设置为 None (默认)然后分析每个帧(即，与 step=1 )。

返回:

pca_space

返回类型:

array, shape (n_frames, n_components)

在 0.19.0 版本发生变更: 现在的转型要求 run() 之前已被调用过，否则将引发 ValueError 都被养大了。

MDAnalysis.analysis.pca.cosine_content(pca_space, i)[源代码]

测量PCA投影的余弦含量。

如果模拟收敛，则主成分分析投影的余弦含量可以作为一个指标。接近1的值表示模拟未收敛。对于低于0.7的值，不能发表任何声明。如果您使用此函数，请引用 [Hess2002] 。

参数:

  • pca_space (array, shape (number of frames, number of components)) -- 待分析的主成分分析空间。

  • i (int) -- 待分析的主成分投影的指标。

返回:

  • 反映PCA中第i个投影的余弦含量的浮点数

  • 太空。输出以0和1为界，其中1表示一致

  • 带余弦，而0表示完全不一致。

引用

[Hess2002] (1,2,3)

Berk Hess. Convergence of sampling in protein simulations. Phys. Rev. E, 65:031910, Mar 2002. doi:10.1103/PhysRevE.65.031910.

MDAnalysis.analysis.pca.rmsip(a, b, n_components=None)[源代码]

计算子空间之间的均方根内积。

仅当组件数量相同时，这才是对称的 a 和 b 。RMSIP有效地衡量向量之间的相关性 a 是对那些 b 。

请举出 [Amadei1999] 和 [Leo-Macias2004] 如果您使用此功能。

参数:

  • a (array, shape (n_components, n_features)) -- 第一个子空间。必须具有与相同数量的要素 b 。如果您使用的是 PCA ，这是转置的 p_components (即 p_components.T )。

  • b (array, shape (n_components, n_features)) -- 第二个子空间。必须具有与相同数量的要素 a 。如果您使用的是 PCA ，这是转置的 p_components (即 p_components.T )。

  • n_components (int or tuple of ints, optional) -- 要为内积计算的组件数。 None 计算所有这些参数。

返回:

所选子空间的均方根内积。0表示它们是相互正交的，1表示它们是相同的。

返回类型:

float

示例

您可以比较同一轨迹的不同间隔之间的RMSIP。例如，要比较前三个主成分中的相似性：

>>> first_interval = pca.PCA(u, select="backbone").run(start=0, stop=25)
>>> second_interval = pca.PCA(u, select="backbone").run(start=25, stop=50)
>>> last_interval = pca.PCA(u, select="backbone").run(start=75)
>>> pca.rmsip(first_interval.results.p_components.T,
...           second_interval.results.p_components.T,
...           n_components=3)
0.38147609331128324
>>> pca.rmsip(first_interval.results.p_components.T,
...           last_interval.results.p_components.T,
...           n_components=3)
0.17478244043688052

在 1.0.0 版本加入.

MDAnalysis.analysis.pca.cumulative_overlap(a, b, i=0, n_components=None)[源代码]

计算子空间中向量的累积重叠。

这是不对称的。累积重叠测量所选向量在 a ，在 b 子空间。

请举出 [Yang2008] 如果您使用此功能。

参数:

  • a (array, shape (n_components, n_features) or vector, length n_features) -- 包含感兴趣向量的第一个子空间。或者，实际的向量。必须具有与相同数量的要素 b 。

  • b (array, shape (n_components, n_features)) -- 第二个子空间。必须具有与相同数量的要素 a 。

  • i (int, optional) -- 待分析的特征向量指标。

  • n_components (int, optional) -- 中的组件数量 b 以计算累计重叠。 None 计算所有这些参数。

返回:

中所选向量的累积重叠 a 发送到 b 子空间。0表示它们是相互正交的，1表示它们是相同的。

返回类型:

float

在 1.0.0 版本加入.