In this example, we'll compute the solvent accessible surface area of one of the residues in our protien accross each frame in a MD trajectory. We're going to use our trustly alanine dipeptide trajectory for this calculation, but in a real-world situtation you'll probably want to use a more interesting peptide or protein, especially one with a Trp residue.
In [ ]:%matplotlib inline from __future__ import print_function import numpy as np import mdtraj as md
We'll use the algorithm from Shrake and Rupley for computing the SASA. Here's the function in MDTraj:
In [ ]:help(md.shrake_rupley)
In [ ]:trajectory = md.load('ala2.h5') sasa = md.shrake_rupley(trajectory) print(trajectory) print('sasa data shape', sasa.shape)
sasa array contains the solvent accessible surface area for each atom in each frame of the trajectory. Let's sum over all of the atoms to get the total SASA from all of the atoms in each frame.
In [ ]:total_sasa = sasa.sum(axis=1) print(total_sasa.shape)
In [ ]:from matplotlib.pylab import * plot(trajectory.time, total_sasa) xlabel('Time [ps]', size=16) ylabel('Total SASA (nm)^2', size=16) show()
In [ ]:def autocorr(x): "Compute an autocorrelation with numpy" x = x - np.mean(x) result = np.correlate(x, x, mode='full') result = result[result.size//2:] return result / result semilogx(trajectory.time, autocorr(total_sasa)) xlabel('Time [ps]', size=16) ylabel('SASA autocorrelation', size=16) show()
In [ ]: