docking_single_ligand.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: center;\">\n",
    "        <img src=\"images/molssi_logo_header.png\" alt=\"MolSSI log\" style=\"height: 250px;\">\n",
    "</div>\n",
    "\n",
    "# Docking with AutoDock Vina\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"> \n",
    "<h2>Overview</h2>\n",
    "<strong>Questions</strong>\n",
    "\n",
    "* How can I dock a ligand using AutoDock Vina?\n",
    "  \n",
    "* How can I visualize docking structures and ligand interactions?\n",
    "\n",
    "<strong>Learning Objectives</strong>\n",
    "\n",
    "* Learn how to define a ligand box for docking.\n",
    "  \n",
    "* Execute ligand docking using AutoDock Vina.\n",
    "\n",
    "* Analyze and interpret docking results using ProLIF.\n",
    "</div>\n",
    "\n",
    "In this notebook, we are going to dock the ligands we generated previously notebook with [PDB entry 2zq2](https://www.rcsb.org/structure/2zq2), a trypsin structure from the cow *Bos taurus* that we retrieved and processed in the previous notebook.\n",
    "\n",
    "Molecular docking simulates the interaction between a small molecule (ligand) and a protein, predicting how these molecules fit together at the molecular level.\n",
    "It is commonly used in fields such as drug development to design molecules to bind to enzymes or proteins.\n",
    "\n",
    "During a docking process, AutoDock Vina evaluates numerous potential orientations and conformations of the ligand within the protein’s binding pocket. \n",
    "Each pose is assigned a score based on how well it fits, using various scoring functions. \n",
    "These functions can be physics-based, empirical, or even be based on machine learning algorithms to predict binding affinity. \n",
    "For a recent review of docking scoring functions, you can see [this publication](https://link.springer.com/article/10.1007/S12539-019-00327-W).\n",
    "\n",
    "In this notebook, we'll be performing docking with the AutoDock Vina Python API using the \"vina\" scoring function.\n",
    "\n",
    "\n",
    "### Libraries for the IQB workshop\n",
    "\n",
    "| Library         | abbreviation | Purpose |\n",
    "|:-------------|:---------:|:------------|\n",
    "| os           | N/A      | operating system functions - handling file paths and directories. |\n",
    "| MDAnalysis     | mda | molecular dynamics library - used for reading/writing files and selecting atoms |\n",
    "| vina | vina | AutoDock Vina software for Python and Jupyter notebooks |\n",
    "| prolif | plf | ProLIF (Protein-Ligand Interaction Fingerprints) - generates interaction fingerprints for complexes made of ligands, protein, DNA or RNA molecules\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing for Docking: Defining a Ligand Box\n",
    "\n",
    "When we dock our ligands to our protein, we will want to define the binding pocket. \n",
    "This is where the software will attempt to bind our ligand.\n",
    "In some cases, people might use software for finding pockets in proteins to define binding sites.\n",
    "However, as we saw in the previous notebook, the structure we retrieved from the PDB already had a ligand bound.\n",
    "To define our binding box, we will take the position of the bound ligand from our original structure.\n",
    "\n",
    "Luckily, MDAnalysis has tools that we can use to measure our molecule and define a binding box. \n",
    "The approach we will take in this notebook is to find the `center_of_geometry` of our ligand to define the center of our binding pocket.\n",
    "Then we will consider the space around the ligand to be our box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the center of the ligand\n",
    "import MDAnalysis as mda\n",
    "\n",
    "original_structure = mda.Universe(\"protein_structures/2zq2.pdb\")\n",
    "ligand_mda = original_structure.select_atoms(\"resname 13U\")\n",
    "\n",
    "# Get the center of the ligand as the \"pocket center\"\n",
    "pocket_center = ## Fill in code for measuring the center of geometry of the ligand.\n",
    "print(pocket_center)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the pocket center, we will define our ligand box.\n",
    "One simple approach to this is to subtract the min and max of the ligand positions in each dimension.\n",
    "In order to allow for ligand flexibility and potential interactions with nearby residues, we will add an additional five angstroms to each side of our box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute min and max coordinates of the ligand\n",
    "# take the ligand box to be the difference between the max and min in each direction.\n",
    "ligand_box = ligand_mda.positions.max(axis=0) - ## Fill in the rest of this expression.\n",
    "ligand_box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pocket_center` and `ligand_box` variables are NumPy arrays.\n",
    "However, AutoDock Vina expects them to be lists.\n",
    "We convert them to lists in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pocket_center = pocket_center.tolist()\n",
    "ligand_box = ## convert ligand_box to list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "<strong>Error in the cell above?</strong>  \n",
    "\n",
    "If you execute the cell above this one more than once, you will see an error occur. \n",
    "This happens because on first execution, `pocket_center` and `ligand_box` are NumPy arrays and they have methods to convert to lists.\n",
    "After you've executed this once, `pocket_center` and `ligand_box` don't have a `tolist` method because they are already lists.\n",
    "\n",
    "If you execute the cell twice and see an error, you can continue with the rest of the notebook because the variables have already been converted.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Docking Ligands with AutoDock Vina\n",
    "\n",
    "Now that we have PDBQT files of our protein and ligand and have defined our docking box, we are ready to perform the actual docking.\n",
    "Before docking, we will make a directory to store our results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a directory to store our results\n",
    "import os\n",
    "\n",
    "pdb_id = \"2zq2\"\n",
    "ligand = \"13U\"\n",
    "\n",
    "## Make a directory called \"docking_results\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will dock using the AutoDock Vina Python API.\n",
    "First, we import `Vina` from `vina`.\n",
    "We start docking with the line `v = Vina(sf_name=\"vina\")`. \n",
    "This creates a docking calculation, `v`, and sets the scoring function to the `vina` scoring function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vina import Vina\n",
    "v = Vina(sf_name=\"vina\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<strong>Scoring Functions in AutoDock Vina</strong>\n",
    "\n",
    "* Vina (`vina`): `vina` is an empirical scoring function. Binding energy is predicted as the sum of pairwise atomic interactions. It includes terms for hydrogen bonds, hydrophobic interactions, and steric interactions. The parameters for this scoring function were empirically derived from fitting data available in the PDBbind database. You can read more in the [original publication](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3041641/), or in the [Vinardo paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4865195/)\n",
    "* Vinardo (`vinardo`): Vinaro stands for \"Vina RaDii Optimized\". It was developed to improve the scoring by adjusting atom radii and reparameterizing the empirical terms based on the PDBBIND 2013 database. You can read more in the [Vinardo paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4865195/).\n",
    "* AutoDock4 (`ad4`):  uses a physics-based model and is the most computationally intensive of the available scores. The `ad4` score requires the definition of a flexible receptor, so it won't work with the PDBQT we have prepared. If you are interested in flexible docking, see [the tutorial from AutoDock Vina](https://autodock-vina.readthedocs.io/en/latest/docking_flexible.html). <strong>If you try to use the `ad4` scoring function on a receptor that was not prepared to be a flexible receptor, your notebook kernel will crash.</strong>\n",
    "\n",
    "</div>\n",
    "\n",
    "Then, we set the files for our ligand and receptor. We will dock just our ideal ligand first. There are two parameters to docking, the `exhaustiveness` and `n_poses`.\n",
    "The exhaustiveness parameter describes the \"exhaustiveness\" of the docking - a higher exhaustiveness means taht more ligand conformations are tried. Exhaustiveness also corresponds to the amount of computational effort used during a docking experiment. The default exhaustiveness value is 8; increasing this to 32 will give a more consistent docking result. \n",
    "\n",
    "In this notebook, we set the exhaustiveness to 5 to improve speed for the workshop. If you were to do a real docking calculation, you should consider increasing this parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.set_receptor(f\"pdbqt/{pdb_id}.pdbqt\")\n",
    "v.set_ligand_from_file(f\"\") ## Fill in the ligand name\n",
    "v.compute_vina_maps(center=, box_size=) ## Fill in the ligand center and box size calculated previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.dock(exhaustiveness=5, n_poses=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the `dock` function, we can write the poses that were calculated to a file.\n",
    "Note that the output format from AutoDock Vina is a PDBQT file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.write_poses(f\"docking_results/{ligand}.pdbqt\", n_poses=5, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the energies of the calculated poses by calling `energies` on the docking calculation variable.\n",
    "According to the Vina documentaiton, the rows correspond to the poses, while columns correspond to different energy types.\n",
    "The types of energies in the columns are `[\"total\", \"inter\", \"intra\", \"torsions\", \"intra best pose\"]`.\n",
    "The number of columns and the types of energies they represent depend on the scoring function you are using.\n",
    "You can see more information in the [docs for AutoDock Vina](https://autodock-vina.readthedocs.io/en/latest/vina.html#vina.vina.Vina.energies)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Fill in code to see the calculated energies for the poses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might wish to save these energies to return to them later. \n",
    "The cell below creates a pandas dataframe and saves the energies as a comma-separated-value (CSV) file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "\n",
    "# These are the columns for the types of energies according to AutoDock Vina docs.\n",
    "column_names = [\"total\", \"inter\", \"intra\", \"torsions\", \"intra best pose\"]\n",
    "\n",
    "df = pd.DataFrame() ## Fill in arguments to the dataframe to create a dataframe with energies.\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the calculated energies from docking to a CSV file\n",
    "df.to_csv(\"\", index=False) ## Fill in file name."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Docking Results\n",
    "\n",
    "After performing the docking simulation and saving the energies, you might wish to visualize the poses. \n",
    "When visualizing results from molecular docking, scientists often visually inspect the 3D docked structure as well as a 2D representation called an interaction map.\n",
    "We can ues a software called ProLIF (Protein-Ligand Interaction Fingerprints) to make and view these maps in the Jupyter notebook.\n",
    "\n",
    "To generate these visualizations, we have to convert our files (again!) to the correct format.\n",
    "\n",
    "In the step above, we wrote the poses to the file `docking_results/13U.pdbqt`. \n",
    "AutoDock Vina only writes in this file, but in order to visualize your results, we need a more standard format.\n",
    "We will use meeko again to convert our poses to an SDF.\n",
    "Note that meeko will only convert pdbqt files if it prepared the input docking files, which is one reason we used it in the previous notebook.\n",
    "\n",
    "Again, we use a command line script to convert out poses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mk_export.py docking_results/13U.pdbqt -s docking_results/13U.sdf # In the original workshop, `mk_export.py` used the option `-o` instead of `-s`. This was changed to `-s` in late 2024."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After converting to SDF, we can again visualize our results with ProLIF.\n",
    "ProLIF requires that molecules be loaded in and has functions to load molecules in several ways.\n",
    "We will use MDAnalysis for loading our proteins to ProLIF and `sdf_supplier` to load the SDFs we converted in the previous step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import prolif as plf\n",
    "import MDAnalysis as mda\n",
    "\n",
    "pdb_id = \"2zq2\"\n",
    "\n",
    "protein = mda.Universe(f\"protein_structures/protein_h.pdb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we load our protein and ligand into ProLIF.\n",
    "The function we use to do this depends on the format of our input data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_plf = plf.Molecule.from_mda(protein)\n",
    "poses_plf = plf.sdf_supplier() # fill in filename here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To analyze the interactions of the ligand and protein we create a molecular fingerprint object.\n",
    "By default, ProLIF will calculate nine types of interactions: 'Hydrophobic', 'HBAcceptor', 'HBDonor', 'Cationic', 'Anionic', 'CationPi', 'PiCation', 'PiStacking', 'VdWContact'\n",
    "However, you could set this to different interactions. You can see more information about the types of interactions in the [ProLIF docs](https://prolif.readthedocs.io/en/latest/source/modules/interaction-fingerprint.html#detecting-interactions-between-residues-prolif-interactions-interactions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = plf.Fingerprint(count=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will run ProLif on our poses.\n",
    "To do this calculation, we pass in our list of poses (`poses_plf`) and our ProLIF protein."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run on your poses\n",
    "fp.run_from_iterable() ## fill in this function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running this analysis, we can visualize the interaction results.\n",
    "We are using the 2D and 3D visualization maps here, but there are [many other types of analysis](https://prolif.readthedocs.io/en/latest/notebooks/docking.html#analysis) that you can perform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pose_index=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fp.plot_lignetwork(poses_plf[pose_index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view = fp.plot_3d(\n",
    "    poses_plf[pose_index], protein_plf, frame=pose_index, display_all=False\n",
    ")\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h3>Exercise</h3>\n",
    "\n",
    "Try docking one of the ligands we modified in the previous notebook. Does it bind better or worse according to the docking score? Are the interactions different for the poses?\n",
    "</div>"
   ]
  },
  {
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}