Step2_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Network Analysis Tutorial - Step 2\n",
    "\n",
    "In the second part of our generalized network analysis tutorial, the user is presented with another jupyter notebook where all the correlation maps calculated in the first part can be translated into molecular visualizations or plots. Here we provide an opportunity for users more comfortable with python programming to tailor the analysis to their specific scientific questions.\n",
    "\n",
    "The Network Analysis Tutorial is part of the work entitled **Generalized correlation-based dynamical network analysis: a new high-performance approach for identifying allosteric communications in molecular dynamics trajectories**, by Marcelo C. R. Melo, Rafael C. Bernardi, Cesar de la Fuente-Nunez, and Zaida Luthey-Schulten. For more information see http://faculty.scs.illinois.edu/schulten/. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:09:43.829281Z",
     "start_time": "2020-06-04T23:09:34.192482Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load the python package\n",
    "import os\n",
    "from dynetan.toolkit import *\n",
    "from dynetan.viz import *\n",
    "from dynetan.proctraj import *\n",
    "from dynetan.gencor import *\n",
    "from dynetan.contact import *\n",
    "from dynetan.datastorage import *\n",
    "\n",
    "#from numpy.linalg import norm\n",
    "from itertools import islice\n",
    "from itertools import combinations\n",
    "from scipy import stats\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "# For visualization\n",
    "import nglview    as nv\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:09:45.804654Z",
     "start_time": "2020-06-04T23:09:45.609293Z"
    }
   },
   "outputs": [],
   "source": [
    "# Allows better plotting inside jupyter notebooks\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:09:46.385284Z",
     "start_time": "2020-06-04T23:09:46.083635Z"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%R install.packages(\"librarian\", repos = \"https://cran.csiro.au/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:09:48.851893Z",
     "start_time": "2020-06-04T23:09:46.501424Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Make sure all the following packages are installed in your R distribution\n",
    "%R librarian::shelf(data.table, ggplot2, ggrepel, gdata, RColorBrewer, colorRamps, rPref, cran_repo = \"https://cran.csiro.au/\")\n",
    "#%R require(data.table)\n",
    "#%R require(ggplot2)\n",
    "#%R require(ggrepel)\n",
    "#%R require(gdata)\n",
    "#%R require(RColorBrewer)\n",
    "#%R require(colorRamps)\n",
    "#%R require(rPref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function: Open an NGL View window displaying the molecular system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_widget():\n",
    "    # http://proteinformatics.charite.de/ngl/doc/#User_manual/Usage/Selection_language\n",
    "    w = nv.show_mdanalysis(workUviz.select_atoms(\"all\"))\n",
    "    w._remote_call(\"setSize\", target=\"Widget\", args=[\"800px\", \"600px\"])\n",
    "    w.parameters = dict(theme='light')\n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function: Set basic visual representations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_widget_representations(widget):\n",
    "    widget.clear_representations()\n",
    "    widget.add_cartoon(\"backbone\")\n",
    "    widget.add_hyperball(\"charged and (not backbone or .CA or (PRO and .N))\")\n",
    "    w.add_representation(\"hyperball\",\n",
    "        selection='(not protein) and not (.1MC1 or .2MC1 or .1MC2 or .2MC2  or .1MC3 or .2MC3)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Files and system definitions (same as in \"ProcessTrajectory\" notebook) \n",
    "\n",
    "Make sure that all the paths are the same as in the **Step 1** notebook. Also, you might want to have analysis resutls saved to a new location, to avoid mixing results of different systems. If you are copying this tutorial to adapt to your own system, be careful to select new paths. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = Change to the relevant work folder\n",
    "systemExample=\"CFTR\"\n",
    "if systemExample is \"UbqCHARMM\":\n",
    "    %cd /home/zharmad/host/projects/Ubq-md\n",
    "elif systemExample is \"UbqAthi\":\n",
    "    %cd /home/zharmad/host/shared-colleague/Ubq-2017\n",
    "elif systemExample is \"periplasmic\":\n",
    "    %cd /home/zharmad/projects/periplasmic/leucine-binding_protein\n",
    "elif systemExample is \"CFTR\":\n",
    "    %cd /home/zharmad/projects/cftr/DyNetAn\n",
    "else:\n",
    "    %cd ..    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:22.898264Z",
     "start_time": "2020-06-04T23:10:22.892936Z"
    }
   },
   "outputs": [],
   "source": [
    "if systemExample is \"UbqCHARMM\":\n",
    "    %cd ..\n",
    "    # Define mutant file IO locations. wt, P67L, E56K, R75Q\n",
    "    state=\"apo1\"\n",
    "    dataDir = \"./dynetan\"\n",
    "    workDir = \"./dynetan/analysis\"\n",
    "    fileNameRoot = \"%s_5x1000\" % state\n",
    "elif systemExample is \"UbqAthi\":\n",
    "    state=\"UbqWT\"\n",
    "    dataDir = \"./dynetan\"\n",
    "    workDir = \"./dynetan/analysis\"\n",
    "    fileNameRoot = \"%s_5x200\" % state\n",
    "elif systemExample is \"peripllasmic\":    \n",
    "    state=\"holo\"\n",
    "    dataDir = \"./dynetan\"\n",
    "    workDir = \"./dynetan/analysis\"\n",
    "    fileNameRoot = \"%s_5x200\" % state\n",
    "elif systemExample is \"CFTR\":\n",
    "    # Define mutant file IO locations. wt, P67L, E56K, R75Q\n",
    "    allele=\"dF508\" ; temperature=\"310K\"\n",
    "    dataDir = \"./results/%s/%s/\" % (allele, temperature)\n",
    "    #Path where results will be written (you may want plots and data files in a new location)\n",
    "    workDir = \"./results/%s/%s/analysis\" % (allele, temperature)\n",
    "    fileNameRoot = \"1to3\"\n",
    "\n",
    "fullPathRoot = os.path.join(dataDir, fileNameRoot)\n",
    "dcdVizFile = fullPathRoot + \"_reducedTraj.dcd\"\n",
    "pdbVizFile = fullPathRoot + \"_reducedTraj.pdb\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:23.725216Z",
     "start_time": "2020-06-04T23:10:23.717418Z"
    }
   },
   "outputs": [],
   "source": [
    "# Creates the workdirectory, and Plots directory.\n",
    "import os\n",
    "if not os.path.exists(workDir):\n",
    "    os.makedirs(workDir)\n",
    "    os.makedirs(os.path.join(workDir, \"Plots\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:09:45.368030Z",
     "start_time": "2020-06-04T23:09:45.364339Z"
    }
   },
   "outputs": [],
   "source": [
    "dnad = DNAdata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:26.333160Z",
     "start_time": "2020-06-04T23:10:25.940523Z"
    }
   },
   "outputs": [],
   "source": [
    "dnad.loadFromFile(fullPathRoot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get some basic information about the graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( nx.info(dnad.nxGraphs[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( dnad.__dict__.keys() )\n",
    "#dnad.nxGraphs[0].number_of_edges()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load reduced trajectory, in case this notebook did not load ALL trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:27.286250Z",
     "start_time": "2020-06-04T23:10:27.123096Z"
    }
   },
   "outputs": [],
   "source": [
    "workUviz = mda.Universe(pdbVizFile, dcdVizFile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are loading the first trajectory of the notebook, we need to create the dnad.nodesAtmSel structure to make selections on the universe based on node indices. This is mostly used for visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:29.005597Z",
     "start_time": "2020-06-04T23:10:29.001742Z"
    }
   },
   "outputs": [],
   "source": [
    "# We add this to the object for ease of access.\n",
    "dnad.nodesAtmSel = workUviz.atoms[ dnad.nodesIxArray ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define some shorthand nodes for later visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:29.211185Z",
     "start_time": "2020-06-04T23:10:29.202288Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define the segID of the Ligand being studied.\n",
    "ligandSegID = \"\"\n",
    "if ligandSegID is not \"\":\n",
    "    nodesSite = getNodeFromSel(\"segid \" + ligandSegID + \" and resid 1503\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "    headNodeSite = getNodeFromSel(\"segid \" + ligandSegID + \" and name PB and resid 1503\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "    print(\"Nodes for active site: %s of %s\" % (headNodeSite, nodesSite))\n",
    "else:\n",
    "    nodesSite = getNodeFromSel(\"protein and resid 15 23 30 43 67 69\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "    headNodeSite = getNodeFromSel(\"protein and resid 43\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "    \n",
    "nodeAll = getNodeFromSel(\"all\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "nodeLabels=nodeAll[::10]\n",
    "\n",
    "nodePathA = getNodeFromSel(\"protein and resid 43\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "nodePathB = getNodeFromSel(\"protein and resid 69\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "#print(nodeLabels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and Visualize Paths and Communities\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show the optimal path between two nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:30.517687Z",
     "start_time": "2020-05-29T04:01:29.627802Z"
    }
   },
   "outputs": [],
   "source": [
    "w = create_widget()\n",
    "set_widget_representations(w)\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optpath = getPath( nodePathA[0], nodePathB[0], dnad.nodesAtmSel, dnad.preds)\n",
    "viewPath(w, optpath, dnad.distsAll, dnad.maxDirectDist, dnad.nodesAtmSel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show the contacts with N highest betweeness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.633150Z",
     "start_time": "2020-05-29T04:01:30.563342Z"
    }
   },
   "outputs": [],
   "source": [
    "w = create_widget()\n",
    "set_widget_representations(w)\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.754040Z",
     "start_time": "2020-05-29T04:01:31.636810Z"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the 100 highest betweennesses (edges most used by optimal paths)\n",
    "w.clear_representations()\n",
    "w.add_cartoon(\"backbone\")\n",
    "w.add_line(\"not protein\")\n",
    "\n",
    "# Retrieve information from one window. k denotes an edge between two nodes. v denotes betweenness.\n",
    "# Retrieve top 5% of such edges.\n",
    "window = 0\n",
    "numTop=int(dnad.nxGraphs[0].number_of_edges()/20)\n",
    "for k,v in islice(dnad.btws[window].items(),numTop):\n",
    "    #print(k,v, dnad.corrMatAll[0, k[0], k[1]])\n",
    "    viewPath(w, getPath(k[0], k[1], dnad.nodesAtmSel, dnad.preds), \n",
    "             dnad.distsAll, dnad.maxDirectDist, dnad.nodesAtmSel, window)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.759732Z",
     "start_time": "2020-05-29T04:01:31.756567Z"
    }
   },
   "outputs": [],
   "source": [
    "# Save the current visualization in a figue.\n",
    "#w.download_image(filename='nglview_frame.png', factor=15, trim=True, antialias=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.854567Z",
     "start_time": "2020-05-29T04:01:31.763392Z"
    }
   },
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual analysis of betweenness. Take mean betweenness of the above consensus edges from all windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the top betweenness from all windows, then identify a consensus subset by intersection.\n",
    "# print(\"= = Dimensions of the correlation matrix (windows, node, node): %i %i %i \" % dnad.corrMatAll.shape )\n",
    "setsBetweenness=[]\n",
    "meanNumEdges=np.mean([dnad.nxGraphs[x].number_of_edges() for x in range(dnad.numWinds) ])\n",
    "numTop=int(meanNumEdges/20)\n",
    "print(\"= = Retrieving top ~5%% (%i of ~%i edges) of node betweenness from all %i windows...\" \\\n",
    "      % (numTop,meanNumEdges, dnad.numWinds) )\n",
    "for window in range(dnad.numWinds):\n",
    "    s=set()\n",
    "    for k,v in islice(dnad.btws[window].items(),numTop):\n",
    "        s.add(k)\n",
    "    setsBetweenness.append(s)\n",
    "\n",
    "overlapsBetweenness=[]\n",
    "for i, j in combinations(setsBetweenness,2):\n",
    "    overlap = len(i.intersection(j))/numTop\n",
    "    #print( \"= = %2.2f%% overlap\" % (100*overlap) )\n",
    "    overlapsBetweenness.append( overlap )\n",
    "print( \"= = Mean pairwise overlap: %2.2f%%\" % (100*np.mean(overlap)) )\n",
    "\n",
    "setOverlapAll=setsBetweenness[0]\n",
    "for i in range(1,dnad.numWinds):\n",
    "    setOverlapAll=setOverlapAll.intersection(setsBetweenness[i])\n",
    "numTopConcensusBtws=len(setOverlapAll)\n",
    "print( \"= = %i edges found in top 5%% of all windows (or %2.2f%% of %i edges).\" % \\\n",
    "      (numTopConcensusBtws, 100*numTopConcensusBtws/numTop, numTop) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a new ordered-dict basedd on the mean value of all windows.\n",
    "tempDict={}\n",
    "tempDictStd={}\n",
    "for k in setOverlapAll:\n",
    "    meanValue=np.mean([ dnad.btws[x][k] for x in range(dnad.numWinds)] )\n",
    "    sigmaValue=np.std([ dnad.btws[x][k] for x in range(dnad.numWinds)] )\n",
    "    tempDict[k]=meanValue\n",
    "    tempDictStd[k]=sigmaValue\n",
    "\n",
    "concensusBtws=OrderedDict()\n",
    "concensusBtwsStd=OrderedDict()\n",
    "for k in sorted(tempDict, key=tempDict.__getitem__, reverse=True):\n",
    "    concensusBtws[k]=tempDict[k]\n",
    "    concensusBtwsStd[k]=tempDictStd[k]\n",
    "    #print( \"...nodes %i ~ %i : %.4f +- %.4f\" % (k[0], k[1], tempDict[k], tempDictStd[k]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get some selections from the top betweenness candidates\n",
    "setBtwsResidues=set()\n",
    "for k,v in islice(concensusBtws.items(),numTopConcensusBtws):\n",
    "#for k,v in islice(concensusBtws.items(),10):\n",
    "    path=getPath(k[0], k[1], dnad.nodesAtmSel, dnad.preds)\n",
    "    l=[]\n",
    "    for node in path:\n",
    "        selectionText=getSelFromNode(node, dnad.nodesAtmSel)\n",
    "        m=selectionText.split()              \n",
    "        #l.append(m[7]+\":\"+m[1]+\":\"+m[4])\n",
    "        l.append(int(m[4]))\n",
    "    #print(l, v)\n",
    "    setBtwsResidues.update(l)\n",
    "print(\" = = Number of distinct residues in set: %i\" % len(setBtwsResidues))\n",
    "print(sorted(setBtwsResidues))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.633150Z",
     "start_time": "2020-05-29T04:01:30.563342Z"
    }
   },
   "outputs": [],
   "source": [
    "w = create_widget()\n",
    "set_widget_representations(w)\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.754040Z",
     "start_time": "2020-05-29T04:01:31.636810Z"
    }
   },
   "outputs": [],
   "source": [
    "# Plot highest betweennesses (edges most used by optimal paths)\n",
    "w.clear_representations()\n",
    "w.add_representation(\"cartoon\", selection='backbone and not 1451-1480')\n",
    "w.add_representation(\"hyperball\",\n",
    "        selection='508 and (not backbone or .CA or (PRO and .N))')\n",
    "w.add_representation(\"hyperball\",\n",
    "        selection='(not protein) and (not POPC)')\n",
    "\n",
    "distOnes=np.ones_like(dnad.distsAll)\n",
    "window=0\n",
    "# Retrieve information from one window. k denotes an edge between two nodes. v denotes betweenness.\n",
    "# Retrieve top 5% of such edges.\n",
    "# The dists argument controls the thickness of cylinders, so put the average correlation coefficient matrix here.\n",
    "# The win arguments seems just to be fed into the distsa argument, so put 0.\n",
    "meanCorMat = np.zeros( (1,dnad.numNodes,dnad.numNodes), dtype=np.float32)\n",
    "meanCorMat[0] = np.mean(dnad.corrMatAll,axis=0)\n",
    "maxMeanCorMat = np.max(meanCorMat)\n",
    "for k,v in islice(concensusBtws.items(),numTopConcensusBtws):\n",
    "    #print(k,v, dnad.corrMatAll[0, k[0], k[1]])\n",
    "    viewPath(w, path=getPath(k[0], k[1], dnad.nodesAtmSel, dnad.preds), \n",
    "             dists=meanCorMat, maxDirectDist=maxMeanCorMat, nodesAtmSel=dnad.nodesAtmSel,\n",
    "             win=0, opacity=0.75, color='blue')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.759732Z",
     "start_time": "2020-05-29T04:01:31.756567Z"
    }
   },
   "outputs": [],
   "source": [
    "# Save the current visualization in a figue.\n",
    "w.download_image(filename='nglview_highestBetweenness.png', factor=4, trim=True, antialias=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:31.854567Z",
     "start_time": "2020-05-29T04:01:31.763392Z"
    }
   },
   "outputs": [],
   "source": [
    " w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyze Communities \n",
    "We keep the communities that have more than 1% of nodes in all windows. Then we group communities across replicas by largest intersection. This is needed because we have no guarantee that the same community will be assigned the same ID in different windows of the same simulation.\n",
    "\n",
    "We fin ally rank the communities by modularity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:35.825349Z",
     "start_time": "2020-06-04T23:10:35.722709Z"
    }
   },
   "outputs": [],
   "source": [
    "import networkx.algorithms.community.quality as nxquality\n",
    "\n",
    "# Creates a list of windows and order them according to graph modularity.\n",
    "windModul = []\n",
    "for window in range(dnad.numWinds):\n",
    "    modul = nxquality.modularity(dnad.nxGraphs[window], \n",
    "                         [ set(nodesList) for nodesList in dnad.nodesComm[window][\"commNodes\"].values()])\n",
    "    windModul.append((window, modul))\n",
    "    \n",
    "windModul.sort(key=lambda x:x[1], reverse=True)\n",
    "\n",
    "# Keep the window with the highest modularity as a reference for community matching\n",
    "refWindow = windModul[0][0]\n",
    "\n",
    "for wind, mod in windModul[:5]:\n",
    "    print( \"Window {} has modularity {:1.4f}.\".format(wind, mod) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:36.127259Z",
     "start_time": "2020-06-04T23:10:36.094895Z"
    }
   },
   "outputs": [],
   "source": [
    "cutoff = max(10, np.ceil(0.01*dnad.numNodes))\n",
    "\n",
    "def matchComm(mCommID, mWindow, refWindow, dnad, cutoff=1):\n",
    "    \"\"\"\n",
    "    Returns the community ID for the reference window that has the largest\n",
    "    intersection with the matching community at the matching window.\n",
    "    Communities at the reference window with less than *cutoff* percent of nodes\n",
    "    are ignored.\n",
    "    \"\"\"\n",
    "    \n",
    "    trgtComm = -1\n",
    "    intersectSize = 0\n",
    "    for commID in dnad.nodesComm[refWindow][\"commOrderSize\"]:\n",
    "        # Skip community if it has less than one percent of the nodes.\n",
    "        commSize = len(dnad.nodesComm[refWindow][\"commNodes\"][commID])\n",
    "        if commSize < cutoff:\n",
    "            continue\n",
    "        \n",
    "        tmpSize = len( set(dnad.nodesComm[refWindow][\"commNodes\"][commID]).intersection( \n",
    "            set(dnad.nodesComm[mWindow][\"commNodes\"][mCommID]) ) )\n",
    "        \n",
    "        # Selects the largets intersection\n",
    "        if intersectSize < tmpSize:\n",
    "            intersectSize = tmpSize\n",
    "            trgtComm = commID\n",
    "    \n",
    "    return trgtComm, intersectSize\n",
    "\n",
    "print(\"Using reference window {0} with highest modularity {1:<1.4}\".format(*windModul[0]))\n",
    "\n",
    "communities = defaultdict(list)\n",
    "for window in range(dnad.numWinds):\n",
    "    for commID in dnad.nodesComm[window][\"commOrderSize\"]:\n",
    "        \n",
    "        # Skip community if it has less than one percent of the nodes.\n",
    "        commSize = len(dnad.nodesComm[window][\"commNodes\"][commID])\n",
    "        if commSize < cutoff:\n",
    "            continue\n",
    "        \n",
    "        matchID, interSize = matchComm(commID, window, refWindow, dnad, cutoff)\n",
    "        \n",
    "        communities[matchID].append( (commID, interSize, window) )\n",
    "        \n",
    "communities = {key:val for (key,val) in communities.items() }\n",
    "communities.keys()\n",
    "\n",
    "# Creates a list of communities ID from the dictionary keys\n",
    "# Orders the keys according to mean intersection size over all windows.\n",
    "tmpList = []\n",
    "for key,val in communities.items():\n",
    "    tmpList.append((key, np.mean([pair[1] for pair in val]), len(val)))\n",
    "tmpList.sort(key=lambda x:x[1], reverse=True)\n",
    "tmpList\n",
    "\n",
    "# Creates a pandas data frame for plotting and analysis\n",
    "commList = []\n",
    "genCommID = 0\n",
    "for key in [x[0] for x in tmpList]:\n",
    "    val = communities[key]\n",
    "    for valList in val:\n",
    "        commList.append( [genCommID, *valList ] )\n",
    "    genCommID += 1\n",
    "\n",
    "commDF = pd.DataFrame(data=commList, columns=[\"genCommID\",\"commID\",\"interSize\",\"Window\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:36.816115Z",
     "start_time": "2020-06-04T23:10:36.794165Z"
    }
   },
   "outputs": [],
   "source": [
    "# Changes \"genCommID\" for communities that are matched to the same community in the reference window.\n",
    "c = commDF.groupby([\"genCommID\",\"Window\"]).cumcount()\n",
    "c *= 0.1\n",
    "commDF[ \"genCommID\" ] += c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:37.185313Z",
     "start_time": "2020-06-04T23:10:37.122868Z"
    }
   },
   "outputs": [],
   "source": [
    "# Creates a NumPy 2D array to organize data and transform it in a pandas DF.\n",
    "# Not pretty but its pynthon...\n",
    "nodeCommNP = np.empty([dnad.numNodes, dnad.numWinds])\n",
    "nodeCommNP.fill(-1)\n",
    "\n",
    "#Group by general community ID\n",
    "grpBy = commDF.groupby(\"genCommID\")\n",
    "for genCommID, group in grpBy:\n",
    "    for winIndx,commID in group[[\"Window\",\"commID\"]].values:\n",
    "        for node in range(dnad.numNodes):\n",
    "            if dnad.nxGraphs[winIndx].nodes[node][\"modularity\"] == commID:\n",
    "                nodeCommNP[node, winIndx] = genCommID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:37.328246Z",
     "start_time": "2020-06-04T23:10:37.314017Z"
    }
   },
   "outputs": [],
   "source": [
    "# Removes nodes that were not classified in a \"big-nough\" (bigger than 1%) cluster in *any* window.\n",
    "nodeCommDF = pd.DataFrame(data=nodeCommNP,columns=[\"Window\"+str(i) for i in range(dnad.numWinds)])\n",
    "nodeCommDF[\"Node\"] = [i for i in range(dnad.numNodes)]\n",
    "nodeCommDF = nodeCommDF[ nodeCommDF.min(1) >= 0]\n",
    "# So we don't get \"blank\"/empty areas in the plot\n",
    "nodeCommDF[\"NodePlot\"] = [i for i in range(len(np.unique(nodeCommDF[\"Node\"])))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:37.789366Z",
     "start_time": "2020-06-04T23:10:37.770084Z"
    }
   },
   "outputs": [],
   "source": [
    "# Checks that target nodes are classified in ALL windows\n",
    "nodeCommDF.loc[ nodeCommDF[\"Node\"].isin(nodesSite) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:38.201629Z",
     "start_time": "2020-06-04T23:10:38.170785Z"
    }
   },
   "outputs": [],
   "source": [
    "# Melts for plotting.\n",
    "nodeCommDFmelt = nodeCommDF.melt(id_vars=[\"Node\",\"NodePlot\"], value_name=\"Cluster\", var_name=\"Window\")\n",
    "# Makes it easier to plot\n",
    "nodeCommDFmelt[\"Cluster\"] = nodeCommDFmelt[\"Cluster\"].astype('category')\n",
    "# Makes it easier to plot\n",
    "for i in range(dnad.numWinds):\n",
    "    nodeCommDFmelt.replace(\"Window\"+str(i),i, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:38.792384Z",
     "start_time": "2020-06-04T23:10:38.762821Z"
    }
   },
   "outputs": [],
   "source": [
    "nodeCommDFmelt.loc[nodeCommDFmelt[\"Node\"].isin(nodesSite)].groupby(\"Node\")[\"Cluster\"].apply(np.unique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:39.108358Z",
     "start_time": "2020-06-04T23:10:39.098285Z"
    }
   },
   "outputs": [],
   "source": [
    "trgtClusters = np.unique( nodeCommDFmelt.loc[nodeCommDFmelt[\"Node\"].isin(nodesSite), \"Cluster\"].values )\n",
    "trgtClusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:39.988933Z",
     "start_time": "2020-06-04T23:10:39.535287Z"
    }
   },
   "outputs": [],
   "source": [
    "# Add readable info to nodes\n",
    "def getTagStr(i):\n",
    "    # Store atom names for residues with multiple nodes\n",
    "    if len(getNodeFromSel( getSelFromNode(i, dnad.nodesAtmSel), dnad.nodesAtmSel, dnad.atomToNode)) > 1:\n",
    "        atmStr = \":\" + dnad.nodesAtmSel.atoms[i].name\n",
    "    else:\n",
    "        atmStr = \"\"\n",
    "        \n",
    "    retStr = dnad.nodesAtmSel.atoms[i].resname.capitalize() + \\\n",
    "            \":\" + str(dnad.nodesAtmSel.atoms[i].resid) + \\\n",
    "            atmStr + \\\n",
    "            \"_\" + dnad.nodesAtmSel.atoms[i].segid\n",
    "            \n",
    "    return retStr\n",
    "\n",
    "nodeCommDFmelt['resid']     = np.vectorize(getTagStr)(nodeCommDFmelt[\"Node\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:40.109538Z",
     "start_time": "2020-06-04T23:10:40.058119Z"
    }
   },
   "outputs": [],
   "source": [
    "# Write data for Ploting (plots from ggplot in R are much better!)\n",
    "nodeCommDFmelt.to_csv(os.path.join(workDir, \"cluster.csv\"),index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:40.336573Z",
     "start_time": "2020-06-04T23:10:40.273435Z"
    }
   },
   "outputs": [],
   "source": [
    "# Get all nodes that make contact with target nodes in any window\n",
    "contactNodes = np.unique( np.where( dnad.corrMatAll[:,nodesSite,:] > 0 )[2] )\n",
    "contactNodesTrgts = list(nodesSite)\n",
    "for node in contactNodes:\n",
    "    if len( set(trgtClusters).intersection( \n",
    "            set(np.unique(nodeCommDFmelt.loc[ nodeCommDFmelt[\"Node\"] == node].Cluster)) ) ) :\n",
    "        contactNodesTrgts.append(node)\n",
    "\n",
    "# Save data to file\n",
    "pd.DataFrame(contactNodesTrgts, columns=[\"contactNodesTrgts\"]).to_csv(\n",
    "    os.path.join(workDir, \"contactNodesTrgts.csv\"),index=False)\n",
    "\n",
    "pd.DataFrame(nodesSite, columns=[\"nodesSite\"]).to_csv(\n",
    "    os.path.join(workDir, \"nodesSite.csv\"),index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare pandas data frame with community data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:10:42.968288Z",
     "start_time": "2020-06-04T23:10:42.962222Z"
    }
   },
   "outputs": [],
   "source": [
    "# Prepares variable names for multi-system comparisons. In this tutorial, we only have one system.\n",
    "system1 = \"all\"\n",
    "cDF = nodeCommDFmelt\n",
    "cDF[\"system\"] = system1\n",
    "refWindow1 = refWindow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:39.531285Z",
     "start_time": "2020-06-04T23:11:39.498328Z"
    }
   },
   "outputs": [],
   "source": [
    "# Loads VMD-compatible color scales to match community colors in R plots, NGLView, and VMD figures.\n",
    "comColorScale = getCommunityColors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:40.648977Z",
     "start_time": "2020-06-04T23:11:40.577426Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,comColorScale -o colorValues,clusterIDs -w 800 -h 450\n",
    "\n",
    "dataPath = file.path(workDir, \"cluster.csv\")\n",
    "\n",
    "dt <- fread(dataPath)\n",
    "clusterIDs = dt[, unique(Cluster)]\n",
    "\n",
    "colourCount = length(unique(dt$Cluster))\n",
    "\n",
    "# We only have 50 availabl colors\n",
    "colourCount <- min(colourCount,50)\n",
    "\n",
    "rgbCodes <- data.table(comColorScale)\n",
    "\n",
    "colorValues <- sapply(seq(colourCount), function(x) rgb(rgbCodes[x, .(R,G,B) ],  maxColorValue = 255) )\n",
    "\n",
    "setorder(dt, Cluster)                      \n",
    "colorValues = setNames(colorValues, dt[, unique(Cluster)])\n",
    "colorValues           \n",
    "\n",
    "print(paste(\"Creating palette for\",colourCount,\"clusters\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:43.424320Z",
     "start_time": "2020-06-04T23:11:43.417006Z"
    }
   },
   "outputs": [],
   "source": [
    "# Builds a color dictionary from the cluster color scale built in R\n",
    "colorValDict = {}\n",
    "colorValDictRGB = {}\n",
    "for key,val in zip(clusterIDs, list(colorValues)):\n",
    "    colorValDict[key] = val\n",
    "\n",
    "for key,val in colorValDict.items():\n",
    "    colorValDictRGB[key] = tuple(int(val.lstrip('#')[i:i+2], 16) for i in (0, 2 ,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View community in the structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Connections to the ligand: \n",
    "First, we show all direct network connections between the ligand and surounding nodes from the active site. The next **two** jupyter notebook cell should be executed to show a sliding selection widget that will allow you to compare the network architectures of the acitve site in different simulation windows.\n",
    "\n",
    "Solid lines indicate connections between nodes in the same community, while dashed lines indicate connections between nodes of different communities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:43:25.251621Z",
     "start_time": "2020-05-29T04:43:24.528123Z"
    }
   },
   "outputs": [],
   "source": [
    "w = nv.show_mdanalysis(workUviz.select_atoms(\"all\"))\n",
    "\n",
    "w._remote_call(\"setSize\", target=\"Widget\", args=[\"800px\", \"600px\"])\n",
    "w.parameters = dict(theme='light')\n",
    "\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:43:27.314671Z",
     "start_time": "2020-05-29T04:43:27.291089Z"
    }
   },
   "outputs": [],
   "source": [
    "windowSlider = widgets.IntSlider(\n",
    "    value=refWindow, min=0, max=dnad.numWinds-1, step=1, description='Window:', readout_format='d',\n",
    "    disabled=False, continuous_update=False, orientation='horizontal', readout=True )\n",
    "\n",
    "def updateCommByTarget(value):\n",
    "    w.clear_representations()\n",
    "    w.add_cartoon(\"backbone\")\n",
    "    #w.add_hyperball()\n",
    "    showCommunityByTarget(w, nodeCommDF, nodesSite, windowSlider.value, dnad.nodesAtmSel, dnad, colorValDict)\n",
    "    \n",
    "windowSlider.observe(updateCommByTarget,  names='value')\n",
    "\n",
    "windowSlider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:43:40.943133Z",
     "start_time": "2020-05-29T04:43:40.935999Z"
    }
   },
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Entire communities:\n",
    "\n",
    "Second, we show all nodes from a chosen community. The following jupyter notebook cell different options to choose which community is shown. The first will select a community ID and show all nodes. The second allows you to select a set of nodes and show them in the structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:35.748681Z",
     "start_time": "2020-05-29T04:01:34.811521Z"
    }
   },
   "outputs": [],
   "source": [
    "w = nv.show_mdanalysis(workUviz.select_atoms(\"all\"))\n",
    "\n",
    "w._remote_call(\"setSize\", target=\"Widget\", args=[\"800px\", \"600px\"])\n",
    "w.parameters = dict(theme='light')\n",
    "\n",
    "opac = 1\n",
    "\n",
    "w.clear_representations()\n",
    "w.add_cartoon(\"protein\", color='lightgray', opacity=opac)\n",
    "#w.add_hyperball(ligandSegID)\n",
    "\n",
    "shapeCounter = [0]\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:36.095702Z",
     "start_time": "2020-05-29T04:01:35.751685Z"
    }
   },
   "outputs": [],
   "source": [
    "# Dysplays all nodes of a cluster and colors them by cluster ID.\n",
    "\n",
    "clusterID = 0\n",
    "showCommunityByID(w, cDF, clusterID, system1, refWindow, \n",
    "                  shapeCounter, dnad.nodesAtmSel, colorValDictRGB, system1, refWindow1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:36.755540Z",
     "start_time": "2020-05-29T04:01:36.098710Z"
    }
   },
   "outputs": [],
   "source": [
    "# Dysplays a given list of nodes and colors them by cluster ID.\n",
    "\n",
    "clusterID = 0\n",
    "nodeList = cDF.loc[ (cDF.system == system1) & (cDF.Window == refWindow1) & (cDF.Cluster == clusterID) ].Node.values\n",
    "\n",
    "showCommunityByNodes(w, cDF, nodeList, system1, refWindow1, shapeCounter, dnad.nodesAtmSel, colorValDictRGB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:44.454763Z",
     "start_time": "2020-05-29T04:01:36.758085Z"
    }
   },
   "outputs": [],
   "source": [
    "# Color ALL nodes by community\n",
    "\n",
    "nodeList = cDF.loc[ (cDF.system == system1) & (cDF.Window == refWindow1) ].Node.values\n",
    "showCommunityByNodes(w, cDF, nodeList, system1, refWindow1, shapeCounter, dnad.nodesAtmSel, colorValDictRGB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:01:44.462432Z",
     "start_time": "2020-05-29T04:01:44.457746Z"
    }
   },
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make more plots: Clustering of nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:51.340335Z",
     "start_time": "2020-06-04T23:11:51.335796Z"
    }
   },
   "outputs": [],
   "source": [
    "plotFilePrefix = \"test_\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:53.451259Z",
     "start_time": "2020-06-04T23:11:52.080711Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,plotFilePrefix -w 1200 -h 500\n",
    "\n",
    "dataPath = file.path(workDir, \"cluster.csv\")\n",
    "plotPath = file.path(workDir, paste0(\"Plots/\",plotFilePrefix,\"Clusters_Node_vs_Window.png\"))\n",
    "\n",
    "dt <- fread(dataPath)\n",
    "dt <- dt[,.(NodePlot,Window,Cluster)]\n",
    "dt <- dt[, Cluster := as.factor(Cluster) ]\n",
    "\n",
    "p <- ggplot(dt) + \n",
    "    geom_raster(aes(x=NodePlot, y=Window, fill=Cluster)) + \n",
    "    scale_fill_manual(values = colorValues) +\n",
    "    labs(x=\"Node\", y=\"Window\") +\n",
    "    theme_bw(base_size=20)\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:11:57.969450Z",
     "start_time": "2020-06-04T23:11:56.766480Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,nodeLabels,plotFilePrefix -w 2000 -h 500\n",
    "\n",
    "dataPath = file.path(workDir, \"cluster.csv\")\n",
    "plotPath = file.path(workDir, paste0(\"Plots/\",plotFilePrefix,\"Clusters_Node_vs_Window_Grouped.png\"))\n",
    "\n",
    "dt <- fread(dataPath)\n",
    "dt <- dt[, Cluster := as.factor(Cluster) ]\n",
    "\n",
    "setorder(dt, Cluster, Window)\n",
    "dt <- dt[, NodePlot := as.factor(NodePlot) ]\n",
    "\n",
    "# Get the actual indices of the nodes in the re-ordered x-axis (grouped by cluster)\n",
    "trgtIndices = data.table(trgt = c(which(dt[Window == refWindow]$Node %in% nodeLabels))  )\n",
    "\n",
    "breaksNodePlot = dt[Window == refWindow,][ Node %in% nodeLabels, NodePlot ]\n",
    "\n",
    "labelsNode = dt[Window == refWindow,][ Node %in% nodeLabels, Node ]\n",
    "\n",
    "# Build base plot\n",
    "p <- ggplot(dt) + \n",
    "    geom_raster(aes(x=NodePlot, y=Window, fill=Cluster)) \n",
    "\n",
    "# Add intercept lines for target nodes\n",
    "p <- p + geom_vline(data=trgtIndices, aes(xintercept=trgt), alpha=0.8, linetype = \"dashed\")\n",
    "p <- p + geom_hline(aes(yintercept=refWindow), alpha=0.9, linetype = \"dashed\")\n",
    "\n",
    "# Finish building plot\n",
    "p <- p + scale_fill_manual(values = colorValues) +\n",
    "    scale_x_discrete(limits=dt[Window == refWindow]$NodePlot, breaks=breaksNodePlot, labels=labelsNode) +\n",
    "    labs(x=\"Nodes\", y=\"Window\") +\n",
    "    theme_bw(base_size=20) +\n",
    "    theme(axis.text.x = element_text(angle = 90))\n",
    "\n",
    "#scale_x_discrete(limits=dt[Window == refWindow]$NodePlot, breaks=breaksNodePlot, labels=labelsNode) +\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:02.374531Z",
     "start_time": "2020-06-04T23:12:01.374346Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,nodesSite,contactNodesTrgts,plotFilePrefix -w 2000 -h 500\n",
    "\n",
    "dataPath = file.path(workDir, \"cluster.csv\")\n",
    "plotPath = file.path(workDir, paste0(\"Plots/\",plotFilePrefix,\"Clusters_Node_vs_Window_Grouped_ActiveSite.png\"))\n",
    "\n",
    "dt <- fread(dataPath)\n",
    "dt <- dt[,.(Node,Window,Cluster,resid)]\n",
    "dt <- dt[, Cluster := as.factor(Cluster) ]\n",
    "dt <- dt[, Node := as.factor(Node) ]\n",
    "setorder(dt, Cluster, Window)\n",
    "\n",
    "dt <- dt[ Node %in% contactNodesTrgts, ]\n",
    "\n",
    "dt <- dt[, resid := sapply(strsplit(resid, \"_\"), '[', 1) ]\n",
    "\n",
    "# Build base plot\n",
    "p <- ggplot(dt) + \n",
    "    geom_raster(aes(x=Node, y=Window, fill=Cluster)) \n",
    "\n",
    "# Get the actual indices of the nodes in the re-ordered x-axis (grouped by cluster)\n",
    "trgtIndices = data.table(trgt = c(which(dt[Window == refWindow]$Node %in% nodesSite))  )\n",
    "\n",
    "# Add intercept lines for target nodes\n",
    "p <- p + geom_vline(data=trgtIndices, aes(xintercept=trgt), alpha=0.9, linetype = \"dashed\")\n",
    "p <- p + geom_hline(aes(yintercept=refWindow), alpha=0.9, linetype = \"dashed\")\n",
    "\n",
    "# Finish building plot\n",
    "p <- p + scale_fill_manual(values = colorValues) +\n",
    "    scale_x_discrete(limits=dt[Window == refWindow]$Node, label=dt[Window == refWindow]$resid) +\n",
    "    labs(x=\"Nodes\", y=\"Window\") +\n",
    "    theme_bw(base_size=20) + \n",
    "    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:12.966612Z",
     "start_time": "2020-06-04T23:12:12.254606Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,headNodeSite,nodesSite,plotFilePrefix,ligandSegID -w 600 -h 450\n",
    "\n",
    "headNodeSite = as.integer(headNodeSite)\n",
    "\n",
    "dataPath = file.path(workDir, \"cluster.csv\")\n",
    "plotPath = file.path(workDir, paste0(\"Plots/\",plotFilePrefix,\"Clusters_Node_vs_Window_Grouped_\",ligandSegID,\".png\"))\n",
    "\n",
    "dt <- fread(dataPath)\n",
    "dt <- dt[,.(Node,Window,Cluster,resid)]\n",
    "dt <- dt[ Node %in% nodesSite, ]\n",
    "dt <- dt[, Cluster := as.factor(Cluster) ]\n",
    "dt <- dt[, Node := as.factor(Node) ]\n",
    "dt <- dt[, Window := as.factor(Window) ]\n",
    "dt <- dt[, resid := sapply(strsplit(resid, \"_\"), '[', 1) ]\n",
    "dt <- dt[, resid := as.factor(resid) ]\n",
    "\n",
    "setorder(dt, Cluster, Node)\n",
    "\n",
    "# Build base plot\n",
    "p <- ggplot(dt) + \n",
    "    geom_raster(aes(x=Node, y=Window, fill=Cluster))\n",
    "\n",
    "# Finish building plot\n",
    "p <- p + scale_fill_manual(values = colorValues) +\n",
    "    scale_x_discrete(limits=as.factor(nodesSite), label=dt[ Node %in% as.factor(nodesSite), unique(resid)]) +\n",
    "    scale_y_discrete(limits=dt[Node == headNodeSite]$Window, labels=NULL) +\n",
    "    labs(x=\"Nodes\", y=\"Window\") +\n",
    "    theme_bw(base_size=20) + \n",
    "    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Combine data for Interface and Active Site analysis\n",
    "\n",
    "Combines interface edges with label information for plots.\n",
    "\n",
    "Creates structure wirh optimal paths and labels for Active Site analysis, shaping information for interaction between ligand and active site residues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:22.083358Z",
     "start_time": "2020-06-04T23:12:21.799500Z"
    }
   },
   "outputs": [],
   "source": [
    "window = 0\n",
    "\n",
    "# Gets all pairs of nodes with non-zero correlations\n",
    "# nonzeroPairs = [(i,j) for i,j in np.asarray(np.where(dnad.corrMatAll[window, :, :] > 0)).T if i < j]\n",
    "nonzeroPairs = [(i,j) for i,j in np.asarray(np.where(dnad.contactMat > 0)).T if i < j]\n",
    "\n",
    "# Combines cartesian distance and correlation\n",
    "# We use the mean cartesian distance (0:mean, 1:SEM, 2:Min, 3:Max)\n",
    "distType = 0\n",
    "# carCorMat = [ [i, j, \n",
    "#                getCartDist(i,j,dnad.numNodes,dnad.nodeDists, distType), \n",
    "#                np.mean(dnad.corrMatAll[window,i,j]) ] \n",
    "#              for i,j in nonzeroPairs ]\n",
    "carCorMat = [ [i, j, \n",
    "               getCartDist(i,j,dnad.numNodes,dnad.nodeDists, 0), \n",
    "               getCartDist(i,j,dnad.numNodes,dnad.nodeDists, 1), \n",
    "               np.mean(dnad.corrMatAll[:,i,j]),\n",
    "               stats.sem(dnad.corrMatAll[:,i,j])] \n",
    "             for i,j in nonzeroPairs if np.mean(dnad.corrMatAll[:,i,j]) > 0 ]\n",
    "\n",
    "carCorMat = pd.DataFrame(carCorMat, columns=[\"i\",\"j\",\"Cart\",\"CartSEM\",\"Corr\",\"CorrSEM\"])\n",
    "\n",
    "def interNode(i,j,dnad):\n",
    "    # Checks if the pair of nodes exists in the \"interNodePairs\" 2D array.\n",
    "    return len( ( dnad.interNodePairs == [i, j] ).all(axis=1).nonzero()[0] )\n",
    "\n",
    "# Adds interface information (true/false)\n",
    "carCorMat[\"EdgeType\"] = carCorMat.apply( lambda x: \"Interface\" if interNode(x[\"i\"], x[\"j\"], dnad) \n",
    "                                        else \"Internal\" , axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:25.321625Z",
     "start_time": "2020-06-04T23:12:24.303655Z"
    }
   },
   "outputs": [],
   "source": [
    "# Initialized data frame with pairs of source node (all nodes in the system) and target node (defined above).\n",
    "\n",
    "dataTmp = []\n",
    "\n",
    "window = 0\n",
    "for window in range(dnad.numWinds):\n",
    "    \n",
    "    df = pd.DataFrame(np.asarray([ (i,trgt) for trgt in nodesSite  for i in range(dnad.numNodes)]),\n",
    "                 columns=[\"node\", \"targets\"])\n",
    "    \n",
    "    df[\"distances\"]  = df.apply( lambda row: dnad.distsAll[window, row[\"node\"], row[\"targets\"]], axis=1)\n",
    "    # Selects the mean distance (0:mean, 1:SEM, 2:Min, 3:Max)\n",
    "    distType = 0 \n",
    "    df[\"cdistances\"] = df.apply( lambda row: getCartDist(row[\"node\"], row[\"targets\"], dnad.numNodes, \n",
    "                                                         dnad.nodeDists, distType), axis=1)\n",
    "    # Selects the standard error of the mean distance (0:mean, 1:SEM, 2:Min, 3:Max)\n",
    "    distType = 1\n",
    "    df[\"cdistSEM\"] = df.apply( lambda row: getCartDist(row[\"node\"], row[\"targets\"], dnad.numNodes, \n",
    "                                                         dnad.nodeDists, distType), axis=1)\n",
    "    df[\"path\"]       = df.apply( lambda row: list(getPath(row[\"node\"], row[\"targets\"], \n",
    "                                                          dnad.nodesAtmSel, dnad.preds, win= window)), axis=1)\n",
    "    df['path_lens']  = df.apply( lambda row: len(row[\"path\"]), axis=1)\n",
    "#     df[\"mincdist\"]   = df.groupby(\"node\")[[\"cdistances\"]].transform(lambda x: np.min(x) )\n",
    "\n",
    "    def getTagStr(i):\n",
    "        return dnad.nodesAtmSel.atoms[i].resname.capitalize() + \":\" + str(dnad.nodesAtmSel.atoms[i].resid) + \\\n",
    "                \"_\" + dnad.nodesAtmSel.atoms[i].segid\n",
    "\n",
    "    df['resid']     = np.vectorize(getTagStr)(df[\"node\"])\n",
    "\n",
    "    ptnIXs = dnad.nodesAtmSel.select_atoms(\"protein\").ix_array\n",
    "    nclIXs = dnad.nodesAtmSel.select_atoms(\"nucleic\").ix_array\n",
    "    def getType(nodeIndx):\n",
    "\n",
    "        if dnad.nodesAtmSel.atoms[nodeIndx].ix in ptnIXs:\n",
    "            return \"Aminoacid\"\n",
    "\n",
    "        elif dnad.nodesAtmSel.atoms[nodeIndx].ix in nclIXs:\n",
    "            return \"Nucleotide\"\n",
    "\n",
    "        else:\n",
    "            return dnad.nodesAtmSel.atoms[nodeIndx].resname\n",
    "\n",
    "    df[\"type\"] = np.vectorize(getType)(df[\"node\"])\n",
    "\n",
    "    df[\"Interface\"] = df[\"node\"].apply( lambda x: x in dnad.contactNodesInter)\n",
    "\n",
    "    # \"Cleans\" all -1 distances betweeen nodes not connected by any path.\n",
    "    df.loc[ df[\"path_lens\"] == 0, \"distances\" ] = 0\n",
    "    \n",
    "    df[\"window\"] = window\n",
    "    \n",
    "    dataTmp.append(df.copy())\n",
    "    del df\n",
    "    \n",
    "df = pd.concat(dataTmp)\n",
    "del dataTmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:26.253576Z",
     "start_time": "2020-06-04T23:12:26.230770Z"
    }
   },
   "outputs": [],
   "source": [
    "# df.head()\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:26.817921Z",
     "start_time": "2020-06-04T23:12:26.791442Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df[ (np.isin(df[\"node\"],getNodeFromSel(\"segid \" + ligandSegID, dnad.nodesAtmSel, dnad.atomToNode))) & (df[\"window\"] == 0) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Add labels to selected residues\n",
    "\n",
    "Using references form the literature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:31.608634Z",
     "start_time": "2020-06-04T23:12:31.590262Z"
    }
   },
   "outputs": [],
   "source": [
    "df[\"Label\"] = 0\n",
    "\n",
    "tags = dict()\n",
    "#\"resname H2U and element N\"\n",
    "tags[\"OMPn\"] = getNodeFromSel(\"resname OMP and name N1\", dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "#\"resname H2U and element P\"\n",
    "tags[\"OMPp\"] = getNodeFromSel(\"resname OMP and name P\", dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "\n",
    "for label, node in tags.items():\n",
    "    if not node in df.node:\n",
    "        print(\"Skiping node not found:\",node)\n",
    "        continue\n",
    "    df.loc[ df.node == node, 'Label'] = label\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:32.687560Z",
     "start_time": "2020-06-04T23:12:32.663679Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df.loc[ (df[\"resid\"] == \"Omp:301_OMP\") & (df[\"window\"] == 0)  ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adds label and type information for all pairs of connected residues in the system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:44.755637Z",
     "start_time": "2020-06-04T23:12:36.374275Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "carCorMat[\"iRes\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"i\"]) & \n",
    "                               (df[\"targets\"] == headNodeSite[0])][\"resid\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"jRes\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"j\"]) & \n",
    "                               (df[\"targets\"] == headNodeSite[0])][\"resid\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"iType\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"i\"]) & \n",
    "                               (df[\"targets\"] == headNodeSite[0])][\"type\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"jType\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"j\"]) & \n",
    "                               (df[\"targets\"] == headNodeSite[0])][\"type\"].iloc[0] , axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:44.776217Z",
     "start_time": "2020-06-04T23:12:44.758762Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "carCorMat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:44.867790Z",
     "start_time": "2020-06-04T23:12:44.779617Z"
    }
   },
   "outputs": [],
   "source": [
    "carCorMat.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:44.968226Z",
     "start_time": "2020-06-04T23:12:44.870894Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a copy only with interface connections, and check that i < j\n",
    "carCorMatInterface = carCorMat[ carCorMat[\"EdgeType\"] == \"Interface\" ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare Intra and Inter segment correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:12:51.480310Z",
     "start_time": "2020-06-04T23:12:48.669368Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i carCorMat,workDir -w 800 -h 500\n",
    "\n",
    "plotPath = file.path(workDir, \"Plots/Interf_Intern_Cart_vs_Corr.png\")\n",
    "\n",
    "cartCutoff = 4.0\n",
    "\n",
    "dt = data.table(carCorMat)\n",
    "dt  = dt[Cart < cartCutoff]\n",
    "\n",
    "p <- ggplot(dt) + \n",
    "    geom_point( aes(x=Cart, y=Corr, color=EdgeType), alpha=0.7, size=3 ) +\n",
    "    geom_smooth( aes(x=Cart, y=Corr, color=EdgeType) ) +\n",
    "    labs(x=\"Mean Cartesian Distance (A)\", y=\"Mean Correlation\", color=\"Edge Type\") +\n",
    "    scale_color_brewer(type='qual', palette=6) +\n",
    "    theme_linedraw(base_size=20) + xlim(c(2.5, cartCutoff))\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare network connectivity, network distance and cartesian distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:13:02.246460Z",
     "start_time": "2020-06-04T23:13:01.883625Z"
    }
   },
   "outputs": [],
   "source": [
    "# Select all nodes involved in at least one \"Interface\" connection.\n",
    "tmpDF = carCorMat[ carCorMat[\"EdgeType\"] == \"Interface\" ]\n",
    "\n",
    "tmpNodeSet = set(tmpDF[\"i\"])\n",
    "tmpNodeSet.update( set(tmpDF[\"j\"]) )\n",
    "\n",
    "nodeContacs = []\n",
    "for node in tmpNodeSet:\n",
    "\n",
    "    tmp = tmpDF.loc[  (tmpDF[\"i\"] == node) | (tmpDF[\"j\"] == node) ][\"Corr\"]\n",
    "    tmp2 = tmpDF.loc[  (tmpDF[\"i\"] == node) | (tmpDF[\"j\"] == node) ][\"Cart\"]\n",
    "    \n",
    "    label = df.loc[ df[\"node\"] == node ][\"resid\"].unique()[0]\n",
    "    \n",
    "    nodeContacs.append( [node, int(tmp.size), \n",
    "                         tmp.mean(), sp.stats.sem( tmp ),\n",
    "                         tmp2.mean(), sp.stats.sem( tmp2 ),\n",
    "                         label] )\n",
    "\n",
    "nodeContacs = pd.DataFrame(nodeContacs, \n",
    "                           columns=[\"Node\", \"NumContacts\", \"MeanCorr\", \"SEMCorr\", \n",
    "                                    \"MeanCart\", \"SEMCart\", \"label\"])\n",
    "nodeContacs = nodeContacs.fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:13:11.069484Z",
     "start_time": "2020-06-04T23:13:10.443508Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i nodeContacs -w 800 -h 600\n",
    "\n",
    "dt <- data.table(nodeContacs)\n",
    "sapply(dt, class)\n",
    "\n",
    "ggplot(dt) + \n",
    "#     geom_point(aes(x=NumContacts, y=MeanCorr, color=MeanCart), size=3) +\n",
    "    geom_point(aes(x=NumContacts, y=MeanCart, color=as.double(MeanCorr)), size=3) +\n",
    "    geom_label_repel(data=dt[NumContacts >= 10 | MeanCart <= 3 | MeanCorr > 0.8 ], \n",
    "                    aes(x=NumContacts, y=MeanCart, label=label)) + \n",
    "    scale_y_log10() +\n",
    "    labs(x=\"Number of Contacts\", \n",
    "         y=\"Mean Cartesian Distance\",\n",
    "         color=\"Mean Correlation\") +\n",
    "    scale_color_gradient(low=\"blue\",high=\"red\") + \n",
    "    theme_linedraw(base_size=20)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Active site connections\n",
    "Here we visualize the residues that make direct connections to the ligand, and display their generalized correlation coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-06-04T23:17:07.975264Z",
     "start_time": "2020-06-04T23:17:06.982510Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i carCorMatInterface,workDir,ligandSegID -w 2000 -h 500\n",
    "\n",
    "plotPath = file.path(workDir, \"Plots/Interface_AAres_vs_Corr_NoLabel.png\")\n",
    "\n",
    "corrCutoff = 0.0\n",
    "\n",
    "dt <- data.table(carCorMatInterface)\n",
    "\n",
    "iResList <- dt[iType == \"Aminoacid\" & jType != \"SOL\"][Corr > corrCutoff][,iRes]\n",
    "dt <- dt[ iRes %in% iResList, ]\n",
    "dt <- dt[jType != \"SOL\", .SD[which.max(Corr)] , by=.(iRes,jRes)]\n",
    "\n",
    "dt <- dt[, iRes := sapply(strsplit(iRes, \"_\"), '[', 1) ]\n",
    "\n",
    "colours <- c(\"Aminoacid\" = \"red\", \n",
    "          \"SOL\" = \"blue\", \n",
    "          \"Nucleotide\" = \"darkorange\")\n",
    "colours <- c(colours, setNames(\"purple\",ligandSegID))\n",
    "\n",
    "shapes <- c(\"Aminoacid\" = 19,\n",
    "          \"SOL\" = 15, \n",
    "          \"Nucleotide\" = 17)\n",
    "shapes <- c(shapes, setNames(19,ligandSegID))\n",
    "\n",
    "p <- ggplot(dt) + \n",
    "    geom_linerange( aes(x=iRes, y=Corr, ymin=Corr-CorrSEM, ymax=Corr+CorrSEM), size=1 )+\n",
    "    geom_point( aes(x = iRes, y = Corr, color = jType, size=Cart) )+\n",
    "    labs(x=\"Aminoacid Residue\", y=\"Mean Correlation\", color=\"Residue Type\", size=\"Mean Distance\") + \n",
    "    scale_colour_manual(name = \"Residue Type\",\n",
    "                     values = colours) +\n",
    "    guides(colour = guide_legend(override.aes = list(size=5))) +\n",
    "    theme_classic(base_size = 20) +\n",
    "    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),\n",
    "          panel.grid.major.x = element_line(color=\"black\", size=0.5, linetype=\"dotted\"))\n",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:08.127424Z",
     "start_time": "2020-05-29T04:02:07.299339Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "w = nv.show_mdanalysis(workUviz.select_atoms(\"all\"))\n",
    "w._remote_call(\"setSize\", target=\"Widget\", args=[\"800px\", \"600px\"])\n",
    "w.parameters = dict(theme='light')\n",
    "w\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:08.882463Z",
     "start_time": "2020-05-29T04:02:08.129756Z"
    }
   },
   "outputs": [],
   "source": [
    "# Paint all connections to the ligand\n",
    "window = 1\n",
    "\n",
    "set_widget_representations(w)\n",
    "\n",
    "tmpDF = carCorMat[(carCorMat[\"i\"].isin(nodesSite)) | (carCorMat[\"j\"].isin(nodesSite))]\n",
    "tmpDF = tmpDF.loc[carCorMat[\"EdgeType\"] == \"Interface\"]\n",
    "\n",
    "for pair in tmpDF[[\"i\",\"j\"]].values:\n",
    "    viewPath(w, getPath(pair[0], pair[1], dnad.nodesAtmSel, dnad.preds), \n",
    "             dnad.distsAll, dnad.maxDirectDist, dnad.nodesAtmSel, window)\n",
    "    w.add_ball_and_stick(getNGLSelFromNode(pair[0], dnad.nodesAtmSel, atom=True))\n",
    "    w.add_ball_and_stick(getNGLSelFromNode(pair[1], dnad.nodesAtmSel, atom=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:08.889747Z",
     "start_time": "2020-05-29T04:02:08.885856Z"
    }
   },
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pretty Figures - Preparing files for VMD\n",
    "\n",
    "The next few steps will define the files that will be created for loading into VMD, a popular molecular visualization software.\n",
    "\n",
    "**Important:**\n",
    "\n",
    "Make sure you have all the files created for VMD. For instance, if you are not interested in allosteric communications and suboptimal paths, just select two arbitrary nodes so that the path files are created correctly. Alternatively, you may want to adapt the tcl scripts provided at the end of the notebook to your particular case. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:12.717797Z",
     "start_time": "2020-05-29T04:02:08.892077Z"
    }
   },
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "# This will be used to look for the maximum and  minimum betweenness value in the graph.\n",
    "# The maximum value will be used ot normalize all betweenness values for better vizualization.\n",
    "# The minimum value will be used in case a betweenness value could not be assigned for a given edge,\n",
    "#    also helping visualization.\n",
    "\n",
    "# Initialize variable with high value.\n",
    "minimumBetweeness = 100\n",
    "# Initialize variable with low value.\n",
    "maximumBetweeness = -1\n",
    "\n",
    "for pair,btw in dnad.btws[winIndx].items():\n",
    "    if btw < minimumBetweeness:\n",
    "            minimumBetweeness = btw\n",
    "    if btw > maximumBetweeness:\n",
    "            maximumBetweeness = btw\n",
    "\n",
    "# Normalize the value.\n",
    "minimumBetweeness /= maximumBetweeness\n",
    "            \n",
    "for winIndx in range(dnad.numWinds):\n",
    "\n",
    "    normCorMat = copy.deepcopy( dnad.corrMatAll[winIndx,:,:] )\n",
    "    normCorMat /= normCorMat.max()\n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create PDB file with the system in the first step of each window, for VMD vizualization.\n",
    "    \n",
    "    pdbVizFile = os.path.join(workDir, \n",
    "                            \"networkData_Structure_window_{}.pdb\".format(winIndx))\n",
    "    \n",
    "    # Calculate number of frames per window.\n",
    "    winLen = int(np.floor(workUviz.trajectory.n_frames/dnad.numWinds))\n",
    "    \n",
    "    # Positions the trajectory at the middle of each window.\n",
    "    workUviz.trajectory[(winIndx+1)*round(winLen/2)]\n",
    "    \n",
    "    with mda.Writer(pdbVizFile, multiframe=False, bonds=\"conect\", n_atoms=workUviz.atoms.n_atoms) as PDB:\n",
    "        PDB.write(workUviz.atoms)\n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create network data file with ALL edges and their normalized weights.\n",
    "    \n",
    "    fileName = os.path.join(workDir, \n",
    "                            \"networkData_AllEdges_window_{}.dat\".format(winIndx))\n",
    "    with open(fileName, \"w\") as outfile:\n",
    "\n",
    "        for pair in np.asarray( np.where( np.triu(normCorMat[:,:]) ) ).T:\n",
    "\n",
    "            node1 = pair[0]\n",
    "            node2 = pair[1]\n",
    "\n",
    "            # Get VMD indices for the atoms\n",
    "            pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "            pdbIndx2 = dnad.nodesAtmSel.atoms[node2].id -1\n",
    "\n",
    "            string = \"{} {} {}\".format(pdbIndx1, pdbIndx2, normCorMat[ node1, node2])\n",
    "            \n",
    "            outfile.write( string + \"\\n\" )\n",
    "    \n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create network data file with ALL NODES, the maximum normalized weight of edges it belongs to,\n",
    "    ### and the community it belongs to.\n",
    "    \n",
    "    fileName = os.path.join(workDir, \n",
    "                            \"networkData_AllNodes_window_{}.dat\".format(winIndx))\n",
    "    with open(fileName, \"w\") as outfile:\n",
    "        \n",
    "        for node1 in range(dnad.numNodes):\n",
    "            \n",
    "            # Get the VMD index for the atom\n",
    "            pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "            \n",
    "            # Get the community the node belongs to\n",
    "            community1 = int(nodeCommNP[node1, winIndx])\n",
    "            \n",
    "            # Find the node to which \"node1\" is connected with highest correlation.\n",
    "            node2 = np.where( normCorMat[node1,:] == normCorMat[node1,:].max() )[0][0]\n",
    "            \n",
    "            # Skip nodes not assigned to any community\n",
    "            if community1 < 0:\n",
    "                continue\n",
    "            \n",
    "            string = \"{} {} {}\".format(pdbIndx1, normCorMat[ node1, node2], community1)\n",
    "            \n",
    "            outfile.write( string + \"\\n\" )\n",
    "    \n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create network data file with INTRA-COMMUNITY edges and their normalized weights.\n",
    "    \n",
    "    fileName = os.path.join(workDir, \n",
    "                            \"networkData_IntraCommunities_window_{}.dat\".format(winIndx))\n",
    "    with open(fileName, \"w\") as outfile:\n",
    "        \n",
    "        for pair in np.asarray( np.where( np.triu(normCorMat[:,:]) ) ).T:\n",
    "            \n",
    "            node1 = pair[0]\n",
    "            node2 = pair[1]\n",
    "\n",
    "            # Checks if both nodes belong to the same community. \n",
    "            # If they don't, skip this edge. We only write intra-community edges in this file!\n",
    "            if nodeCommNP[node1, winIndx] != nodeCommNP[node2, winIndx] :\n",
    "                continue\n",
    "            \n",
    "            # If both nodes do not belong to any community (assigned to community -1), also skip the edge.\n",
    "            if nodeCommNP[node1, winIndx] < 0:\n",
    "                continue\n",
    "            \n",
    "            community1 = int(nodeCommNP[node1, winIndx])\n",
    "            \n",
    "            # Get VMD indices for the atoms\n",
    "            pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "            pdbIndx2 = dnad.nodesAtmSel.atoms[node2].id -1\n",
    "\n",
    "            string = \"{} {} {} {}\".format(pdbIndx1, pdbIndx2, normCorMat[ node1, node2], community1)\n",
    "            \n",
    "            outfile.write( string + \"\\n\" )\n",
    "    \n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create network data file with INTER-COMMUNITY edges and their normalized weights.\n",
    "    \n",
    "    fileName = os.path.join(workDir, \n",
    "                                \"networkData_InterCommunities_window_{}.dat\".format(winIndx))\n",
    "    with open(fileName, \"w\") as outfile:\n",
    "        \n",
    "        for pair in np.asarray( np.where( np.triu(normCorMat[:,:]) ) ).T:\n",
    "            \n",
    "            node1 = pair[0]\n",
    "            node2 = pair[1]\n",
    "\n",
    "            # Checks if both nodes belong to the same community. \n",
    "            # If they don't, skip this edge. We only write intra-community edges in this file!\n",
    "            if nodeCommNP[node1, winIndx] == nodeCommNP[node2, winIndx] :\n",
    "                continue\n",
    "            \n",
    "            # If either node does not belong to any community (assigned to community -1), also skip the edge.\n",
    "            if (nodeCommNP[node1, winIndx] < 0) or (nodeCommNP[node2, winIndx] < 0):\n",
    "                continue\n",
    "            \n",
    "            community1 = int(nodeCommNP[node1, winIndx])\n",
    "            community2 = int(nodeCommNP[node2, winIndx])\n",
    "            \n",
    "            # Get VMD indices for the atoms\n",
    "            # VMD uses a 0-based index, so we subtract 1 from the PDB index\n",
    "            pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "            pdbIndx2 = dnad.nodesAtmSel.atoms[node2].id -1\n",
    "\n",
    "            string = \"{} {} {} {} {}\".format(pdbIndx1, pdbIndx2, \n",
    "                                             normCorMat[ node1, node2], community1, community2)\n",
    "            \n",
    "            outfile.write( string + \"\\n\" )\n",
    "    \n",
    "    \n",
    "    ##########################################################################################\n",
    "    ### Create file with edges listed by betweeness value (highest to lowest).\n",
    "    \n",
    "    fileName = os.path.join(workDir, \n",
    "                            \"networkData_Betweenness_window_{}.dat\".format(winIndx))\n",
    "    with open(fileName, \"w\") as outfile:\n",
    "        \n",
    "        for pair,btw in dnad.btws[winIndx].items():\n",
    "            \n",
    "            node1 = pair[0]\n",
    "            node2 = pair[1]\n",
    "\n",
    "            # If either node does not belong to any community (assigned to community -1), also skip the edge.\n",
    "            if (nodeCommNP[node1, winIndx] < 0) or (nodeCommNP[node2, winIndx] < 0):\n",
    "                continue\n",
    "            \n",
    "            community1 = int(nodeCommNP[node1, winIndx])\n",
    "            community2 = int(nodeCommNP[node2, winIndx])\n",
    "            \n",
    "            # Get VMD indices for the atoms\n",
    "            # VMD uses a 0-based index, so we subtract 1 from the PDB index\n",
    "            pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "            pdbIndx2 = dnad.nodesAtmSel.atoms[node2].id -1\n",
    "\n",
    "            string = \"{} {} {} {} {} {}\".format(pdbIndx1, pdbIndx2, \n",
    "                                                normCorMat[ node1, node2], btw/maximumBetweeness, \n",
    "                                                community1, community2)\n",
    "            \n",
    "            outfile.write( string + \"\\n\" )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write Optimal and Sub-Optimal Paths\n",
    "\n",
    "Using the convenience functions \"getNodeFromSel\" and \"getSelFromNode\", one can easily probe the system and determine the relationship between node in the network graph and the atoms and residues they represent in the actual system.\n",
    "\n",
    "- getSelFromNode\n",
    "- getNodeFromSel\n",
    "\n",
    "See examples below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:12.730219Z",
     "start_time": "2020-05-29T04:02:12.720388Z"
    }
   },
   "outputs": [],
   "source": [
    "#srcNode = getNodeFromSel(\"resname PHE and resid 508\",dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "#srcNode = getNodeFromSel(\"resname TRP and resid 496\",dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "#print(\"Source node:\", srcNode)\n",
    "\n",
    "#trgNode = getNodeFromSel(\"resname ARG and resid 352\",dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "#print(\"Target node:\", trgNode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create a list of important paths:**\n",
    "\n",
    "Once you have chosen the nodes that define each path of interest, create a list in the cell below with the indices of the source and target nodes.\n",
    "\n",
    "To make sure that the VMD scripts will run without the need to be adapted, the user must select at least one pair of nodes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:12.959829Z",
     "start_time": "2020-05-29T04:02:12.852749Z"
    }
   },
   "outputs": [],
   "source": [
    "# For example, to write the paths between node 0 (Valine 11) and nodes 58 and 60, the following\n",
    "#   list must be created:  \n",
    "# nodesForPaths = [ [0,58], [0,60] ]\n",
    "#\n",
    "nodesForPaths = [ [*nodePathA,*nodePathB] ]\n",
    "nodesForPaths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-29T04:02:13.365287Z",
     "start_time": "2020-05-29T04:02:12.962583Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Determine how many extra sub-optimal paths will be written.\n",
    "numSuboptimalPaths = 5\n",
    "\n",
    "pathListFile = open(os.path.join(workDir, \"paths.list\"), \"w\")\n",
    "\n",
    "for srcNode, trgNode in nodesForPaths:\n",
    "    \n",
    "    tmpList = getSelFromNode(srcNode,dnad.nodesAtmSel, atom=True).split()\n",
    "    srcNodeSel = \"\".join([tmpList[1],tmpList[4],tmpList[10]])\n",
    "\n",
    "    tmpList = getSelFromNode(trgNode,dnad.nodesAtmSel, atom=True).split()\n",
    "    trgNodeSel = \"\".join([tmpList[1],tmpList[4],tmpList[10]])\n",
    "    \n",
    "    # Adds the path suffix to the file\n",
    "    pathListFile.write(\"_{}_{}\\n\".format(srcNodeSel, trgNodeSel))\n",
    "    \n",
    "    for winIndx in range(dnad.numWinds):\n",
    "\n",
    "        normCorMat = copy.deepcopy( dnad.corrMatAll[winIndx,:,:] )\n",
    "        normCorMat /= normCorMat.max()\n",
    "\n",
    "        ##########################################################################################\n",
    "        ### Create file with edges listed by betweeness value (highest to lowest).\n",
    "        \n",
    "        # File name is created based on selections, not node index, for readability.\n",
    "        \n",
    "        \n",
    "        \n",
    "        fileName = os.path.join(workDir, \n",
    "                                    \"networkData_Paths_window_{}_{}_{}.dat\".format(winIndx, \n",
    "                                                                                srcNodeSel, trgNodeSel))\n",
    "        with open(fileName, \"w\") as outfile:\n",
    "\n",
    "            allPaths = []\n",
    "\n",
    "            # Reconstructs the optimal path from Floyd-Warshall algorithm\n",
    "            pathFW = nx.reconstruct_path(srcNode, trgNode, dnad.preds[winIndx])\n",
    "\n",
    "            allPaths.append(pathFW)\n",
    "\n",
    "            # Behind the scenes, use Dijkstra algorithm to find sub-optimal paths\n",
    "            for pathSO in islice(nx.shortest_simple_paths(dnad.nxGraphs[0], \n",
    "                                                srcNode, trgNode, weight=\"dist\"), 1, numSuboptimalPaths + 1):\n",
    "                allPaths.append(pathSO)\n",
    "\n",
    "            # Create a counter of number of paths that go though each edge, among all (sub-)optimal path(s).\n",
    "            pathCounter = defaultdict(int)\n",
    "            for pathIndx, pathIter in enumerate(allPaths):\n",
    "                # Iterate over edges in the path\n",
    "                for i in range(len(pathIter)-1):\n",
    "\n",
    "                    node1 = pathIter[i]\n",
    "                    node2 = pathIter[i+1]\n",
    "\n",
    "                    pathCounter[(node1, node2)] += 1\n",
    "\n",
    "            # Normalize the count\n",
    "            maxCount = np.max(list(pathCounter.values()))\n",
    "            for pair, count in pathCounter.items():\n",
    "                pathCounter[pair] = count/maxCount\n",
    "\n",
    "            for pathIndx, pathIter in enumerate(allPaths):\n",
    "                # Iterate over edges in the path\n",
    "                for i in range(len(pathIter)-1):\n",
    "\n",
    "                    node1 = pathIter[i]\n",
    "                    node2 = pathIter[i+1]\n",
    "\n",
    "                    # Get the community each node belongs to\n",
    "                    community1 = int(nodeCommNP[node1, winIndx])\n",
    "                    community2 = int(nodeCommNP[node2, winIndx])\n",
    "\n",
    "                    # If either node does not belong to any community (assigned to community -1), \n",
    "                    #     also skip the edge.\n",
    "                    if (community1 < 0) or (community2 < 0):\n",
    "                        continue\n",
    "\n",
    "                    # Get the betweeness value\n",
    "                    try:\n",
    "                        btw = dnad.btws[winIndx][( node1, node2)]\n",
    "                    except:\n",
    "                        # If one could not be calculated (very few paths going though this edge)\n",
    "                        # set an arbitrarily low value.\n",
    "                        btw = minimumBetweeness\n",
    "\n",
    "                    # Get VMD indices for the atoms\n",
    "                    # VMD uses a 0-based index, so we subtract 1 from the PDB index\n",
    "                    pdbIndx1 = dnad.nodesAtmSel.atoms[node1].id -1\n",
    "                    pdbIndx2 = dnad.nodesAtmSel.atoms[node2].id -1\n",
    "\n",
    "                    string = \"{} {} {} {} {} {}\".format(pdbIndx1, pdbIndx2, \n",
    "                                                     normCorMat[ node1, node2], \n",
    "                                                     btw/maximumBetweeness, pathCounter[(node1, node2)], \n",
    "                                                     pathIndx)\n",
    "\n",
    "                    outfile.write( string + \"\\n\" )\n",
    "\n",
    "pathListFile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write TCL files to load all data we just wrote down and visualize in VMD.\n",
    "prepTclViz(\"networkData\", str(dnad.numWinds), ligandSegID, workDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write PDB files for the last frame of each window rather than the first.\n",
    "selection = workUviz.select_atoms(\"all\")\n",
    "print(\"Loaded trajectory is %i frames divided into %i windows.\" % (workUviz.trajectory.n_frames, dnad.numWinds))\n",
    "frameList=np.arange(workUviz.trajectory.n_frames/dnad.numWinds-1,\\\n",
    "          workUviz.trajectory.n_frames,\\\n",
    "          workUviz.trajectory.n_frames/dnad.numWinds,\\\n",
    "          dtype=int)\n",
    "for i,frame in enumerate(frameList):\n",
    "    fileName=\"%s/networkData_Structure_window_%i_last.pdb\" % (workDir,i)\n",
    "    selection.write( fileName, frames=workUviz.trajectory[[frame]] )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rendering high-quality images with VMD - Step 3\n",
    "\n",
    "In the third step, the user can load files produced by the Jupyter notebooks into VMD. An easy to load script will handle all the work, creating a simple graphical  user  interface  (GUI)  that can be used to  easily  render publication-quality  images.   These  renderings  can  represent many aspects of the biomolecular system.\n",
    "\n",
    "There are *two options* to load the system into VMD and open the GUI:\n",
    "\n",
    "**1. From inside VMD.**\n",
    "\n",
    "    After opening VMD, in the VMD Main window, go to: Extensions > Tk Console. \n",
    "    \n",
    "    In the Tk Console, type the following command to navigate to the folder where the Analysis Results were saved (same as in the beginning of this tutorial Step 2): \n",
    "    \n",
    "`cd << path to folder >>`\n",
    "    \n",
    "    In the results folder, type the following command to load the results and the GUI: \n",
    "    \n",
    "`source network_view_2.tcl`\n",
    "\n",
    "**2. When loading VMD.**\n",
    "\n",
    "    If you choose to load VMD from a terminal (command line) window, navigate to the folder where the Analysis Results were saved (same as in the beginning of this tutorial Step 2): \n",
    "    \n",
    "`cd << path to folder >>`\n",
    "    \n",
    "    Load VMD with the following command: \n",
    "    \n",
    "`vmd -e network_view_2.tcl`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to use the Network View 2 GUI\n",
    "\n",
    "The Network View 2.0 GUI was created in a way that allows for easy interaction without deep background knowledge on VMD. If you are an expert VMD user, you can still change anything in the representation, but the GUI might erase your changes when loading some of the specialized features.\n",
    "\n",
    "The GUI allows user to visualize and render all the properties of the network presented in the Step 2. For instance:\n",
    "\n",
    "**1.** To visualize the communities just click on \"All Communities\". You can navigate diferent windows by clicking in the window \"step\" at the top-left corner. \n",
    "\n",
    "**2.** To load the betweenness just click on \"Betweenness\". You can, at the same time, visualize the communities of the protein nodes, by clicking on \"Show/Hide Colors\" in the \"Color Protein by Communities\" tab (only works for proteins).\n",
    "\n",
    "**3.** The \"Representations\" tab allows the user to show or hide parts of the structure.\n",
    "\n",
    "**4.** If you want to render the network in a higher or lower resolution, the \"Network Drawing Resultion\" tab provides options. Note that after selecting a new resolution, you must load the network representation again.\n",
    "\n",
    "**5.** Three options are available for quickly start rendering. The first two are GPU-only, and therefore depend on your computer having compatible GPU hardware. The third option uses CPU and should run in most computers. If you want to use a different rendering option, in the VMD Main window, go to: File > Render. Then, select the desired rendering option. **Attention** the quick menu for rendering will always save files with the same name. If you want to render multiple figures, and keep all of them, rename the just-rendered images to avoid overwriting the files. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-16T13:58:37.462224Z",
     "start_time": "2020-04-16T13:58:37.453579Z"
    }
   },
   "source": [
    "# ---- The End ----"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}