Aligning a trajectory to a reference¶
We use align.AlignTraj to align a trajectory to a reference frame and write it to a file.
Last executed: Feb 06, 2020 with MDAnalysis 0.20.2-dev0
Last updated: February 2020
Minimum version of MDAnalysis: 0.17.0
Packages required:
Optional packages for molecular visualisation:
Note
MDAnalysis implements RMSD calculation using the fast QCP algorithm ([The05]) and a rotation matrix R that minimises the RMSD ([LAT09]). Please cite ([The05]) and ([LAT09]) when using the MDAnalysis.analysis.align module in published work.
[1]:
import MDAnalysis as mda
from MDAnalysis.analysis import align
from MDAnalysis.tests.datafiles import CRD, PSF, DCD, DCD2
import nglview as nv
Loading files¶
The test files we will be working with here are trajectories of a adenylate kinase (AdK), a phosophotransferase enzyme. ([BDPW09]) The trajectories sample a transition from a closed to an open conformation.
[2]:
adk_open = mda.Universe(CRD, DCD2)
adk_closed = mda.Universe(PSF, DCD)
Currently, the proteins are not aligned to each other. The difference becomes obvious when the closed conformation is compared to the open. Below, we set adk_open to the last frame and see the relative positions of each protein in a merged Universe.
[3]:
adk_open.trajectory[-1] # last frame
merged = mda.Merge(adk_open.atoms, adk_closed.atoms)
nv.show_mdanalysis(merged)
Aligning a trajectory with AlignTraj¶
While align.alignto aligns structures, or a frame of a trajectory, align.AlignTraj efficiently aligns an entire trajectory to a reference.
[4]:
align.AlignTraj(adk_closed, adk_open,
select='name CA',
filename='aligned.dcd').run()
[4]:
<MDAnalysis.analysis.align.AlignTraj at 0x1132587f0>
[5]:
aligned = mda.Universe(PSF, 'aligned.dcd')
aligned.segments.segids = ['Aligned']
adk_open.segments.segids = ['Open']
[6]:
merged2 = mda.Merge(aligned.atoms, adk_open.atoms)
nv.show_mdanalysis(merged2)
Copying coordinates into a new Universe¶
MDAnalysis.Merge does not automatically load coordinates for a trajectory. We can do this ourselves. Below, we copy the coordinates of the 98 frames in the aligned universe.
[7]:
from MDAnalysis.analysis.base import AnalysisFromFunction
import numpy as np
from MDAnalysis.coordinates.memory import MemoryReader
def copy_coords(ag):
return ag.positions.copy()
aligned_coords = AnalysisFromFunction(copy_coords,
aligned.atoms).run().results
print(aligned_coords.shape)
(98, 3341, 3)
To contrast, we will keep the open conformation static.
[8]:
adk_coords = adk_open.atoms.positions.copy()
adk_coords.shape
[8]:
(3341, 3)
Because there are 98 frames of the aligned Universe, we copy the coordinates of the adk_open positions and stack them.
[9]:
adk_traj_coords = np.stack([adk_coords] * 98)
adk_traj_coords.shape
[9]:
(98, 3341, 3)
We join aligned_coords and adk_traj_coords on the second axis with np.hstack and load the coordinates into memory into the merged2 Universe.
[10]:
merged_coords = np.hstack([aligned_coords,
adk_traj_coords])
merged2.load_new(merged_coords, format=MemoryReader)
m2_view = nv.show_mdanalysis(merged2)
m2_view
Online notebooks do not show the molecule trajectory, but here you can use nglview.contrib.movie.MovieMaker to make a gif of the trajectory. This requires you to install moviepy.
[11]:
from nglview.contrib.movie import MovieMaker
movie = MovieMaker(m2_view, output='merged.gif', in_memory=True)
movie.make()
Writing trajectories to a file¶
We can also save this new trajectory to a file.
[12]:
with mda.Writer('aligned.xyz', merged2.atoms.n_atoms) as w:
for ts in merged2.trajectory:
w.write(merged2.atoms)
References¶
[1] Oliver Beckstein, Elizabeth J. Denning, Juan R. Perilla, and Thomas B. Woolf. Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of Open↔Closed Transitions. Journal of Molecular Biology, 394(1):160–176, November 2009. 00107. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164, doi:10.1016/j.jmb.2009.09.009.
[2] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. Proceedings of the 15th Python in Science Conference, pages 98–105, 2016. 00152. URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.
[3] Pu Liu, Dimitris K. Agrafiotis, and Douglas L. Theobald. Fast determination of the optimal rotational matrix for macromolecular superpositions. Journal of Computational Chemistry, pages n/a–n/a, 2009. URL: http://doi.wiley.com/10.1002/jcc.21439, doi:10.1002/jcc.21439.
[4] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry, 32(10):2319–2327, July 2011. 00778. URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.
[5] Hai Nguyen, David A Case, and Alexander S Rose. NGLview–interactive molecular graphics for Jupyter notebooks. Bioinformatics, 34(7):1241–1242, April 2018. 00024. URL: https://academic.oup.com/bioinformatics/article/34/7/1241/4721781, doi:10.1093/bioinformatics/btx789.
[6] Douglas L. Theobald. Rapid calculation of RMSDs using a quaternion-based characteristic polynomial. Acta Crystallographica Section A Foundations of Crystallography, 61(4):478–480, July 2005. 00127. URL: http://scripts.iucr.org/cgi-bin/paper?S0108767305015266, doi:10.1107/S0108767305015266.