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: June 25, 2020 with MDAnalysis 1.0.0

Minimum version of MDAnalysis: 1.0.0

Packages required:

• MDAnalysis ([MADWB11], [GLB+16])

• MDAnalysisTests

Optional packages for molecular visualisation:

• nglview ([NCR18])

See also

• Aligning a trajectory to a frame from itself

• Aligning a structure to another

• RMSD

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.

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.

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_closed and 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:

1. Load the trajectory into memory in the first place

adk_closed = mda.Universe(PSF, DCD, in_memory=True)

2. Select in_memory=True during AlignTraj (below)

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 aligned universe.

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.

adk_coords = adk_open.atoms.positions.copy()
adk_coords.shape

(3341, 3)

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)

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.

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.

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)

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.