{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecule substructural alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are cases when we want to align molecules which only share \n",
    "a fraction of their structure, like ligands in a congeneric series, where\n",
    "a common scaffold or core is common to all ligands in the series.\n",
    "We can perform this alignment using HTMD function: ```maximalSubstructureAlignment``` \n",
    "\n",
    "In this example, we will use two ligands which bind to the Beta Adrenergic receptor.\n",
    "They share a fraction of their scaffold and bind in the same pocket."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.ui import *\n",
    "from moleculekit.tools.graphalignment import maximalSubstructureAlignment\n",
    "from htmd.home import home\n",
    "\n",
    "path = home(dataDir='test-molecule-graphalignment')\n",
    "\n",
    "# Load the molecule which will be used as reference\n",
    "reference_ligand = Molecule(os.path.join(path, 'ref_lig.pdb'))\n",
    "reference_ligand.view(style='Lines', name='Reference')\n",
    "\n",
    "# Load other ligand from this congeneric series\n",
    "ligand_to_align = Molecule(os.path.join(path, 'lig2align.pdb'))\n",
    "ligand_to_align.view(style='Licorice', name='ToAlign')\n",
    "\n",
    "# Align\n",
    "lig_aligned = maximalSubstructureAlignment(reference_ligand,ligand_to_align)\n",
    "\n",
    "# View the aligned molecule\n",
    "lig_aligned.view(style='Licorice', name='Aligned')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detailed explanation\n",
    "\n",
    "First, we load the molecule which will be used as reference, cocrystallized in PDB entry 2Y02 with resname WHJ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-12 13:08:39,766 - moleculekit.molecule - INFO - Removed 370 atoms. 4694 atoms remaining in the molecule.\n",
      "2024-06-12 13:08:39,797 - moleculekit.molecule - INFO - Removed 4640 atoms. 54 atoms remaining in the molecule.\n",
      "2024-06-12 13:08:39,798 - moleculekit.molecule - INFO - Removed 27 atoms. 27 atoms remaining in the molecule.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
       "       17, 18, 19, 20, 21, 22, 23, 24, 25, 26], dtype=int32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reference_crystal = Molecule('2Y02')\n",
    "reference_crystal.filter('protein or resname WHJ')\n",
    "reference_ligand = reference_crystal.copy()\n",
    "# Extracts the ligand from the protein\n",
    "reference_ligand.filter('resname WHJ')\n",
    "# There are two ligands with the same resname, so we select one of them\n",
    "reference_ligand.filter('residue 1') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's extract the other ligand from this congeneric series, cocrystallized in PDB entry 2Y03 with resname 5FW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-12 13:08:42,724 - moleculekit.molecule - INFO - Removed 4832 atoms. 30 atoms remaining in the molecule.\n",
      "2024-06-12 13:08:42,726 - moleculekit.molecule - INFO - Removed 15 atoms. 15 atoms remaining in the molecule.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29],\n",
       "      dtype=int32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ligand_to_align = Molecule('2Y03')\n",
    "ligand_to_align.filter('resname 5FW')\n",
    "# Again, there are two residues with the same name, and\n",
    "# residue 1 happens to be already aligned with the reference ligand, so we use residue 0\n",
    "ligand_to_align.filter('residue 0') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at both molecules as they are now. The reference ligand is displayed with narrower bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_ligand.view(style='Lines',name='Reference')\n",
    "ligand_to_align.view(style='Licorice',name='ToAlign')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we align the extracted ligand to the crystal of the molecule we are using as reference. Then, we can see the result. Again, the reference ligand is displayed with narrower bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "lig_aligned = maximalSubstructureAlignment(reference_ligand,ligand_to_align)\n",
    "lig_aligned.view(style='Licorice',name='Aligned')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can see how the aligned ligand looks together with the protein."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_crystal.view(sel='protein',style='NewCartoon',name='Protein')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}