Aligning a trajectory to a reference
We use align.AlignTraj to align a trajectory to a frame in a reference trajectory and write it to a file.
Last updated: December 2022
Minimum version of MDAnalysis: 2.0.0
Packages required:
Optional packages for molecular visualisation:
Throughout this tutorial we will include cells for visualising Universes with the NGLView library. However, these will be commented out, and we will show the expected images generated instead of the interactive widgets.
See also
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
import warnings
# suppress some MDAnalysis warnings when writing PDB files
warnings.filterwarnings('ignore')
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)
[4]:
# merged_view = nv.show_mdanalysis(merged)
# merged_view
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.
[5]:
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)
[6]:
# merged1_view = nv.show_mdanalysis(merged1)
# merged1_view
As you can see, the adk_closed and adk_open trajectories still look the same. However, when we load our aligned trajectory from aligned.dcd, we can see them superposed:
[7]:
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)
[8]:
# merged2_view = nv.show_mdanalysis(merged2)
# merged2_view
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)
Select
in_memory=TrueduringAlignTraj(below)
[9]:
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)
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.
[10]:
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['timeseries'].shape)
(98, 3341, 3)
To contrast, we will keep the open conformation static.
[11]:
adk_coords = adk_open.atoms.positions.copy()
adk_coords.shape
[11]:
(3341, 3)
Because there are 98 frames of the aligned Universe, we copy the coordinates of the adk_open positions and stack them.
[12]:
adk_traj_coords = np.stack([adk_coords] * 98)
adk_traj_coords.shape
[12]:
(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.
[13]:
merged_coords = np.hstack([aligned_coords['timeseries'],
adk_traj_coords])
merged2.load_new(merged_coords, format=MemoryReader)
[13]:
<Universe with 6682 atoms>
[14]:
# 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.
[15]:
# from nglview.contrib.movie import MovieMaker
# movie = MovieMaker(
# m2_view,
# step=4, # only render every 4th frame
# output='merged.gif',
# render_params={"factor": 3}, # set to 4 for higher quality
# )
# movie.make()
Writing trajectories to a file
Finally, we can also save this new trajectory to a file.
[16]:
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.