From 15f0a714e209cb9f23b4c16cefe03ca548bb542c Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Thu, 17 May 2012 16:12:31 -0400
Subject: [PATCH] Added functions to calculate the distance matrix from a
trajectory from pairwise distances (d-RMSD).
---
modules/mol/alg/pymod/trajectory_analysis.py | 65 ++++++++++++++++++++
1 file changed, 65 insertions(+)
diff --git a/modules/mol/alg/pymod/trajectory_analysis.py b/modules/mol/alg/pymod/trajectory_analysis.py
index c23603205..f71e42740 100644
--- a/modules/mol/alg/pymod/trajectory_analysis.py
+++ b/modules/mol/alg/pymod/trajectory_analysis.py
@@ -82,3 +82,68 @@ def RMSD_Matrix_From_Traj(t,sele,first=0,last=-1):
raise
+def PairwiseDistancesFromTraj(t,sele,first=0,last=-1,seq_sep=1):
+ """
+ This function calculates the distances between any pair of atoms in the
+ EntityView sele with sequence separation larger than seq_sep from a trajectory t.
+ It return a matrix containing one line for each atom pair and N columns, where
+ N is the length of the trajectory.
+ Its inputs are:
+ t : the trajectory (CoordGroupHandle)
+ sele : the EntityView used to determine the atom pairs
+ first=0 : the first frame of t to be used
+ last=-1 : the last frame of t to be used
+ seq_sep=1 : The minimal sequence separation between
+ Returns a numpy NpairsxNframes matrix.
+ """
+ try:
+ import numpy as npy
+ if last==-1:last=t.GetFrameCount()
+ n_frames=last-first
+ n_var=0
+ for i,a1 in enumerate(sele.atoms):
+ for j,a2 in enumerate(sele.atoms):
+ if not j-i<seq_sep:n_var+=1
+ #n_var=sele.GetAtomCount()
+ #n_var=(n_var-1)*(n_var)/2.
+ dist_matrix=npy.zeros(n_frames*n_var)
+ dist_matrix=dist_matrix.reshape(n_var,n_frames)
+ k=0
+ for i,a1 in enumerate(sele.atoms):
+ for j,a2 in enumerate(sele.atoms):
+ if j-i<seq_sep:continue
+ dist_matrix[k]=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
+ k+=1
+ return dist_matrix
+ except ImportError:
+ LogError("Function needs numpy, but I could not import it.")
+ raise
+
+def DistanceMatrixFromPairwiseDistances(distances,p=2):
+ """
+ This function calculates an distance matrix M(NxN) from
+ the pairwise distances matrix D(MxN), where N is the number
+ of frames in the trajectory and M the number of atom pairs.
+ M[i,j] is the distance between frame i and frame j
+ calculated as a p-norm of the differences in distances
+ from the two frames (distance-RMSD for p=2).
+ Inputs:
+ distances : a pairwise distance matrix as obtained from PairwiseDistancesFromTraj()
+ Returns a numpy NxN matrix, where N is the number of frames.
+ """
+ try:
+ import numpy as npy
+ n1=distances.shape[0]
+ n2=distances.shape[1]
+ dist_mat=npy.identity(n2)
+ for i in range(n2):
+ for j in range(n2):
+ if j<=i:continue
+ d=(((abs(distances[:,i]-distances[:,j])**p).sum())/float(n1))**(1./p)
+ dist_mat[i,j]=d
+ dist_mat[j,i]=d
+ return dist_mat
+ except ImportError:
+ LogError("Function needs numpy, but I could not import it.")
+ raise
+
--
GitLab