Intermediate Level MDAnalysis hacks

MDAnalysis comes with a number of existing analysis code in the MDAnalysis.analysis module and example scripts (see also the Examples on the MDAnalysis wiki).

RMSD

As an example we will use the MDAnalysis.analysis.rms.rmsd() function from the MDAnalysis.analysis.rms module. It computes the coordinate root mean square distance between two sets of coordinates. For example for the AdK trajectory the backbone RMSD between first and last frame is

>>> u = Universe(PSF,DCD)
>>> bb = u.selectAtoms('backbone')
>>> A = bb.positions  # coordinates of first frame
>>> u.trajectory[-1]      # forward to last frame
>>> B = bb.positions  # coordinates of last frame
>>> rmsd(A,B)
6.8342494129169804

Superposition of structure

In order to superimpose two structures in a way that minimizes the RMSD we have functions in the MDAnalysis.analysis.align module.

The example uses files provided as part of the MDAnalysis test suite (in the variables PSF, DCD, and PDB_small). For all further examples execute first

>>> from MDAnalysis import Universe
>>> from MDAnalysis.analysis.align import *
>>> from MDAnalysis.tests.datafiles import PSF, DCD, PDB_small

In the simplest case, we can simply calculate the C-alpha RMSD between two structures, using rmsd():

>>> ref = Universe(PDB_small)
>>> mobile = Universe(PSF,DCD)
>>> rmsd(mobile.atoms.CA.positions, ref.atoms.CA.positions)
18.858259026820352

Note that in this example translations have not been removed. In order to look at the pure rotation one needs to superimpose the centres of mass (or geometry) first:

>>> ref0 =  ref.atoms.CA.positions - ref.atoms.CA.centerOfMass()
>>> mobile0 =  mobile.atoms.CA.positions - mobile.atoms.CA.centerOfMass()
>>> rmsd(mobile0, ref0)
 6.8093965864717951

The rotation matrix that superimposes mobile on ref while minimizing the CA-RMSD is obtained with the rotation_matrix() function

>>> R, rmsd = rotation_matrix(mobile0, ref0)
>>> print rmsd
6.8093965864717951
>>> print R
[[ 0.14514539 -0.27259113  0.95111876]
 [ 0.88652593  0.46267112 -0.00268642]
 [-0.43932289  0.84358136  0.30881368]]

Putting all this together one can superimpose all of mobile onto ref:

>>> mobile.atoms.translate(-mobile.atoms.CA.centerOfMass())
>>> mobile.atoms.rotate(R)
>>> mobile.atoms.translate(ref.atoms.CA.centerOfMass())
>>> mobile.atoms.write("mobile_on_ref.pdb")

Exercise 5

Use the above in order to investigate how rigid the CORE, NMP, and LID domains are during the transition: Compute time series of the CA RMSD of each domain relative to its own starting structure, when superimposed on the starting structure.

• You will need to make a copy of the starting reference coordinates that are needed for the shifts, e.g.

  NMP = u.selectAtoms("resid 30:59")
  u.trajectory[0]   # make sure to be on initial frame
  ref_com = NMP.selectAtoms("name CA").centerOfMass()
  ref0 = NMP.positions - ref_com
  

  which is then used instead of ref.atoms.CA.centerOfMass() (which would change for each time step).

• I suggest writing a function that does the superposition for a given time step, reference, and mobile AtomGroup to make the code more manageable (or use MDAnalysis.analysis.align.alignto())

Possible solution

The code contains a function superpose() and rmsd(). The latter is marginally faster because we only need the calculated RMSD and not the full rotation matrix. (We are calling the lower-level function MDAnalysis.core.qcprot.CalcRMSDRotationalMatrix() directly, which has somewhat non-intuitive calling conventions). superpose() also does the superposition of the mobile group to the references (but alignto() is actually a more flexible tool for doing this). Otherwise it is mostly book-keeping, which is solved by organizing everything in dictionaries with keys “CORE”, “NMP”, “LID”.

import numpy as np
from MDAnalysis.analysis.align import rotation_matrix
from MDAnalysis.core.qcprot import CalcRMSDRotationalMatrix

def superpose(mobile, xref0, xref_com=None):
    """Superpose the AtomGroup *mobile* onto the coordinates *xref0* centered at the orgin.

    The original center of mass of the reference group *xref_com* must
    be supplied or the superposition is done at the origin of the
    coordinate system.
    """
    # 995 us
    xref_com = xref_com if xref_com is not None else np.array([0., 0., 0.])
    xmobile0 = mobile.positions - mobile.centerOfMass()
    R, rmsd = rotation_matrix(xmobile0, xref0)
    mobile.rotate(R)
    mobile.translate(xref_com)
    return rmsd

def rmsd(mobile, xref0):
    """Calculate optimal RMSD for AtomGroup *mobile* onto the coordinates *xref0* centered at the orgin.

    The coordinates are not changed. No mass weighting.
    """
    # 738 us
    xmobile0 = mobile.positions - mobile.centerOfMass()
    return CalcRMSDRotationalMatrix(xref0.T.astype(np.float64), xmobile0.T.astype(np.float64), mobile.numberOfAtoms(), None, None)


if __name__ == "__main__":
    import MDAnalysis
    import matplotlib
    import matplotlib.pyplot as plt

    # load AdK DIMS trajectory
    from MDAnalysis.tests.datafiles import PSF, DCD
    u = MDAnalysis.Universe(PSF, DCD)

    # one AtomGroup per domain
    domains = {
        'CORE': u.selectAtoms("(resid 1:29 or resid 60:121 or resid 160:214) and name CA"),
        'LID': u.selectAtoms("resid 122-159 and name CA"),
        'NMP': u.selectAtoms("resid 30-59 and name CA"),
        }
    colors = {'CORE': 'black', 'NMP': 'blue', 'LID': 'red'}

    u.trajectory[0]   # rewind trajectory
    xref0 = dict((name, g.positions - g.centerOfMass()) for name,g in domains.iteritems())

    nframes = len(u.trajectory)
    results = dict((name, np.zeros((nframes, 2), dtype=np.float64)) for name in domains)

    for iframe,ts in enumerate(u.trajectory):
        for name, g in domains.iteritems():
            results[name][iframe, :] = u.trajectory.time, rmsd(g, xref0[name])


    # plot
    fig = plt.figure(figsize=(5,5))
    ax = fig.add_subplot(111)
    for name in "CORE", "NMP", "LID":
        data = results[name]
        ax.plot(data[:,0], data[:,1], linestyle="-", color=colors[name], lw=2, label=name)
    ax.legend(loc="best")
    ax.set_xlabel(r"time  $t$ (ps)")
    ax.set_ylabel(r"C$_\alpha$ RMSD from $t=0$, $\rho_{\mathrm{C}_\alpha}$ ($\AA$)")

    for ext in ('svg', 'pdf', 'png'):
        fig.savefig("AdK_domain_rigidity.{0}".format(ext))