python - 使用MDA分析的PCA(python3.7)

Question

我开始在计算化学领域工作，并被要求对分子动力学的一些轨迹进行主成分分析。我被告知要使用 MDAnalysis 包，因此我在他们的页面上找到了一个教程，并尝试遵循它(当然我包括了我自己的输入)以查看它是否有效。我从来没有做过像这个广告这样的分析，我也是 python 编码的新手。我附上了受教程启发的代码。但它对我不起作用，它引发了许多错误，其中一个错误是它无法接受我的输入(拓扑是 PDB 文件，坐标是 XTC 文件)，但这些格式是在支持的格式中列出的，或者其他错误是“PCA 类”未定义。我没有从其他人那里找到太多关于使用 MDA 分析处理 PCA 的信息，因此我希望在这里我能找到曾经做过类似事情的人，并且可以帮助我。我已经尝试过相关的 subreddits，但没有结果。

from __future__ import division, absolute_import
import MDAnalysis as mda
import MDAnalysis.analysis.pca as pca

from six.moves import range
import warnings

import numpy as np
import scipy.integrate

from MDAnalysis import Universe
from MDAnalysis.analysis.align import _fit_to
from MDAnalysis.lib.log import ProgressMeter


u = mda.Universe("L22trial.pdb", "L22trial.xtc") 


PCA = mda.analysis.pca.PCA
class PCA():
    pca = PCA(u, select='backbone').run()
    pca_space =  pca.transform(u.select_atoms('backbone'))

    def __init__(self, universe, select='all', align=False, mean=None,
                 n_components=None, **kwargs):
            super(PCA, self).__init__(universe.trajectory, **kwargs)
            self._u = universe

            self.align = align
            self._calculated = False
            self.n_components = n_components
            self._select = select
            self._mean = mean

    def _prepare(self):
        self._u.trajectory[self.start]
        self._reference = self._u.select_atoms(self._select)
        self._atoms = self._u.select_atoms(self._select)
        self._n_atoms = self._atoms.n_atoms

        if self._mean is None:
            self.mean = np.zeros(self._n_atoms*3)
            self._calc_mean = True
        else:
            self.mean = self._mean.positions
            self._calc_mean = False

        if self.n_frames == 1:
            raise ValueError('No covariance information can be gathered from a single trajectory  frame.\n')

        n_dim = self._n_atoms * 3
        self.cov = np.zeros((n_dim, n_dim))
        self._ref_atom_positions = self._reference.positions
        self._ref_cog = self._reference.center_of_geometry()
        self._ref_atom_positions -= self._ref_cog

        if self._calc_mean:
            interval = int(self.n_frames // 100)
            interval = interval if interval > 0 else 1
            format = ("Mean Calculation Step %(step)5d/%(numsteps)d [%(percentage)5.1f%%]")
            mean_pm = ProgressMeter(self.n_frames if self.n_frames else 1, interval=interval, verbose=self._verbose, format=format)
            for i, ts in enumerate(self._u.trajectory[self.start:self.stop:self.step]):
                if self.align:
                    mobile_cog = self._atoms.center_of_geometry()
                    mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, self._ref_atom_positions, self._atoms, mobile_com=mobile_cog, ref_com=self._ref_cog)
                else:
                    self.mean += self._atoms.positions.ravel()
                mean_pm.echo(i)
            self.mean /= self.n_frames

        self.mean_atoms = self._atoms
        self.mean_atoms.positions = self._atoms.positions

    def _single_frame(self):
        if self.align:
            mobile_cog = self._atoms.center_of_geometry()
            mobile_atoms, old_rmsd = _fit_to(self._atoms.positions, self._ref_atom_positions, self._atoms, mobile_com=mobile_cog, ref_com=self._ref_cog)
            x = mobile_atoms.positions.ravel()
        else:
            x = self._atoms.positions.ravel()
        x -= self.mean
        self.cov += np.dot(x[:, np.newaxis], x[:, np.newaxis].T)

    def _conclude(self):
        self.cov /= self.n_frames - 1
        e_vals, e_vects = np.linalg.eig(self.cov)
        sort_idx = np.argsort(e_vals)[::-1]
        self.variance = e_vals[sort_idx]
        self.variance = self.variance[:self.n_components]
        self.p_components = e_vects[:self.n_components, sort_idx]
        self.cumulated_variance = (np.cumsum(self.variance) / np.sum(self.variance))

        self._calculated = True

    def transform(self, atomgroup, n_components=None, start=None, stop=None, step=None):
        if not self._calculated:
            raise ValueError('Call run() on the PCA before using transform')
        if isinstance(atomgroup, Universe):
            atomgroup = atomgroup.atoms
        if(self._n_atoms != atomgroup.n_atoms):
            raise ValueError('PCA has been fit for {} atoms. Your atomgroup has {} atoms'.format(self._n_atoms, atomgroup.n_atoms))
        if not (self._atoms.types == atomgroup.types).all():
            warnings.warn('Atom types do not match with types used to fit PCA')

        traj = atomgroup.universe.trajectory
        start, stop, step = traj.check_slice_indices(start, stop, step)
        n_frames = len(range(start, stop, step))

        dim = (n_components if n_components is not None else self.p_components.shape[1])
        dot = np.zeros((n_frames, dim))

        for i, ts in enumerate(traj[start:stop:step]):
            xyz = atomgroup.positions.ravel() - self.mean
            dot[i] = np.dot(xyz, self.p_components[:, :n_components])

        return dot

def cosine_content(pca_space, i):
    t = np.arange(len(pca_space))
    T = len(pca_space)
    cos = np.cos(np.pi * t * (i + 1) / T)
    return ((2.0 / T) * (scipy.integrate.simps(cos*pca_space[:, i])) ** 2 /
            scipy.integrate.simps(pca_space[:, i] ** 2))

Answer 1

看来您复制并粘贴了 PCA 类本身。我的猜测是你不需要这样做(我从未使用过该模块，所以这只是一个猜测)。文档( https://www.mdanalysis.org/docs/documentation_pages/analysis/pca.html )似乎表明您唯一需要做的就是以下

    import MDAnalysis as mda
    import MDAnalysis.analysis.pca as pca

    u = mda.Universe("L22trial.pdb", "L22trial.xtc") 

    mypca = pca.PCA(u, select='backbone').run()
    pca_space =  mypca.transform(u.select_atoms('backbone'))

如果您收到错误消息“没有名为‘MDAnalysis.analysis.pca.PCA’的模块；‘MDAnalysis.analysis.pca’不是一个包”，则说明其含义如下:-)。这意味着您的计算机上没有名为 MDAnalysis 的软件包。要解决此问题，您需要使用 pip install 命令进行安装，如果您使用 conda 包管理器，则需要使用 conda 进行安装。请参阅此链接https://www.mdanalysis.org/pages/installation_quick_start/

查看链接https://www.mdanalysis.org/docs/_modules/MDAnalysis/analysis/pca.html您从中受到启发，它证实了我的第一个猜测，我认为我的答案应该允许您使用该软件包。