docs/notebooks/quickstart.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quickstart\n",
    "\n",
    "This is a very short guide on how to use ProLIF to generate an interaction fingerprint for a ligand-protein complex from an MD simulation.\n",
    "\n",
    "Let's start by importing MDAnalysis and ProLIF to read our input files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "import prolif as plf\n",
    "# load trajectory\n",
    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
    "# create selections for the ligand and protein\n",
    "lig = u.atoms.select_atoms(\"resname LIG\")\n",
    "prot = u.atoms.select_atoms(\"protein\")\n",
    "lig, prot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MDAnalysis should automatically recognize the file type that you're using from its extension. Click [here](https://userguide.mdanalysis.org/stable/quickstart.html) to learn more about loading files with MDAnalysis, and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.\n",
    "\n",
    "Next, lets make sure that our ligand was correctly read by MDAnalysis.\n",
    "\n",
    "⚠ This next step is crucial if you're loading a structure from a file that doesn't explicitely contain bond orders and formal charges. MDAnalysis will infer those from the atoms connectivity, which requires all atoms including hydrogens to be present in the input file.\n",
    "\n",
    "ProLIF molecules are built on top of RDKit and are compatible with its drawing code. Let's have a quick look at our ligand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "from rdkit.Chem import Draw\n",
    "# create a molecule from the MDAnalysis selection\n",
    "lmol = plf.Molecule.from_mda(lig)\n",
    "# cleanup before drawing\n",
    "mol = Chem.RemoveHs(lmol)\n",
    "mol.RemoveAllConformers()\n",
    "Draw.MolToImage(mol, size=(400,200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do the same for the residues in the protein (I'll only show the first 20 to keep the notebook short):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmol = plf.Molecule.from_mda(prot)\n",
    "frags = []\n",
    "# to show all residues, simply use `for res in pmol:`\n",
    "for i in range(20):\n",
    "    res = pmol[i]\n",
    "    mol = Chem.RemoveHs(res)\n",
    "    mol.RemoveAllConformers()\n",
    "    frags.append(mol)\n",
    "Draw.MolsToGridImage(frags,\n",
    "                     legends=[str(res.resid) for res in pmol], \n",
    "                     subImgSize=(200, 140),\n",
    "                     molsPerRow=4,\n",
    "                     maxMols=prot.n_residues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Everything looks good, we can now compute a fingerprint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use default interactions\n",
    "fp = plf.Fingerprint()\n",
    "# run on a slice of frames from begining to end with a step of 10\n",
    "fp.run(u.trajectory[::10], lig, prot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. Alternatively, you can pass a list of residues like so:\n",
    "\n",
    "```python\n",
    "fp.run(..., residues=[\"TYR38.A\", \"ASP129.A\"])\n",
    "```\n",
    "Or simply use `fp.run(..., residues=\"all\")` to use all residues in the `prot` selection.\n",
    "\n",
    "The best way to access our results is to export our interaction fingerprints to a Pandas DataFrame. By default, the resulting DataFrame only keeps track of residues and interaction types that were seen in at least one of the frames in your trajectory. If needed, you can access the full results with `fp.to_dataframe(drop_empty=False)`."
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop the ligand residue column since there's only a single ligand residue\n",
    "df = df.droplevel(\"ligand\", axis=1)\n",
    "df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show all pi-stacking interactions\n",
    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# show all interactions with a specific protein residue\n",
    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).head(5)\n",
    "# or more simply\n",
    "df[\"ASP129.A\"].head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a simple example to plot the interactions over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "\n",
    "# reorganize data\n",
    "data = df.reset_index()\n",
    "data = pd.melt(data, id_vars=[\"Frame\"], var_name=[\"residue\",\"interaction\"])\n",
    "data = data[data[\"value\"] != False]\n",
    "data.reset_index(inplace=True, drop=True)\n",
    "\n",
    "# plot\n",
    "sns.set_theme(font_scale=.8, style=\"white\", context=\"talk\")\n",
    "g = sns.catplot(\n",
    "    data=data, x=\"interaction\", y=\"Frame\", hue=\"interaction\", col=\"residue\",\n",
    "    hue_order=[\"Hydrophobic\", \"HBDonor\", \"HBAcceptor\", \"PiStacking\", \"CationPi\", \"Cationic\"],\n",
    "    height=3, aspect=0.2, jitter=0, sharex=False, marker=\"_\", s=8, linewidth=3.5,\n",
    ")\n",
    "g.set_titles(\"{col_name}\")\n",
    "g.set(xticks=[], ylim=(-.5, data.Frame.max()+1))\n",
    "g.set_xticklabels([])\n",
    "g.set_xlabels(\"\")\n",
    "g.fig.subplots_adjust(wspace=0)\n",
    "g.add_legend()\n",
    "g.despine(bottom=True)\n",
    "for ax in g.axes.flat:\n",
    "    ax.invert_yaxis()\n",
    "    ax.set_title(ax.get_title(), pad=15, rotation=60, ha=\"center\", va=\"baseline\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the occurence of each interaction on the trajectory\n",
    "occ = df.mean()\n",
    "# restrict to the frequent ones\n",
    "occ.loc[occ > 0.3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# regroup all interactions together and do the same\n",
    "g = (df.groupby(level=[\"protein\"], axis=1)\n",
    "       .sum()\n",
    "       .astype(bool)\n",
    "       .mean())\n",
    "g.loc[g > 0.3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also compute a Tanimoto similarity between each frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import DataStructs\n",
    "bvs = fp.to_bitvectors()\n",
    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bvs[0], bvs)\n",
    "tanimoto_sims"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.5 64-bit ('prolif': conda)",
   "language": "python",
   "name": "python38564bitprolifconda6f88be71e4394ccb9ad64335b88224ec"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}