{
 "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.molecule.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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.molecule.graphalignment import maximalSubstructureAlignment\n",
    "\n",
    "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_crystal.view(sel='protein',style='NewCartoon',name='Protein')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}