Free_Cloud_Docking_2D_3D.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "py35-paddle1.2.0"
    },
    "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.7.4"
    },
    "toc-autonumbering": true,
    "toc-showmarkdowntxt": false,
    "colab": {
      "name": "Free_Cloud_Docking.ipynb",
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/quantaosun/pl3/blob/main/Free_Cloud_Docking_2D_3D.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-VAb-v6K0EOy"
      },
      "source": [
        "# Programmatic Molecular Docking \n",
        "\n",
        "### Based on Vina-derived Smina. Obabel is used to calculate the protein charges. The docking grid is defined as the centre of the native ligand.\n",
        "\n",
        "- This notebook is written by https://github.com/quantaosun\n",
        "- what you get is ***the docking score + docked complex + a 2D image + a 3D interaction pymol session file*** \n",
        "- the docked complex and 2D image are readily to be used, The 3D pymol session need to be downloaded to a local pymol for further polish if you want a publication quality"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Install conda"
      ],
      "metadata": {
        "id": "nuGUIDxbW_IW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Conda installation\n",
        "!pip install -q condacolab\n",
        "import condacolab\n",
        "condacolab.install()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GCIuO0Lz_vMW",
        "outputId": "f110e540-5899-4846-a53c-da9d85e5f3ea"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "⏬ Downloading https://github.com/jaimergp/miniforge/releases/latest/download/Mambaforge-colab-Linux-x86_64.sh...\n",
            "📦 Installing...\n",
            "📌 Adjusting configuration...\n",
            "🩹 Patching environment...\n",
            "⏲ Done in 0:00:20\n",
            "🔁 Restarting kernel...\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title ### **Import Google Drive**\n",
        "#@markdown Click in the \"Run\" buttom to make your Google Drive accessible.\n",
        "from google.colab import drive\n",
        "\n",
        "drive.flush_and_unmount()\n",
        "drive.mount('/content/drive', force_remount=True)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kMhdx0Ms1fO9",
        "outputId": "3582ba17-7dc7-4562-a380-8a64c7a430db"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Drive not mounted, so nothing to flush and unmount.\n",
            "Mounted at /content/drive\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Please provide PDB ID and ligand structure you want to dock with"
      ],
      "metadata": {
        "id": "6KKs2IBRXDUh"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "s5B8su8nB0mS",
        "outputId": "18641569-f370-41db-aad4-fe8dd9628b6e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "/content/drive/MyDrive/3rd_review/20230320\n",
            "--2023-03-20 11:46:33--  https://files.rcsb.org/download/4YLK.pdb\n",
            "Resolving files.rcsb.org (files.rcsb.org)... 132.249.210.234\n",
            "Connecting to files.rcsb.org (files.rcsb.org)|132.249.210.234|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: unspecified [application/octet-stream]\n",
            "Saving to: ‘4YLK.pdb’\n",
            "\n",
            "4YLK.pdb                [  <=>               ] 579.10K  2.24MB/s    in 0.3s    \n",
            "\n",
            "2023-03-20 11:46:33 (2.24 MB/s) - ‘4YLK.pdb’ saved [593001]\n",
            "\n"
          ]
        }
      ],
      "source": [
        "#@title Input expected from the user\n",
        "path = \".\" #@param {type:\"string\"}\n",
        "%cd '{path}'\n",
        "PDB_ID = \"3HTB\" #@param {type:\"string\"}\n",
        "pdb = PDB_ID + \".pdb\"\n",
        "! wget https://files.rcsb.org/download/$pdb\n",
        "#@title Ligand three-letter name\n",
        "Ligand = \"JZ4\" #@param {type:\"string\"}\n",
        "smiles = \"c1ccccc1\" #@param {type:\"string\"}\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!conda --version"
      ],
      "metadata": {
        "id": "eeexDciMilD1",
        "outputId": "e7ab28f5-b45e-4505-dd68-568d1d84daa0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "conda 22.11.1\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Install dependencies"
      ],
      "metadata": {
        "id": "WwE9IjM3XMd9"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "jupyter": {
          "source_hidden": true
        },
        "id": "QxFkBBWq0EO1"
      },
      "source": [
        "#@title Install independencies\n",
        "!conda install -c rdkit rdkit -y 2>&1 1>/dev/null\n",
        "!conda install -q -y -c openbabel openbabel 2>&1 1>/dev/null\n",
        "!pip install py3Dmol 2>&1 1>/dev/null\n",
        "!conda install -c conda-forge prolif 2>&1 1>/dev/null\n",
        "!conda install -c conda-forge mdanalysis # for 2D interaction visualization\n",
        "!wget --no-check-certificate https://sourceforge.net/projects/smina/files/smina.static/download -O smina.static 2>&1 1>/dev/null\n",
        "!chmod u+x smina.static\n",
        "#!./smina.static -h\n",
        "import sys\n",
        "import os\n",
        "from rdkit import Chem\n",
        "from rdkit.Chem import AllChem\n",
        "from rdkit.Chem import Draw\n",
        "from rdkit.Chem.Draw import IPythonConsole\n",
        "from rdkit.Chem import DataStructs\n",
        "from rdkit.Chem import RDConfig\n",
        "from rdkit.Chem import rdBase\n",
        "print(rdBase.rdkitVersion)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!conda install conda=22.9.0"
      ],
      "metadata": {
        "id": "GPgYknTDlrKE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!conda --version"
      ],
      "metadata": {
        "id": "b_bOLV8glttN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!conda install -c conda-forge pymol-open-source"
      ],
      "metadata": {
        "id": "6DuIXbFGlKL6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Prepare receptor"
      ],
      "metadata": {
        "id": "AIUJdJuLRzs4"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nbHM8p20i7pD",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "#@title Write out a script for pymol\n",
        "com_file = open('split.pml','w')\n",
        "com_file.write('''\n",
        "load 3HTB.pdb\n",
        "remove resn HOH\n",
        "h_add elem O or elem N\n",
        "select 3HTB-JZ4, resn JZ4 #Create a selection called 1OYT-FSN from the ligand\n",
        "select 3HTB-receptor, 3HTB and not 3HTB-JZ4 #Select all that is not the ligand\n",
        "save 3HTB-JZ4.pdb, 3HTB-JZ4\n",
        "save 3HTB-receptor.pdb, 3HTB-receptor    \n",
        "''')\n",
        "com_file.close()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Correct the script with right PDB ID and ligand name\n",
        "!sed -i 's/3HTB/{PDB_ID}/g' split.pml\n",
        "!sed 's/JZ4/{Ligand}/g' split.pml\n",
        "!sed -i 's/JZ4/{Ligand}/g' split.pml"
      ],
      "metadata": {
        "cellView": "form",
        "id": "Gpjf85fDeTf5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "scOeFAqYi7pD"
      },
      "outputs": [],
      "source": [
        "#@title Run pymol to split protein and native ligand\n",
        "%cd '{path}'\n",
        "! pymol -c split.pml"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Variable assignment"
      ],
      "metadata": {
        "id": "BFtBn4YxVjYk"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "Protein_PDB_file_name = PDB_ID + '-receptor.pdb'\n",
        "ligand_name = PDB_ID + '-' + Ligand + '.pdb'"
      ],
      "metadata": {
        "id": "DTfZKqsqUUDV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# PDBFixer to prepare the receptor"
      ],
      "metadata": {
        "id": "pZKX9tmTypnT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!conda install -c conda-forge openmm -y"
      ],
      "metadata": {
        "id": "9bfwlEoWx-Ft"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!conda install -c conda-forge pdbfixer"
      ],
      "metadata": {
        "id": "lZ-B8goUx-Jd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pdbfixer --help"
      ],
      "metadata": {
        "id": "EQ6wF3F_x-pB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "QcuSUgrX0ldg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pdbfixer  {Protein_PDB_file_name} --ph=7  --replace-nonstandard --add-residues "
      ],
      "metadata": {
        "id": "gld65kMi0Ka8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!cat output.pdb > receptor.pdb"
      ],
      "metadata": {
        "id": "8POZF8t_16K1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Caculating charges for receptor"
      ],
      "metadata": {
        "id": "Yq55r4UdpoJg"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "jupyter": {
          "source_hidden": true
        },
        "id": "y0IUuw0h0EO7"
      },
      "source": [
        "#@title Target and native ligand preparation before docking\n",
        "#!obabel {Protein_PDB_file_name} -xr -O receptor.pdbqt\n",
        "!obabel receptor.pdb -xr -O receptor.pdbqt\n",
        "!obabel {ligand_name} -O native_ligand.pdbqt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Prepare ligand"
      ],
      "metadata": {
        "id": "Nev8kDowRvSK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Convert SMILE strings to SDF format\n",
        "!echo \"{smiles}\" > ligand2.smi\n",
        "!obabel ligand2.smi -osdf -O small.sdf --gen3d"
      ],
      "metadata": {
        "cellView": "form",
        "id": "BYYzix2oCgF2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Define small ligand path and number\n",
        "# File locations\n",
        "sdfFilePath = 'small.sdf' #@param {type:\"string\"}\n",
        "ConfoutputFilePath = 'small_conformation.sdf' #@param {type:\"string\"}\n",
        "\n",
        "inputMols = [x for x in Chem.SDMolSupplier(sdfFilePath,removeHs=False)]\n",
        "# Assign atomic chirality based on the structures:\n",
        "len(inputMols) # Check how many strucures"
      ],
      "metadata": {
        "cellView": "form",
        "id": "2MrRF4lkCjzT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "sKfSM19Jwj-f"
      },
      "outputs": [],
      "source": [
        "#@title SDF conformation generation for small molecule\n",
        "\n",
        "import multiprocessing\n",
        "\n",
        "# Download this from http://pypi.python.org/pypi/futures\n",
        "from concurrent import futures\n",
        "\n",
        "# conda install progressbar\n",
        "import progressbar\n",
        "\n",
        "#Find number cores available, leave two or system might be unresponsive\n",
        "numcores = multiprocessing.cpu_count()\n",
        "max_workers = numcores #-1\n",
        "\n",
        "#Knowledge based torsion generator http://pubs.acs.org/doi/abs/10.1021/acs.jcim.5b00654\n",
        "# This function is called in the subprocess.\n",
        "# The parameters (molecule and number of conformers) are passed via a Python\n",
        "\n",
        "ps = AllChem.ETKDG()\n",
        "ps.pruneRmsThresh=0.5\n",
        "ps.numThreads=0\n",
        "#Edit for number of confs desired eg n = 5\n",
        "n=5\n",
        "def generateconformations(m, n, name):\n",
        "    m = Chem.AddHs(m)\n",
        "    ids=AllChem.EmbedMultipleConfs(m, n, ps)\n",
        "    for id in ids:\n",
        "        AllChem.UFFOptimizeMolecule(m, confId=id)\n",
        "    # EmbedMultipleConfs returns a Boost-wrapped type which\n",
        "    # cannot be pickled. Convert it to a Python list, which can.\n",
        "    return m, list(ids), name\n",
        "\n",
        "smi_input_file, sdf_output_file = sys.argv[1:3]\n",
        "\n",
        "\n",
        "writer = Chem.SDWriter(ConfoutputFilePath)\n",
        "# suppl = [x for x in Chem.SDMolSupplier(sdfFilePath,removeHs=False)]\n",
        "#suppl = Chem.SmilesMolSupplier(smi_input_file, titleLine=False)\n",
        "\n",
        "# for mol in suppl:\n",
        "#     print(mol.GetPropsAsDict(includePrivate=True).get('_Name'))\n",
        "\n",
        "with futures.ProcessPoolExecutor(max_workers=max_workers) as executor:\n",
        "    # Submit a set of asynchronous jobs\n",
        "    jobs = []\n",
        "    for mol in inputMols:\n",
        "        if mol:\n",
        "            name = mol.GetProp('_Name')\n",
        "            job = executor.submit(generateconformations, mol, n, name)\n",
        "            jobs.append(job)\n",
        "\n",
        "    widgets = [\"Generating conformations; \", progressbar.Percentage(), \" \",\n",
        "               progressbar.ETA(), \" \", progressbar.Bar()]\n",
        "    pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(jobs))\n",
        "    for job in pbar(futures.as_completed(jobs)):\n",
        "        mol, ids, name = job.result()\n",
        "        mol.SetProp('_Name', name)\n",
        "        for id in ids:\n",
        "            writer.write(mol, confId=id)\n",
        "writer.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Define Docking inputs variables"
      ],
      "metadata": {
        "id": "0W-wd1QGRjXb"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "6N-veZedu3LD"
      },
      "outputs": [],
      "source": [
        "#@title 👋 Input names needed for Dock1\n",
        "ProteinForDocking = \"receptor.pdbqt\" #@param {type:\"string\"}\n",
        "LigandFromProtein = \"native_ligand.pdbqt\" #@param {type:\"string\"}\n",
        "Docked_out_FilePath = \"Docked.pdb\" #@param {type:\"string\"}\n",
        "FlexibleDockedFilePath = \"FlexibleDocked.sdf.gz\" #@param {type:\"string\"}"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Docking "
      ],
      "metadata": {
        "id": "eBRDz0YhaSEu"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title  Normal Docking (docking box is 8Å greater than the ligand specified. The –exhaustiveness is 8 )\n",
        "!./smina.static --seed 0 --autobox_ligand '{LigandFromProtein}' -r '{ProteinForDocking}' -l '{ConfoutputFilePath}'  -o '{Docked_out_FilePath}'"
      ],
      "metadata": {
        "cellView": "form",
        "id": "HpT1N-2sGFtN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# ⚠️⚠️(Alternatively, if you are patient enough, try this slow but mor accurate docking),otherwise just skip this"
      ],
      "metadata": {
        "id": "QjV2o4OragOp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title ⚠️Flexible docking, use with caution, it take super long time。 You can just skip unless you really want a relative accurate docking.\n",
        "#!./smina.static --seed 0 --autobox_ligand '{LigandFromProtein}' -r '{ProteinForDocking}' -l '{ConfoutputFilePath}' --flexdist_ligand '{LigandFromProtein}' --flexdist 5 -o '{Docked_out_FilePath}'"
      ],
      "metadata": {
        "id": "RC0b4bFbHF16",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Convert docked ligand to sdf"
      ],
      "metadata": {
        "id": "Ydb-TlS7WbAl"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eAwnSc3n0EO7"
      },
      "source": [
        "!obabel Docked.pdb -O Docked.sdf"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Visulization of Docked result\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "rC1idNnsnlkk"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Define display style\n",
        "#@title Define receptor name\n",
        "TargetProtein = Protein_PDB_file_name\n",
        "def DrawDocking(protein,ligand):\n",
        "     complex_pl = Chem.MolToPDBBlock(protein)\n",
        "     docked_pdb = Chem.MolToPDBBlock(ligand)\n",
        "     viewer = py3Dmol.view(width=800,height=800)\n",
        "     viewer.addModel(complex_pl,'pdb')\n",
        "     viewer.addModel(docked_pdb)\n",
        "     prot = {'resn': [\"DMS\", \"UNL\", \"SO4\", \"LIG\", \"HOH\", \"Cl\"], 'invert': 1}  #define prot as all except list\n",
        "     viewer.setStyle(prot,{'cartoon': {'colorscheme':'ssPyMol'}}) # Colour by secondary structure\n",
        "     Lig_373 = {'resn' : 'LIG'} #original ligand in pdb file\n",
        "     MyLig = {'resn':'UNL'} #ligand to be added from docking\n",
        "     viewer.addSurface(py3Dmol.VDW,{'opacity':0.7, 'color': 'white'}, prot)\n",
        "     viewer.setStyle({'resi': '132'}, {'stick': {'colorscheme': 'blueCarbon'}})\n",
        "     viewer.setStyle({'resi': '147'}, {'stick': {'colorscheme': 'blueCarbon'}})\n",
        "     viewer.setStyle({'resi': '311'}, {'stick': {'colorscheme': 'blueCarbon'}})\n",
        "     viewer.setStyle(Lig_373,{'stick':{'colorscheme': 'whiteCarbon','radius':.1}}) \n",
        "     viewer.setStyle(MyLig,{'stick':{'colorscheme' : 'greenCarbon'}})\n",
        "     viewer.zoomTo(MyLig)\n",
        "     return viewer "
      ],
      "metadata": {
        "id": "QG8-WFP0zAJ8",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "jupyter": {
          "source_hidden": true
        },
        "cellView": "form",
        "id": "7MvHqxEU0EO7"
      },
      "source": [
        "#@title Import visualization software\n",
        "\n",
        "import py3Dmol\n",
        "\n",
        "def drawit2(m,confId=-1):\n",
        "    mb = Chem.MolToMolBlock(m,confId=confId)\n",
        "    p = py3Dmol.view(width=400, height=400)\n",
        "    p.addModel(mb,'sdf')\n",
        "    p.setStyle({'stick':{}})\n",
        "    p.setBackgroundColor('0xeeeeee')\n",
        "    p.zoomTo()\n",
        "    return p\n",
        "\n",
        "def DrawComplex(protein,ligand):\n",
        "    complex_pl = Chem.MolToPDBBlock(Chem.CombineMols(protein,ligand))\n",
        "    #complex_mol=Chem.CombineMols(receptor,mols[-1])\n",
        "    view = py3Dmol.view(width=600,height=600)\n",
        "    view.addModel(complex_pl,'pdb')\n",
        "    #view.addModel(Chem.MolToMolBlock(mols[0]),'sdf')\n",
        "    chA = {'chain':['H','L','I']}\n",
        "    chB = {'resn':'UNL'}\n",
        "    view.setStyle(chA,{'cartoon': {'color':'spectrum'}})\n",
        "    #view.setStyle(chA,{'lines': {}})\n",
        "    view.addSurface(py3Dmol.VDW,{'opacity':0.8}, chB)\n",
        "    view.setStyle(chB,{'stick':{}})\n",
        "    view.zoomTo()\n",
        "    return view   \n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title View docked ligand\n",
        "#pose = \"2\" #@param [\"0\", \"1\", \"2\", \"3\", \"4\", \"5\", \"6\", \"7\", \"8\"]\n",
        "mols = [m for m in Chem.SDMolSupplier('Docked.sdf') if m is not None]\n",
        "#drawit2(mols[1])"
      ],
      "metadata": {
        "cellView": "form",
        "id": "8ng6uFiBN_M-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Show docked top1 (zoom out)\n",
        "receptor = Chem.MolFromPDBFile(TargetProtein)\n",
        "DrawComplex(receptor,mols[0])"
      ],
      "metadata": {
        "id": "vWtjk5XVzBf3",
        "cellView": "form"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Show the docked pocket (zoom in)\n",
        "receptor = Chem.MolFromPDBFile(TargetProtein)\n",
        "DrawDocking(receptor,mols[0])"
      ],
      "metadata": {
        "cellView": "form",
        "id": "wa4NnC4MRKEI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Export the receptor and fix it with PDBfixer,\n",
        "Before receptor is merged, we want to make sure it is well fixed"
      ],
      "metadata": {
        "id": "9NMNSGndkqPV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# merge the docked pose and the receptor"
      ],
      "metadata": {
        "id": "9lqDyA6125cN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#!conda install -c conda-forge openmm -y 2>&1 1>/dev/null\n",
        "#!conda install -c conda-forge pdbfixer -y 2>&1 1>/dev/null"
      ],
      "metadata": {
        "id": "NRP2lxoX5s_x"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@ Title To make sure the correct charge for MDanalysis, fix it again\n",
        "#!pdbfixer  {Protein_PDB_file_name} --ph=7  --replace-nonstandard --add-residues "
      ],
      "metadata": {
        "id": "ftXwt7HI59l7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Export the top1 docked ligand and add H.\n",
        "To achieve this, we convert pdb to mol2, then back to pdb. "
      ],
      "metadata": {
        "id": "EJsvj08Z7ybG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "m=Chem.MolFromPDBFile('Docked.pdb')\n",
        "Draw.MolToImage(m)"
      ],
      "metadata": {
        "id": "vigO8Iaw7DCT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!obabel -ipdb Docked.pdb -omol2 -O Docked.mol2 -m -h"
      ],
      "metadata": {
        "id": "zMzqY4Mh7Oyx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "m=Chem.MolFromMol2File('Docked1.mol2',sanitize=False)\n",
        "Draw.MolToImage(m)"
      ],
      "metadata": {
        "id": "HxiI6Tg27TMZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!obabel -imol2 Docked1.mol2 -opdb -O Docked1.pdb -h"
      ],
      "metadata": {
        "id": "tgyLnj_F7TPA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "m=Chem.MolFromPDBFile('Docked1.pdb')\n",
        "Draw.MolToImage(m)"
      ],
      "metadata": {
        "id": "os-vj1Y47TRr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Merge the docked TOP1 and the receptor"
      ],
      "metadata": {
        "id": "2MpY-l_68j4t"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Write out a script for merge\n",
        "com_file = open('merge.pml','w')\n",
        "com_file.write('''\n",
        "load output.pdb\n",
        "load Docked1.pdb\n",
        "h_add elem O or elem N\n",
        "save complex_prepared.pdb, all   \n",
        "''')\n",
        "com_file.close()"
      ],
      "metadata": {
        "cellView": "form",
        "id": "H3vXpVrn25Zm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "!pymol -c merge.pml # output.pdb is the default output name"
      ],
      "metadata": {
        "id": "LQ-IQIdj5gLv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Install dependencies for 2D interaction"
      ],
      "metadata": {
        "id": "m3PBsVv08yRh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!conda install -c conda-forge mdanalysis 2>&1 1>/dev/null\n",
        "!conda install -c conda-forge prolif 2>&1 1>/dev/null"
      ],
      "metadata": {
        "id": "5IjWpXX6hXXG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import MDAnalysis as mda\n",
        "import prolif as plf"
      ],
      "metadata": {
        "id": "sBgJvEoVil_O"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Define the merged complex as a MDAnalysis object"
      ],
      "metadata": {
        "id": "2128SVp385Us"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "# load topology\n",
        "u = mda.Universe('complex_prepared.pdb')\n",
        "lig = u.select_atoms(\"resname UNL\")\n",
        "prot = u.select_atoms(\"protein\")"
      ],
      "metadata": {
        "id": "-aZJe8mAimF4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Clean up the ligand before 2D"
      ],
      "metadata": {
        "id": "idXK66oB9j-W"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from rdkit import Chem\n",
        "from rdkit.Chem import Draw\n",
        "\n",
        "# create a molecule from the MDAnalysis selection\n",
        "lmol = plf.Molecule.from_mda(lig)\n",
        "# cleanup before drawing\n",
        "mol = Chem.RemoveHs(lmol)\n",
        "mol.RemoveAllConformers()\n",
        "Draw.MolToImage(mol, size=(400, 200))"
      ],
      "metadata": {
        "id": "JNcvWmG8kEmT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Draw the 2D interaction"
      ],
      "metadata": {
        "id": "z7EmdPEg9pC0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from prolif.plotting.network import LigNetwork\n",
        "\n",
        "fp = plf.Fingerprint()\n",
        "fp.run(u.trajectory[::10], lig, prot)\n",
        "df = fp.to_dataframe(return_atoms=True)\n",
        "\n",
        "net = LigNetwork.from_ifp(\n",
        "    df,\n",
        "    lmol,\n",
        "    # replace with `kind=\"frame\", frame=0` for the other depiction\n",
        "    kind=\"aggregate\",\n",
        "    threshold=0.3,\n",
        "    rotation=270,\n",
        ")\n",
        "net.display()"
      ],
      "metadata": {
        "id": "4aVukh4EimKF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Towards 3D interaction images\n",
        "\n",
        "The following section will help you obtain a publication-quality 3D protein-ligand interaction image.\n",
        "- A Pymol pre_pymol.pse file will be generated \n",
        "- The pse file shows the residues within 5 angstroms of ligand\n",
        "- You are supposed to download the pse to finish the drawing of the interaction lines like H bond dash lines and hydrophobic interaction lines since these lines are pretty challenging to process programmatically due to it varying too much from protein to protein"
      ],
      "metadata": {
        "id": "fhdrBDiodyFM"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 2. 👋 Write pymol script\n",
        "com_file = open('image.pml','w')\n",
        "com_file.write('''\n",
        "load output.pdb, tmp\n",
        "load Docked1.pdb, docked\n",
        "#select the surrounding then set it to an object \n",
        "# find 5 angstroms around ligand\n",
        "select resn UNL\n",
        "select res, byres sele around 5\n",
        "#set the ligand to an lig object\n",
        "select lig, resn UNL\n",
        "# Select all to hide\n",
        "select all\n",
        "# hide everythong\n",
        "hide\n",
        "# show the ligand and surroundings\n",
        "show sticks, lig|res\n",
        "# hide the H atom on the sticks\n",
        "hide sticks, h.\n",
        "# zoon on lig and res, background faded\n",
        "zoom lig|res\n",
        "# Color ligand carbon atom yellow\n",
        "color yellow, lig&name C*\n",
        "# Color ligand carbon atom yellow\n",
        "color cyan, res&name C*\n",
        "# Set label and font size\n",
        "# note you still have to click \"residues\" form L dropdown menue to show the label\n",
        "set label_size, 28\n",
        "set label_font_id, 8\n",
        "set label_color, grey\n",
        "# We also define some hydrophobic interaction style\n",
        "set sphere_scale, 0.35\n",
        "set sphere_transparency, 0.3\n",
        "set sphere_color, cyan\n",
        "# Define some hydrogen bond. interaciton style\n",
        "set dash_length, 0.3\n",
        "set dash_radius, 0.08\n",
        "set dash_gap, 0.3\n",
        "# set backgroud white\n",
        "bg_color white\n",
        "show cartoon\n",
        "set cartoon_transparency, 0.3\n",
        "bg_color white\n",
        "set ray_trace_mode, 1\n",
        "png binding_pocket.png\n",
        "save pre_pymol.pse\n",
        "''')\n",
        "com_file.close()\n"
      ],
      "metadata": {
        "id": "WmWteOvBSna5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Generate a 3D pymol session"
      ],
      "metadata": {
        "id": "xVd6aKTRdlrl"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_b3oglqfNwOw"
      },
      "source": [
        "#@title 4.👋 Generate a pymol session pse file that 70% quality of publication, please download to polish it more.\n",
        "!pymol -c image.pml\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Please download pre_pymol.pse to continue making the image with your local Pymol"
      ],
      "metadata": {
        "id": "cRD0_uLiaf22"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# References\n",
        "\n",
        "All the dependencies used are open soured, they are \n",
        "\n",
        "*   rdkit\n",
        "*   openbabel\n",
        "*   py3Dmol \n",
        "*   pymol-open-source \n",
        "*   smina.static\n",
        "*   MDAnalysis\n",
        "*   ProLif\n",
        "*   Python 3.8\n",
        "*   conda\n",
        "*   colab conda\n",
        "\n",
        "- This notebook has referred [Is life worth living](https://www.cheminformania.com/ligand-docking-with-smina/)\n",
        "- Thi notebook has referred [Jupyter Notebook for docking either locally or using Colab](https://www.macinchem.org/reviews/JupyterDocking/jupyterdocking.php)\n",
        " "
      ],
      "metadata": {
        "id": "CoByyqgKY_h1"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Post MD simulatio for top 3 poses\n",
        "Please copy Docked1.pdb, Docked2.pdb and Docked3.pdb to a folder called ligands."
      ],
      "metadata": {
        "id": "TVr7HJ6zzMqK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import os\n",
        "import glob\n",
        "import numpy as np\n",
        "from simtk import unit, openmm\n",
        "from simtk.openmm import app\n",
        "from openff.toolkit.topology import Molecule, Topology\n",
        "from openff.toolkit.typing.engines.smirnoff import ForceField\n",
        "from openmmforcefields.generators import SystemGenerator\n",
        "\n",
        "\n",
        "# Define file paths\n",
        "protein_file = \"Output.pdb\"\n",
        "ligand_dir = \"ligands\"\n",
        "ligand_files = glob.glob(os.path.join(ligand_dir, \"*.pdb\"))\n",
        "output_dir = \"output\"\n",
        "\n",
        "# Load the protein and create a topology\n",
        "pdb = app.PDBFile(protein_file)\n",
        "forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3p.xml')\n",
        "modeller = app.Modeller(pdb.topology, pdb.positions)\n",
        "topology = modeller.getTopology()\n",
        "\n",
        "# Create a SystemGenerator\n",
        "ff = ForceField(\"openff-1.3.0.offxml\")\n",
        "system_generator = SystemGenerator(\n",
        "    forcefields=[\"amber14-all.xml\", \"amber14/tip3p.xml\"],\n",
        "    small_molecule_forcefield=ff,\n",
        "    molecules=[],\n",
        "    cache=None\n",
        ")\n",
        "\n",
        "# Loop over the ligand files\n",
        "for i, ligand_file in enumerate(ligand_files):\n",
        "    print(f\"Processing ligand file {i+1}/{len(ligand_files)}: {os.path.basename(ligand_file)}\")\n",
        "    \n",
        "    # Create a Molecule object for the ligand\n",
        "    ligand_mol = Molecule.from_file(ligand_file)\n",
        "\n",
        "    # Generate the topology and system for the ligand\n",
        "    ligand_top = Topology.from_molecules(ligand_mol)\n",
        "    ligand_sys = system_generator.create_system(ligand_top)\n",
        "\n",
        "    # Add the ligand to the modeller and merge topologies\n",
        "    modeller.add(ligand_mol.to_topology().to_openmm())\n",
        "    topology = modeller.getTopology()\n",
        "\n",
        "    # Create a combined system and set up simulation\n",
        "    system = forcefield.createSystem(\n",
        "        topology,\n",
        "        nonbondedMethod=app.PME,\n",
        "        nonbondedCutoff=1.0*unit.nanometers,\n",
        "        constraints=app.HBonds,\n",
        "        rigidWater=True,\n",
        "        ewaldErrorTolerance=0.0005\n",
        "    )\n",
        "    system.addForce(ligand_sys.getForce(0))\n",
        "    integrator = openmm.LangevinIntegrator(\n",
        "        300*unit.kelvin, 1.0/unit.picosecond, 2.0*unit.femtoseconds\n",
        "    )\n",
        "    platform = openmm.Platform.getPlatformByName(\"CUDA\")\n",
        "    properties = {'CudaPrecision': 'mixed'}\n",
        "    simulation = app.Simulation(topology, system, integrator, platform, properties)\n",
        "    simulation.context.setPositions(modeller.positions)\n",
        "\n",
        "    # Run simulation\n",
        "    simulation.reporters.append(app.DCDReporter(os.path.join(output_dir, f\"ligand{i+1}.dcd\"), 250))\n",
        "    simulation.reporters.append(app.StateDataReporter(os.path.join(output_dir, f\"ligand{i+1}.csv\"), 250, step=True, \n",
        "                                                       time=True, potentialEnergy=True, kineticEnergy=True, \n",
        "                                                       totalEnergy=True, temperature=True))\n",
        "    simulation.reporters.append(app.CheckpointReporter(os.path.join(output_dir, f\"ligand{i+1}.chk\"), 250))\n",
        "\n",
        "    print(\"Running simulation...\")\n",
        "    simulation.step(5000)\n",
        "    print(\"Done!\")\n"
      ],
      "metadata": {
        "id": "twRm2VkwzLac"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}