4.10.2. Principal Component Analysis (PCA) — MDAnalysis.analysis.pca
- Authors:
John Detlefs
- Year:
2016
- Copyright:
GNU Public License v3
Added in version 0.16.0.
This module contains the linear dimensions reduction method Principal Component Analysis (PCA). PCA sorts a simulation into 3N directions of descending variance, with N being the number of atoms. These directions are called the principal components. The dimensions to be analyzed are reduced by only looking at a few projections of the first principal components. To learn how to run a Principal Component Analysis, please refer to the PCA Tutorial.
The PCA problem is solved by solving the eigenvalue problem of the covariance matrix, a \(3N \times 3N\) matrix where the element \((i, j)\) is the covariance between coordinates \(i\) and \(j\). The principal components are the eigenvectors of this matrix.
For each eigenvector, its eigenvalue is the variance that the eigenvector
explains. Stored in PCA.results.cumulated_variance
, a ratio for each
number of eigenvectors up to index \(i\) is provided to quickly find out
how many principal components are needed to explain the amount of variance
reflected by those \(i\) eigenvectors. For most data,
PCA.results.cumulated_variance
will be approximately equal to one for some \(n\) that is significantly
smaller than the total number of components. These are the components of
interest given by Principal Component Analysis.
From here, we can project a trajectory onto these principal components and attempt to retrieve some structure from our high dimensional data.
For a basic introduction to the module, the PCA Tutorial shows how to perform Principal Component Analysis.
4.10.2.1. PCA Tutorial
The example uses files provided as part of the MDAnalysis test suite
(in the variables PSF
and
DCD
). This tutorial shows how to use the
PCA class.
First load all modules and test data:
import MDAnalysis as mda
import MDAnalysis.analysis.pca as pca
from MDAnalysis.tests.datafiles import PSF, DCD
Given a universe containing trajectory data we can perform Principal Component
Analysis by using the class PCA
and retrieving the principal
components.:
u = mda.Universe(PSF, DCD)
PSF_pca = pca.PCA(u, select='backbone')
PSF_pca.run()
Inspect the components to determine the principal components you would like
to retain. The choice is arbitrary, but I will stop when 95 percent of the
variance is explained by the components. This cumulated variance by the
components is conveniently stored in the one-dimensional array attribute
PCA.results.cumulated_variance
. The value at the ith index of
PCA.results.cumulated_variance
is the sum of the variances from 0 to
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)
From here, inspection of the pca_space
and conclusions to be drawn from the
data are left to the user.
4.10.2.2. Classes and Functions
- class MDAnalysis.analysis.pca.PCA(universe, select='all', align=False, mean=None, n_components=None, **kwargs)[source]
Principal component analysis on an MD trajectory.
After initializing and calling method with a universe or an atom group, principal components ordering the atom coordinate data by decreasing variance will be available for analysis. As an example::
pca = PCA(universe, select='backbone').run() pca_space = pca.transform(universe.select_atoms('backbone'), 3)
generates the principal components of the backbone of the atomgroup and then transforms those atomgroup coordinates by the direction of those variances. Please refer to the PCA Tutorial for more detailed instructions. When using mean selections, the first frame of the selected trajectory slice is used as a reference.
- Parameters:
universe (Universe) – Universe
select (string, optional) – A valid selection statement for choosing a subset of atoms from the atomgroup.
align (boolean, optional) – If True, the trajectory will be aligned to a reference structure.
mean (array_like, optional) – Optional reference positions to be be used as the mean of the covariance matrix.
n_components (int, optional) – The number of principal components to be saved, default saves all principal components
verbose (bool (optional)) – Show detailed progress of the calculation if set to
True
.
- results.p_components
Principal components of the feature space, representing the directions of maximum variance in the data. The column vector p_components[:, i] is the eigenvector corresponding to the variance[i].
Added in version 2.0.0.
- Type:
array, (n_atoms * 3, n_components)
- p_components
Alias to the
results.p_components
.Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use
results.p_components
instead.- Type:
array, (n_atoms * 3, n_components)
- results.variance
Raw variance explained by each eigenvector of the covariance matrix.
Added in version 2.0.0.
- Type:
array (n_components, )
- variance
Alias to the
results.variance
.Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use
results.variance
instead.- Type:
array (n_components, )
- results.cumulated_variance
Percentage of variance explained by the selected components and the sum of the components preceding it. If a subset of components is not chosen then all components are stored and the cumulated variance will converge to 1.
Added in version 2.0.0.
- Type:
array, (n_components, )
- cumulated_variance
Alias to the
results.cumulated_variance
.Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use
results.cumulated_variance
instead.- Type:
array, (n_components, )
Notes
Computation can be sped up by supplying precalculated mean positions.
Changed in version 0.19.0: The start frame is used when performing selections and calculating mean positions. Previously the 0th frame was always used.
Changed in version 1.0.0:
n_components
now limits the correct axis ofp_components
.cumulated_variance
now accurately represents the contribution of each principal component and does not change whenn_components
is given. Ifn_components
is not None or is less than the number ofp_components
,cumulated_variance
will not sum to 1.align=True
now correctly aligns the trajectory and computes the correct means and covariance matrix.Changed in version 2.0.0:
mean_atoms
removed, as this did not reliably contain the mean positions.mean
input now accepts coordinate arrays instead of atomgroup.p_components
,variance
andcumulated_variance
are now stored in aMDAnalysis.analysis.base.Results
instance.Changed in version 2.8.0:
self.run()
can now appropriately useframes
parameter (bug described by #4425 and fixed by #4423). Previously, behaviour was to manually iterate throughself._trajectory
, which would incorrectly handle cases where theframe
argument was passed.- cumulative_overlap(other, i=0, n_components=None)[source]
Compute the cumulative overlap of a vector in a subspace.
This is not symmetric. The cumulative overlap measures the overlap of the chosen vector in this instance, in the
other
subspace.Please cite [Yang2008] if you use this function.
- Parameters:
- Returns:
Cumulative overlap of the chosen vector in this instance to the
other
subspace. 0 indicates that they are mutually orthogonal, whereas 1 indicates that they are identical.- Return type:
See also
Added in version 1.0.0.
- classmethod get_supported_backends()
Tuple with backends supported by the core library for a given class. User can pass either one of these values as
backend=...
torun()
method, or a custom object that hasapply
method (see documentation forrun()
):‘serial’: no parallelization
‘multiprocessing’: parallelization using multiprocessing.Pool
‘dask’: parallelization using dask.delayed.compute(). Requires installation of mdanalysis[dask]
If you want to add your own backend to an existing class, pass a
backends.BackendBase
subclass (see its documentation to learn how to implement it properly), and specifyunsupported_backend=True
.- Returns:
names of built-in backends that can be used in
run(backend=...)()
- Return type:
Added in version 2.8.0.
- property parallelizable
Boolean mark showing that a given class can be parallelizable with split-apply-combine procedure. Namely, if we can safely distribute
_single_frame()
to multiple workers and then combine them with a proper_conclude()
call. If set toFalse
, no backends except forserial
are supported.Note
If you want to check parallelizability of the whole class, without explicitly creating an instance of the class, see
_analysis_algorithm_is_parallelizable
. Note that you setting it to other value will break things if the algorithm behind the analysis is not trivially parallelizable.- Returns:
if a given
AnalysisBase
subclass instance is parallelizable with split-apply-combine, or not- Return type:
Added in version 2.8.0.
- project_single_frame(components=None, group=None, anchor=None)[source]
Computes a function to project structures onto selected PCs
Applies Inverse-PCA transform to the PCA atomgroup. Optionally, calculates one displacement vector per residue to extrapolate the transform to atoms not in the PCA atomgroup.
- Parameters:
components (int, array, optional) – Components to be projected onto. The default
None
maps onto all components.group (AtomGroup, optional) – The AtomGroup containing atoms to be projected. The projection applies to whole residues in
group
. The atoms in the PCA class are not affected by this argument. The defaultNone
does not extrapolate the projection to non-PCA atoms.anchor (string, optional) – The string to select the PCA atom whose displacement vector is applied to non-PCA atoms in a residue. The
anchor
selection is applied togroup
.The resulting atomselection must have exactly one PCA atom in each residue ofgroup
. The defaultNone
does not extrapolate the projection to non-PCA atoms.
- Returns:
The resulting function f(ts) takes as input a
Timestep
ts, and returns ts with the projected structureWarning
The transformation function takes a
Timestep
as input because this is required for Trajectory transformations (“on-the-fly” transformations). However, the inverse-PCA transformation is applied on the atoms of the Universe that was used for the PCA. It is expected that the ts is from the same Universe but this is currently not checked.- Return type:
function
Notes
When the PCA class is run for an atomgroup, the principal components are cached. The trajectory can then be projected onto one or more of these principal components. Since the principal components are sorted in the order of decreasing explained variance, the first few components capture the essential molecular motion.
If N is the number of atoms in the PCA group, each component has the length 3N. A PCA score \(w\_i\), along component \(u\_i\), is calculated for a set of coordinates \((r(t))\) of the same atoms. The PCA scores are then used to transform the structure, \((r(t))\) at a timestep, back to the original space.
\[\begin{split}w_{i}(t) = ({\textbf r}(t) - \bar{{\textbf r}}) \cdot {\textbf u}_i \\ {\textbf r'}(t) = (w_{i}(t) \cdot {\textbf u}_i^T) + \bar{{\textbf r}}\end{split}\]For each residue, the projection can be extended to atoms that were not part of PCA by applying the displacement vector of a PCA atom to all the atoms in the residue. This could be useful to preserve the bond distance between a PCA atom and other non-PCA atoms in a residue.
If there are r residues and n non-PCA atoms in total, the displacement vector has the size 3r. This needs to be broadcasted to a size 3n. An extrapolation trick is used to shape the array, since going over each residue for each frame can be expensive. Non-PCA atoms’ displacement vector is calculated with fancy indexing on the anchors’ displacement vector. index_extrapolate saves which atoms belong to which anchors. If there are two non-PCA atoms in the first anchor’s residue and three in the second anchor’s residue, index_extrapolate is [0, 0, 1, 1, 1]
Examples
Run PCA class before using this function. For backbone PCA, run:
pca = PCA(universe, select='backbone').run()
Obtain a transformation function to project the backbone trajectory onto the first principal component:
project = pca.project_single_frame(components=0)
To project onto the first two components, run:
project = pca.project_single_frame(components=[0,1])
Alternatively, the transformation can be applied to PCA atoms and extrapolated to other atoms according to the CA atom’s translation in each residue:
all = u.select_atoms('all') project = pca.project_single_frame(components=0, group=all, anchor='name CA')
Finally, apply the transformation function to a timestep:
project(u.trajectory.ts)
or apply the projection to the universe:
u.trajectory.add_transformations(project)
Added in version 2.2.0.
- rmsip(other, n_components=None)[source]
Compute the root mean square inner product between subspaces.
This is only symmetric if the number of components is the same for both instances. The RMSIP effectively measures how correlated the vectors of this instance are to those of
other
.Please cite [Amadei1999] and [Leo-Macias2004] if you use this function.
- Parameters:
- Returns:
Root mean square inner product of the selected subspaces. 0 indicates that they are mutually orthogonal, whereas 1 indicates that they are identical.
- Return type:
Examples
You can compare the RMSIP between different intervals of the same trajectory. For example, to compare similarity within the top three principal components:
>>> 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) >>> round(first_interval.rmsip(second_interval, n_components=3), 6) 0.381476 >>> round(first_interval.rmsip(last_interval, n_components=3), 6) 0.174782
See also
Added in version 1.0.0.
- run(start: int | None = None, stop: int | None = None, step: int | None = None, frames: Iterable | None = None, verbose: bool | None = None, n_workers: int | None = None, n_parts: int | None = None, backend: str | BackendBase | None = None, *, unsupported_backend: bool = False, progressbar_kwargs=None)
Perform the calculation
- Parameters:
start (int, optional) – start frame of analysis
stop (int, optional) – stop frame of analysis
step (int, optional) – number of frames to skip between each analysed frame
frames (array_like, optional) –
array of integers or booleans to slice trajectory;
frames
can only be used instead ofstart
,stop
, andstep
. Setting bothframes
and at least one ofstart
,stop
,step
to a non-default value will raise aValueError
.Added in version 2.2.0.
verbose (bool, optional) – Turn on verbosity
progressbar_kwargs (dict, optional) – ProgressBar keywords with custom parameters regarding progress bar position, etc; see
MDAnalysis.lib.log.ProgressBar
for full list. Available only forbackend='serial'
backend (Union[str, BackendBase], optional) –
By default, performs calculations in a serial fashion. Otherwise, user can choose a backend:
str
is matched to a builtin backend (one ofserial
,multiprocessing
anddask
), or aMDAnalysis.analysis.results.BackendBase
subclass.Added in version 2.8.0.
n_workers (int) –
positive integer with number of workers (processes, in case of built-in backends) to split the work between
Added in version 2.8.0.
n_parts (int, optional) –
number of parts to split computations across. Can be more than number of workers.
Added in version 2.8.0.
unsupported_backend (bool, optional) –
if you want to run your custom backend on a parallelizable class that has not been tested by developers, by default False
Added in version 2.8.0.
Changed in version 2.2.0: Added ability to analyze arbitrary frames by passing a list of frame indices in the frames keyword argument.
Changed in version 2.5.0: Add progressbar_kwargs parameter, allowing to modify description, position etc of tqdm progressbars
Changed in version 2.8.0: Introduced
backend
,n_workers
,n_parts
andunsupported_backend
keywords, and refactored the method logic to support parallelizable execution.
- transform(atomgroup, n_components=None, start=None, stop=None, step=None, verbose=False)[source]
Apply the dimensionality reduction on a trajectory
- Parameters:
atomgroup (AtomGroup or Universe) – The AtomGroup or Universe containing atoms to be PCA transformed.
n_components (int, optional) – The number of components to be projected onto. The default
None
maps onto all components.start (int, optional) – The frame to start on for the PCA transform. The default
None
becomes 0, the first frame index.stop (int, optional) – Frame index to stop PCA transform. The default
None
becomes the total number of frames in the trajectory. Iteration stops before this frame number, which means that the trajectory would be read until the end.step (int, optional) – Include every step frames in the PCA transform. If set to
None
(the default) then every frame is analyzed (i.e., same asstep=1
).verbose (bool, optional) –
verbose = True
option displays a progress bar for the iterations of transform.verbose = False
disables the progress bar, just returns the pca_space array when the calculations are finished.
- Returns:
pca_space
- Return type:
array, shape (n_frames, n_components)
Changed in version 0.19.0: Transform now requires that
run()
has been called before, otherwise aValueError
is raised.Changed in version 2.8.0: Transform now has shows a tqdm progressbar, which can be toggled on with
verbose = True
, or off withverbose = False
- MDAnalysis.analysis.pca.cosine_content(pca_space, i)[source]
Measure the cosine content of the PCA projection.
The cosine content of pca projections can be used as an indicator if a simulation is converged. Values close to 1 are an indicator that the simulation isn’t converged. For values below 0.7 no statement can be made. If you use this function please cite [1].
- Parameters:
pca_space (array, shape (number of frames, number of components)) – The PCA space to be analyzed.
i (int) – The index of the pca_component projection to be analyzed.
- Returns:
A float reflecting the cosine content of the ith projection in the PCA
space. The output is bounded by 0 and 1, with 1 reflecting an agreement
with cosine while 0 reflects complete disagreement.
References
- MDAnalysis.analysis.pca.rmsip(a, b, n_components=None)[source]
Compute the root mean square inner product between subspaces.
This is only symmetric if the number of components is the same for
a
andb
. The RMSIP effectively measures how correlated the vectors ofa
are to those ofb
.Please cite [Amadei1999] and [Leo-Macias2004] if you use this function.
- Parameters:
a (array, shape (n_components, n_features)) – The first subspace. Must have the same number of features as
b
. If you are using the results ofPCA
, this is the TRANSPOSE ofp_components
(i.e.p_components.T
).b (array, shape (n_components, n_features)) – The second subspace. Must have the same number of features as
a
. If you are using the results ofPCA
, this is the TRANSPOSE ofp_components
(i.e.p_components.T
).n_components (int or tuple of ints, optional) – number of components to compute for the inner products.
None
computes all of them.
- Returns:
Root mean square inner product of the selected subspaces. 0 indicates that they are mutually orthogonal, whereas 1 indicates that they are identical.
- Return type:
Examples
You can compare the RMSIP between different intervals of the same trajectory. For example, to compare similarity within the top three principal components:
>>> 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) >>> round(pca.rmsip(first_interval.results.p_components.T, ... second_interval.results.p_components.T, ... n_components=3), 6) 0.381476 >>> round(pca.rmsip(first_interval.results.p_components.T, ... last_interval.results.p_components.T, ... n_components=3), 6) 0.174782
Added in version 1.0.0.
- MDAnalysis.analysis.pca.cumulative_overlap(a, b, i=0, n_components=None)[source]
Compute the cumulative overlap of a vector in a subspace.
This is not symmetric. The cumulative overlap measures the overlap of the chosen vector in
a
, in theb
subspace.Please cite [Yang2008] if you use this function.
- Parameters:
a (array, shape (n_components, n_features) or vector, length n_features) – The first subspace containing the vector of interest. Alternatively, the actual vector. Must have the same number of features as
b
.b (array, shape (n_components, n_features)) – The second subspace. Must have the same number of features as
a
.i (int, optional) – The index of eigenvector to be analysed.
n_components (int, optional) – number of components in
b
to compute for the cumulative overlap.None
computes all of them.
- Returns:
Cumulative overlap of the chosen vector in
a
to theb
subspace. 0 indicates that they are mutually orthogonal, whereas 1 indicates that they are identical.- Return type:
Added in version 1.0.0.