PCA_analysis_of_MD_simulations-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Part 1: Comparison of trajectories of wild-type neuraminidase and of the I233R/H275Y (IRHY) double mutant.\n",
    "The data for this part of the workshop comes from:\n",
    "[Long time scale GPU dynamics reveal the mechanism of drug resistance of the dual mutant I223R/H275Y neuraminidase from H1N1-2009 influenza virus.](https://www.ncbi.nlm.nih.gov/pubmed/22574858)\n",
    "Woods CJ, Malaisree M, Pattarapongdilok N, Sompornpisut P, Hannongbua S, Mulholland AJ.\n",
    "Biochemistry. 2012 May 29;51(21):4364-75. doi: 10.1021/bi300561n.\n",
    "\n",
    "You have been provided with two trajectory files (AMBER binpos format): `wt_ca.binpos` and `irhy_ca.binpos`. The 2400 frames in each trajectory file are spaced every 200ps from 20ns to 500ns. For computational simplicity, the files have been stripped down to just the coordinates of the C-alpha atoms.\n",
    "\n",
    "Let's begin by loading the two trajectories into MDTraj trajectory objects, joining them together into a single trajectory, then visualising the dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mdtraj as mdt\n",
    "import nglview as nv\n",
    "\n",
    "# Load the data for the wt and irhy simulations:\n",
    "t_wt = mdt.load('data/wt_ca.binpos', top='data/wt_ca.pdb')\n",
    "t_irhy = mdt.load('data/irhy_ca.binpos', top='data/irhy_ca.pdb')\n",
    "# Combine the two sets of trajectory data into one trajectory:\n",
    "trajdata = t_wt.join(t_irhy, check_topology=False)\n",
    "view = nv.show_mdtraj(trajdata)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "Notice that half way through the movie the protein jumps - this marks the transition from viewing the dynamics of the wild-type protein to viewing the dynamics of the mutant - the two simulations happen to have been set up in different parts of coordinate space.\n",
    "\n",
    "---\n",
    "\n",
    "### Part 1A: RMSD Analysis ###\n",
    "To begin with, we will plot the RMSD of each snaphot in each file relative to the first.\n",
    "\n",
    "Run the following cell. The python code here loads the two trajectories, defines a function to calculate and then plot RMSDs, and then applies the function to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams.update({'font.size': 15}) #This sets a better default label size for plots\n",
    "\n",
    "# Create labels for the two datasets:\n",
    "datanames = ['wildtype', 'irhy mutant']\n",
    "\n",
    "# define the plotting function:\n",
    "def plot_rmsd(traj, datanames):\n",
    "    \"\"\"\n",
    "    This function takes a MDTraj trajectory and a list of data names and produces an rmsd plot.\n",
    "    \n",
    "    \"\"\"\n",
    "    traj.superpose(traj[0]) # least squares fits each snapshot to the first.\n",
    "    frames_per_set = len(traj) // len(datanames) # we assume each trajectory file is the same length.\n",
    "    for i in range(len(datanames)):\n",
    "        # The next two lines do the rmsd calculation:\n",
    "        diff = traj.xyz[i * frames_per_set : (i + 1) * frames_per_set] - traj.xyz[0]\n",
    "        rmsd = np.sqrt((diff * diff).sum(axis=2).mean(axis=1))\n",
    "        plt.plot(rmsd, label=datanames[i]) # plot the line for this dataset on the graph.\n",
    "    plt.xlabel('Frame number')\n",
    "    plt.ylabel('RMSD (nm.)')\n",
    "    plt.legend(loc='lower right')\n",
    "\n",
    "# now use the plotting function:    \n",
    "plot_rmsd(trajdata, datanames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### What do you conclude from these graphs? Do the simulations look well-equilibrated?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's calculate the root-mean-square fluctuations (RMSFs) of the atoms in each trajectory. This will tell us which parts of the structure are most mobile, and which are most rigid. The most flexible parts of the system are likely to be those most difficult to equilibrate and to sample well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def plot_rmsf(traj):\n",
    "    \"\"\"\n",
    "    Plots the root mean square fluctuations of the atoms in a MDTraj trajectory.\n",
    "    \n",
    "    \"\"\"\n",
    "    diff = traj.xyz - traj.xyz.mean(axis=0)\n",
    "    rmsf = np.sqrt((diff * diff).sum(axis=2).mean(axis=0))\n",
    "    plt.xlabel('atom number')\n",
    "    plt.ylabel('RMSF (nm.)')\n",
    "    plt.plot(rmsf)\n",
    "    \n",
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(121)\n",
    "frames_per_set = len(trajdata) // len(datanames)\n",
    "plot_rmsf(trajdata[:frames_per_set]) # the first half of the cofasu has the wt data.\n",
    "plt.title(datanames[0])\n",
    "plt.subplot(122)\n",
    "plot_rmsf(trajdata[frames_per_set:]) # the second half of the cofasu has the irhy data.\n",
    "plt.title(datanames[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The left-hand plot is for the wild-type trajectory, the right-hand one is for the irhy mutant. It's clear that in both cases the C-terminus of the protein is exceptionally dynamic, and probably the least well-sampled. Since this region is well away from the oseltamivir binding site, let's repeat the RMSD analysis, leaving residues 370 onwards out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Create an MDTraj trajectory for just a selection of the atoms in the system:\n",
    "selection = trajdata.topology.select('resid 1 to 370')\n",
    "seldata = mdt.Trajectory(trajdata.xyz[:, selection], trajdata.topology.subset(selection))\n",
    "\n",
    "plot_rmsd(seldata, datanames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Are the results as you expect? If you discount the particularly flexible C-terminus of the protein, do these simulations otherwise look well equilibrated?\n",
    "\n",
    "**_EXERCISE_**:\n",
    "\n",
    "**The RMSF plot shows that there are also quite dynamic regions towards the N-terminus of the protein. What happens if you leave these out of the RMSD analysis too? Try editing the cell above so the selection is `resid 100 to 300` and repeat the analysis. Experiment with different selections.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Part 1B: PCA Analysis\n",
    "Now let's move from analysis by RMSD to analysis by PCA.\n",
    "\n",
    "The code below does a PCA analysis on both trajectory sets combined, and then plots the projection of each trajectory onto the common PC1/PC2 subspace:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "from mdplus import pca\n",
    "# First define a plotting function:\n",
    "def plot_pca(scores, datanames, highlight=None):\n",
    "    \"\"\"\n",
    "    Plots the projection of each trajectory in the cofasu in the PC1/PC2 subspace.\n",
    "    \n",
    "    If highlight is a number, this dataset is plotted in red against all others in grey.\n",
    "    \n",
    "    \"\"\"\n",
    "    p1 = scores[:,0] # the projection of each snapshot along the first principal component\n",
    "    p2 = scores[:,1] # the projection along the second.\n",
    "\n",
    "    frames_per_rep = len(p1) // len(datanames) # number of frames (snapshots) in each dataset - assumed equal length\n",
    "    for i in range(len(datanames)):\n",
    "        start = i * frames_per_rep\n",
    "        end = (i + 1) * frames_per_rep\n",
    "        if highlight is None: # each dataset is plotted with a different colour\n",
    "            plt.plot(p1[start:end], p2[start:end], label=datanames[i]) \n",
    "            plt.text(p1[start], p2[start], 'start')\n",
    "            plt.text(p1[end-1], p2[end-1], 'end')\n",
    "        else:\n",
    "            if i != highlight:\n",
    "                plt.plot(p1[start:end], p2[start:end], color='grey')\n",
    "    if highlight is not None:\n",
    "        start = highlight * frames_per_rep\n",
    "        end = (highlight + 1) * frames_per_rep\n",
    "        plt.plot(p1[start:end], p2[start:end], color='red', label=datanames[highlight])\n",
    "        plt.text(p1[start], p2[start], 'start')\n",
    "        plt.text(p1[end-1], p2[end-1], 'end')\n",
    "\n",
    "    plt.xlabel('PC1')\n",
    "    plt.ylabel('PC2')\n",
    "    plt.legend()\n",
    "\n",
    "# Now use it:\n",
    "selection = trajdata.topology.select('resid 1 to 370')\n",
    "seldata = trajdata.atom_slice(selection)\n",
    "p = pca.PCA()\n",
    "scores = p.fit_transform(seldata.xyz)\n",
    "plot_pca(scores, datanames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The left-hand plot shows the track in PC1/PC2 space of each trajectory. Both start in the bottom right hand corner but then move off in separate directions.\n",
    "\n",
    "\n",
    "**_EXERCISE_**:\n",
    "\n",
    "**As you did for the RMSD analysis, experiment with leaving out the most dynamic (N- and C-terminal) regions of the proteins from the PCA analysis, by editing the atom selections.**\n",
    "\n",
    "#### Compare the RMSD and PCA plots. Consider how clearly they identify the similarities and differences between the two trajectories, and illustrate sampling and convergence.\n",
    "\n",
    "[Next: Part 2: Analysis of replicate simulations of dynein light chain](PCA_analysis_of_MD_simulations-2.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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
}