notebooks/Analysis_2q7q.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Most Beautiful Disulfide Bond in the World\n",
    "Eric G. Suchanek, PhD 2/24/24\n",
    "\n",
    "In this notebook illustrate some of the features of proteusPy by analyzing the lowest energy disulfide bond in the RCSB protein structure database. If you are not familiar with `proteusPy` you can find the API at: https://suchanek.githubio.com/proteusPy.html\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pyvista as pv\n",
    "from pyvista import set_plot_theme\n",
    "from proteusPy import Disulfide, DisulfideList, Load_PDB_SS\n",
    "\n",
    "# pyvista setup for notebooks\n",
    "pv.set_jupyter_backend(\"trame\")\n",
    "\n",
    "set_plot_theme(\"default\")\n",
    "LIGHT = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the RCSB Disulfide Database\n",
    "We load the database and get its properties as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--> DisulfideLoader: Downloading Disulfide Database from Drive...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading...\n",
      "From (original): https://drive.google.com/uc?id=1igF-sppLPaNsBaUS7nkb13vtOGZZmsFp\n",
      "From (redirected): https://drive.google.com/uc?id=1igF-sppLPaNsBaUS7nkb13vtOGZZmsFp&confirm=t&uuid=95441470-1539-40af-8463-8236bf404286\n",
      "To: /Users/egs/miniforge3/envs/proteusPy/lib/python3.11/site-packages/proteusPy/data/PDB_SS_ALL_LOADER.pkl\n",
      "100%|████████████████████████████████████████████████████████████████| 340M/340M [00:33<00:00, 10.1MB/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--> DisulfideLoader: Downloading Disulfide Subset Database from Drive...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading...\n",
      "From: https://drive.google.com/uc?id=1puy9pxrClFks0KN9q5PPV_ONKvL-hg33\n",
      "To: /Users/egs/miniforge3/envs/proteusPy/lib/python3.11/site-packages/proteusPy/data/PDB_SS_SUBSET_LOADER.pkl\n",
      "100%|██████████████████████████████████████████████████████████████| 9.64M/9.64M [00:01<00:00, 8.45MB/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-> load_PDB_SS(): Reading /Users/egs/miniforge3/envs/proteusPy/lib/python3.11/site-packages/proteusPy/data/PDB_SS_ALL_LOADER.pkl... \n",
      "-> load_PDB_SS(): Done reading /Users/egs/miniforge3/envs/proteusPy/lib/python3.11/site-packages/proteusPy/data/PDB_SS_ALL_LOADER.pkl... \n",
      "    =========== RCSB Disulfide Database Summary ==============\n",
      "       =========== Built: 2024-03-16 17:12:24 ==============\n",
      "PDB IDs present:                    35818\n",
      "Disulfides loaded:                  120494\n",
      "Average structure resolution:       2.34 Å\n",
      "Lowest Energy Disulfide:            2q7q_75D_140D\n",
      "Highest Energy Disulfide:           1toz_456A_467A\n",
      "Cα distance cutoff:                 8.00 Å\n",
      "Total RAM Used:                     30.72 GB.\n",
      "    ================= proteusPy: 0.93.0 =======================\n"
     ]
    }
   ],
   "source": [
    "PDB_SS = Load_PDB_SS(verbose=True)\n",
    "PDB_SS.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see from the statistics above that disulfide 2q7q_75D_140D has the lowest energy, so let's extract it from the database and display it. \n",
    "\n",
    "A few notes about the display window. You might need to click into the window to refresh it. Click drag to rotate the structures, mousewheel to zoom. The window titles display several parameters about the disulfide bonds including their approximate torsional energy, their Ca-Ca distance, and the *torsion length*. \n",
    "\n",
    "The latter parameter is formally, the Euclidean length of the sidechain dihedral angle when treated as a five-dimensional vector. This sounds all mathy and complicated, but in essence it gives a measure of how 'long' that five dimensional vector is. This is used by the package to compare individual structures and gauge their structural similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Disulfide 2q7q_75D_140D, Source: 2q7q, Resolution: 1.6 Å \n",
      "Χ1-Χ5: -59.36°, -59.28°, -83.66°, -59.82° -59.91°, -25.17°, 0.49 kcal/mol \n",
      "Cα Distance: 5.50 Å \n",
      "Torsion length: 145.62 deg>\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4860cc7de8e44ae88daead1d86678e61",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value='<iframe src=\"http://localhost:64666/index.html?ui=P_0x28b6d0090_0&reconnect=auto\" class=\"pyvista…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ssmin, ssmax = PDB_SS.SSList.minmax_energy\n",
    "ssmin_energy = ssmin.energy\n",
    "\n",
    "best_ss = PDB_SS[\"2q7q_75D_140D\"]\n",
    "best_dihedrals = best_ss.dihedrals\n",
    "best_ss.pprint()\n",
    "best_ss.display(style=\"sb\", light=LIGHT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that, gentle reader, is it. *The most beautiful disulfide bond in the world*! Look at it. This is the lowest energy structure in the entire database. The sidechain dihdedral angles (Χ1-Χ5: -59.36°, -59.28°, -83.66°, -59.82° -59.91°), and the estimated energy, (0.49 kcal/mol). \n",
    "\n",
    "How does this compare to an analytical (modelled) minimum? We can use the ``minimize`` module from ``scipy`` to check. We know from chemistry that a reasonable guess for a low energy conformation would have the dihedral angles: (Χ1-Χ5: -60.00°, -60.00°, -90.00°, -60.00° -60.00°, 0.60 kcal/mol). Let's run this through scipy and compute a minimum energy conformation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Modeled minimum energy: 0.489 kcal/mol for conformation: ['-60.000', '-60.000', '-83.048', '-60.000', '-60.000']\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "from proteusPy.Disulfide import Disulfide_Energy_Function\n",
    "\n",
    "initial_guess = [\n",
    "    -60.0,\n",
    "    -60.0,\n",
    "    -90.0,\n",
    "    -60.0,\n",
    "    -60.0,\n",
    "]  # initial guess for chi1, chi2, chi3, chi4, chi5\n",
    "result = minimize(Disulfide_Energy_Function, initial_guess, method=\"Nelder-Mead\")\n",
    "minimum_energy = result.fun\n",
    "min_conf = result.x\n",
    "print(\n",
    "    f'Modeled minimum energy: {minimum_energy:.3f} kcal/mol for conformation: {[f\"{x:.3f}\" for x in min_conf]}'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the computed minimum energy structure is *0.489* kcal/mol, estimated. The difference from the actual conformation is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Modeled - actual energy difference is: -0.002797862530647066 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "diff = minimum_energy - ssmin_energy\n",
    "print(f\"Modeled - actual energy difference is: {diff} kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The real structure is actually at a *lower* energy than the predicted! What an amazing Disulfide bond! Let's build a model for its conformation. We do that by creating an empty disulfide and then using the `Disulfide.build_yourself` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "modelled_min = Disulfide(\"model\")\n",
    "modelled_min.dihedrals = min_conf\n",
    "modelled_min.build_yourself()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now make a ``DisulfideList`` list and put the real structure and modelled structure into it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "minmax = DisulfideList([modelled_min, ssmin], \"minmax\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, display them in a common reference frame:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "minmax.display_overlay()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two structures overlap with an overall RMS error of 2.14 A. Not bad considering the modeled structure is made using idealized bond lengths and angles. This serves to remind us that nature is the ultimate modeling engine, and that we still have much to learn."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "* *Application of Artificial Intelligence in Protein Design* - Doctoral Dissertation, EG Suchanek, 1987, Johns Hopkins Medical School\n",
    "* https://doi.org/10.1021/bi00368a023\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "proteusPy",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}