This is a julia framework to perform mode coupling analysis as described in this paper.
It was designed to provide an infrastructure for possible implementations of other analyses on protein dynamics.
This is a type that keeps a MD trajectory for a protein. Some pre-processing is necessary on the original trajectory.
-
Use your preffered MD trajectory tool (I use vmd)
-
We will only use alpha carbons (name CA) so you should generate a trimmed down version of your trajectory. This will also result in a much smaller file that is easier to process.
-
You should also align the trajectory. in VMD, this is done by RMSD trajectory tool
-
Save the aligned trajectory of CA atoms in xyz format (traj.xyz). This is a plain text format that is easy to parse.
-
Use preprocess.awk script included here
-
It parses xyz format trajectory file and generates a matrix formatted plaintext.
-
Use it as
preprocess.awk traj-aligned-calphas.xyz > R.dat
-
Now you have a matrix with each column representing a degree of freedom in
R.dat
In the end, R is a matrix where every three consecutive column represents x, y, z coordinates of an alpha carbon respectively and every row is a snapshot in time.
R is aligned which means all translation and rotation is removed from the system. This can be later confirmed by the number of zero eigenvalues (it should be 6).
You can read the R matrix from a file using the readdlm function in julia.
R=readdlm("R.dat");
You can then initialize a protein instance like so
p=protein(R);
or you can optionally provide an index
p=protein(R,index);
where index is a vector of integers that map each carbon apha in the R matrix to an actual residue number in the protein. This feature is very useful when you have shifts and gaps in the protein.
length(index) should be size(R,2)/3 which is the number of residues in the protein.
this will give you a protein with
- p.R : shifted coordinates with mean=0 for each dof.
- p.mean : mean for each dof. used for shifting
- p.stdPerRes : std dev. for each residue (in cartesian distance)
- p.covariance : covariance matrix for the protein (calculated for p.R)
I have implemented this in hope that it can also be used for some other analysis other than the mode coupling analysis that we did.
This type represents a whole mode-coupling analysis on a protein
. It is initialized identically as a protein
(or it can take a protein as an input).