Aligning a trajectory to a reference¶
align.AlignTraj to align a trajectory to a frame in a reference trajectory and write it to a file.
Last updated: June 25, 2020 with MDAnalysis 1.0.0
Minimum version of MDAnalysis: 1.0.0
Optional packages for molecular visualisation:
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.
import MDAnalysis as mda from MDAnalysis.analysis import align from MDAnalysis.tests.datafiles import CRD, PSF, DCD, DCD2 import nglview as nv
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.
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.
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 (API docs) efficiently aligns an entire trajectory to a reference. Unlike most other analysis modules,
AlignTraj allows you to write the output of the analysis to a file. This
is because when
Universes are created by loading from a file, changes to frame-by-frame (dynamic) information do not persist when the frame is changed. If the trajectory is not written to a file, or pulled into memory (below),
AlignTraj will have no effect.
align.AlignTraj(adk_closed, # trajectory to align adk_open, # reference select='name CA', # selection of atoms to align filename='aligned.dcd', # file to write the trajectory to match_atoms=True, # whether to match atoms based on mass ).run() # merge adk_closed and adk_open into the same universe merged1 = mda.Merge(adk_closed.atoms, adk_open.atoms) nv.show_mdanalysis(merged1)
As you can see, the
adk_open trajectories still look the same. However, when we load our aligned trajectory from
aligned.dcd, we can see them superposed:
aligned = mda.Universe(PSF, 'aligned.dcd') aligned.segments.segids = ['Aligned'] # rename our segments adk_open.segments.segids = ['Open'] # so they're coloured differently merged2 = mda.Merge(aligned.atoms, adk_open.atoms) nv.show_mdanalysis(merged2)
If you don’t want to write a file, you can also choose to load the entire trajectory into memory. (This is not always feasible depending on how large your trajectory is, and how much memory your device has, in which case it is much more efficient to write an aligned trajectory to a file as above). You can accomplish this in one of two ways:
Load the trajectory into memory in the first place
adk_closed = mda.Universe(PSF, DCD, in_memory=True)
align.AlignTraj(adk_closed, # trajectory to align adk_open, # reference select='name CA', # selection of atoms to align filename='aligned.dcd', # file to write the trajectory to match_atoms=True, # whether to match atoms based on mass in_memory=True ).run() # merge adk_closed and adk_open into the same universe merged3 = mda.Merge(adk_closed.atoms, adk_open.atoms) nv.show_mdanalysis(merged3)
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
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)
/Users/lily/anaconda3/envs/mda-user-guide/lib/python3.7/site-packages/MDAnalysis/analysis/base.py:282: DeprecationWarning: The structure of the `results` array will change in MDAnalysis version 2.0. "MDAnalysis version 2.0.", category=DeprecationWarning
To contrast, we will keep the open conformation static.
adk_coords = adk_open.atoms.positions.copy() adk_coords.shape
Because there are 98 frames of the
aligned Universe, we copy the coordinates of the
adk_open positions and stack them.
adk_traj_coords = np.stack([adk_coords] * 98) adk_traj_coords.shape
(98, 3341, 3)
adk_traj_coords on the second axis with
np.hstack and load the coordinates into memory into the
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
from nglview.contrib.movie import MovieMaker movie = MovieMaker(m2_view, output='merged.gif', in_memory=True) movie.make()
You have to install moviepy, imageio and ffmeg pip install moviepy==0.2.2.11 pip install imageio==1.6
Writing trajectories to a file¶
We can also save this new trajectory to a file.
with mda.Writer('aligned.xyz', merged2.atoms.n_atoms) as w: for ts in merged2.trajectory: w.write(merged2.atoms)
 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.
 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.
 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.
 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.
 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.
 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.