{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. Analysing MD data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Important note\n",
    "**This part of the tutorial was originally developed with the assumption that the reader had already gone through the Molecular Dynamics (MD) OxCompBio workshop. It makes some assumptions of prior knowledge about molecular dynamics and the trajectories generated from such simulations.**\n",
    "\n",
    "**We therefore advise readers to first complete the MD tutorial and then return to this section afterwards.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As previously discussed, there are many things for which python can be used. One such example is the analysis of Molecular Dynamics (MD) simulations.\n",
    "\n",
    "There exist several MD analysis libraries, including (but not limited to):\n",
    "- [MDAnalysis](https://www.mdanalysis.org/)\n",
    "- [MDTraj](http://mdtraj.org/1.9.3/)\n",
    "- [PyTraj](https://amber-md.github.io/pytraj/latest/index.html)\n",
    "- [LOOS](http://grossfieldlab.github.io/loos/)\n",
    "\n",
    "In this section, we will look at how we can use the MDAnalysis library to process simulation structures and trajectories. We will also look at how we can interface with the [NGLView](https://amber-md.github.io/pytraj/latest/index.html) to visualise our simulations.\n",
    "\n",
    "This tutorial will aim to put together all you have learned so far into one final exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.1 The MDAnalysis Universe module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To work with structures and trajectories in MDAnalysis, you must first create an MDAnalysis `Universe`.\n",
    "\n",
    "This is an object that contains all the necessary information about your system, including atomic coordinates, atom types, box vectors (if available), etc...\n",
    "\n",
    "A `Universe` is a 'class' object, which stores all the information above your system, including;\n",
    "- Atomic coordinates (i.e. positions).\n",
    "- Atom names.\n",
    "- Atom types.\n",
    "- Bond information (if available).\n",
    "\n",
    ".. and so on. Having access to this information is very useful for analysis, such as if you are trying to measure distances or dihedrals. \n",
    "\n",
    "When we load a trajectory (as we will do later in the section), the `Universe` also stores information like the current frame number, etc. \n",
    "\n",
    "First we import the MDAnalysis `Universe` module, we will also import `numpy` for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis import Universe\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You then load the structure file into a `Universe` instance named `alanine` with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alanine = Universe('datafiles/ALA.pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note 1**: We could have also used `import MDAnalysis`, in which case we would have to use `MDAnalysis.Universe` in place of `Universe`. As previously discussed, the `from ... import ...` syntax allows us to import only a sub-module from a library, rather than the full thing, which can be convenient in various circumstances.\n",
    "\n",
    "**Note 2**: On loading the PDB file you will get a warning that element information is absent or missing. This is normal and occurs because the input PDB file we use does not have element information for all the atoms (which is expected in a normal PDB file). Therefore MDAnalysis warns us that there are no elements present for us to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check the different properties of your universe using MDAnalysis. The general syntax for this is `universename.attributename`; some properties (attributes) are themselves classes (e.g. `alanine.atoms` groups all the information about our atoms) and will have their own attributes. For example, if you want to check the names of the atoms of your structure, you can type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alanine.atoms.names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The atom \"name\" is essentially a type of tag that defines the atoms in your system. Quite often this will be used to label similar atoms, for example in the PDB convention alpha carbons are always named \"CA\". When analysing MD simulations, the naming scheme will vary depending on the force field used (see MD tutorial)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### How would you find the values of the Phi and Psi dihedral angles of alanine dipeptide, in degrees or radians?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: the Phi dihedral is calculated from atoms CLP, NL, CA, and CRP, and the Psi dihedral is calculated using atoms NR, CRP, CA, and NL.\n",
    "\n",
    "For this we can make use of the built-in functions associated with the Universe class (or 'methods'). We can use `select_atoms(selection)` to isolate a group of atoms that match `selection`; the syntax of the selections strings is generally the same as for VMD (which you may have use in the MD tutorial). \n",
    "\n",
    "`alanine.select_atoms(selection)` will give us an 'atom group' class instance. We can then use this instance's method `dihedral` to turn four atoms into a 'dihedral' instance, and finally using `value()` on this instance will return the (current) value of that dihedral angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(alanine))\n",
    "print(type(alanine.select_atoms('name CLP NL CA CRP')))\n",
    "print(type(alanine.select_atoms('name CLP NL CA CRP').dihedral))\n",
    "\n",
    "phi = alanine.select_atoms('name CLP NL CA CRP').dihedral.value()\n",
    "psi = alanine.select_atoms('name NR CRP CA NL').dihedral.value()\n",
    "print(f'The Phi dihedral is {phi:.2f} in degrees and {numpy.deg2rad(phi):.2f} in radians')\n",
    "print(f'The Phi dihedral is {psi:.2f} in degrees and {numpy.deg2rad(psi):.2f} in radians')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### How would you plot the change in Phi and Psi over a trajectory?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to do this, let's break the process down into small steps.\n",
    "\n",
    "You will first need to load the trajectory. The format to do this is:\n",
    "\n",
    "universe_name = Universe(PDBfile, TRJfile)\n",
    "\n",
    "In this case we have provided you a trajectory in the `datafiles` folder, the PDB file is called 'ALA.pdb' and the TRJ file is called 'ALA.xtc'. So we can call them as `./datafiles/ALA.pdb` and `./datafiles/ALA.xtc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 12.1.1: Create a Universe called 'ala_trajectory' using the PDB file 'ALA.pdb' \n",
    "# and the trajectory file 'ALA.xtc'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have a Universe with your trajectory, we can access another feature of the Universe: trajectory data.\n",
    "\n",
    "First we should look at what the current trajectory frame number is.\n",
    "\n",
    "The frame number is stored under `universename.trajectory.frame` and the total number of frames is found under `universename.trajectory.n_frames`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 12.1.2: Print the current frame number and the total number of frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `Universe.trajectory` acts kind of like a list storing each frame in a simulation and various information about it, so we can loop through it the same way we did for lists in previous sections. As we go through each frame, the coordinates of each atom, stored in Universe.atoms, will be updated.\n",
    "\n",
    "We can therefore quite a quick `for` loop that prints the frame number for each frame in the trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# When using MDanalysis, we usually refer to each frame as a 'timestep' or 'ts'\n",
    "# Here at each loop instance, we allocate the current frame to 'ts'\n",
    "for ts in ala_trajectory.trajectory:\n",
    "    # the comma at the end here will stop it printing on a new line every time,\n",
    "    # so this doesn't take up too much space\n",
    "    print(ts.frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, use the same logic as Exercise 11.1 to print the phi and psi dihedral angles at each frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 12.1.3: Print the frame number and at least one of the dihedrals.\n",
    "# First create the phi and psi selections\n",
    "\n",
    "# We then loop over the the trajectory\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now that you can access both the frame number and the phi and psi dihedrals, it's time to plot them.\n",
    "\n",
    "First, we will need to import a plotting library. In this case we import pyplot from matplotlib.\n",
    "\n",
    "We then create three empty lists (`frames`, `all_phi`, and `all_psi`) to hold our trajectory data.\n",
    "\n",
    "Then we fill each of those lists with the relevant values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We declare matplotlib inline to make sure it plots properly\n",
    "%matplotlib inline\n",
    "# We need pyplot from matplotlib to plot the dihedrals\n",
    "from matplotlib import pyplot\n",
    "\n",
    "# Create three empty lists\n",
    "frames = []\n",
    "all_phi = []\n",
    "all_psi = []\n",
    "\n",
    "# iterate through the trajectory\n",
    "for frame in ala_trajectory.trajectory:\n",
    "    # calculate phi and psi\n",
    "    phi = ala_trajectory.select_atoms('name CLP or name NL or name CA or name CRP').dihedral.value()\n",
    "    psi = ala_trajectory.select_atoms('name NR or name CRP or name CA or name NL').dihedral.value()\n",
    "    # append frame number, phi, and psi to the lists\n",
    "    frames.append(frame.frame)\n",
    "    all_phi.append(phi)\n",
    "    all_psi.append(psi)\n",
    "\n",
    "# We then plot our data\n",
    "pyplot.plot(frames, all_psi)\n",
    "pyplot.plot(frames, all_phi)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, you can use the matplotlib.pyplot module to plot `x` and `y` values in the following manner:\n",
    "\n",
    "```python\n",
    "pyplot.plot(x_value, y_value)\n",
    "```\n",
    "\n",
    "We then use `pyplot.show()` to make the plot appear.\n",
    "\n",
    "**Note:** You can adjust your matplotlib pyplot plot using various different options such as; `pyplot.title()`, `pyplot.xlim()`, `pyplot.ylim()`, `pyplot.xlabel()`, and `pyplot.ylabel()`.\n",
    "\n",
    "For example, `pyplot.title()` will add a title.\n",
    "\n",
    "Check the `help()` of these functions for more details!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 12.2 Analysing a protein simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the MD tutorial you will do some molecular dynamics simulations of a HIV-1 protease protein using gromacs.\n",
    "Let's now look at how you could use MDAnalysis, and an associated python visulisation library NGLview, to analyse an MD trajectory of this system.\n",
    "\n",
    "_Note: depending on whether you have already done the MD tutorial, you may find that you have already done the analysis shown here. You may still find this section useful as we provide a few extra details which explain how the MD trajectory analysis works._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualising the trajectory using nglview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nglview library is a python widget for visualising simulation trajectories, achieving a similar task to the VMD program that you may also use in the MD tutorial. One of the interesting advantages of nglview is that it interfaces directly with analysis packages such as MDAnalysis and runs within jupyter notebooks.\n",
    "\n",
    "Let's see how we can use nglview to visualise an MDAnalysis universe object.\n",
    "\n",
    "First we need to create a universe (let's call it `protein`) from the simulation output files \"pre_md.pdb\" and \"md_cent.xtc\" (which are present in the `datafiles` directory, hence will be passed to Universe as `./datafiles/pre_md.pdb` and `./datafiles/md_cent.xtc`).\n",
    "\n",
    "_Note 1: We have pre-aligned the trajectory to the first frame for you so as to remove any motions related to translation._\n",
    "\n",
    "_Note 2: When loading the trajectory you will get a warning that MDAnalysis needs to reload offsets, this is normal in this scenario._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise 12.2.1: Let's load a universe named protein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's load nglview and use it's `show_mdanalysis` function to load the MDAnalysis universe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview\n",
    "protein_view = nglview.show_mdanalysis(protein)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default this pre-sets the nglview to show the protein in the cartoon representation. Let's add a few options to colour the protein by secondary structure, show water oxygens and change the background colour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's update the cartoon representation to colour the protein by secondary structure\n",
    "protein_view.update_cartoon(color='sstruc')\n",
    "\n",
    "# We then add a transparent hyperball representation of the water oxygens \n",
    "#(play with the opacity value, see what you get)\n",
    "protein_view.add_hyperball('SOL and not hydrogen', opacity=0.4)\n",
    "\n",
    "# Let's change the display a little bit\n",
    "protein_view.parameters = dict(camera_type='orthographic', clip_dist=0)\n",
    "\n",
    "# Set the background colour to black\n",
    "protein_view.background = 'black'\n",
    "\n",
    "# Call protein_view to visualise the trajectory\n",
    "protein_view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The nglview output can be controlled in the following way:\n",
    "\n",
    "- play / pause button: play the trajectory \n",
    "- double click window: enter or exit full screen mode \n",
    "- left mouse button: rotate system \n",
    "- middle mouse wheel: zoom in/out \n",
    "- right mouse button: translate system \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can be seen from the trajectory, the HIV-1 protease structure does indeed move, but by how much? In the next section we will see how we can use MDAnalysis to quantify backbone fluctuations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the root-mean-square deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to gain a quantitative description of how the HIV-1 protease moves in our simulation we can calculate the root-mean-square deviation (RMSD) of the protein backbone.\n",
    "\n",
    "The RMSD gives us an idea of how 'stable' our protein is when compared to our starting, static, structure. The lower the RMSD is the, more stable we can say our protein is. \n",
    "\n",
    "The RMSD as a function of time, $\\rho (t)$, can be defined by the following equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\\\\n",
    "\\rho (t) = \\sqrt{\\frac{1}{N}\\sum^N_{i=1}w_i\\big(\\mathbf{x}_i(t) - \\mathbf{x}^{\\text{ref}}_i\\big)^2}\n",
    "\\end{equation}\n",
    "\n",
    "Luckily MDAnalysis has its own built-in function to calcualte this, we can import it like we did before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.analysis.rms import RMSD as rmsd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to calculate the RMSD for every frame in our trajectory we will need:\n",
    "\n",
    "- A reference structure\n",
    "- A universe object\n",
    "- A selection of atoms\n",
    "\n",
    "In our case the reference structure will be the HIV-1 protease structure in the first frame.\n",
    "\n",
    "Our universe object will be the 'protein' object we defined above.\n",
    "\n",
    "For our selection we will use the backbone atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref = Universe('datafiles/pre_md.pdb', 'datafiles/md_cent.xtc')\n",
    "\n",
    "# Set the ref trajectory to the first frame\n",
    "ref.trajectory[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the way that GROMACS post processes the trajectory file we need to edit it slightly before running our RMSD.\n",
    "\n",
    "This is done by aligning all frames to the reference structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.analysis import align\n",
    "\n",
    "# Create the MD simulation universe\n",
    "protein = Universe('datafiles/pre_md.pdb', 'datafiles/md_cent.xtc')\n",
    "# Call the MDAnalysis align function to align the MD simulation unvierse to the reference (first frame) universe\n",
    "align_strucs = align.AlignTraj(protein, ref, select=\"backbone\", weights=\"mass\", in_memory=True, verbose=True)\n",
    "\n",
    "R = align_strucs.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will have noticed that running this function stores it in the variable 'R', we can now access the RMSD values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmsd_data = R.rmsd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'd like to visualise how the RMSD changes over time and this can be done in the same way you did in Excercise 11.1.5.\n",
    "\n",
    "Take a look at the 'rmsd_data' variable (it's a numpy array) and try plotting it below.\n",
    "\n",
    "You will need to access 'rmsd_data' (a numpy array) in order to plot both the time and the RMSD as a line plot.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Excercise 12.2.2: Plot the RMSD data for the HIV-1 protease system. \n",
    "\n",
    "# Make sure to add appropriate axis titles.\n",
    "# What does the RMSD tell you about the protein?\n",
    "# [If you have time] What happens when you calculate the RMSD using more atoms (i.e. not just the backbone)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the root-mean-square fluctuation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To look at how each residue flucuates over it's average postion we can use the closely related measurement of root-mean-square fluctuation (RMSF).\n",
    "\n",
    "The RMSF for an atom, $\\rho_i$, is given by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho_i = \\sqrt{\\sum^N_{i=1} \\big\\langle(\\mathbf{x}_i - \\langle \\mathbf{x}_i \\rangle )^2 \\big\\rangle }\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.analysis.rms import RMSF as rmsf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset the trajectory to the first frame\n",
    "protein.trajectory[0]\n",
    "\n",
    "# We will need to select the alpha Carbons only\n",
    "calphas = protein.select_atoms(\"name CA\")\n",
    "rmsf_calc = rmsf(calphas, verbose=True).run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Excercise 12.2.3: Plot the RMSF data for the HIV-1 protease system. \n",
    "    # Tip, in order to plot the resids you will need to access them through the rmsf_calc object\n",
    "\n",
    "# Make sure to add appropriate axis titles.\n",
    "# What parts of the protein have a high RMSF, can you locate these on the protein structure?\n",
    "# [If you have time] What happens when you calculate the RMSF using more atoms (i.e. not just the backbone)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}