notebooks/basic_molecular_docking.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HlvKcskyTzu4"
      },
      "source": [
        "# **Basic Molecular Docking**\n",
        "\n",
        "This <font color='orange'>**LABO**</font>**DOCK** notebook is designed for conducting basic molecular docking procedures using **Autodock Vina 1.2.5** and performing binding interaction analysis with **PLIP 2.3.0**.\n",
        "\n",
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/RyanZR/labodock/blob/main/notebooks/basic_molecular_docking.ipynb)\n",
        "[![version](https://img.shields.io/badge/version-2.0.0-blue)](https://github.com/RyanZR/labodock/tree/main/notebooks)\n",
        "[![DOI](https://img.shields.io/badge/DOI-10.5281/zenodo.8246977-blue)](https://doi.org/10.5281/zenodo.8246977)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FDIu2L1jfx5_"
      },
      "source": [
        "---\n",
        "---\n",
        "# **01 | Setting Up Environment**\n",
        "\n",
        "Firstly, we install all the necessary libraries and packages."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zr-anaLWgaHq"
      },
      "source": [
        "The main packages that will be installed include:\n",
        "\n",
        "+ AutoDock Vina 1.2.5 (https://vina.scripps.edu/)\n",
        "+ Condacolab 0.1.7 (https://github.com/conda-incubator/condacolab)\n",
        "+ Meeko 0.5.0 (https://github.com/forlilab/Meeko)\n",
        "+ OpenBabel 3.1.1 (https://github.com/openbabel/openbabel)\n",
        "+ PLIP 2.3.0 (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index)\n",
        "+ Py3Dmol 2.0.3 (https://pypi.org/project/py3Dmol/)\n",
        "+ Rdkit 2022.9.5 (https://github.com/rdkit/rdkit)\n",
        "+ Spyrmsd 0.6.0 (https://github.com/RMeli/spyrmsd)\n",
        "\n",
        "All the above libraries and packages are currently supported by **Python 3.10**."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dr15mv6Ner0h",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Install packages and dependencies**\n",
        "# @markdown Thanks to **`mamba`**, the installation takes **less than 2 mins**. \\\n",
        "# @markdown It will **restart** the kernel (session).\n",
        "\n",
        "import sys\n",
        "import time\n",
        "import contextlib\n",
        "\n",
        "with open('/content/labodock_install.log', 'a') as inpt:\n",
        "    with contextlib.redirect_stdout(inpt):\n",
        "\n",
        "        # -- Start installation --\n",
        "        start = time.time()\n",
        "        !rm -r /content/sample_data\n",
        "        !wget https://github.com/ccsb-scripps/AutoDock-Vina/releases/download/v1.2.5/vina_1.2.5_linux_x86_64 -O vina\n",
        "        !chmod u+x vina\n",
        "\n",
        "        !pip install py3Dmol==2.0.3\n",
        "        !pip install rdkit-pypi==2022.9.5\n",
        "        !pip install meeko==0.5.0\n",
        "        !pip install condacolab==0.1.7\n",
        "\n",
        "        import condacolab\n",
        "        condacolab.install_mambaforge()\n",
        "        !mamba install -c conda-forge spyrmsd=0.6.0 openbabel=3.1.1 plip=2.3.0\n",
        "        end = time.time()\n",
        "        # -- End installation --\n",
        "\n",
        "        print(f'+ Time elapsed: ' + time.strftime('%Mm %Ss', time.gmtime(end - start)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3X-mmyEnfimS",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Import Python modules**\n",
        "# @markdown This will import the necessary Python modules.\n",
        "\n",
        "# Internal Modules\n",
        "import os\n",
        "import glob\n",
        "import time\n",
        "import shutil\n",
        "\n",
        "# External modules\n",
        "import py3Dmol\n",
        "import plip\n",
        "\n",
        "# Data-related\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Google Colab-related\n",
        "from google.colab import drive, files\n",
        "from tqdm.notebook import tqdm\n",
        "\n",
        "# Docking-related\n",
        "from rdkit import Chem, RDLogger\n",
        "from rdkit.Chem import rdFMCS, AllChem, Draw\n",
        "from spyrmsd import io, rmsd\n",
        "\n",
        "# Binding interaction-related\n",
        "from plip.exchange.report import BindingSiteReport\n",
        "from plip.structure.preparation import PDBComplex\n",
        "\n",
        "print(f'+ Import completed')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PST20Bo-k4W5",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Create folders**\n",
        "# @markdown Enter a **< Job Name >** to create a working directory containing:\n",
        "# @markdown + **`PROTEIN`** : Stores the target protein\n",
        "# @markdown + **`LIGAND`** : Stores the ligand to be docked\n",
        "# @markdown + **`NATIVE`** : Stores the native ligand of target protein\n",
        "# @markdown + **`DOCK`** : Stores the docking pose and scoring data\n",
        "# @markdown + **`INTERACTION`** : Stores the binding interaction data\n",
        "\n",
        "Job_name = '7KBS'   # @param {type: 'string'}\n",
        "\n",
        "invalid_chars = '^<>/\\\\{}[]~`$ '\n",
        "assert Job_name, 'Do not leave this blank.'\n",
        "assert not set(invalid_chars).intersection(Job_name), 'Disallowed characters.'\n",
        "\n",
        "DIR = os.getcwd()\n",
        "WRK_DIR = os.path.join(DIR, Job_name)\n",
        "PRT_FLD = os.path.join(WRK_DIR, 'PROTEIN')\n",
        "LIG_FLD = os.path.join(WRK_DIR, 'LIGAND')\n",
        "NTV_FLD = os.path.join(WRK_DIR, 'NATIVE')\n",
        "DCK_FLD = os.path.join(WRK_DIR, 'DOCKING')\n",
        "INT_FLD = os.path.join(WRK_DIR, 'INTERACTION')\n",
        "\n",
        "folders = [WRK_DIR, PRT_FLD, LIG_FLD, NTV_FLD, DCK_FLD, INT_FLD]\n",
        "\n",
        "for folder in folders:\n",
        "    if os.path.exists(folder):\n",
        "        print(f'+ Folder existed: {folder}')\n",
        "    else:\n",
        "        os.mkdir(folder)\n",
        "        print(f'+ Folder created: {folder}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "peLezlh0pJ05",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Set up utilities**\n",
        "# @markdown This creates important custom functions and methods for later\n",
        "# @markdown docking and binding interaction study.\n",
        "\n",
        "%alias vina /content/vina\n",
        "\n",
        "#############################################\n",
        "# Suppress Warnings\n",
        "\n",
        "RDLogger.DisableLog('rdApp.warning')\n",
        "\n",
        "#############################################\n",
        "# Grid Box Calculation Methods\n",
        "\n",
        "class GridBox:\n",
        "\n",
        "    ranges = tuple[list[float], list[float], list[float]]\n",
        "    coords = tuple[float, float, float]\n",
        "    center_bxsize = tuple[tuple[float, float, float], tuple[float, float, float]]\n",
        "\n",
        "    def __init__(self, inpt_file: str) -> None:\n",
        "        self.inpt = open(inpt_file, 'r')\n",
        "        self.data = self.inpt.read()\n",
        "        self.cmol = Chem.MolFromPDBBlock(self.data)\n",
        "        self.conf = self.cmol.GetConformer()\n",
        "        self.ntom = self.cmol.GetNumAtoms()\n",
        "        self.inpt.close()\n",
        "\n",
        "    def update_gridbox(self, mol_block: str) -> None:\n",
        "        self.cmol = Chem.MolFromPDBBlock(mol_block)\n",
        "        self.conf = self.cmol.GetConformer()\n",
        "        self.ntom = self.cmol.GetNumAtoms()\n",
        "\n",
        "    def compute_coords(self) -> ranges:\n",
        "        x_coord = [self.conf.GetAtomPosition(c).x for c in range(self.ntom)]\n",
        "        y_coord = [self.conf.GetAtomPosition(c).y for c in range(self.ntom)]\n",
        "        z_coord = [self.conf.GetAtomPosition(c).z for c in range(self.ntom)]\n",
        "        return x_coord, y_coord, z_coord\n",
        "\n",
        "    def compute_ranges(self) -> ranges:\n",
        "        x, y, z = self.compute_coords()\n",
        "        x_range = [min(x), max(x)]\n",
        "        y_range = [min(y), max(y)]\n",
        "        z_range = [min(z), max(z)]\n",
        "        return x_range, y_range, z_range\n",
        "\n",
        "    def compute_center(self, use_range: bool = True) -> coords:\n",
        "        x, y, z = self.compute_ranges() if use_range else self.compute_coords()\n",
        "        x_center = round(np.mean(x), 3)\n",
        "        y_center = round(np.mean(y), 3)\n",
        "        z_center = round(np.mean(z), 3)\n",
        "        return x_center, y_center, z_center\n",
        "\n",
        "    def generate_res_molblock(self, residues_list: list[str]) -> str:\n",
        "        res_lines = [line for line in self.data.split('\\n')\n",
        "                     if line[22:26].lstrip() in residues_list\n",
        "                     and 'END' not in line]\n",
        "        res_block = '\\n'.join(res_lines)\n",
        "        return res_block\n",
        "\n",
        "    def labox(self, scale: float = 2.0) -> coords:\n",
        "        xr, yr, zr = self.compute_ranges()\n",
        "        center = self.compute_center()\n",
        "        bxsize = (round(abs(xr[0] - xr[1]) * scale, 3),\n",
        "                  round(abs(yr[0] - yr[1]) * scale, 3),\n",
        "                  round(abs(zr[0] - zr[1]) * scale, 3))\n",
        "        return center, bxsize\n",
        "\n",
        "    def eboxsize(self, gy_box_ratio: float = 0.23, modified: bool = False) -> center_bxsize:\n",
        "        xc, yc, zc = self.compute_coords()\n",
        "        center = self.compute_center(modified)\n",
        "        distsq = [(x-center[0])**2 + (y-center[1])**2 + (z-center[2])**2\n",
        "                  for x, y, z in zip(xc, yc, zc)]\n",
        "        bxsize = (round(np.sqrt(sum(distsq) / len(xc)) / gy_box_ratio, 3),) * 3\n",
        "        return center, bxsize\n",
        "\n",
        "    def autodock_grid(self) -> center_bxsize:\n",
        "        xr, yr, zr = self.compute_ranges()\n",
        "        center = self.compute_center()\n",
        "        bxsize = (22.5, 22.5, 22.5)\n",
        "        return center, bxsize\n",
        "\n",
        "    def defined_by_res(self, residue_number: str, scale: float = 1.25) -> center_bxsize:\n",
        "        res_list = residue_number.replace(',', ' ').split()\n",
        "        res_block = self.generate_res_molblock(res_list)\n",
        "        self.update_gridbox(res_block)\n",
        "        return self.labox(scale=scale)\n",
        "\n",
        "#############################################\n",
        "# RMSD Calculation Methods\n",
        "\n",
        "class ComputeRMSD:\n",
        "\n",
        "    def __init__(self) -> None:\n",
        "        self.MCS_mol = None\n",
        "        self.MCS_png = None\n",
        "\n",
        "    def load_molecule(self, inpt_file: str, remove_Hs: bool = True) -> tuple:\n",
        "        molecule = io.loadmol(inpt_file)\n",
        "        molecule.strip() if remove_Hs else None\n",
        "        name = os.path.basename(inpt_file).split('.')[0]\n",
        "        coor = molecule.coordinates\n",
        "        anum = molecule.atomicnums\n",
        "        mtrx = molecule.adjacency_matrix\n",
        "        cmol = Chem.MolFromPDBFile(inpt_file)\n",
        "        return name, coor, anum, mtrx, cmol\n",
        "\n",
        "    def mol_to_png(self, mol: object) -> object:\n",
        "        legend = 'Maximum Common Substructure'\n",
        "        png = Draw.MolToImage(mol, legend=legend)\n",
        "        return png\n",
        "\n",
        "    def find_MCS(self, ref: tuple, lig: tuple) -> object:\n",
        "        if self.MCS_mol is None:\n",
        "            MCS_obj = rdFMCS.FindMCS([ref[4], lig[4]])\n",
        "            MCS_mol = Chem.MolFromSmarts(MCS_obj.smartsString)\n",
        "            MCS_png = self.mol_to_png(MCS_mol)\n",
        "            self.MCS_mol = MCS_mol\n",
        "            self.MCS_png = MCS_png\n",
        "        return self.MCS_mol\n",
        "\n",
        "    def hung_RMSD(self, ref: tuple, lig: tuple) -> float:\n",
        "        try:\n",
        "            hRMSD = round(rmsd.hrmsd(ref[1], lig[1], ref[2], lig[2]), 3)\n",
        "        except:\n",
        "            hRMSD = 'ERROR'\n",
        "        return hRMSD\n",
        "\n",
        "    def symm_RMSD(self, ref: tuple, lig: tuple, minimise: bool = False) -> float:\n",
        "        try:\n",
        "            sRMSD = round(rmsd.symmrmsd(ref[1], lig[1], ref[2], lig[2], ref[3], lig[3], minimize=minimise), 3)\n",
        "        except:\n",
        "            sRMSD = 'ERROR'\n",
        "        return sRMSD\n",
        "\n",
        "    def labo_RMSD(self, ref: tuple, lig: tuple) -> float:\n",
        "        mol_substr = self.find_MCS(ref, lig)\n",
        "        ref_substr = ref[4].GetSubstructMatch(mol_substr)\n",
        "        lig_substr = lig[4].GetSubstructMatch(mol_substr)\n",
        "\n",
        "        distsq = []\n",
        "        for ref_atom, lig_atom in zip(ref_substr, lig_substr):\n",
        "            ref_pos = ref[4].GetConformer().GetAtomPosition(ref_atom)\n",
        "            lig_pos = lig[4].GetConformer().GetAtomPosition(lig_atom)\n",
        "            ref_coord = np.array((ref_pos.x, ref_pos.y, ref_pos.z))\n",
        "            lig_coord = np.array((lig_pos.x, lig_pos.y, lig_pos.z))\n",
        "            coo_dist = np.linalg.norm(ref_coord - lig_coord)\n",
        "            distsq.append(coo_dist ** 2)\n",
        "\n",
        "        try:\n",
        "            lRMSD = round(np.sqrt(sum(distsq)/len(distsq)), 3)\n",
        "        except:\n",
        "            lRMSD = 'ERROR'\n",
        "        return lRMSD\n",
        "\n",
        "    def rmsd_report(self,\n",
        "                    ref: tuple,\n",
        "                    lig: tuple,\n",
        "                    lRMSD: bool = True,\n",
        "                    hRMSD: bool = True,\n",
        "                    sRMSD: bool = True\n",
        "                    ) -> dict[str: list[float]]:\n",
        "        report = {}\n",
        "        report['NAME'] = [lig[0]]\n",
        "        report['LABO_RMSD'] = [self.labo_RMSD(ref, lig)] if lRMSD else None\n",
        "        report['HUNG_RMSD'] = [self.hung_RMSD(ref, lig)] if hRMSD else None\n",
        "        report['SYMM_RMSD'] = [self.symm_RMSD(ref, lig)] if sRMSD else None\n",
        "        report = {k: v for k, v in report.items() if v is not None}\n",
        "        return report\n",
        "\n",
        "#############################################\n",
        "# AA Consntant and Bond Colour Dictionary\n",
        "\n",
        "# Kyte and Doolittle Hydropathy Scale (1982)\n",
        "AA_HB = {'ALA':  1.8, 'ARG': -4.5, 'ASN': -3.5, 'ASP': -3.5, 'CYS':  2.5,\n",
        "         'GLN': -3.5, 'GLU': -3.5, 'GLY': -0.4, 'HIS': -3.2, 'ILE':  4.5,\n",
        "         'LEU':  3.8, 'LYS': -3.9, 'MET':  1.9, 'PHE':  2.8, 'PRO': -1.6,\n",
        "         'SER': -0.8, 'THR': -0.7, 'TRP': -0.9, 'TYR': -1.3, 'VAL':  4.2}\n",
        "\n",
        "# University of Calgary PI Scale\n",
        "AA_PI = {'ALA':  6.0, 'ARG': 10.76, 'ASN': 5.41, 'ASP': 2.77, 'CYS': 5.07,\n",
        "         'GLN': 5.65, 'GLU':  3.22, 'GLY': 5.97, 'HIS': 7.59, 'ILE': 6.02,\n",
        "         'LEU': 5.98, 'LYS':  9.74, 'MET': 5.74, 'PHE': 5.48, 'PRO':  6.3,\n",
        "         'SEC': 5.68, 'SER':  5.68, 'THR':  5.6, 'TRP': 5.89, 'TYR': 5.66,\n",
        "         'VAL': 5.96}\n",
        "\n",
        "BOND_COL = {'HYDROPHOBIC': ['0x59e382', 'GREEN'],\n",
        "            'HBOND': ['0x59bee3', 'LIGHT BLUE'],\n",
        "            'WATERBRIDGE': ['0x4c4cff', 'BLUE'],\n",
        "            'SALTBRIDGE': ['0xefd033', 'YELLOW'],\n",
        "            'PISTACKING': ['0xb559e3', 'PURPLE'],\n",
        "            'PICATION': ['0xe359d8', 'VIOLET'],\n",
        "            'HALOGEN': ['0x59bee3', 'LIGHT BLUE'],\n",
        "            'METAL':['0xe35959', 'ORANGE']}\n",
        "\n",
        "#############################################\n",
        "# AA-to-Colour Converter Function\n",
        "\n",
        "def sequential_gradient(value: float,\n",
        "                        min_value: float,\n",
        "                        max_value: float,\n",
        "                        targ_colour: str = '00ff00',\n",
        "                        interpolation: float = 0.0\n",
        "                        ) -> str:\n",
        "    norm_val = (value - min_value) / (max_value - min_value)\n",
        "\n",
        "    rgb = tuple(int(targ_colour[d:d+2], 16) for d in (0, 2, 4))\n",
        "    r = int(255 - (255 - rgb[0]) * (1 - interpolation) * norm_val)\n",
        "    g = int(255 - (255 - rgb[1]) * (1 - interpolation) * norm_val)\n",
        "    b = int(255 - (255 - rgb[2]) * (1 - interpolation) * norm_val)\n",
        "\n",
        "    hex_code = f'#{r:02x}{g:02x}{b:02x}'\n",
        "    return hex_code\n",
        "\n",
        "def diverging_gradient(value: float,\n",
        "                       min_value: float,\n",
        "                       max_value: float,\n",
        "                       base_colour: str = 'ff0000',\n",
        "                       targ_colour: str = '0000ff',\n",
        "                       interpolation: float = 0.3\n",
        "                       ) -> str:\n",
        "    norm_val = (value - min_value) / (max_value - min_value)\n",
        "\n",
        "    white = (255, 255, 255)\n",
        "    rgb_A = tuple(int(base_colour[d:d+2], 16) for d in (0, 2, 4))\n",
        "    rgb_B = tuple(int(targ_colour[d:d+2], 16) for d in (0, 2, 4))\n",
        "\n",
        "    if norm_val < 0.5 - interpolation / 2:\n",
        "        factor = norm_val / (0.5 - interpolation / 2)\n",
        "        r = int(rgb_A[0] + (white[0] - rgb_A[0]) * factor)\n",
        "        g = int(rgb_A[1] + (white[1] - rgb_A[1]) * factor)\n",
        "        b = int(rgb_A[2] + (white[2] - rgb_A[2]) * factor)\n",
        "    elif norm_val > 0.5 + interpolation / 2:\n",
        "        factor = (norm_val - 0.5 - interpolation / 2) / (0.5 - interpolation / 2)\n",
        "        r = int(white[0] + (rgb_B[0] - white[0]) * factor)\n",
        "        g = int(white[1] + (rgb_B[1] - white[1]) * factor)\n",
        "        b = int(white[2] + (rgb_B[2] - white[2]) * factor)\n",
        "    else:\n",
        "        r, g, b = white\n",
        "\n",
        "    hex_code = f'#{r:02x}{g:02x}{b:02x}'\n",
        "    return hex_code\n",
        "\n",
        "def a2c_converter(aa_map: dict, grad_func: 'function') -> dict:\n",
        "    min_value = min(aa_map.values())\n",
        "    max_value = max(aa_map.values())\n",
        "    aa_dict = {aa: grad_func(value, min_value, max_value)\n",
        "               for aa, value in aa_map.items()}\n",
        "    return aa_dict\n",
        "\n",
        "#############################################\n",
        "# Built-in Styling Function\n",
        "\n",
        "def builtin_style(style: str, opacity: float = 1.0) -> dict:\n",
        "    match style:\n",
        "        case _ if any(kw in style for kw in ('Carbon', 'chain', 'ssJmol', 'ssPyMol')):\n",
        "            style_dict = {'colorscheme': style}\n",
        "        case 'hydrophobicity':\n",
        "            style_dict = {'colorscheme': {\n",
        "                'prop': 'resn', 'map': a2c_converter(AA_HB, sequential_gradient)}}\n",
        "        case 'isoelectric points':\n",
        "            style_dict = {'colorscheme': {\n",
        "                'prop': 'resn', 'map': a2c_converter(AA_PI, diverging_gradient)}}\n",
        "        case 'b factor':\n",
        "            style_dict = {'colorscheme': {\n",
        "                'prop': 'b', 'gradient': 'rwb', 'min': 90, 'max': 50}}\n",
        "        case _:\n",
        "            style_dict = {'color': style}\n",
        "\n",
        "    style_dict.update({'opacity': opacity, 'singleBonds': False})\n",
        "    return style_dict\n",
        "\n",
        "#############################################\n",
        "# Built-in Colour Scale Function\n",
        "\n",
        "def colour_scale(aa_map: dict, grad_func: 'function') -> None:\n",
        "    min_value = min(aa_map.values())\n",
        "    max_value = max(aa_map.values())\n",
        "\n",
        "    linear_values = np.linspace(min_value, max_value, 100)\n",
        "    colours = [grad_func(value, min_value, max_value)\n",
        "               for value in linear_values]\n",
        "\n",
        "    fig, ax = plt.subplots(figsize=(4.85, 0.25))\n",
        "    norm_value = plt.Normalize(min_value, max_value)\n",
        "    colour_map = plt.cm.colors.ListedColormap(colours)\n",
        "    scalar_map = plt.cm.ScalarMappable(norm_value, colour_map)\n",
        "    scalar_map.set_array([])\n",
        "\n",
        "    cscale = plt.colorbar(scalar_map, ax, orientation='horizontal')\n",
        "    cscale.set_ticks([min_value, max_value])\n",
        "\n",
        "def show_cscale(rept_info: dict, surf_info: dict) -> None:\n",
        "\n",
        "    def cs_selector() -> str:\n",
        "        if any(surf_info):\n",
        "            style = [*surf_info.values()][0]\n",
        "        elif any(rept_info):\n",
        "            style = [*rept_info.values()][0]\n",
        "        else:\n",
        "            style = None\n",
        "        return style\n",
        "\n",
        "    def cs_display(style: str):\n",
        "        if style == 'hydrophobicity':\n",
        "            label_title(style, 'Less', 'More')\n",
        "            colour_scale(AA_HB, sequential_gradient)\n",
        "        elif style == 'isoelectric points':\n",
        "            label_title(style, 'Acid', 'Base')\n",
        "            colour_scale(AA_PI, diverging_gradient)\n",
        "        else:\n",
        "            pass\n",
        "\n",
        "    def label_title(text: str, min: str, max: str) -> None:\n",
        "        print(f'-' * 55)\n",
        "        print(f'{min}{text.upper():^47}{max}')\n",
        "        print(f'-' * 55)\n",
        "\n",
        "    cs_display(cs_selector())\n",
        "\n",
        "#############################################\n",
        "# Other Functions\n",
        "\n",
        "def extract_config(inpt_file: str) -> tuple:\n",
        "    with open(inpt_file, 'r') as inpt:\n",
        "        data = [line.split() for line in inpt.readlines()]\n",
        "    center = (float(data[0][2]), float(data[1][2]), float(data[2][2]))\n",
        "    bxsize = (float(data[4][2]), float(data[5][2]), float(data[6][2]))\n",
        "    return center, bxsize\n",
        "\n",
        "def interaction_dict(inpt_file: str, interactions: str = '', usage: str = 'view' or 'lbsp') -> dict:\n",
        "\n",
        "    usg_map = {'lbsp': 0, 'view': 1}\n",
        "\n",
        "    def filter_df(int_df: pd.DataFrame, interactions: list = []) -> pd.DataFrame:\n",
        "        int_df = int_df[int_df['BOND'].isin(interactions)] if interactions else int_df\n",
        "        return int_df\n",
        "\n",
        "    def s2f_dict(item: dict) -> dict:\n",
        "        return {key: tuple(float(val) for val in value[1:-1].split(','))\n",
        "                for key, value in item.items()}\n",
        "\n",
        "    def b2c_dict(item: dict) -> dict:\n",
        "        return {key: BOND_COL[val][usg_map[usage]] for key, val in item.items()}\n",
        "\n",
        "    intrxn = interactions.replace(',', ' ').split()\n",
        "    inter_df = pd.read_csv(inpt_file)\n",
        "    int_dict = filter_df(inter_df, intrxn).to_dict()\n",
        "    int_dict['LIGCOO'] = s2f_dict(int_dict['LIGCOO'])\n",
        "    int_dict['PROTCOO'] = s2f_dict(int_dict['PROTCOO'])\n",
        "    int_dict['COLOR'] = b2c_dict(int_dict['BOND'])\n",
        "\n",
        "    return int_dict\n",
        "\n",
        "def find_midpoint(coords: list) -> tuple[float, float, float]:\n",
        "    return tuple(round(coord, 3) for coord in np.mean(coords, axis=0))\n",
        "\n",
        "#############################################\n",
        "# LaboSpace Viewer\n",
        "\n",
        "class LaboSpace:\n",
        "\n",
        "    residue_style = {\n",
        "        'stick':\n",
        "         {'colorscheme': 'orangeCarbon', 'radius': 0.15}}\n",
        "    residue_label = {\n",
        "        'alignment': 'bottomLeft',\n",
        "        'showBackground': False,\n",
        "        'inFront': True,\n",
        "        'fontSize': 14,\n",
        "        'fontColor': '0x000000',\n",
        "        'screenOffset': {'x': 25, 'y': 25}}\n",
        "    atom_label = {\n",
        "        'alignment': 'bottomLeft',\n",
        "        'showBackground': False,\n",
        "        'inFront': True,\n",
        "        'fontSize': 14,\n",
        "        'fontColor': '0x000000',\n",
        "        'screenOffset': {'x': 10, 'y': 10}}\n",
        "\n",
        "    def __init__(self, vw: int = 500, vh: int = 500) -> None:\n",
        "        self.mview = py3Dmol.view(width=vw, height=vh)\n",
        "        self.count = -1\n",
        "        self.residues = []\n",
        "\n",
        "    def read_moldata(self, inpt_file: str) -> str:\n",
        "        inpt = open(inpt_file, 'r')\n",
        "        data = inpt.read()\n",
        "        inpt.close()\n",
        "        return data\n",
        "\n",
        "    def load_receptor(self, inpt_file: str) -> object:\n",
        "        data = self.read_moldata(inpt_file)\n",
        "        self.mview.addModel(data, 'pdb')\n",
        "        self.count += 1\n",
        "        return self\n",
        "\n",
        "    def load_ligand(self, inpt_file: str) -> object:\n",
        "        data = self.read_moldata(inpt_file)\n",
        "        self.mview.addModel(data)\n",
        "        self.count += 1\n",
        "        return self\n",
        "\n",
        "    def set_style(self,\n",
        "                  show_represent: bool = True,\n",
        "                  represent_type: str = 'cartoon',\n",
        "                  represent_style: dict = {}\n",
        "                  ) -> object:\n",
        "        if show_represent:\n",
        "            self.mview.setStyle(\n",
        "                {'model': self.count},\n",
        "                {represent_type: represent_style})\n",
        "        else:\n",
        "            self.mview.setStyle(\n",
        "                {'model': self.count},\n",
        "                {})\n",
        "        return self\n",
        "\n",
        "    def add_style(self,\n",
        "                  show_represent: bool = True,\n",
        "                  represent_style: dict = {}\n",
        "                  ) -> object:\n",
        "        if show_represent:\n",
        "            self.mview.addStyle(\n",
        "                {'model': self.count},\n",
        "                represent_style)\n",
        "        return self\n",
        "\n",
        "    def add_residues(self,\n",
        "                     show_residues: bool = True,\n",
        "                     residue_number: str = ''\n",
        "                     ) -> object:\n",
        "        if show_residues and residue_number:\n",
        "            res = residue_number.replace(',', ' ').split()\n",
        "            self.residues.extend(list(set(res)))\n",
        "            self.mview.addStyle(\n",
        "                {'and': [{'model': self.count}, {'resi': self.residues}]},\n",
        "                self.residue_style)\n",
        "            self.mview.addResLabels(\n",
        "                {'and': [{'model': self.count}, {'resi': self.residues}]},\n",
        "                self.residue_label)\n",
        "        return self\n",
        "\n",
        "    def add_surface(self,\n",
        "                    show_surface: bool = True,\n",
        "                    surface_type: str = 'SES',\n",
        "                    surface_style: dict = {}\n",
        "                    ) -> object:\n",
        "        if show_surface:\n",
        "            self.mview.addSurface(\n",
        "                surface_type,\n",
        "                surface_style,\n",
        "                {'model': self.count})\n",
        "        return self\n",
        "\n",
        "    def add_gridbox(self,\n",
        "                    show_gridbox: bool,\n",
        "                    center: list[float],\n",
        "                    bxsize: list[float]\n",
        "                    ) -> object:\n",
        "        if show_gridbox:\n",
        "            bxi, byi, bzi = center\n",
        "            bxf, byf, bzf = bxsize\n",
        "            self.mview.addBox({\n",
        "                'center': {'x': bxi, 'y': byi, 'z': bzi},\n",
        "                'dimensions': {'w': bxf, 'h': byf, 'd': bzf},\n",
        "                'color': 'skyBlue',\n",
        "                'opacity': 0.6})\n",
        "            self.mview.addLabel(\n",
        "                f'center: {bxi:>8}, {byi:>8}, {bzi:>8}',\n",
        "                {'showBackground': False,\n",
        "                 'fontSize': 14,\n",
        "                 'fontColor': '0x000000',\n",
        "                 'useScreen': True,\n",
        "                 'screenOffset': {'x': 15, 'y': 0}})\n",
        "            self.mview.addLabel(\n",
        "                f'bxsize: {bxf:>8}, {byf:>8}, {bzf:>8}',\n",
        "                {'showBackground': False,\n",
        "                 'fontSize': 14,\n",
        "                 'fontColor': '0x000000',\n",
        "                 'useScreen': True,\n",
        "                 'screenOffset': {'x': 15, 'y': -20}})\n",
        "        return self\n",
        "\n",
        "    def add_interaction(self,\n",
        "                        interaction_file: str,\n",
        "                        show_interaction: bool = True,\n",
        "                        select_interaction: list = []\n",
        "                        ) -> object:\n",
        "        if show_interaction:\n",
        "            int_dict = interaction_dict(interaction_file, select_interaction, 'lbsp')\n",
        "            dist = int_dict['DIST'].values()\n",
        "            bond = int_dict['BOND'].values()\n",
        "            resn = int_dict['RESNR'].values()\n",
        "            ligcoo = int_dict['LIGCOO'].values()\n",
        "            prtcoo = int_dict['PROTCOO'].values()\n",
        "            color = int_dict['COLOR'].values()\n",
        "\n",
        "            int_res = list(set(resn) - set(self.residues))\n",
        "            self.residues.extend(int_res)\n",
        "            self.mview.addStyle(\n",
        "                {'and': [{'model': 0}, {'resi': int_res}]},\n",
        "                self.residue_style)\n",
        "            self.mview.addResLabels(\n",
        "                {'and': [{'model': 0}, {'resi': int_res}]},\n",
        "                self.residue_label)\n",
        "\n",
        "            for dis, col, lig, prt in zip(dist, color, ligcoo, prtcoo):\n",
        "                mid = find_midpoint([lig, prt])\n",
        "                self.mview.addCylinder(\n",
        "                    {'start': {'x': lig[0], 'y': lig[1], 'z': lig[2]},\n",
        "                     'end': {'x': prt[0], 'y': prt[1], 'z': prt[2]},\n",
        "                     'radius': 0.05,\n",
        "                     'fromCap': 1,\n",
        "                     'toCap': 1,\n",
        "                     'color': col,\n",
        "                     'dashed': True})\n",
        "                self.mview.addLabel(\n",
        "                    str(dis) + ' Å',\n",
        "                    {'position': {'x': mid[0], 'y': mid[1], 'z': mid[2]},\n",
        "                     'alignment': 'bottomLeft',\n",
        "                     'inFront': False,\n",
        "                     'backgroundColor': col,\n",
        "                     'fontSize': 10,\n",
        "                     'screenOffset': {'x': 10, 'y': 10}})\n",
        "        return self\n",
        "\n",
        "    def label_atoms(self, show_label: bool = False) -> object:\n",
        "        # WARNING: Avoid applying on protein !!!\n",
        "        if show_label:\n",
        "            self.mview.addPropertyLabels(\n",
        "                'atom',\n",
        "                {'model': self.count},\n",
        "                self.atom_label)\n",
        "        return self\n",
        "\n",
        "    def view_space(self,\n",
        "                   zoom_model: int = -1,\n",
        "                   slab_view: bool = False,\n",
        "                   slab_model: int = -1,\n",
        "                   background_colour: str = '0xFFFFFF'\n",
        "                   ) -> None:\n",
        "        self.mview.setBackgroundColor(background_colour)\n",
        "        self.mview.setProjection('orthographic')\n",
        "        self.mview.zoomTo({'model': zoom_model})\n",
        "        self.mview.fitSlab({'model': slab_model}) if slab_view else None\n",
        "        self.mview.show()\n",
        "\n",
        "print(f'+ Methods and functions successfully built')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F4tmLSywoHH2"
      },
      "source": [
        "---\n",
        "---\n",
        "# **02 | Preparing the Receptor**\n",
        "\n",
        "In this step, we lay the groundwork for molecular docking study by preparing the protein of interest. We will obtain the protein's crystal structure in **`pdb`** format from the [RCSB Protein Data Bank (PDB)](https://www.rcsb.org). Then, we clean the structure by extraction method and seperate its existing subunits into individual files. Finally, we parameterise our target strucutre to generate the **`pdbqt`** file required for docking."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "WIxoT1jKpD0u",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Download protein PDB file**\n",
        "# @markdown Enter the **< Accession ID >** to download the protein of interest.\n",
        "\n",
        "Accession_ID = '7KBS' # @param {type: 'string'}\n",
        "\n",
        "PDB_pdb = Accession_ID.upper() + '.pdb'\n",
        "PDB_pdb_pFile = os.path.join(PRT_FLD, PDB_pdb)\n",
        "Accession_ID_url = 'http://files.rcsb.org/download/' + PDB_pdb\n",
        "\n",
        "!wget {Accession_ID_url} -O {PDB_pdb_pFile} -q\n",
        "print(f'+ PDB downloaded: {PDB_pdb} > PROTEIN folder')\n",
        "print(f'+ RCSB PDB link: https://www.rcsb.org/structure/{Accession_ID}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0qVIuNFpeKFp",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Clean protein structure and extract subunits**\n",
        "# @markdown This generates a **clean** protein structure containing all the subunit(s). \\\n",
        "# @markdown The subunits are then extracted into individual **`pdb`** files. \\\n",
        "# @markdown *Note: `Subunits` and `Chains` are interchangeable to describe the\n",
        "# @markdown multiple polypeptides that made up the protein.*\n",
        "\n",
        "def oupt_parse(inpt_file: str) -> tuple:\n",
        "    name = os.path.basename(inpt_file).split('.')[0]\n",
        "    path = os.path.dirname(inpt_file)\n",
        "    dirn = os.path.basename(path)\n",
        "    return name, dirn, path\n",
        "\n",
        "def ter_line(asn: str, resn: str, chid: str, resi: str) -> str:\n",
        "    return f'TER   {asn}      {resn} {chid}{resi:>4} {\" \"*54}\\n'\n",
        "\n",
        "def extract_entity(inpt_file: str, oupt_file: str, keywords: list) -> None:\n",
        "    oupt_name, oupt_dirn, oupt_path = oupt_parse(oupt_file)\n",
        "    extract = []\n",
        "    headers = ['ATOM', 'HETATM']\n",
        "    with open(inpt_file, 'r') as inpt:\n",
        "        for line in inpt:\n",
        "            record = line.split()[0]\n",
        "            if all(keyword in line for keyword in keywords) and \\\n",
        "               any(header in record for header in headers):\n",
        "                extract.append(line)\n",
        "\n",
        "    assert extract, 'Expected \\'keywords\\' in PDB lines; found 0 lines'\n",
        "\n",
        "    with open(oupt_file, 'w') as oupt:\n",
        "        for line in extract:\n",
        "            oupt.write(line)\n",
        "    print(f'+ Entity extracted: {oupt_name}.pdb > {oupt_dirn} folder')\n",
        "\n",
        "def extract_chains(inpt_file: str) -> None:\n",
        "    inpt_name, oupt_dirn, oupt_path = oupt_parse(inpt_file)\n",
        "    with open(inpt_file, 'r') as inpt:\n",
        "        data = inpt.readlines()\n",
        "\n",
        "    chid = sorted(set(line[21:22] for line in data))\n",
        "    chid_leng = len(chid)\n",
        "    chid_list = ', '.join(chid)\n",
        "    print(f'+ Chains detected: {chid_leng} ({chid_list})')\n",
        "\n",
        "    for id in chid:\n",
        "        oupt_name = inpt_name + '_' + id\n",
        "        oupt_file = os.path.join(oupt_path, oupt_name + '.pdb')\n",
        "        with open(oupt_file, 'w') as oupt:\n",
        "            for line in data:\n",
        "                if line[21:22] in id:\n",
        "                    oupt.write(line)\n",
        "                    asn = f'{int(line[6:11])+1:>5}'\n",
        "                    resn = line[17:20]\n",
        "                    resi = line[22:26]\n",
        "            oupt.write(ter_line(asn, resn, id, resi))\n",
        "            oupt.write('END')\n",
        "        print(f'+ Chain extracted: {oupt_name}.pdb > {oupt_dirn} folder')\n",
        "\n",
        "PDB_prot = Accession_ID + '_prot'\n",
        "PDB_prot_pdb = PDB_prot + '.pdb'\n",
        "PDB_prot_pdb_pFile = os.path.join(PRT_FLD, PDB_prot_pdb)\n",
        "\n",
        "extract_entity(PDB_pdb_pFile, PDB_prot_pdb_pFile, ['ATOM'])\n",
        "extract_chains(PDB_prot_pdb_pFile)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GIHOfCmahoK7",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Visualise 3D protein structure**\n",
        "# @markdown Enter the **< Protein >** to be viewed.\n",
        "\n",
        "View = '7KBS_prot_A' # @param {type: 'string'}\n",
        "Representation_type = 'cartoon' # @param ['cartoon', 'cross', 'line', 'sphere', 'stick']\n",
        "Representation_style = 'white' # @param ['chain', 'white', 'whiteCarbon', 'ssJmol', 'ssPyMol', 'b factor', 'hydrophobicity', 'isoelectric points']\n",
        "Representation_opacity = 1 # @param {type: 'slider', min:0, max:1, step:0.1}\n",
        "Residue_number = '' # @param {type: 'string'}\n",
        "Surface_type = 'SES' # @param ['VDW', 'SAS', 'SES', 'MS']\n",
        "Surface_style = 'hydrophobicity' # @param ['chain', 'white', 'whiteCarbon', 'ssJmol', 'ssPyMol', 'b factor', 'hydrophobicity', 'isoelectric points']\n",
        "Surface_opacity = 1 # @param {type: 'slider', min:0, max:1, step:0.1}\n",
        "Show_representation = True # @param {type: 'boolean'}\n",
        "Show_residue = False # @param {type: 'boolean'}\n",
        "Show_surface = True # @param {type: 'boolean'}\n",
        "\n",
        "PROT_view_pFile = os.path.join(PRT_FLD, View + '.pdb')\n",
        "\n",
        "LBSP = LaboSpace(960, 640)\n",
        "LBSP.load_receptor(PROT_view_pFile)\\\n",
        "    .set_style(\n",
        "        show_represent=Show_representation,\n",
        "        represent_type=Representation_type,\n",
        "        represent_style=builtin_style(\n",
        "            style=Representation_style,\n",
        "            opacity=Representation_opacity))\\\n",
        "    .add_residues(\n",
        "        show_residues=Show_residue,\n",
        "        residue_number=Residue_number)\\\n",
        "    .add_surface(\n",
        "        show_surface=Show_surface,\n",
        "        surface_type=Surface_type,\n",
        "        surface_style=builtin_style(\n",
        "            style=Surface_style,\n",
        "            opacity=Surface_opacity))\n",
        "LBSP.view_space()\n",
        "\n",
        "show_cscale(\n",
        "    {Show_representation: Representation_style},\n",
        "    {Show_surface: Surface_style})"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SDRJdQ7bkB6z",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Parameterise protein with Gasteiger charges**\n",
        "# @markdown Enter the **< Target Protein >** to be parameterised.\\\n",
        "# @markdown This generate a **`protein.pdbqt`** file after:\n",
        "# @markdown + Addition of Gasteiger Partial Charge\n",
        "# @markdown + Addition of polar hydrogens\n",
        "# @markdown + Removal of non-polar hydrogens\n",
        "\n",
        "Target_protein = '7KBS_prot_A' # @param {type: 'string'}\n",
        "\n",
        "PROT_pdb = Target_protein + '.pdb'\n",
        "PROT_pdbqt = Target_protein + '.pdbqt'\n",
        "PROT_pdb_pFile = os.path.join(PRT_FLD, PROT_pdb)\n",
        "PROT_pdb_dFile = os.path.join(DCK_FLD, PROT_pdb)\n",
        "PROT_pdbqt_dFile = os.path.join(DCK_FLD, PROT_pdbqt)\n",
        "\n",
        "!obabel {PROT_pdb_pFile} -xr -O {PROT_pdbqt_dFile} -h --partialcharge gasteiger > /dev/null 2>&1\n",
        "!mk_prepare_receptor.py --pdbqt {PROT_pdbqt_dFile} -o {PROT_pdbqt_dFile} --skip_gpf > /dev/null 2>&1\n",
        "print(f'+ Parameterisation: {PROT_pdb} > {PROT_pdbqt}')\n",
        "\n",
        "shutil.copy(PROT_pdb_pFile, PROT_pdb_dFile)\n",
        "print(f'+ {PROT_pdbqt} > DOCKING folder')\n",
        "print(f'+ {PROT_pdb} > DOCKING folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JdKt0KAL7QRG"
      },
      "source": [
        "---\n",
        "---\n",
        "# **03 | Preparing the Native Ligand** (optional)\n",
        "\n",
        "Next, we retrieve the native ligand from the **`pdb`** file to serve as a reference for later comparison. Noted that **this step is optional** as not all **`pdb`** files contain native ligand. However, we highly recommend proceeding with this step for an insightful visual (3D inspection) and qualtitative (RMSD) comparison of binding pose."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LHAremEDIR7p",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Retrieve native ligand and extract subunits**\n",
        "# @markdown Enter the **< Keyword >** assigned for query \\\n",
        "# @markdown This generates a **clean** native ligand structure containing multiple\n",
        "# @markdown binding poses corresponding to each protein subunit(s). \\\n",
        "# @markdown The binding poses are then extracted and corrected into individual\n",
        "# @markdown **`pdb`** files.\n",
        "\n",
        "Keyword = 'RAL' # @param {type: 'string'}\n",
        "\n",
        "def get_molblock(keyword: str) -> str:\n",
        "    url_path = 'http://files.rcsb.org/ligands/download/' + keyword + '_model.sdf'\n",
        "    sdf_file = os.path.join(NTV_FLD, keyword + '.sdf')\n",
        "    os.system(f'wget {url_path} -O {sdf_file} -q')\n",
        "    molblock = [mol for mol in  Chem.SDMolSupplier(sdf_file) if mol is not None][0]\n",
        "    os.remove(sdf_file)\n",
        "    return molblock\n",
        "\n",
        "def correct_bond_order(inpt_list: list, temp: Chem.rdchem.Mol) -> None:\n",
        "    for inpt_file in inpt_list:\n",
        "        targ = AllChem.MolFromPDBFile(inpt_file)\n",
        "        cmol = AllChem.AssignBondOrdersFromTemplate(temp, targ)\n",
        "        pdbb = Chem.MolToPDBBlock(cmol, flavor=4)\n",
        "        with open(inpt_file, 'w') as oupt:\n",
        "            oupt.write(pdbb)\n",
        "\n",
        "true_id = Keyword[-3:] if len(Keyword) > 3 else Keyword\n",
        "print(f'+ RCSB PDB link: https://www.rcsb.org/ligand/{true_id}')\n",
        "\n",
        "ntv_pdb = Keyword.upper() + '.pdb'\n",
        "ntv_pdb_nFile = os.path.join(NTV_FLD, ntv_pdb)\n",
        "extract_entity(PDB_pdb_pFile, ntv_pdb_nFile, [Keyword, 'HETATM'])\n",
        "extract_chains(ntv_pdb_nFile)\n",
        "\n",
        "ntv_nFiles = sorted(glob.glob(NTV_FLD + '/' + Keyword + '_*.pdb'))\n",
        "ntv_smiles = get_molblock(true_id)\n",
        "correct_bond_order(ntv_nFiles, ntv_smiles)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vqfYDslyFY-N",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Visualise 3D native ligand**\n",
        "# @markdown Enter the **< Native Ligand >** to be viewed.\n",
        "\n",
        "View = 'RAL_A' # @param {type: 'string'}\n",
        "Representation_type = 'stick' # @param ['line', 'sphere', 'stick']\n",
        "Show_atom_labels = False # @param {type: 'boolean'}\n",
        "\n",
        "NTV_view_nFile = os.path.join(NTV_FLD, View + '.pdb')\n",
        "\n",
        "LBSP = LaboSpace(960, 640)\n",
        "LBSP.load_ligand(NTV_view_nFile)\\\n",
        "    .set_style(\n",
        "        show_represent=True,\n",
        "        represent_type=Representation_type,\n",
        "        represent_style={'colorscheme': 'lightGreyCarbon'})\\\n",
        "    .label_atoms(\n",
        "        show_label=Show_atom_labels)\n",
        "LBSP.view_space()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IlAlwF3rUfet",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Select native ligand**\n",
        "# @markdown Enter the **< Target Native Ligand >** for later comparison.\n",
        "\n",
        "Target_native_ligand = 'RAL_A' # @param {type: 'string'}\n",
        "\n",
        "NTV_pdb = Target_native_ligand + '.pdb'\n",
        "NTV_pdb_nFile = os.path.join(NTV_FLD, NTV_pdb)\n",
        "NTV_pdb_dFile = os.path.join(DCK_FLD, NTV_pdb)\n",
        "\n",
        "shutil.copy(NTV_pdb_nFile, NTV_pdb_dFile)\n",
        "print(f'+ {NTV_pdb} > DOCKING folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_tluzw-AJrEa"
      },
      "source": [
        "---\n",
        "---\n",
        "# **04 | Preparing the Ligand**\n",
        "\n",
        "Now, we start to prepare our ligand of interest. We will be using its SMILES notation to generate a 3D energetically minimised ligand **`sd`** file based on specifed force field and steepest gradient alogrithm. We also parameterise our target ligand to generate the **`pdbqt`** file required for docking."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "T2jsTHBFbY1w",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Provide ligand information**\n",
        "# @markdown Enter the **< ID >** and its **< SMILES >** notation.\n",
        "\n",
        "ID = 'RLXF' # @param {type: 'string'}\n",
        "SMILES = 'C1CCN(CC1)CCOC2=CC=C(C=C2)C(=O)C3=C(SC4=C3C=CC(=C4)O)C5=CC=C(C=C5)O' # @param {type: 'string'}\n",
        "\n",
        "assert ID.strip(), 'Expected \\'ID\\'; found \\'none\\''\n",
        "assert SMILES.strip(), 'Expected \\'SMILES\\'; found \\'none\\''\n",
        "\n",
        "print(f'+ Ligand \\'{ID}\\' with \\'{SMILES}\\' ready for SDF generation')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8623P-E6J1Fg",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Generate ligand SD file**\n",
        "# @markdown Select **< Force Field >** to set MM method for energy minimisation. \\\n",
        "# @markdown Select **< Convergence Criteria >** to set minimum iteration difference required before stopping. **[ Default is 0.00001 ]** \\\n",
        "# @markdown Select **< Maximum Steps >** to set maximum iterations allowed for minimisation. **[ Default is 10000 ]**\n",
        "\n",
        "Force_field = 'UFF' # @param ['GAFF', 'Ghemical', 'MMFF94', 'MMFF94s', 'UFF']\n",
        "Convergence_criteria = '0.00001' # @param ['0.1', '0.01','0.001', '0.0001', '0.00001', '0.000001', '0.0000001']\n",
        "Maximum_steps = 10000 # @param {type: 'slider', min:1000, max:100000, step:1000}\n",
        "\n",
        "def check_convergence(log_file: str) -> None:\n",
        "    with open(log_file, 'r') as inpt:\n",
        "        data = inpt.read()\n",
        "    if 'CONVERGED' in data:\n",
        "        verb = 'has'\n",
        "        step = data.split('\\n')[-4].split()[0]\n",
        "        step_str = f'({step} steps)'\n",
        "    else:\n",
        "        verb = 'has not'\n",
        "        step_str = ''\n",
        "    print(f'+ Steepest gradient {verb} converged {step_str}')\n",
        "\n",
        "LIG_ID = ID\n",
        "LIG_sdf = LIG_ID + '.sdf'\n",
        "LIG_obmin_log = LIG_ID + '_obmin.log'\n",
        "LIG_sdf_lFile = os.path.join(LIG_FLD, LIG_sdf)\n",
        "LIG_obmin_log_lFile = os.path.join(LIG_FLD, LIG_obmin_log)\n",
        "\n",
        "print(f'+ Selected {Force_field} for energy minimisation '\n",
        "      f'up to {Convergence_criteria} iteration difference or '\n",
        "      f'at most {Maximum_steps:,} steps')\n",
        "\n",
        "!obabel -:{'\"'+SMILES+'\"'} -O {LIG_sdf_lFile} --title {LIG_ID} --gen3d \\\n",
        "--best --minimize --ff {Force_field} --steps {Maximum_steps} --sd \\\n",
        "--crit {Convergence_criteria} --log &> {LIG_obmin_log_lFile}\n",
        "\n",
        "check_convergence(LIG_obmin_log_lFile)\n",
        "print(f'+ {LIG_sdf} > LIGAND folder')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UtBrkZ-Pg1hQ",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Visualise 3D ligand**\n",
        "# @markdown Enter the **< Ligand >** to be viewed.\n",
        "\n",
        "View = 'RLXF' # @param {type: 'string'}\n",
        "Representation_type = 'stick' # @param ['line', 'sphere', 'stick']\n",
        "Show_atom_labels = False # @param {type: 'boolean'}\n",
        "\n",
        "LIG_view_lFile = os.path.join(LIG_FLD, View + '.sdf')\n",
        "\n",
        "LBSP = LaboSpace(960, 640)\n",
        "LBSP.load_ligand(LIG_view_lFile)\\\n",
        "    .set_style(\n",
        "        show_represent=True,\n",
        "        represent_type=Representation_type,\n",
        "        represent_style={'colorScheme': 'lightGreyCarbon'})\\\n",
        "    .label_atoms(\n",
        "        show_label=Show_atom_labels)\n",
        "LBSP.view_space()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Qgp6BDDymQtX",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Parameterise ligand with Gasteiger charges**\n",
        "# @markdown This generates a **`ligand.pdbqt`** file after:\n",
        "# @markdown + Addition of Gasteiger Partial Charge\n",
        "# @markdown + Removal of non-polar hydrogens\n",
        "\n",
        "def pdbqt_add_chid(inpt_file: str) -> None:\n",
        "    with open(inpt_file, 'r') as inpt:\n",
        "        data = inpt.read()\n",
        "    new_data = data.replace('  UNL  ', '  UNL A')\n",
        "    with open(inpt_file, 'w') as oupt:\n",
        "        oupt.write(new_data)\n",
        "\n",
        "LIG_pdbqt = ID + '.pdbqt'\n",
        "LIG_dFLD = os.path.join(DCK_FLD, ID)\n",
        "LIG_pdbqt_dFFile = os.path.join(LIG_dFLD, LIG_pdbqt)\n",
        "\n",
        "os.makedirs(LIG_dFLD, exist_ok=True)\n",
        "!mk_prepare_ligand.py -i {LIG_sdf_lFile} -o {LIG_pdbqt_dFFile} > /dev/null 2>&1\n",
        "pdbqt_add_chid(LIG_pdbqt_dFFile)\n",
        "\n",
        "print(f'+ Parameterisation: {LIG_sdf} > {LIG_pdbqt}')\n",
        "print(f'+ {LIG_pdbqt} > DOCKING folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LMf2ZMMIHrxq"
      },
      "source": [
        "---\n",
        "---\n",
        "# **05 | Setting Up Grid Box**\n",
        "\n",
        "Here, we define a chemical search space with the use of grid box, which often centered within the binding, active or allosteric site of the target protein."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sZ8FYg49IDQR"
      },
      "source": [
        "Built-in methods for defining grid box include:\n",
        "\n",
        "+ **`LaBOX`** : Use the mean of XYZ extreme values from native ligand. [Link](https://github.com/RyanZR/LaBOX)\n",
        "+ **`eBoxSize`** : Use the radius of gyration of native ligand. [Link](https://github.com/michal-brylinski/eboxsize)\n",
        "+ **`eBoxSize-Mod`** : Similar to eBoxSize, but box center is computed using LaBOX method.\n",
        "+ **`Autodock-Grid`** : 22.5 × 22.5 × 22.5 Å\n",
        "+ **`Defined-by-Res`** : Use when important residues involved in binding interaction are known.\n",
        "+ **`Manual-Mode`** : Use the sliders below to define grid box."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cQooGBRdHrez",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Place grid box at binding site**\n",
        "# @markdown Select **< Method >** to set grid box at the center of binding site.\n",
        "\n",
        "Residue_number = '' # @param {type: 'string'}\n",
        "Method = 'LaBOX' # @param ['LaBOX', 'eBoxSize', 'eBoxSize-Mod', 'Autodock-Grid', 'Defined-by-Res', 'Manual-Mode']\n",
        "Show_residues = True # @param {type: 'boolean'}\n",
        "\n",
        "# @markdown **Manual Mode**\n",
        "X = 0 # @param {type: 'number'}\n",
        "Y = 0 # @param {type: 'number'}\n",
        "Z = 0 # @param {type: 'number'}\n",
        "Width = 10 # @param {type: 'number'}\n",
        "Height = 10 # @param {type: 'number'}\n",
        "Depth = 10 # @param {type: 'number'}\n",
        "Scale = 1.0 # @param {type: 'number'}\n",
        "\n",
        "LBSP = LaboSpace(960, 640)\n",
        "LBSP_view = True\n",
        "\n",
        "try:\n",
        "    if Method == 'Defined-by-Res':\n",
        "        GB = GridBox(PROT_pdb_dFile)\n",
        "        center, bxsize = GB.defined_by_res(Residue_number)\n",
        "    if Method == 'Manual-Mode':\n",
        "        center = (X, Y, Z)\n",
        "        bxsize = (round(Width, 3),\n",
        "                  round(Height, 3),\n",
        "                  round(Depth, 3))\n",
        "except Exception as excp:\n",
        "    print(f'+ {excp}')\n",
        "    print(f'+ Error can be:\\n'\n",
        "           '  - Expected \\'Residue_number\\'; found \\'none\\'\\n'\n",
        "           '  - Expected \\'protein PDB file\\'; found \\'none\\'\\n'\n",
        "           '  - Due to atomic error in \\'protein PDB file\\'')\n",
        "    LBSP_view = False\n",
        "else:\n",
        "    LBSP.load_receptor(PROT_pdb_dFile)\\\n",
        "        .set_style(\n",
        "            show_represent=True,\n",
        "            represent_type='cartoon',\n",
        "            represent_style={'color': 'white'})\\\n",
        "        .add_residues(\n",
        "            show_residues=Show_residues,\n",
        "            residue_number=Residue_number)\n",
        "\n",
        "try:\n",
        "    if Method == 'LaBOX':\n",
        "        GB = GridBox(NTV_pdb_nFile)\n",
        "        center, bxsize = GB.labox()\n",
        "    if Method == 'eBoxSize':\n",
        "        GB = GridBox(NTV_pdb_nFile)\n",
        "        center, bxsize = GB.eboxsize()\n",
        "    if Method == 'eBoxSize-Mod':\n",
        "        GB = GridBox(NTV_pdb_nFile)\n",
        "        center, bxsize = GB.eboxsize(modified=True)\n",
        "    if Method == 'Autodock-Grid':\n",
        "        GB = GridBox(NTV_pdb_nFile)\n",
        "        center, bxsize = GB.autodock_grid()\n",
        "except Exception as excp:\n",
        "    print(f'+ {excp}')\n",
        "    print(f'+ Expected \\'native ligand PDB file\\'; found \\'none\\'')\n",
        "    print(f'+ Upload file to \\'{NTV_FLD}\\' and run \\'Select native ligand\\' cell to load it.')\n",
        "    LBSP_view = False\n",
        "else:\n",
        "    try:\n",
        "        LBSP.load_ligand(NTV_pdb_nFile)\\\n",
        "            .set_style(\n",
        "                show_represent=True,\n",
        "                represent_type='stick',\n",
        "                represent_style={'colorScheme': 'greyCarbon'})\n",
        "    except:\n",
        "        pass\n",
        "\n",
        "if LBSP_view:\n",
        "    LBSP.add_gridbox(show_gridbox=True, center=center, bxsize=bxsize)\n",
        "    LBSP.view_space()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zD34rhfQrR4t",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Generate docking config file**\n",
        "# @markdown This generates **`config_file`** for AutoDock Vina.\n",
        "\n",
        "cfg_name = 'config_file'\n",
        "cfg_dFile = os.path.join(DCK_FLD, cfg_name)\n",
        "\n",
        "with open(cfg_dFile, 'w') as cfg:\n",
        "    cfg.write(f'center_x = {center[0]}\\n')\n",
        "    cfg.write(f'center_y = {center[1]}\\n')\n",
        "    cfg.write(f'center_z = {center[2]}\\n')\n",
        "    cfg.write(f'\\n')\n",
        "    cfg.write(f'size_x = {bxsize[0]}\\n')\n",
        "    cfg.write(f'size_y = {bxsize[1]}\\n')\n",
        "    cfg.write(f'size_z = {bxsize[2]}\\n')\n",
        "\n",
        "print(f'+ {cfg_name} > DOCKING folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oviol2hQtGZF"
      },
      "source": [
        "---\n",
        "---\n",
        "# **06 | Performing Molecular Docking**\n",
        "\n",
        "We now proceed with performing molecular docking using Autodock Vina for our target ligand. The duration can vary, mainly influenced by factors including **number of rotatable bonds** and **level of exhaustiveness** selected. Do be patient as the process may take a few minutes to complete."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "prePleU-tEeK",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Run AutoDock Vina**\n",
        "# @markdown Select **< Exhaustiveness >** for molecular docking. **[ Default is 16 ]** \\\n",
        "# @markdown This generates a **`ligand_output.pdbqt`** containing different binding poses predicted.\n",
        "\n",
        "Exhaustiveness = '16' # @param [4, 8, 16, 32, 64, 128, 256]\n",
        "\n",
        "cpu_cores = os.cpu_count()\n",
        "LIG_oupt_log = LIG_ID + '_output.log'\n",
        "LIG_oupt_pdbqt = LIG_ID + '_output.pdbqt'\n",
        "LIG_oupt_log_dFFile = os.path.join(LIG_dFLD, LIG_oupt_log)\n",
        "LIG_oupt_pdbqt_dFFile = os.path.join(LIG_dFLD, LIG_oupt_pdbqt)\n",
        "\n",
        "# -- Start docking --\n",
        "start = time.time()\n",
        "%vina --receptor {PROT_pdbqt_dFile} --ligand {LIG_pdbqt_dFFile} \\\n",
        "--out {LIG_oupt_pdbqt_dFFile} --config {cfg_dFile} --cpu {cpu_cores} \\\n",
        "--exhaustiveness {Exhaustiveness} --verbosity 2 | tee {LIG_oupt_log_dFFile}\n",
        "end = time.time()\n",
        "# -- End docking --\n",
        "\n",
        "print(f'')\n",
        "print(f'+ Docking completed')\n",
        "print(f'+ Time elapsed: ' + time.strftime('%Mm %Ss', time.gmtime(end - start)))\n",
        "print(f'+ {LIG_oupt_pdbqt} > DOCKING folder')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZYwqvY5zFR9u",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Process output file**\n",
        "# @markdown This generates **`ligand_[n].pdb`** for each predicted binding poses.\n",
        "\n",
        "LIG_dash_pdb = LIG_ID + '_.pdb'\n",
        "LIG_dash_pdb_dFFile = os.path.join(LIG_dFLD, LIG_dash_pdb)\n",
        "\n",
        "with open(LIG_oupt_pdbqt_dFFile, 'r') as oupt:\n",
        "    data = oupt.read()\n",
        "    pose = data.count('MODEL')\n",
        "\n",
        "!obabel {LIG_oupt_pdbqt_dFFile} -O {LIG_dash_pdb_dFFile} -m > /dev/null 2>&1\n",
        "print(f'+ {pose} {LIG_ID}_[n].pdb > DOCKING folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RRCearBPB_49"
      },
      "source": [
        "---\n",
        "---\n",
        "# **07 | Profiling Binding Interactions**\n",
        "\n",
        "After docking, we use PLIP to determine the non-covalent binding interactions between the docked ligand and target protein. We will prepare the required complex **`pdb`** file containing docked ligand pose and target protein for PLIP to profile interactions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TjZGfa8TlTLK",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Generate complex PDB file**\n",
        "# @markdown This merges the target protein with each docked ligand poses to\n",
        "# @markdown produce individual **`cmpx.pdb`** file.\n",
        "\n",
        "def generate_cmpx_pdb(inpt_prot: str, inpt_pose: str, oupt_cmpx: str) -> None:\n",
        "\n",
        "    def write_line(line: str, keywords: list, oupt_file: str) -> None:\n",
        "        header = line.split()[0]\n",
        "        if header in keywords:\n",
        "            oupt_file.write(line)\n",
        "\n",
        "    def cmpx_writer() -> None:\n",
        "        with open(oupt_cmpx, 'w') as oupt_file, \\\n",
        "             open(inpt_prot, 'r') as prot_file, \\\n",
        "             open(inpt_pose, 'r') as pose_file:\n",
        "            for prot_line in prot_file:\n",
        "                write_line(prot_line, prot_headers, oupt_file)\n",
        "            for pose_line in pose_file:\n",
        "                write_line(pose_line, pose_headers, oupt_file)\n",
        "\n",
        "    prot_headers = ['ATOM', 'CONECT', 'TER']\n",
        "    pose_headers = ['ATOM', 'CONECT', 'END']\n",
        "    cmpx_writer()\n",
        "\n",
        "LIG_iFLD = os.path.join(INT_FLD, LIG_ID)\n",
        "os.makedirs(LIG_iFLD, exist_ok=True)\n",
        "LIG_pdb_dFFiles = sorted(glob.glob(LIG_dFLD + '/*.pdb'))\n",
        "\n",
        "count = 0\n",
        "for pose_file in tqdm(LIG_pdb_dFFiles):\n",
        "    pose_name = os.path.basename(pose_file).split('.')[0]\n",
        "    cmpx_pdb = pose_name + '_cmpx.pdb'\n",
        "    cmpx_pdb_iFile = os.path.join(LIG_iFLD, cmpx_pdb)\n",
        "    generate_cmpx_pdb(PROT_pdb_dFile, pose_file, cmpx_pdb_iFile)\n",
        "    count += 1\n",
        "\n",
        "print(f'+ {count} {LIG_ID}_[n]_cmpx.pdb > INTERACTION folder')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "__4AV4OM0p79",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Run PLIP**\n",
        "# @markdown This generates an **`interpt.csv`** containing all detected binding\n",
        "# @markdown interaction for each individual **`cmpx.pdb`** file. \\\n",
        "# @markdown Supported: **`HYDROPHOBIC`**, **`HBOND`**, **`WATERBRIDGE`**,\n",
        "# @markdown **`SALTBRIDGE`**, **`PISTACKING`**, **`PICATION`**, **`HALOGEN`**, **`METAL`**. \\\n",
        "\n",
        "def interaction_profiler(inpt_cmpx: str, oupt_csv: str) -> None:\n",
        "\n",
        "    int_bonds = ['HYDROPHOBIC', 'HBOND', 'WATERBRIDGE', 'SALTBRIDGE',\n",
        "                 'PISTACKING', 'PICATION', 'HALOGEN', 'METAL']\n",
        "\n",
        "    def BSR(inpt_cmpx: str) -> object:\n",
        "        cmpx_mol = PDBComplex()\n",
        "        cmpx_mol.load_pdb(inpt_cmpx)\n",
        "        cmpx_lig = [lig for lig in cmpx_mol.ligands if lig.hetid == 'UNL'][0]\n",
        "        cmpx_mol.characterize_complex(cmpx_lig)\n",
        "        cmpx_int = cmpx_mol.interaction_sets['UNL:A:1']\n",
        "        cmpx_rpt = BindingSiteReport(cmpx_int)\n",
        "        return cmpx_rpt\n",
        "\n",
        "    def BSR_feat(bsr: object) -> list:\n",
        "        return [list(getattr(bsr, bond.lower() + '_features')) + ['BOND']\n",
        "                for bond in int_bonds]\n",
        "\n",
        "    def BSR_info(bsr: object) -> list:\n",
        "        return [[list(info) + [bond] for info in getattr(bsr, bond.lower() + '_info')]\n",
        "                for bond in int_bonds]\n",
        "\n",
        "    def replace_column(df: pd.DataFrame, col_A: str, cols: list) -> pd.DataFrame:\n",
        "        for col in cols:\n",
        "            if col in df.columns:\n",
        "                df[col_A] = df[col_A].fillna(df[col])\n",
        "            else:\n",
        "                pass\n",
        "        return df\n",
        "\n",
        "    def BSR_dataframe(bsr_feat: list, bsr_info: list) -> pd.DataFrame:\n",
        "        bsr_df = []\n",
        "        for feat, info in zip(bsr_feat, bsr_info):\n",
        "            if info:\n",
        "                df = pd.DataFrame(info, columns=feat)\n",
        "            else:\n",
        "                df = pd.DataFrame(columns=bsr_feat[0])\n",
        "            bsr_df.append(df)\n",
        "        BSR_df = pd.concat(bsr_df, ignore_index=True)\n",
        "        BSR_df = replace_column(BSR_df, 'DIST', ['DIST_D-A', 'CENTDIST'])\n",
        "        return BSR_df\n",
        "\n",
        "    cmpx_bsr = BSR(inpt_cmpx)\n",
        "    bsr_feat = BSR_feat(cmpx_bsr)\n",
        "    bsr_info = BSR_info(cmpx_bsr)\n",
        "    bsr_data = BSR_dataframe(bsr_feat, bsr_info)\n",
        "    bsr_data.to_csv(oupt_csv, index=False)\n",
        "\n",
        "cmpx_pdb_iFFiles = sorted(glob.glob(LIG_iFLD + '/*_cmpx.pdb'))\n",
        "\n",
        "count = 0\n",
        "for cmpx_file in tqdm(cmpx_pdb_iFFiles):\n",
        "    cmpx_name = os.path.basename(cmpx_file).split('.')[0][:-5]\n",
        "    inter_csv = cmpx_name + '_interpt.csv'\n",
        "    inter_csv_iFile = os.path.join(LIG_iFLD, inter_csv)\n",
        "    interaction_profiler(cmpx_file, inter_csv_iFile)\n",
        "    count += 1\n",
        "\n",
        "print(f'+ {count} {LIG_ID}_[n]_interpt.csv > INTERACTION folder')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xIfJWFo30Xf6"
      },
      "source": [
        "---\n",
        "---\n",
        "# **08 | Analysing Docking Results**\n",
        "\n",
        "We process all the docking and binding internaction results generated into for further evaluation. We will compute the RMSD between native ligand and docked ligands using three different calculation methods. We will also visualise the docking pose in three-dimensional chemical space."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "V8k-pIWDvBY2",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Generate docking report**\n",
        "# @markdown This clusters the docking score from the docked ligand pose **`pdb`**\n",
        "# @markdown file and computes the RMSD with different 3 methods. \\\n",
        "# @markdown The report is then saved into a **`ligand_dockrpt.csv`** file.\n",
        "\n",
        "def vina_report(inpt_file: str) -> dict:\n",
        "\n",
        "    def vn_scores(inpt_file: str) -> float:\n",
        "        with open(inpt_file, 'r') as inpt_data:\n",
        "            for line in inpt_data:\n",
        "                if 'REMARK VINA RESULT' in line:\n",
        "                    dock_sc = float(line.split()[3])\n",
        "                    rmsd_lb = float(line.split()[4])\n",
        "                    rmsd_ub = float(line.split()[5])\n",
        "        return dock_sc, rmsd_lb, rmsd_ub\n",
        "\n",
        "    def vn_report(name: str, scores: tuple) -> dict:\n",
        "        return {'NAME': [name], 'DOCK_SC': [scores[0]],\n",
        "                'RMSD_LB': [scores[1]], 'RMSD_UB': [scores[2]]}\n",
        "\n",
        "    inpt_name = os.path.basename(inpt_file).split('.')[0]\n",
        "    scores = vn_scores(inpt_file)\n",
        "    report = vn_report(inpt_name, scores)\n",
        "    return report\n",
        "\n",
        "LIG_MCS_png = LIG_ID + '_MCS.png'\n",
        "LIG_dockrpt_csv = LIG_ID + '_dockrpt.csv'\n",
        "LIG_MCS_png_dFFile = os.path.join(LIG_dFLD, LIG_MCS_png)\n",
        "LIG_dockrpt_csv_dFFile = os.path.join(LIG_dFLD, LIG_dockrpt_csv)\n",
        "\n",
        "RMSD = ComputeRMSD()\n",
        "dock_rpt = pd.DataFrame()\n",
        "\n",
        "try:\n",
        "    ref_info = RMSD.load_molecule(NTV_pdb_dFile)\n",
        "except:\n",
        "    print('+ Expected \\'native ligand PDB file\\' for RMSD calculation; found \\'none\\'.')\n",
        "    print('+ Docking report generation proceeded without RMSD calculation.')\n",
        "    for lig_file in LIG_pdb_dFFiles:\n",
        "        vina_rpt = pd.DataFrame.from_dict(vina_report(lig_file))\n",
        "        dock_rpt = pd.concat([dock_rpt, vina_rpt], ignore_index=True)\n",
        "else:\n",
        "    for lig_file in LIG_pdb_dFFiles:\n",
        "        lig_info = RMSD.load_molecule(lig_file)\n",
        "        vina_rpt = pd.DataFrame.from_dict(vina_report(lig_file))\n",
        "        rmsd_rpt = pd.DataFrame.from_dict(RMSD.rmsd_report(ref_info, lig_info)).drop(['NAME'], axis=1)\n",
        "        pose_rpt = pd.concat([vina_rpt, rmsd_rpt], axis=1)\n",
        "        dock_rpt = pd.concat([dock_rpt, pose_rpt], ignore_index=True)\n",
        "    RMSD.MCS_png.save(LIG_MCS_png_dFFile)\n",
        "\n",
        "dock_rpt.to_csv(LIG_dockrpt_csv_dFFile, index=False)\n",
        "print(f'+ {LIG_ID}_dockrpt.csv > DOCKING folder')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7gD-K72KsLY6",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Show docking report**\n",
        "# @markdown Enter the **< Ligand Name >** to view docking score and RMSD.\n",
        "# @markdown + **`DOCK_SC`** : Autodock Vina Docking Score\n",
        "# @markdown + **`RMSD_LB`** : RMSD Lower Bound (vs Best Pose)\n",
        "# @markdown + **`RMSD_UB`** : RMSD Upper Bound (vs Best Pose)\n",
        "# @markdown + **`LABO_RMSD`** : RMSD by Common Substructure (vs Native Ligand)\n",
        "# @markdown + **`HUNG_RMSD`** : RMSD by Hungarian Algorithm (vs Native Ligand)\n",
        "# @markdown + **`SYMM_RMSD`** : Symmetrical-corrected RMSD (vs Native Ligand)\n",
        "\n",
        "Ligand_name = 'RLXF' # @param {type : 'string'}\n",
        "\n",
        "LIG_dockrpt_csv = Ligand_name + '_dockrpt.csv'\n",
        "LIG_dockrpt_csv_dFile = os.path.join(LIG_dFLD, LIG_dockrpt_csv)\n",
        "LIG_dockrpt_df = pd.read_csv(LIG_dockrpt_csv_dFile)\n",
        "LIG_dockrpt_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "f4wYT6u-kDZN",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Show interaction profile**\n",
        "# @markdown Enter the **< Pose Name >** to view the binding interactions.\n",
        "\n",
        "Pose_name = 'RLXF_1' # @param {type : 'string'}\n",
        "\n",
        "LIG_interpt_csv = Pose_name + '_interpt.csv'\n",
        "LIG_interpt_csv_iFile = os.path.join(LIG_iFLD, LIG_interpt_csv)\n",
        "LIG_interpt_dict = interaction_dict(LIG_interpt_csv_iFile, usage='view')\n",
        "LIG_interpt_df = pd.DataFrame.from_dict(LIG_interpt_dict)\n",
        "LIG_interpt_df = LIG_interpt_df[['RESTYPE', 'RESNR', 'DIST', 'BOND', 'COLOR']]\n",
        "LIG_interpt_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vjkiHT29vDMs",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Visualise 3D docking pose**\n",
        "# @markdown **PROTEIN MODEL** \\\n",
        "# @markdown Enter the **< Protein >** to be viewed.\n",
        "\n",
        "View_protein = '7KBS_prot_A' # @param {type: 'string'}\n",
        "Protein_type = 'cartoon' # @param ['cartoon', 'cross', 'line', 'sphere', 'stick']\n",
        "Protein_style = 'white' # @param ['chain', 'white', 'whiteCarbon', 'ssJmol', 'ssPyMol', 'b factor', 'hydrophobicity', 'isoelectric points']\n",
        "Protein_opacity = 0.5 # @param {type: 'slider', min:0, max:1, step:0.1}\n",
        "Residue_number = '' # @param {type: 'string'}\n",
        "Surface_type = 'SES' # @param ['VDW', 'SAS', 'SES', 'MS']\n",
        "Surface_style = 'isoelectric points' # @param ['chain', 'white', 'whiteCarbon', 'ssJmol', 'ssPyMol', 'b factor', 'hydrophobicity', 'isoelectric points']\n",
        "Surface_opacity = 0.5 # @param {type: 'slider', min:0, max:1, step:0.1}\n",
        "Show_protein = True # @param {type: 'boolean'}\n",
        "Show_residue = False # @param {type: 'boolean'}\n",
        "Show_surface = False # @param {type: 'boolean'}\n",
        "Show_gridbox = False # @param {type: 'boolean'}\n",
        "\n",
        "# @markdown **LIGAND MODEL** \\\n",
        "# @markdown Enter the **< Ligand >** to be viewed.\n",
        "\n",
        "View_native_ligand = 'RAL_A' # @param {type: 'string'}\n",
        "Native_ligand_style = 'stick' # @param ['cross', 'line', 'sphere', 'stick']\n",
        "Show_native_ligand = True # @param {type: 'boolean'}\n",
        "View_docked_ligand = 'RLXF_1' # @param {type: 'string'}\n",
        "Docked_ligand_style = 'stick' # @param ['cross', 'line', 'sphere', 'stick']\n",
        "Show_docked_ligand = True # @param {type: 'boolean'}\n",
        "\n",
        "# @markdown **INTERACTION MODEL** \\\n",
        "# @markdown Enter the **< Interaction Type >** to be viewed. \\\n",
        "# @markdown Select or combine from **`HYDROPHOBIC`**, **`HBOND`**, **`WATERBRIDGE`**,\n",
        "# @markdown **`SALTBRIDGE`**, **`PISTACKING`**, **`PICATION`**, **`HALOGEN`**, **`METAL`**. \\\n",
        "# @markdown *Note: All interactions are selected if not provided.*\n",
        "\n",
        "Interaction_type = \"HBOND, WATERBRIDGE, SALTBRIDGE, PISTACKING, PICATION, HALOGEN, METAL\" # @param {type: 'string'}\n",
        "Show_interaction = True # @param {type: 'boolean'}\n",
        "\n",
        "# @markdown **OTHER OPTIONS** \\\n",
        "# @markdown Miscellaneous visualisation settings.\n",
        "\n",
        "Slab_view = False # @param {type: 'boolean'}\n",
        "\n",
        "dLIG_dFLD = os.path.join(DCK_FLD, View_docked_ligand[:-2])\n",
        "iLIG_iFLD = os.path.join(INT_FLD, View_docked_ligand[:-2])\n",
        "PROT_view_dFile = os.path.join(DCK_FLD, View_protein + '.pdb')\n",
        "NTV_view_dFile = os.path.join(DCK_FLD, View_native_ligand + '.pdb')\n",
        "LIG_view_dFFile = os.path.join(dLIG_dFLD, View_docked_ligand + '.pdb')\n",
        "INT_view_iFFile = os.path.join(iLIG_iFLD, View_docked_ligand + '_interpt.csv')\n",
        "cfg_center, cfg_bxsize = extract_config(cfg_dFile)\n",
        "\n",
        "LBSP = LaboSpace(960, 640)\n",
        "LBSP.load_receptor(PROT_view_dFile)\\\n",
        "    .set_style(\n",
        "        show_represent=Show_protein,\n",
        "        represent_type=Protein_type,\n",
        "        represent_style=builtin_style(\n",
        "            style=Protein_style,\n",
        "            opacity=Protein_opacity))\\\n",
        "    .add_residues(\n",
        "        show_residues=Show_residue,\n",
        "        residue_number=Residue_number)\\\n",
        "    .add_surface(\n",
        "        show_surface=Show_surface,\n",
        "        surface_type=Surface_type,\n",
        "        surface_style=builtin_style(\n",
        "            style=Surface_style,\n",
        "            opacity=Surface_opacity))\n",
        "\n",
        "try:\n",
        "    LBSP.load_ligand(NTV_view_dFile)\\\n",
        "        .set_style(\n",
        "            show_represent=Show_native_ligand,\n",
        "            represent_type=Native_ligand_style,\n",
        "            represent_style={'color': 'grey'})\n",
        "except:\n",
        "    pass\n",
        "\n",
        "LBSP.load_ligand(LIG_view_dFFile)\\\n",
        "    .set_style(\n",
        "        show_represent=Show_docked_ligand,\n",
        "        represent_type=Docked_ligand_style,\n",
        "        represent_style={'colorscheme': 'salmonCarbon'})\n",
        "LBSP.add_interaction(\n",
        "    interaction_file=INT_view_iFFile,\n",
        "    show_interaction=Show_interaction,\n",
        "    select_interaction=Interaction_type)\n",
        "LBSP.add_gridbox(\n",
        "    show_gridbox=Show_gridbox,\n",
        "    center=cfg_center,\n",
        "    bxsize=cfg_bxsize)\n",
        "LBSP.view_space(\n",
        "    slab_view=Slab_view)\n",
        "\n",
        "show_cscale({Show_protein: Protein_style}, {Show_surface: Surface_style})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FpDiF-oz-I0h"
      },
      "source": [
        "---\n",
        "---\n",
        "# **09 | Save to Google Drive**\n",
        "\n",
        "Finally, we save all the generated data into Google Drive. A copy of the current working directory will be stored in our specified destination in Google Drive."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qQGwh9nh-IWQ",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Mount Google Drive**\n",
        "# @markdown This flush and mount Google Drive.\n",
        "\n",
        "drive.flush_and_unmount()\n",
        "drive.mount('/content/drive', force_remount=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1pvLVkGS-iDh",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title **Store data**\n",
        "# @markdown Enter the **< Destination >** for storage. \\\n",
        "# @markdown This copies all generated data from current working directory into\n",
        "# @markdown specified destination.\n",
        "\n",
        "Destination = '/content/drive/MyDrive/Docking' # @param {type: 'string'}\n",
        "\n",
        "DST_FLD = os.path.join(Destination, Job_name)\n",
        "shutil.copytree(WRK_DIR, DST_FLD)\n",
        "\n",
        "print(f'+ {Job_name} stored in {DST_FLD}')"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "FDIu2L1jfx5_",
        "F4tmLSywoHH2",
        "JdKt0KAL7QRG",
        "_tluzw-AJrEa",
        "LMf2ZMMIHrxq",
        "oviol2hQtGZF",
        "RRCearBPB_49",
        "xIfJWFo30Xf6",
        "FpDiF-oz-I0h"
      ],
      "provenance": [],
      "authorship_tag": "ABX9TyM9L5AclANZ9fYbngyyVxSe"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}