# 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”.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 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))