analysis_md.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "afd82513",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Trajectory Analysis, Visualization and Manipulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8930ef64",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Loading trajectory data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ccae86c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "\n",
    "# Load data\n",
    "simulation = mda.Universe('adk_equilibrium\\\\adk4AKE.psf',\n",
    "                          'adk_equilibrium\\\\1ake_007-nowater-core-dt240ps.dcd', \n",
    "                          in_memory=True)\n",
    "\n",
    "# Inspect the simulation\n",
    "print(len(simulation.atoms), len(simulation.trajectory))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dd44113",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Atom selections\n",
    "CA_group = simulation.select_atoms(\"name CA\")\n",
    "print(len(CA_group.atoms))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b27bf94",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Atomgroups read directly from the frame being processed\n",
    "for frame in simulation.trajectory[:10]:\n",
    "    print(CA_group.atoms.positions[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98003408",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# We are now on the first frame\n",
    "simulation.trajectory[0]\n",
    "print(\"Frame %s\" % simulation.trajectory.frame)\n",
    "print(CA_group.atoms.positions[0])\n",
    "\n",
    "# This moves the simulation to frame 4\n",
    "simulation.trajectory[4]\n",
    "print(\"Frame %s\" % simulation.trajectory.frame)\n",
    "\n",
    "print(CA_group.atoms.positions[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3fb97a9",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Root Mean Square Displacements (RMSD)\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Root Mean Square Deviation (RMSD) evaluates the similarity between proteins or molecular structures over time. It quantifies the average displacement between corresponding atoms in two structures, often a simulated and a reference structure. First it aligns the structures to minimize the distances and then computes the RMSD.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccacec2c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename='RMSD.png', width=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a040305",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "RMSD gives an idea of broad conformational changes and is often used to observe the presence of structure rearrengements or displacements. It can also be used to have a first pick into the convergence of a simulation. RMSD profiles tend to increase at the beggining and then reach a plateau. \n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "With Gromacs it can be computed using the gmx rms command:\n",
    "</div>\n",
    "<br/>\n",
    "gmx rms -s file.tpr -f file.xtc -n file.ndx -o output.xvg    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46becc77",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from MDAnalysis.analysis import rms, align\n",
    "import pandas as pd\n",
    "\n",
    "# Align trajectory to itself\n",
    "aligner = align.AlignTraj(simulation, reference=simulation, select=\"backbone\", in_memory=True)\n",
    "\n",
    "# RMSD can be calculated from atom selection or from trajectories by giving selection atoms\n",
    "rmsd_ca = rms.RMSD(CA_group, reference=CA_group, ref_frame=0).run()\n",
    "rmsd_ca2 = rms.RMSD(simulation, reference=simulation, select=\"name CA\", ref_frame=0).run()\n",
    "# It's important to call the method .run() to obtain the results\n",
    "\n",
    "# check that every point matches\n",
    "False in (rmsd_ca.results.rmsd == rmsd_ca2.results.rmsd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7568cd4",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Load results into a dataframe (Table)\n",
    "rmsd_df = pd.DataFrame(rmsd_ca.results.rmsd, \n",
    "                       columns=['Frame', 'Time','CA Backbone'])\n",
    "\n",
    "rmsd_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a4dce18",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Add columns to the table\n",
    "rmsd_df[\"Time (ns)\"] = rmsd_df[\"Time\"]/1000\n",
    "rmsd_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8de793a",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import plotly.express as px\n",
    "\n",
    "# Plotting\n",
    "figure = px.line(rmsd_df, x=\"Time (ns)\", y=\"CA Backbone\")\n",
    "figure.update_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fafe4154",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute estimates from the data\n",
    "print(rmsd_df[\"CA Backbone\"][1:].mean())\n",
    "#print(rmsd_df[\"CA Backbone\"][rmsd_df[\"Time (ns)\"] > 20].mean())\n",
    "equilibrium_frame = int(min(rmsd_df[\"Frame\"][rmsd_df[\"Time (ns)\"] > 20]))\n",
    "print(equilibrium_frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f27b3732",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Root Mean Square Displacements (RMSD) pitfalls\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "RMSD is usually not very informative. Simulations contain a high number of atoms which means that, by doing the square mean, relevant changes may get trapped into the noise. One solution could be to select subsets of atoms or to weight differently certain atoms.\n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Another limitation is that RMSD suffers from degeneracy. Different conformation can have the same RMSD respect the reference frame. One can compute pairwise RMSD (expensive) or use other metrics.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a2f95e9",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div style=\"text-align: justify\"> \n",
    "For simulations that have not reached equilibrium special care should be used when analyzing RMSD profiles. Ideally one should discard all the non equilibirum frames.\n",
    "</div>\n",
    "</br>\n",
    "<div style=\"text-align: justify\"> \n",
    "RMSD differences that appear relatively early are usually due to the equilibration procedure. It is Important to use replicates and to careful inspect the simulations\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baa6ded8",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Root Mean Square Fluctuation (RMSF)\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Root Mean Square Fluctuation (RMSF) is a metric used to quantify the flexibility or mobility of atoms within a biomolecular system, such as a protein or a nucleic acid. It provides insight into how much individual atoms deviate from their average positions during the simulation.\n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "The first step is to align the simulation and compute an average structure. Then all the frames are aligned to the average positions and finally the RMSF can be computed.\n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "With Gromacs it can be computed with the gmx rmsf command.\n",
    "</div>\n",
    "<br/>\n",
    "gmx rmsf -s file.tpr -f file.xtc -n file.ndx -o output.xvg    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "545fe1a7",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Image(filename='RMSF.png', width=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6771d3cc",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# First aling and compute the average structure\n",
    "\n",
    "average = align.AverageStructure(simulation, reference=simulation, select=\"name CA\").run(start=equilibrium_frame)\n",
    "ref = average.results.universe\n",
    "\n",
    "aligner = align.AlignTraj(simulation, ref, select='name CA', in_memory=True).run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0c9d9a0",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Compute RMSF\n",
    "rmsf_ca = rms.RMSF(CA_group).run(start=equilibrium_frame)\n",
    "\n",
    "rmsf_df = pd.DataFrame(rmsf_ca.results.rmsf, \n",
    "                       columns=['Fluctuations'])\n",
    "\n",
    "# Format data\n",
    "rmsf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e896360e",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "rmsf_df[\"Atom number\"] = CA_group.indices\n",
    "rmsf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88baa0aa",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "figure = px.line(rmsf_df, x=\"Atom number\", y=\"Fluctuations\")\n",
    "figure.update_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cf4c7f1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Visualization\n",
    "# Add fluctuations on B-factor information\n",
    "# Since gromacs files don't have B-factors we have to initialize the attirbute\n",
    "simulation.add_TopologyAttr('tempfactors')\n",
    "\n",
    "# Color all the atoms of each residue based on the fluctuations of its CA\n",
    "protein = simulation.select_atoms('protein') # select protein atoms\n",
    "for residue, r_value in zip(protein.residues, rmsf_ca.results.rmsf):\n",
    "    residue.atoms.tempfactors = r_value\n",
    "\n",
    "# Save pdb with b-factors as colors\n",
    "simulation.atoms.write('rmsf_tempfactors.pdb')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bd1850d",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Nglview to have a short pic of the results\n",
    "import nglview as nv\n",
    "\n",
    "view = nv.show_mdanalysis(simulation.atoms)\n",
    "view.update_representation(color_scheme='bfactor')\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b7b8df3",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Root Mean Square Fluctuation (RMSF) pitfalls\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "RMSF can require long timescales to properly converge, specially on highly mobile parts whose energy landscape contains deep minimas.</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Since RMSF provides a single measure for the whole trajectory, parts of the protein that behave in a \"multimodal way\" (conformation with low fluctuations and a conformation with high fluctuations) may not be described properly. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6d7865f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Radius of gyration (R<u>g</u>)\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "The radius of gyration (R<u>g</u>) of a structure provides a metric of its compactness. It is a useful metric to detect possible folding and unfolding events.</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "To compute the radius of gyration first it is necessary to compute the center of mass of the selected atoms and the distance of all the atoms from the center of mass (r<u>i</u>). Then the R<u>g</u> can be computed as follows: \n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "With Gromacs it can be computed with the gmx gyrate command\n",
    "</div>\n",
    "<br/>\n",
    "gmx gyrate -s file.tpr -f file.xtc -n file.ndx -o output.xvg    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f580e17f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Image(filename='gyration.png', width=500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56a89290",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Data holder for the radius of gyration and frames\n",
    "gyration = []\n",
    "time = []\n",
    "\n",
    "for i, ts in enumerate(simulation.trajectory[equilibrium_frame:]):\n",
    "    time.append(simulation.trajectory.time/1000)\n",
    "    gyration.append(simulation.atoms.radius_of_gyration())\n",
    "\n",
    "gyration_df = pd.DataFrame({\"Time (ns)\": time, \"Gyration\": gyration})\n",
    "gyration_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1869cb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "figure = px.line(gyration_df, x=\"Time (ns)\", y=\"Gyration\")\n",
    "figure.update_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b956b95",
   "metadata": {},
   "source": [
    "## Principal component analysis or dimensionality reduction\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Principal Component Analysis (PCA) is a dimensionality reduction technique commonly used in various fields, including molecular dynamics (MD) simulations, to simplify complex data while retaining its most important features. In the context of MD simulations, PCA is used to analyze and visualize the high-dimensional trajectory data of a molecular system.\n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "PCA begins by calculating a covariance matrix from the centered data, revealing how different atom positions vary together.  Then an Eigenvalue decomposition of the covariance matrix is performed which yields eigenvectors (representing variance directions) and eigenvalues (indicating variance amounts). Then the principal components can be determined and sorted based on the ammount of varaince they explain. These define a new coordinate system capturing significant variability in the original data.\n",
    "By selecting a subset of principal components explaining most variance (usually 2 or 3), data dimensionality is reduced while preserving critical information.\n",
    "Finally, the original data is projected creating transformed values. This simplifies visualization and represents the MD trajectory in reduced dimensions.</div>\n",
    "<br/> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a573fe61",
   "metadata": {},
   "outputs": [],
   "source": [
    "from MDAnalysis.analysis import pca\n",
    "\n",
    "# Align the trajectory\n",
    "aligner = align.AlignTraj(simulation, simulation, select='name CA',\n",
    "                          in_memory=True).run(start=equilibrium_frame)\n",
    "\n",
    "# Compute PCA\n",
    "pc = pca.PCA(simulation, select='name CA',\n",
    "             align=True, mean=None,\n",
    "             n_components=None).run(start=equilibrium_frame)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34b5651b",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Inspect variance explained by the first component\n",
    "print(pc.results.variance[0])\n",
    "\n",
    "# Accumulative variance explained by the first 3 components\n",
    "pc.results.cumulated_variance[3]\n",
    "\n",
    "explained_variance_df = pd.DataFrame({\"Eigenvector\": [i+1 for i in range(10)],\n",
    "                                      \"Explained Variance\": pc.results.cumulated_variance[:10]})\n",
    "\n",
    "# plot explained variance of 10 first components\n",
    "figure = px.line(explained_variance_df, x=\"Eigenvector\", y=\"Explained Variance\")\n",
    "figure.update_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63bab93f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Transform selection to project them to the selected components\n",
    "transformed = pc.transform(CA_group, n_components=3)\n",
    "# Create Dataframe\n",
    "pca_df = pd.DataFrame(transformed[equilibrium_frame:],\n",
    "                  columns=['PC{}'.format(i+1) for i in range(3)])\n",
    "# add time colum\n",
    "pca_df[\"Time (ns)\"] = (pca_df.index + equilibrium_frame) * simulation.trajectory.dt/1000 \n",
    "\n",
    "pca_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caad28ff",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# 3d PCA plot\n",
    "fig = px.scatter_3d(pca_df, x=\"PC1\", y=\"PC2\", z=\"PC3\", color=\"Time (ns)\")\n",
    "fig.update_traces(**{\"marker.size\": 3})\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfcca83",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Principal component analysis or dimensionality reduction\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Although PCA are really good exploratory analysis one should treat them with care. Due to the sampling timescales limitations that current computational resources have; it is practically impossible to have fully converged simulations. This means that a PCA analysis will be able to capture those differences, even if in reality they are not meaningful. When using PCA analysis is crucial to use replicates and to carefully asses convergence.\n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Another common mistake occurs during the analysis itself. It is important, in order to compare PCA's, that the simulations have been projected on the same space. This can be achieved by concatenating the simulations first, running the analysis and then colour them based on the simulation that they provide from. </div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "Because of this, it is of crucial importance that the atoms of the simulations being concatenated match and that there is no difference in the initial setup and orientation between conditions. Ideally the different conditions should have started from simulations that are as similar as possible to avoid the PCA from picking differences originating from the setup. </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42d5b6ac",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from MDAnalysis.core.universe import Merge\n",
    "from MDAnalysis.coordinates.memory import MemoryReader\n",
    "import numpy as np\n",
    "\n",
    "# Example of how to create a concatenated simulation of points of interest\n",
    "\n",
    "# Create a universe object\n",
    "combined_simulation = Merge(simulation.select_atoms(\"name CA\"))\n",
    "\n",
    "# Extract coordinates from the simulations\n",
    "coordinates1 = np.array([CA_group.positions.copy() for frame in simulation.trajectory[equilibrium_frame:]])\n",
    "coordinates2 = np.array([CA_group.positions.copy() for frame in simulation.trajectory[equilibrium_frame:]])\n",
    "\n",
    "# Load data into the combined simulation\n",
    "combined_simulation.load_new(np.vstack([coordinates1, coordinates2]), format=MemoryReader)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f029fd06",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other analysis techniques\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "The techniques covered here encompass some of the standard analysis that one commonly uses when analysing MD simulations. However, there are plenty of other analysis available. Which analysis technique is more descriptive greatly depends on the particular study case. Finding a measurement that correctly describes and discriminates the process of interest can be challenging. \n",
    "</div>\n",
    "<br/>\n",
    "<div style=\"text-align: justify\"> \n",
    "The experimental setup used can also greatly influence the results one may obtain. In addition to that, it is also possible to design more complex MD simulations that open the door to more complex analysis, such as enhanced sampling techniques and free energy calculation.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "base",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}