PCA_analysis_of_MD_simulations-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Part 2: Analysis of replicate simulations of dynein light chain ##\n",
    "In this section of the workshop we move from a comparative analysis of MD trajectories of two different proteins, to the comparative analysis of multple independent simulations of the same protein.\n",
    "\n",
    "You have ben provided with eight trajectory files (Amber netcdf format) from eight independent 25 nanosecond simulations of the dyenin light chain LC8 (PDB code 1rhw) - a small protein of 89 amino acids. In this case the trajectory files contain all protein atoms, not just the C-alpha atoms, and are sampled every 10 ps (2500 frames per trajectory). You will examine convergence and sampling in this ensemble. \n",
    "\n",
    "### Part 2A: RMSD Analysis\n",
    "To begin with you can generate RMSD plots for each of the eight trajectories.\n",
    "\n",
    "Run the code in the following cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "from utilities import *\n",
    "import numpy as np\n",
    "import mdtraj as mdt\n",
    "from mdplus import pca\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#=====================\n",
    "# You may play with the variables in this section later on in the tutorial:\n",
    "n_reps = 8\n",
    "#n_reps = 6\n",
    "atom_selection = 'all'\n",
    "#atom_selection = 'name CA'\n",
    "#atom_selection = 'name CA and resid 7 to 89'\n",
    "first_frame = 1\n",
    "#first_frame = 100\n",
    "last_frame = 2500\n",
    "#======================\n",
    "\n",
    "trajectory_files = ['data/rep{}/1rhw.md1.nc'.format(i + 1) for i in range(n_reps)]\n",
    "replicate_ids = ['data/rep{}'.format(i + 1) for i in range(n_reps)]\n",
    "topology_file = 'data/1rhw_prot.pdb'\n",
    "trajlist = []\n",
    "for tf in trajectory_files:\n",
    "    traj = mdt.load(tf, top=topology_file)\n",
    "    sel = traj.topology.select(atom_selection)\n",
    "    trajlist.append(mdt.Trajectory(traj.xyz[first_frame - 1:last_frame, sel], topology=traj.topology.subset(sel)))\n",
    "dynein_data = trajlist[0].join(trajlist[1:])\n",
    "dynein_data.superpose(dynein_data[0])\n",
    "\n",
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(121)\n",
    "plot_rmsd(dynein_data, replicate_ids)\n",
    "\n",
    "plt.subplot(122)\n",
    "plot_rmsf(dynein_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The RMSD plots for each replicate should appear broadly similar, rising quickly at the start but then more-or-less plateauing. The RMSF plot on the right shows which atoms are most mobile.\n",
    "\n",
    "These plots calculated the RMSDs over all atoms - try repeating the calculation, but just including C-alpha atoms, as you did for the neuraminidase example before. To do this, just uncomment (remove the '#' symbol from) the line in the cell above that reads `#atom_selection = 'name CA'` and comment out (add a '#' to the start of) the line above that currently reads `atom_selection = 'all'`, then press shift-enter again to re-run the calculation. \n",
    "\n",
    "#### Does much change? Do the curves plateau at a significantly lower RMSD?\n",
    "\n",
    "\n",
    "The RMSF plot should make you think that it might be better to leave out the first few (maybe 6?) amino acids as well. By adding/removing '#'s in the cell above, set the value of `atom_selection` to be `'resid 7 to 89'`, and recalculate.\n",
    "\n",
    "#### Does that improve things?\n",
    "\n",
    "\n",
    "We might guess that the simulations include a preliminary relaxation stage, when the RMSD rises quickly. What happens if you discard the first, say, 100 snapshots from each trajectory file? To do this uncomment the line above that reads `#first_frame = 100`, comment out the line above (that reads `first_frame  = 1`) and recalculate. \n",
    "\n",
    "#### Does that help?\n",
    "\n",
    "\n",
    "#### Experiment with leaving out more of the initial simulation data.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Part 2B: PCA Analysis\n",
    "\n",
    "Let's now see if PCA reveals anything more. The code below does just the same as for the neuraminidase example you studied earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "p = pca.PCA()\n",
    "scores = p.fit_transform(dynein_data.xyz)\n",
    "\n",
    "plt.figure(figsize=(15,15))\n",
    "for dataset in range(8):\n",
    "    plt.subplot(3, 3, dataset + 1)\n",
    "    plot_pca(scores, replicate_ids, highlight=dataset)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Each panel shows the path of a single replicate trajectory, over the background of the rest. What do you see about the behaviour of individual replicates? Are there \"outliers\"? You may be suspicious about replicate 8. If you have time, you could go back to the start of this section, and by changing the value of `n_reps`, repeat the analysis without it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Hopefully this part of the workshop has provided you with some tools that you can apply to your own simulation studies. You will have seen the limitations of conventional RMSD analysis for looking at equilibration and convergence, and how PCA can help. You will see how data analysis is frequently an interactive, iterative, process. It has probably also made it very obvious that getting a well-equilibrated, well-sampled dataset from MD simulation from which reliable and reproducible conclusions can be drawn is not a trivial issue!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}