Trajectory analysis: MDAnalysis (Part II)
In this lecture, we will perform some analyses on a trajectory of the T4 lysozyme. We will do this using both MDAnalysis and Gromacs, plus VMD to check the hydrogen bonds. You can find the necessary material here.
Analysis with Gromacs
Download the folder.
The output of an MD calculation performed using Gromacs is usually the trajectory in .xtc format. As we saw in lecture n.2, the trajectory is usually processed before starting the data analysis.
Indeed, we remove the effects of the Periodic Boundary Conditions and the diffusive dynamics by using the trjconv command:
gmx trjconv -s md_tpr.tpr -f md_trajectory.xtc -o output_trajectory_noPBC.xtc -pbc mol -center
a further usefull step is to remove the rotations and translation of the protein with respect to the initial condidions:
gmx trjconv -s md_tpr.tpr -f output_trajectory_noPBC.xtc -o output_trajectory_fitted.xtc -fit rot+trans
Now we can start the data analysis using the fitted trajectory output_trajectory_fitted.xtc that we can rename in traj_fitted.xtc.
Apply the commands above to process the 3 trajectories in the folder
Pay attention to specify the name of the system in the output file names! Keep the input names also in the output files.
RMSD
The RMSD is defined as following:
\[RMSD = \sqrt{\frac{1}{N} \sum_{i=1}^N (\vec{r}_i - \vec{r}_{ref})^2 }\]The RMSD can be computed by using the program gmx rms.
It compares two structures by computing the root mean square deviation (RMSD). Each structure from a trajectory (-f) is compared to a reference structure.
The reference structure is taken from the structure file (-s)
gmx rms -f traj_fitted.xtc -s tpr_file.tpr -o output_file.xvg
option: -n index_file.ndx
Calculate the RMSD of the two trajectories (about p53 protein) and plot the results.
We can plot the data using gnuplot (type it in a new shell), and then:
>plot “file.dat” w lp
>replot “file2.dat” w lp
To quit from gnuplot, type q and Enter.
Which information can you extrapolate about the systems from the plots?
Create an index file and calculate the RMSD for a subset of residues
Suggestion:
to create the index file start by using the command gmx make_ndx
gmx make_ndx -f tpr_file.tpr -o index_file.ndx
RMSF
The RMSF is defined as following:
\[RMSF = \sqrt{\frac{1}{T} \sum_{t=1}^T |\vec{r}_i (t) - \langle \vec{r}_{i} \rangle |^2}\]The RMSF (Root mean square fluctuation) can be computed by the program gmx rmsf.
gmx rmsf computes the root mean square fluctuation (RMSF, i.e. standard deviation) of atomic positions in the trajectory (supplied with -f) after (optionally) fitting to a reference frame (supplied with -s).
gmx rmsf -f trajectory_file.xtc -s tpr_file.tpr -o output_file.xvg
option: -n index_file.ndx
Calculate the RMSF of the trajectory (p53_wt) and plot the results.
Which information can you extrapolate about the system from the plots?
Radius of Gyration
The radius of gyration is calculated as:
\[R_g^2 = \frac{\sum_{i=1}^{N} m_i (\vec{r}_i - \vec{r}_{CM})^2 }{\sum_{i=1}^N m_i}\]It is a quantity that quantifies you how compact your protein is.
gmx gyrate -f trajectory_file.xtc -s tpr_file.tpr -o output_file.xvg
Which processes can be monitored by calculating the radius of gyration in time?
Calculate the radius of gyration for the trajectory (GB1), and plot the results over time.
Calculate the RMSD for the trajectory (GB1), and plot the results over time.
How we can interpret these two results? Visualize the trajectory (GB1) in VMD.
Chimera
Now we can try to download another useful software for molecular visualization: Chimera.
Visualize the initial structure of each trajectory (in .pdb format)
In order to be able to extract the first frame of the trajectory and to save it in .pdb format we have to use the gmx trjconv command:
gmx trjconv -f traj_fitted.xtc -s tpr_file.tpr -o first_frame.pdb -dump 0
Pay attention to specify the name of the system in the output file names! keep the input names also for the output files.