How to do covariance analysis

From ChengLab
Revision as of 21:06, 9 April 2019 by Kevin (talk | contribs) (From scratch)
Jump to: navigation, search

From scratch

The source code can be found here

  1. Get coordinates using MDTraj
  2. Align the structures
  3. Subtract by mean
  4. Divided by SDs
  5. Combine residual fluctuations (3Nx3N > NxN)

The key computation is to compute the covariances:

X = sel_traj.reshape(traj.n_frames, sel.shape[0] * 3)
X -= X.mean(axis=0)
X /= X.std(axis=0)
emp_cov = np.dot(X.T, X)/traj.n_frames

and to combine intra-residue covariances:

for ii in range(0, sel.shape[0]):
   for jj in range(0, sel.shape[0]):
       emp_cov_residue[ii][jj] = ( emp_cov[ii*3][jj*3] + emp_cov[ii*3+1][jj*3+1] + emp_cov[ii*3+2][jj*3+2] ) / 3