使用MDAanalysis的PCA（python3.7）

Question

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

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'))

如果出现错误消息“No module named”MDAnalysis.analysis.pca分析.PCA'；'MDAnalysis.analysis.pca分析“不是一个包裹”，意思是：-）。 这意味着您的计算机上没有名为MDAnalysis的包。要解决这个问题，您需要使用pip install命令安装，如果使用conda package manager，则需要使用conda命令安装。查看此链接https://www.mdanalysis.org/pages/installation_quick_start/

看看你从中得到灵感的链接https://www.mdanalysis.org/docs/_modules/MDAnalysis/analysis/pca.html，它证实了我的第一个猜测，我认为我的答案应该允许你使用这个包。你知道吗