docs/notebooks/md-ligand-protein.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ligand-protein MD\n",
    "\n",
    "This tutorial showcases how to use ProLIF to generate an interaction fingerprint for a ligand-protein complex from an MD simulation.\n",
    "\n",
    "ProLIF uses MDAnalysis to process MD simulations, as it supports many [file formats](https://userguide.mdanalysis.org/stable/formats/index.html) and is inter-operable with RDKit (which ProLIF uses under the hood to implement the different types of interactions). To learn more on how to use MDAnalysis, you can find their user guide [here](https://userguide.mdanalysis.org/stable/index.html).\n",
    "\n",
    "We will use 3 objects from the MDAnalysis library:\n",
    "\n",
    "- The `Universe` which bundles the atoms and bonds of your system with the coordinates in your trajectory.\n",
    "- The `AtomGroup` which is a collection of atoms that you can define by applying a [selection](https://userguide.mdanalysis.org/stable/selections.html) on the `Universe`.\n",
    "- The trajectory (or most often a subset of it) to know which frames to process.\n",
    "\n",
    ":::{important}\n",
    "For convenience, the topology and trajectory files for this tutorial are included with the ProLIF installation, and you can access the path to both files through `prolif.datafiles.TOP` and `prolif.datafiles.TRAJ` respectively. **Remember to switch these with the actual paths to your inputs outside of this tutorial.**\n",
    ":::\n",
    "\n",
    ":::{tip}\n",
    "At the top of the page you can find links to either download this notebook or run it in\n",
    "Google Colab. You can install the dependencies for the tutorials with the command:\n",
    "\n",
    "```shell\n",
    "pip install prolif[tutorials]\n",
    "```\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "Let's start by importing MDAnalysis and ProLIF to read our tutorial files, and create selections for the ligand and protein:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "import prolif as plf\n",
    "\n",
    "# load topology and trajectory\n",
    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
    "\n",
    "# create selections for the ligand and protein\n",
    "ligand_selection = u.select_atoms(\"resname LIG\")\n",
    "protein_selection = u.select_atoms(\"protein\")\n",
    "ligand_selection, protein_selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MDAnalysis should automatically recognize the file type that you're using from its extension. \n",
    "\n",
    ":::{note}\n",
    "Click [here](https://userguide.mdanalysis.org/stable/examples/quickstart.html) to learn more about loading files with MDAnalysis,\n",
    "and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case the protein is of reasonable size, but if you're working with a very large system, to save some time and memory it may be wise to restrict the protein selection to a sphere around the ligand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_selection = u.select_atoms(\n",
    "    \"protein and byres around 20.0 group ligand\", ligand=ligand_selection\n",
    ")\n",
    "protein_selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on your system, it may be of interest to also include water molecules in the protein selection. There are none in this tutorial example but something like this could be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_selection = u.select_atoms(\n",
    "    \"(protein or resname WAT) and byres around 20.0 group ligand\",\n",
    "    ligand=ligand_selection,\n",
    ")\n",
    "protein_selection"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, lets make sure that our ligand was correctly read by MDAnalysis.\n",
    "\n",
    ":::{important}\n",
    "This next step is crucial if you're loading a structure from a file that doesn't\n",
    "explicitely contain bond orders and formal charges, such as a PDB file or most MD\n",
    "trajectory files. MDAnalysis will infer those from the atoms connectivity, which\n",
    "**requires all atoms including hydrogens to be present in the input file**.\n",
    ":::\n",
    "\n",
    "Since ProLIF molecules are built on top of RDKit, we can use RDKit functions to display molecules. Let's have a quick look at our ligand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a molecule from the MDAnalysis selection\n",
    "ligand_mol = plf.Molecule.from_mda(ligand_selection)\n",
    "# display\n",
    "plf.display_residues(ligand_mol, size=(400, 200))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same for the residues in the protein (only showing the first 20 to keep the notebook short):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_mol = plf.Molecule.from_mda(protein_selection)\n",
    "# remove the `slice(20)` part to show all residues\n",
    "plf.display_residues(protein_mol, slice(20))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{admonition} Troubleshooting\n",
    "\n",
    "If one of the two molecules was not processed correctly, it might be because your input file does not contain bond information. If that's the case, please add this snippet right after creating your selections:\n",
    "\n",
    "```python\n",
    "ligand_selection.guess_bonds()\n",
    "protein_selection.guess_bonds()\n",
    "```\n",
    "\n",
    "In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of MDAnalysis molecules to RDKit. Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
    "\n",
    "```python\n",
    "protein_selection.guess_bonds(vdwradii={\"H\": 1.05, \"O\": 1.48})\n",
    "```\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fingerprint generation\n",
    "\n",
    "Everything looks good, we can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
    "You can list all interactions that are available with the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plf.Fingerprint.list_available()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair.\n",
    "\n",
    "If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. when there are two atoms in the ligand that make an HBond-donor interaction with residue X), use `plf.Fingerprint(count=True)`.\n",
    "This is also quite useful for visualization purposes as you can then display the atom pair that has the shortest distance which will look more accurate.\n",
    "This fingerprint type is however a bit slower to compute.\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use default interactions\n",
    "fp = plf.Fingerprint()\n",
    "# run on a slice of the trajectory frames: from begining to end with a step of 10\n",
    "fp.run(u.trajectory[::10], ligand_selection, protein_selection)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. You can modify the 6.0 Å cutoff by specifying `plf.Fingerprint(vicinity_cutoff=7.0)`, but this is only useful if you decide to change the distance parameters for an interaction class (see in the advanced section of the tutorials).\n",
    "\n",
    "Alternatively, you can pass a list of residues like so:\n",
    "\n",
    "```python\n",
    "fp.run(<other parameters>, residues=[\"TYR38.A\", \"ASP129.A\"])\n",
    "```\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp.to_pickle(\"fingerprint.pkl\")\n",
    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "Once the execution is done, you can access the results through `fp.ifp` which is a nested dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame_number = 0\n",
    "ligand_residue = \"LIG1.G\"\n",
    "protein_residue = \"ASP129.A\"\n",
    "\n",
    "fp.ifp[frame_number][(ligand_residue, protein_residue)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this contains all the details about the different interactions that were detected, it's not the easiest thing to digest.\n",
    "\n",
    "The best way to analyse our results is to export the interaction fingerprint to a Pandas DataFrame. You can read more about pandas in their\n",
    "[user_guide](https://pandas.pydata.org/docs/user_guide/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = fp.to_dataframe()\n",
    "# show only the 10 first frames\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some common pandas snippets to extract useful information from the fingerprint table.\n",
    "\n",
    ":::{important}\n",
    "Make sure to remove the `.head(5)` at the end of the commands to display the results for all the frames.\n",
    ":::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hide an interaction type (Hydrophobic)\n",
    "df.drop(\"Hydrophobic\", level=\"interaction\", axis=1).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show only one protein residue (ASP129.A)\n",
    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show only an interaction type (PiStacking)\n",
    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# percentage of the trajectory where each interaction is present\n",
    "(df.mean().sort_values(ascending=False).to_frame(name=\"%\").T * 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# same but we regroup all interaction types\n",
    "(\n",
    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
    "    .sum()\n",
    "    .T.astype(bool)\n",
    "    .mean()\n",
    "    .sort_values(ascending=False)\n",
    "    .to_frame(name=\"%\")\n",
    "    .T\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# percentage of the trajectory where PiStacking interactions are present, by residue\n",
    "(\n",
    "    df.xs(\"PiStacking\", level=\"interaction\", axis=1)\n",
    "    .mean()\n",
    "    .sort_values(ascending=False)\n",
    "    .to_frame(name=\"%\")\n",
    "    .T\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# percentage of the trajectory where interactions with SER212 occur, by interaction type\n",
    "(\n",
    "    df.xs(\"SER212.A\", level=\"protein\", axis=1)\n",
    "    .mean()\n",
    "    .sort_values(ascending=False)\n",
    "    .to_frame(name=\"%\")\n",
    "    .T\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# percentage of the trajectory where each interaction type is present\n",
    "(\n",
    "    df.T.groupby(level=\"interaction\")\n",
    "    .sum()\n",
    "    .T.astype(bool)\n",
    "    .mean()\n",
    "    .sort_values(ascending=False)\n",
    "    .to_frame(name=\"%\")\n",
    "    .T\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10 residues most frequently interacting with the ligand\n",
    "(\n",
    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
    "    .sum()\n",
    "    .T.astype(bool)\n",
    "    .mean()\n",
    "    .sort_values(ascending=False)\n",
    "    .head(10)\n",
    "    .to_frame(\"%\")\n",
    "    .T\n",
    "    * 100\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can compute a Tanimoto similarity between frames:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tanimoto similarity between the first frame and the rest\n",
    "from rdkit import DataStructs\n",
    "\n",
    "bitvectors = fp.to_bitvectors()\n",
    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bitvectors[0], bitvectors)\n",
    "tanimoto_sims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compare binding modes in your trajectory, it's possible to compute the entire similarity matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Tanimoto similarity matrix\n",
    "bitvectors = fp.to_bitvectors()\n",
    "similarity_matrix = []\n",
    "for bv in bitvectors:\n",
    "    similarity_matrix.append(DataStructs.BulkTanimotoSimilarity(bv, bitvectors))\n",
    "similarity_matrix = pd.DataFrame(similarity_matrix, index=df.index, columns=df.index)\n",
    "\n",
    "# display heatmap\n",
    "fig, ax = plt.subplots(figsize=(3, 3), dpi=200)\n",
    "colormap = sns.diverging_palette(\n",
    "    300, 145, s=90, l=80, sep=30, center=\"dark\", as_cmap=True\n",
    ")\n",
    "sns.heatmap(\n",
    "    similarity_matrix,\n",
    "    ax=ax,\n",
    "    square=True,\n",
    "    cmap=colormap,\n",
    "    vmin=0,\n",
    "    vmax=1,\n",
    "    center=0.5,\n",
    "    xticklabels=5,\n",
    "    yticklabels=5,\n",
    ")\n",
    "ax.invert_yaxis()\n",
    "plt.yticks(rotation=\"horizontal\")\n",
    "fig.patch.set_facecolor(\"white\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation\n",
    "\n",
    "There are a few different options builtin when it comes to visualisation.\n",
    "\n",
    "You can start by plotting the interactions over time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib ipympl\n",
    "\n",
    "fp.plot_barcode()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{tip}\n",
    "If you uncomment `%matplotlib ipympl` at the top of the above cell, you should be able to see an interactive version of the plot.\n",
    ":::\n",
    "\n",
    "You can also display the interactions in a 2D interactive diagram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = fp.plot_lignetwork(ligand_mol)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This diagram is interactive, you can:\n",
    "- zoom and pan,\n",
    "- move residues around,\n",
    "- click on the legend to display or hide types of residues or interactions,\n",
    "- hover an interaction line to display the distance.\n",
    "\n",
    ":::{note}\n",
    "After arranging the residues to your liking, you can save the plot as a PNG image with:\n",
    "```python\n",
    "view.save_png()\n",
    "```\n",
    "Note that this only works in notebooks and cannot be used in regular Python scripts.\n",
    ":::\n",
    "\n",
    "You can generate 2 types of diagram with this function, controlled by the `kind` argument:\n",
    "- `frame`: shows a single specific frame (specified with `frame`, corresponds to the frame number in the simulation).\n",
    "- `aggregate` (default): the interactions from all frames are grouped and displayed. An optional `threshold` parameter controls the minimum frequency required for an interaction to be displayed (default `0.3`, meaning that interactions occuring in less than 30% of frames will be hidden). The width of interactions is linked to the frequency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some more examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = fp.plot_lignetwork(ligand_mol, threshold=0.0)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up to now we've only been using the default fingerprint generator, but you can optionally enable the `count` parameter to enumerate all occurences of an interaction (the default fingerprint generator will stop at the first occurence), and then display all of them by specifying `display_all=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp_count = plf.Fingerprint(count=True)\n",
    "fp_count.run(u.trajectory[0:1], ligand_selection, protein_selection)\n",
    "view = fp_count.plot_lignetwork(ligand_mol, kind=\"frame\", frame=0, display_all=True)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also visualize this information in 3D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = 0\n",
    "# seek specific frame\n",
    "u.trajectory[frame]\n",
    "ligand_mol = plf.Molecule.from_mda(ligand_selection)\n",
    "protein_mol = plf.Molecule.from_mda(protein_selection)\n",
    "# display\n",
    "view = fp_count.plot_3d(ligand_mol, protein_mol, frame=frame, display_all=False)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the lignetwork plot, you can hover atoms and interactions to display more information.\n",
    "\n",
    "The advantage of using a count fingerprint in that case is that it will automatically select the interaction occurence with the shortest distance for a more intuitive visualization.\n",
    "\n",
    "Once you're satisfied with the orientation, you can export the view as a PNG image with:\n",
    "\n",
    "```python\n",
    "view.save_png()\n",
    "```\n",
    "Note that this only works in notebooks and cannot be used in regular Python scripts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also compare two different frames on the same view, and the protein residues that have\n",
    "different interactions in the other frame or are missing will be highlighted in magenta:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from prolif.plotting.complex3d import Complex3D\n",
    "\n",
    "frame = 0\n",
    "u.trajectory[frame]\n",
    "ligand_mol = plf.Molecule.from_mda(ligand_selection)\n",
    "protein_mol = plf.Molecule.from_mda(protein_selection)\n",
    "comp3D = Complex3D.from_fingerprint(fp, ligand_mol, protein_mol, frame=frame)\n",
    "\n",
    "frame = 120\n",
    "u.trajectory[frame]\n",
    "ligand_mol = plf.Molecule.from_mda(ligand_selection)\n",
    "protein_mol = plf.Molecule.from_mda(protein_selection)\n",
    "other_comp3D = Complex3D.from_fingerprint(fp, ligand_mol, protein_mol, frame=frame)\n",
    "\n",
    "view = comp3D.compare(other_comp3D)\n",
    "view"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
  },
  "kernelspec": {
   "display_name": "Python 3.8.5 ('prolif')",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}