00001 """
00002 **This Module requires numpy**
00003
00004 This module contains functions to analyze trajectories, mainly
00005 similiraty measures baed on RMSDS and pairwise distances.
00006
00007 Author: Niklaus Johner (niklaus.johner@unibas.ch)
00008 """
00009
00010 import ost.mol.alg
00011 import ost.geom
00012 from ost import LogError
00013 import os
00014
00015 def smooth(vec,n):
00016
00017
00018
00019 try:
00020 vec2=vec.copy()
00021 except:
00022 vec2=vec[:]
00023 for i in range(n):
00024 v=0.0
00025 count=1.0
00026 v+=vec[i]
00027 for j in range(n):
00028 count+=1
00029 v+=vec[i+j+1]
00030 for j in range(i):
00031 count+=1
00032 v+=vec[i-(j+1)]
00033 vec2[i]=v/float(count)
00034 for i in range(1,n+1):
00035 v=0.0
00036 count=1.0
00037 v+=vec[-i]
00038 for j in range(n):
00039 count+=1
00040 v+=vec[-(i+j+1)]
00041 for j in range(i-1):
00042 count+=1
00043 v+=vec[-i+j+1]
00044 vec2[-i]=v/float(count)
00045 for i in range(n,len(vec2)-n):
00046 v=vec[i]
00047 for j in range(n):
00048 v+=vec[i+j+1]
00049 v+=vec[i-j-1]
00050 vec2[i]=v/float(2.*n+1.)
00051 return vec2
00052
00053
00054 """
00055 From here on the module needs numpy
00056 """
00057
00058 def RMSD_Matrix_From_Traj(t,sele,first=0,last=-1,align=True,align_sele=None):
00059 """
00060 This function calculates a matrix M such that M[i,j] is the
00061 RMSD (calculated on **sele**) between frames i and j of the trajectory **t**
00062 aligned on sele.
00063
00064 :param t: the trajectory
00065 :param sele: the selection used for alignment and RMSD calculation
00066 :param first: the first frame of t to be used
00067 :param last: the last frame of t to be used
00068 :type t: :class:`~ost.mol.CoordGroupHandle`
00069 :type sele: :class:`~ost.mol.EntityView`
00070 :type first: :class:`int`
00071 :type last: :class:`int`
00072
00073 :return: Returns a numpy N\ :subscript:`frames`\ xN\ :subscript:`frames` matrix,
00074 where N\ :subscript:`frames` is the number of frames.
00075 """
00076 if not align_sele:align_sele=sele
00077 try:
00078 import numpy as npy
00079 if last==-1:last=t.GetFrameCount()
00080 n_frames=last-first
00081 rmsd_matrix=npy.identity(n_frames)
00082 for i in range(n_frames):
00083 if align:
00084 t=ost.mol.alg.SuperposeFrames(t,align_sele,begin=first,end=last,ref=i)
00085 eh=t.GetEntity()
00086 t.CopyFrame(i)
00087 rmsd_matrix[i,:]=ost.mol.alg.AnalyzeRMSD(t,sele,sele)
00088 if i==0:
00089 last=last-first
00090 first=0
00091 return rmsd_matrix
00092 except ImportError:
00093 LogError("Function needs numpy, but I could not import it.")
00094 raise
00095
00096
00097 def PairwiseDistancesFromTraj(t,sele,first=0,last=-1,seq_sep=1):
00098 """
00099 This function calculates the distances between any pair of atoms in **sele**
00100 with sequence separation larger than **seq_sep** from a trajectory **t**.
00101 It return a matrix containing one line for each atom pair and N\ :subscript:`frames` columns, where
00102 N\ :subscript:`frames` is the number of frames in the trajectory.
00103
00104 :param t: the trajectory
00105 :param sele: the selection used to determine the atom pairs
00106 :param first: the first frame of t to be used
00107 :param last: the last frame of t to be used
00108 :param seq_sep: The minimal sequence separation between atom pairs
00109 :type t: :class:`~ost.mol.CoordGroupHandle`
00110 :type sele: :class:`~ost.mol.EntityView`
00111 :type first: :class:`int`
00112 :type last: :class:`int`
00113 :type seq_sep: :class:`int`
00114
00115 :return: a numpy N\ :subscript:`pairs`\ xN\ :subscript:`frames` matrix.
00116 """
00117 try:
00118 import numpy as npy
00119 if last==-1:last=t.GetFrameCount()
00120 n_frames=last-first
00121 n_var=0
00122 for i,a1 in enumerate(sele.atoms):
00123 for j,a2 in enumerate(sele.atoms):
00124 if not j-i<seq_sep:n_var+=1
00125
00126
00127 dist_matrix=npy.zeros(n_frames*n_var)
00128 dist_matrix=dist_matrix.reshape(n_var,n_frames)
00129 k=0
00130 for i,a1 in enumerate(sele.atoms):
00131 for j,a2 in enumerate(sele.atoms):
00132 if j-i<seq_sep:continue
00133 dist_matrix[k]=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
00134 k+=1
00135 return dist_matrix
00136 except ImportError:
00137 LogError("Function needs numpy, but I could not import it.")
00138 raise
00139
00140 def DistanceMatrixFromPairwiseDistances(distances,p=2):
00141 """
00142 This function calculates an distance matrix M(N\ :subscript:`frames`\ xN\ :subscript:`frames`\ ) from
00143 the pairwise distances matrix D(N\ :subscript:`pairs`\ xN\ :subscript:`frames`\ ), where
00144 N\ :subscript:`frames` is the number of frames in the trajectory
00145 and N\ :subscript:`pairs` the number of atom pairs.
00146 M[i,j] is the distance between frame i and frame j
00147 calculated as a p-norm of the differences in distances
00148 from the two frames (distance-RMSD for p=2).
00149
00150 :param distances: a pairwise distance matrix as obtained from
00151 :py:func:`~mol.alg.trajectory_analysis.PairwiseDistancesFromTraj`
00152 :param p: exponent used for the p-norm.
00153
00154 :return: a numpy N\ :subscript:`frames`\ xN\ :subscript:`frames` matrix, where N\ :subscript:`frames`
00155 is the number of frames.
00156 """
00157 try:
00158 import numpy as npy
00159 n1=distances.shape[0]
00160 n2=distances.shape[1]
00161 dist_mat=npy.identity(n2)
00162 for i in range(n2):
00163 for j in range(n2):
00164 if j<=i:continue
00165 d=(((abs(distances[:,i]-distances[:,j])**p).sum())/float(n1))**(1./p)
00166 dist_mat[i,j]=d
00167 dist_mat[j,i]=d
00168 return dist_mat
00169 except ImportError:
00170 LogError("Function needs numpy, but I could not import it.")
00171 raise
00172
00173 def DistRMSDFromTraj(t,sele,ref_sele,radius=7.0,average=False,seq_sep=4,first=0,last=-1):
00174 """
00175 This function calculates the distance RMSD from a trajectory.
00176 The distances selected for the calculation are all the distances
00177 between pair of atoms from residues that are at least **seq_sep** apart
00178 in the sequence and that are smaller than **radius** in **ref_sel**.
00179 The number and order of atoms in **ref_sele** and **sele** should be the same.
00180
00181 :param t: the trajectory
00182 :param sele: the selection used to calculate the distance RMSD
00183 :param ref_sele: the reference selection used to determine the atom pairs and reference distances
00184 :param radius: the upper limit of distances in ref_sele considered for the calculation
00185 :param seq_sep: the minimal sequence separation between atom pairs considered for the calculation
00186 :param average: use the average distance in the trajectory as reference instead of the distance obtained from ref_sele
00187 :param first: the first frame of t to be used
00188 :param last: the last frame of t to be used
00189
00190 :type t: :class:`~ost.mol.CoordGroupHandle`
00191 :type sele: :class:`~ost.mol.EntityView`
00192 :type ref_sele: :class:`~ost.mol.EntityView`
00193 :type radius: :class:`float`
00194 :type average: :class:`bool`
00195 :type first: :class:`int`
00196 :type last: :class:`int`
00197 :type seq_sep: :class:`int`
00198
00199 :return: a numpy vecor dist_rmsd(N\ :subscript:`frames`).
00200 """
00201 if not sele.GetAtomCount()==ref_sele.GetAtomCount():
00202 print 'Not same number of atoms in the two views'
00203 return
00204 try:
00205 import numpy as npy
00206 if last==-1:last=t.GetFrameCount()
00207 n_frames=last-first
00208 dist_rmsd=npy.zeros(n_frames)
00209 pair_count=0.0
00210 for i,a1 in enumerate(ref_sele.atoms):
00211 for j,a2 in enumerate(ref_sele.atoms):
00212 if j<=i:continue
00213 r1=a1.GetResidue()
00214 c1=r1.GetChain()
00215 r2=a2.GetResidue()
00216 c2=r2.GetChain()
00217 if c1==c2 and abs(r2.GetNumber().num-r1.GetNumber().num)<seq_sep:continue
00218 d=ost.geom.Distance(a1.pos,a2.pos)
00219 if d<radius:
00220 a3=sele.atoms[i]
00221 a4=sele.atoms[j]
00222 d_traj=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a3.GetHandle(),a4.GetHandle())[first:last]
00223 if average:d=npy.mean(d_traj)
00224 for k,el in enumerate(d_traj):
00225 dist_rmsd[k]+=(el-d)**2.0
00226 pair_count+=1.0
00227 return (dist_rmsd/float(pair_count))**0.5
00228 except ImportError:
00229 LogError("Function needs numpy, but I could not import it.")
00230 raise
00231
00232 def AverageDistanceMatrixFromTraj(t,sele,first=0,last=-1):
00233 """
00234 This function calcultes the distance between each pair of atoms
00235 in **sele**, averaged over the trajectory **t**.
00236
00237 :param t: the trajectory
00238 :param sele: the selection used to determine the atom pairs
00239 :param first: the first frame of t to be used
00240 :param last: the last frame of t to be used
00241 :type t: :class:`~ost.mol.CoordGroupHandle`
00242 :type sele: :class:`~ost.mol.EntityView`
00243 :type first: :class:`int`
00244 :type last: :class:`int`
00245
00246 :return: a numpy N\ :subscript:`pairs`\ xN\ :subscript:`pairs` matrix, where N\ :subscript:`pairs`
00247 is the number of atom pairs in **sele**.
00248 """
00249 try:
00250 import numpy as npy
00251 except ImportError:
00252 LogError("Function needs numpy, but I could not import it.")
00253 raise
00254 n_atoms=sele.GetAtomCount()
00255 M=npy.zeros([n_atoms,n_atoms])
00256 for i,a1 in enumerate(sele.atoms):
00257 for j,a2 in enumerate(sele.atoms):
00258 if j>i:continue
00259 d=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
00260 M[i,j]=npy.mean(d)
00261 M[j,i]=npy.mean(d)
00262 return M
00263
00264 def AnalyzeDistanceFluctuationMatrix(t,sele,first=0,last=-1):
00265 try:
00266 import numpy as npy
00267 except ImportError:
00268 LogError("Function needs numpy, but I could not import it.")
00269 raise
00270 n_atoms=sele.GetAtomCount()
00271 M=npy.zeros([n_atoms,n_atoms])
00272 for i,a1 in enumerate(sele.atoms):
00273 for j,a2 in enumerate(sele.atoms):
00274 if i>j:continue
00275 d=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last]
00276 M[j,i]=npy.std(d)
00277 M[i,j]=npy.std(d)
00278 return M
00279
00280 def IterativeSuperposition(t,sele,threshold=1.0,initial_sele=None,iterations=5,ref_frame=0):
00281 if initial_sele:current_sele=initial_sele
00282 else: current_sele=sele
00283 for i in range(iterations):
00284 t=ost.mol.alg.SuperposeFrames(t,current_sele,ref=ref_frame)
00285 al=[a for a in sele.atoms if ost.mol.alg.AnalyzeRMSF(t,ost.mol.CreateViewFromAtoms([a]))<threshold]
00286 if len(al)==0:return
00287 current_sele=ost.mol.CreateViewFromAtoms(al)
00288 return current_sele
00289