Analyzing (frataxin) molecular dynamics trajectories¶
2/5/2021¶
Additional material:
MDAnalysis¶
Python package for reading and analyzing molecular dynamics trajectories.
A simulation ("universe") is defined by its topology file (prmtop) and trajectory (dcd or nc)
In [5]:
import MDAnalysis
traj1 = MDAnalysis.Universe('wt/wt_amber.prmtop', 'wt/wt_amber_1.dcd')
traj2 = MDAnalysis.Universe('wt/wt_amber.prmtop', 'wt/wt_amber_2.dcd')
traj3 = MDAnalysis.Universe('wt/wt_amber.prmtop', 'wt/wt_amber_3.dcd')
traj = MDAnalysis.Universe('wt/wt_amber.prmtop', ['wt/wt_amber_1.dcd','wt/wt_amber_2.dcd','wt/wt_amber_3.dcd'])
I have downsampled the original simulation and aligned the backbone.
cpptraj script:
parm wt_amber.prmtop
trajin wt_amber_1_md3.nc 1 last 10
autoimage
rms first @CA
trajout wt_amber_1.dcd dcd noboxIn [3]:
traj1
Out[3]:
In [4]:
traj1.trajectory
Out[4]:
RMSD¶
Root mean squared deviation (RMSD) $$\sqrt{\frac{\sum_i^n(x_i^a-x_i^b)^2+(y_i^a-y_i^b)^2+(z_i^a-z_i^b)^2}{n}}$$
In [7]:
from MDAnalysis.analysis.rms import rmsd #pull in rmsd function
prot = traj.select_atoms("protein") #only care about protein, not water/ions
refcoord = prot.positions #save initial coordinates
rmsds = [rmsd(refcoord,prot.positions) for ts in traj.trajectory] #coordinates implicitly update as you iterate
In [8]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
In [10]:
plt.plot(rmsds)
Out[10]:
What does this tell us?
In [11]:
traj.trajectory[-1] #set to last frame
lastcoord = prot.positions
lastrmsds = [rmsd(lastcoord,prot.positions) for ts in traj.trajectory]
In [18]:
plt.plot(rmsds,label='first')
plt.plot(lastrmsds,label='last')
plt.legend();
In [37]:
from MDAnalysis.analysis.align import AlignTraj
#align to loop
align = AlignTraj(traj, traj, select='backbone',filename='aligned.dcd')
align.run()
Out[37]:
In [39]:
traj = MDAnalysis.Universe('wt/wt_amber.prmtop','aligned.dcd')
prot = traj.select_atoms("protein") #only care about protein, not water/ions
refcoord = prot.positions #save initial coordinates
rmsds = [rmsd(refcoord,prot.positions) for ts in traj.trajectory]
traj.trajectory[-1] #set to last frame
lastcoord = prot.positions
lastrmsds = [rmsd(lastcoord,prot.positions) for ts in traj.trajectory]
plt.plot(rmsds,label='first')
plt.plot(lastrmsds,label='last')
plt.legend();
In [40]:
#loading all the coordinates into memory takes about 1GB
coords = [prot.positions for ts in traj.trajectory]
In [ ]:
import multiprocessing
n = len(coords)
def row_rmsds(i):
return [rmsd(coords[i],coords[j]) for j in range(n)]
pool = multiprocessing.Pool()
rmat = pool.map(row_rmsds,range(n))
Making the above twice as fast is left as an exercise for the reader.
In [42]:
import seaborn as sns
rmat = np.array(rmat)
sns.heatmap(rmat,square=True,xticklabels=n,yticklabels=n,cmap='YlGnBu',cbar_kws={'label':'RMSD'},vmin=0);
What does this tell us?
PCA¶
In [43]:
from MDAnalysis.analysis import pca
pc = pca.PCA(traj, select='backbone').run()
In [50]:
plt.plot(pc.cumulated_variance)
Out[50]:
In [55]:
backbone = traj.select_atoms('backbone')
backbone
Out[55]:
In [54]:
plt.plot(pc.cumulated_variance[:100])
Out[54]:
In [56]:
transformed = pc.transform(backbone, n_components=10)
In [57]:
traj.trajectory
Out[57]:
In [59]:
import pandas as pd
df = pd.DataFrame(transformed,
columns=['PC{}'.format(i+1) for i in range(10)])
df['sim'] = df.index//1000
df['Time (ps)'] = df.index * u.trajectory.dt
df
Out[59]:
In [76]:
g = sns.pairplot(df,vars=('PC1','PC2','PC3','PC4','PC5'),hue='sim')