docking_preparation.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 Preparation\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"> \n",
    "<h2>Overview</h2>\n",
    "    \n",
    "<strong>Questions</strong>\n",
    "\n",
    "* How do I prepare structure files for docking?\n",
    "\n",
    "<strong>Learning Objectives</strong>\n",
    "\n",
    "* Use the RCSB Search API to find a protein structure for docking.\n",
    "* Prepare a protein and ligand structures for input into AutoDock Vina.\n",
    "</div>\n",
    "\n",
    "In this notebook, we will get our structure files ready for docking.\n",
    "We will need to create a special file format called pdbqt, which is used by AutoDock Vina.\n",
    "The PDBQT format is similar to the PDB format, but it includes additional information such as atomic charges.\n",
    "\n",
    "We will need to complete a few steps:\n",
    "\n",
    "1. Download a trypsin protein structure with our ligand of interest docked from the PDB.\n",
    "2. Isolate the protein in the structure (strip any water molecules and ligand).\n",
    "3. Add hydrogens to the protein and clean up the structure.\n",
    "4. Create a PDBQT file for our protein.\n",
    "5. Create PDBQT files for our ligands.\n",
    "### Python Libraries for the IQB workshop\n",
    "\n",
    "We'll be using many of the same libraries that we used in previous notebooks.\n",
    "However, this time we will add one more library to our list - [MDAnalysis](https://www.mdanalysis.org/).\n",
    "MDAnalysis is a Python library that was written for analyzing molecular dynamics (MD) simulations.\n",
    "However, it has reading and writing capabilities for many molecular file formats, as well as a selection language for isolating particular parts of molecules.\n",
    "\n",
    "| Library         | abbreviation | Purpose |\n",
    "|:-------------|:---------:|:------------|\n",
    "| rcsbsearchapi | N/A      | functions for searching the Protein Data Bank based on the mmCIF dictionary |\n",
    "| requests     | N/A  | access web URLs - used with APIs for databases |\n",
    "| os           | N/A      | operating system functions - handling file paths and directories. |\n",
    "| nglview      | nv       | for viewing molecular structures |\n",
    "| MDAnalysis     | mda | molecular dynamics library - used for reading/writing files and selecting atoms |\n",
    "\n",
    "### Command Line Software Used in this Notebook\n",
    "\n",
    "We'll also be adding a few command line scripts and utilities to this notebook.\n",
    "Usually, these would be executed in the terminal, but we'll run them from the Jupyter interface.\n",
    "These will be used to prepare our structures for docking calculations\n",
    "\n",
    "| Software         | Purpose |\n",
    "|:-------------|:---------|\n",
    "| pdb2pqr      | adding hydrogens and missing atoms to protein, adjusting for pH |\n",
    "| meeko        | preparing ligands for docking |\n",
    "\n",
    "## Finding and saving a protein structure\n",
    "\n",
    "To find an appropriate protein for our docking study, we can again use the PDB Search API, similar to how we used it in the first notebook.\n",
    "This time, we will use the same EC number, but the second part of our query will be our ligand name.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rcsbsearchapi.search import TextQuery\n",
    "from rcsbsearchapi import rcsb_attributes as attrs\n",
    "\n",
    "ECnumber = \"3.4.21.4\"    \n",
    "\n",
    "q1 = attrs.rcsb_polymer_entity.rcsb_ec_lineage.id == ECnumber  # looking for trypsins\n",
    "q2 = ## Fill in text query here\n",
    "\n",
    "query = q1 & q2               # combining the two queries into one\n",
    "\n",
    "results = list(query())\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdb_id = ## Get the PDB ID from the list and convert it to lowercase\n",
    "print(pdb_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have identified the PDB ID of a trypsin with our ligand of interest bound. \n",
    "Our PDB ID is `2zq2`, and we saved it in a variable in the cell above.\n",
    "\n",
    "We will now create a directory called `protein_structures` to save this in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os # for making directories\n",
    "import requests\n",
    "\n",
    "# make a directory for pdb files\n",
    "protein_directory = \"protein_structures\"\n",
    "\n",
    "## fill in function to make directories\n",
    "\n",
    "pdb_request = requests.get(f\"https://files.rcsb.org/download/{pdb_id}.pdb\")\n",
    "pdb_request.status_code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f\"{protein_directory}/{pdb_id}.pdb\", \"w+\") as f:\n",
    "    ## fill in write command"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the protein strucure\n",
    "\n",
    "Before we start to really work with our molecule, let's investigate the structure.\n",
    "We will use a library called MDAnalysis to first process our PDB, then visualize it with a library called NGLView.\n",
    "\n",
    "[MDAnalysis](https://www.mdanalysis.org/) is a Python library that is used to process molecular dynamics trajectories and other 3D strucure molecular files.\n",
    "The core object for MDAnalysis is a \"Universe\" and it corresponds to a molecular system.\n",
    "We can load a PDB file into MDAnalysis, then do things like measure distances in our structure or isolate particular parts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "\n",
    "# Load into MDA universe\n",
    "u = \n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After loading our PDB, we can see that we have an MDAnalysis \"universe\" (or molecular system) that contains 2024 atoms. \n",
    "We can inspect this structure visually using a library called NGLView.\n",
    "NGLView is a molecular visualizer made to work on the web and Jupyter notebooks.\n",
    "If you're interested in learning about NGLView, you can see a [video from MolSSI's first Crash Course with the PDB](https://www.youtube.com/watch?v=6QHWhycMuXc). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview as nv\n",
    "## create and show an NGLView from an MDAnalysis universe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This view looks a bit messy, and we likely want to isolate the protein and ligand for viewing.\n",
    "MDAnalysis has a human readable [selection syntax](https://docs.mdanalysis.org/stable/documentation_pages/selections.html)\n",
    "that allows us to isolate parts of our structure. We will take our MDAnalysis Universe (the variable `u`) and use the `select_atoms` function.\n",
    "Inside this function, we will fill in what we want to select.\n",
    "\n",
    "We will create separate variables for the protein and ligand. We can select all protein residues in MDAnalysis using the word \"protein\" in the `select_atoms` function. Then, we will select our ligand using `resname 13U`. This corresponds to the residue name in the PDB we downloaded.\n",
    "We can also select waters in the structure by using `\"resname HOH\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select protein atoms\n",
    "protein = u.select_atoms(\"protein\")\n",
    "ligand = u.select_atoms(\"resname 13U\")\n",
    "water = ## fill in selection for water.\n",
    "\n",
    "water"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After selecting parts of our structure, we can add them individually to an NGLView view.\n",
    "In the cell below, we visualize the protein's surface area colored by hydrophobicity.\n",
    "Waters from the crystal structure are in spacefill representation, and we add the ligand in a ball and stick representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "view = nv.show_mdanalysis(protein)\n",
    "view.clear_representations()\n",
    "## Add representation to view protein surface colored by hydrophobicity\n",
    "lig_view = view.add_component(ligand)\n",
    "lig_view.add_representation(\"ball+stick\")\n",
    "water_view = view.add_component(water)\n",
    "water_view.add_representation(\"spacefill\")\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you rotate this structure so that you are looking at the bottom, you will be able to see our `13U` ligand bound.\n",
    "Upon viewing this structure, you will notice that our ligand seems to appear twice. \n",
    "If you open the PDB file to investigate, you will see the following in the ligand section:\n",
    "\n",
    "```\n",
    "HETATM 1673  C14A13U A 501      18.144  -9.216  12.088  0.61 24.22           C  \n",
    "ANISOU 1673  C14A13U A 501     1755   4793   2654   1752    148   1233       C  \n",
    "HETATM 1674  C14B13U A 501      18.147  -8.840  11.672  0.39 24.46           C  \n",
    "ANISOU 1674  C14B13U A 501     2583   4283   2430   1765    353   1279       C  \n",
    "HETATM 1675  O32A13U A 501      18.209  -8.355  11.186  0.61 24.38           O  \n",
    "ANISOU 1675  O32A13U A 501     2354   5394   1514   2217    238    919       O\n",
    "```\n",
    "\n",
    "This PDB structure provides [**alternate locations**](https://proteopedia.org/wiki/index.php/Alternate_locations) for each ligand atom. \n",
    "These occur when the experimental data supports multiple positions for the same atom.\n",
    "In excerpt above, you will see C14A13U and C14B13U. These are alternate locations of the same atom. \n",
    "Alternate locations can also occur in the protein with some residues.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<strong>Selecting alternate locations using MDAnalysis</strong>  \n",
    "    \n",
    "By checking the [documentation page](https://docs.mdanalysis.org/stable/documentation_pages/selections.html) for MDAnalysis selections, we can see that MDAnalysis is prepared for this scenario. We will want to use the `altloc` keyword. This keyword is described as:\n",
    "\n",
    "> altLoc alternative-location\n",
    "\n",
    "> a selection for atoms where alternative locations are available, which is often the case with high-resolution crystal structures e.g. resid 4 and resname ALA and altLoc B selects only the atoms of ALA-4 that have an altLoc B record.\n",
    "\n",
    "If you wanted to use MDAnalysis to select for a particular ligand location, you could use:\n",
    "\n",
    "```python\n",
    "ligand_A = u.select_atoms(\"resname 13U and altLoc A\")\n",
    "ligand_B = u.select_atoms(\"resname 13U and altLoc B\")\n",
    "```\n",
    "</div>\n",
    "\n",
    "To perform a docking calculation, we will have to isolate the protein.\n",
    "This starting structure for our protein contains extra molecules like ligands and water.\n",
    "You will also notice from examining our visualization that our structure does not include hydrogen atoms.\n",
    "If you were to examine the PDB file, you would also see that there are some missing atoms and some of our protein residues have alternate locations marked just like thie ligand.\n",
    "\n",
    "For docking, we will want to remove all of these extra molecules and only keep the protein.\n",
    "We will also want to remove any alternate locations of residues.\n",
    "We will use MDAnalysis to remove these extra molecules and save our starting protein structure as a new file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write protein to new PDB file\n",
    "protein.write(f\"{protein_directory}/protein_{pdb_id}.pdb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "<strong>Protein Charge</strong>  \n",
    "\n",
    "After saving the protein in the cell above, you may see a warning about information for formal charges not being set in the protein. \n",
    "This warning appears because MDAnalysis did not find specific formal charge data in the PDB file and used a default value instead. \n",
    "This is not a concern for us because we will adjust the protonation states of different residues using PDB2PQR in the next steps. \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixing the protein structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've isolated our protein, we will want to ensure that we've correctly added hydrogen and fixed any missing atoms.\n",
    "\n",
    "For fixing our protein, we will use a specialized program called PDB2PQR that is made for working with biomolecules like proteins.\n",
    "The advantage of using PDB2PQR is that it will check our protein for missing atoms and multiple occupancy in the protein, and it will pick positions and add missing atoms.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<strong>More complicated fixes: PDBFixer</strong>  \n",
    "\n",
    "Another popular software for fixing PDB files is called [PDBFixer](https://github.com/openmm/pdbfixer). PDBFixer is an open-source tool developed by the OpenMM team and is designed to fix common problems in Protein Data Bank (PDB) files before they are used in molecular simulations. It can be used to remove unwanted molecules like water, add missing heavy atoms to incomplete residues, add hydrogen atoms where needed.\n",
    "\n",
    "PDBFixer can be especially useful when there are missing loops or residues. In this workshop, our protein is not missing any residues, but many proteins from the PDB might require more adjustment.\n",
    "\n",
    "To see an example of preparing proteins with PDBFixer, see this [recent YouTube video](https://www.youtube.com/watch?v=pwfKE6wPaMg) posted by the Open Forcefield Initiative. In this video, the presenter first uses PDBFixer, then PDB2PQR to adjust protonation.\n",
    "\n",
    "</div>\n",
    "\n",
    "We will use the command-line interface of this PDB2PQR. This means that you would usually type the command below into your terminal\n",
    "You can run command line commands in the Jupyter notebook by putting a `!` in front of the command. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! pdb2pqr --pdb-output=protein_structures/protein_h.pdb --pH=7.4 protein_structures/protein_2zq2.pdb protein_structures/protein_2zq2.pqr --whitespace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving a protein PDBQT File"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PDB2PQR program outputs two files, a PDB file and a PQR file. The PDB file is similar to PDB files we have worked with before, except that it contains hydrogens.\n",
    "The PQR file is another molecular file format that is similar to a PDB, but contains information about atomic radii and atomic charges.\n",
    "\n",
    "For use with AutoDock Vina, we need our protein file to be in the \"PDBQT\" format. PDBQT is a specialized file format used by AutoDock Vina and other AutoDock tools. Like the PQR format, the PDBQT format can also contain partial charges. We will load our PQR file and use MDAnalysis to write a PDBQT file.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "<strong>What information do we need in the PDBQT file?</strong>\n",
    "\n",
    "We'll be using AutoDock Vina with the \"vina\" scoring function (this will be explained in more detail in the next notebook). The vina scoring function doesn't use charges to dock, so we could have also used the PDB file without charges to convert to a PDBQT file. However, some scoring functions do use partial charges.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a directory for pdb files\n",
    "pdbqt_directory = ## fill in the name of the directory to write PDBQT files to\n",
    "os.makedirs(pdbqt_directory, exist_ok=True)\n",
    "\n",
    "u = mda.Universe(f\"{protein_directory}/protein_{pdb_id}.pqr\")\n",
    "u.atoms.write(f\"{pdbqt_directory}/{pdb_id}.pdbqt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PDBQT file generated by MDAnalysis includes two lines at the start of the structure that AutoDock Vina doesn't accept. \n",
    "These lines start with \"TITLE\" and \"CRYST1\". To resolve this, the following cell replaces these lines with \"REMARK\", which is acceptable to AutoDock Vina.\n",
    "\n",
    "You might have also just chosen to use a different software to write your PDBQT. \n",
    "[OpenBabel](https://openbabel.org/index.html) is a popular choice. However, we are using MDAnalysis here for consistency with the rest of the workshop and to limit the number of libraries we are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the just-written PDBQT file, replace text, and write back\n",
    "with open(f\"{pdbqt_directory}/{pdb_id}.pdbqt\", 'r') as file:\n",
    "    file_content = file.read()\n",
    "\n",
    "# Replace 'TITLE' and 'CRYST1' with 'REMARK'\n",
    "file_content = file_content.replace('TITLE', 'REMARK').replace('CRYST1', 'REMARK')\n",
    "\n",
    "# Write the modified content back to the file\n",
    "with open(f\"{pdbqt_directory}/{pdb_id}.pdbqt\", 'w') as file:\n",
    "    file.write(file_content)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ligand Preparation\n",
    "\n",
    "When preparing small molecule PDBQT files, you could have also chosen to use MDAnalysis or other tools.\n",
    "However, we are going to use a special program for small molecules and docking called [meeko](https://github.com/forlilab/Meeko).\n",
    "We choose to use meeko for our ligands because it will allow us to more easily visualize our results later.\n",
    "Note that when using meeko, ligands should have hydrogens added already.\n",
    "\n",
    "We are using the command line for meeko, similar to PDB2PQR. \n",
    "You could also choose to use the Python API for this, but the command line is simpler for common tasks like converting an SDF to a PDBQT.\n",
    "\n",
    "In the cell below, we execute a command that converts our ligands that in we prepared in the `molecule_manipulation` notebook to a PDBQT file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use meeko to prepare small molecules - using meeko helps us visualize them later.\n",
    "! mk_prepare_ligand.py -i ligands_to_dock/13U.sdf -o pdbqt/13U.pdbqt\n",
    "! mk_prepare_ligand.py -i ligands_to_dock/13U_modified_methyl.sdf -o pdbqt/13U_modified_methyl.pdbqt\n",
    "# Fill in the command for converting the N modified ligand."
   ]
  },
  {
   "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
}