Tutorial-Step_2-AnalysisAndPlots.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. This notebook is intended to provide a starting point for interactive data analysis. Here, we provide an opportunity for users more comfortable with python programming to tailor the analysis to their specific scientific questions.\n",
    "\n",
    "We make extensive use of **Pandas** data frames to analyze results from network data and community detection. We also use **NGL View** to create interactive visualizations of optimal paths and network nodes. Finally, the R-Python bridge provided by **rpy2** is also used to exemplify the creation of plots with the `ggplot2` package, matching colors between plots and VMD visualizations.\n",
    "\n",
    "The last part of this tutorial, **Step 3**, creates the necessary files to visualize and render publication-quality figures using VMD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:01:47.284555Z",
     "start_time": "2023-02-23T03:01:47.280120Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load the python package\n",
    "import dynetan as dna\n",
    "\n",
    "from dynetan.toolkit import getNodeFromSel, getSelFromNode, getNGLSelFromNode\n",
    "from dynetan.toolkit import getPath, getCartDist\n",
    "from dynetan.viz import getCommunityColors, viewPath, prepTclViz\n",
    "from dynetan.viz import showCommunityByID, showCommunityByNodes, showCommunityByTarget\n",
    "\n",
    "#from numpy.linalg import norm\n",
    "from itertools import islice\n",
    "from collections import defaultdict\n",
    "from scipy import stats\n",
    "\n",
    "import os\n",
    "import MDAnalysis as mda\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import networkx as nx\n",
    "\n",
    "# For visualization\n",
    "import nglview    as nv\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:25.114547Z",
     "start_time": "2023-02-23T00:26:25.112710Z"
    }
   },
   "outputs": [],
   "source": [
    "dnad = dna.datastorage.DNAdata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:25.531355Z",
     "start_time": "2023-02-23T00:26:25.115259Z"
    }
   },
   "outputs": [],
   "source": [
    "# Allows better plotting inside jupyter notebooks\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:25.949834Z",
     "start_time": "2023-02-23T00:26:25.532595Z"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.637465Z",
     "start_time": "2023-02-23T00:26:25.950851Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Make sure all the following packages are installed in your R distribution\n",
    "\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": [
    "## Files and system definitions\n",
    "\n",
    "Ensure all the paths are the same as in the **Step 1** notebook (named `Tutorial-Step_1-ProcessTrajectory.ipynb`). Also, you might want to have analysis results 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": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.640430Z",
     "start_time": "2023-02-23T00:26:26.638304Z"
    }
   },
   "outputs": [],
   "source": [
    "dataDir = \"./TutorialResults/\"\n",
    "\n",
    "# Path where results will be written (you may want plots and data files in a new location)\n",
    "workDir = \"./TutorialResults/AnalysisResults\"\n",
    "\n",
    "fileNameRoot = \"dnaData\"\n",
    "fullPathRoot = os.path.join(dataDir, fileNameRoot)\n",
    "\n",
    "# Define the segID of the Ligand being studied.\n",
    "ligandSegID = \"OMP\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.643479Z",
     "start_time": "2023-02-23T00:26:26.641357Z"
    }
   },
   "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": "2023-02-23T00:26:26.669860Z",
     "start_time": "2023-02-23T00:26:26.645059Z"
    }
   },
   "outputs": [],
   "source": [
    "dnad.loadFromFile(fullPathRoot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we chose to write a reduced trajectory of our system at the end of Step 1, we can recreate the dnad.nodesAtmSel structure to easily select parts of our system based on node indices. This is mainly used for visualizations but can power future analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.698502Z",
     "start_time": "2023-02-23T00:26:26.670740Z"
    }
   },
   "outputs": [],
   "source": [
    "dcdVizFile = fullPathRoot + \"_reducedTraj.dcd\"\n",
    "pdbVizFile = fullPathRoot + \"_reducedTraj.pdb\"\n",
    "\n",
    "workUviz = mda.Universe(pdbVizFile, dcdVizFile)\n",
    "\n",
    "# We add this to the object for ease of access.\n",
    "dnad.nodesAtmSel = workUviz.atoms[ dnad.nodesIxArray ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create variables to easily access the nodes associated with the OMP ligand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.702904Z",
     "start_time": "2023-02-23T00:26:26.699997Z"
    }
   },
   "outputs": [],
   "source": [
    "trgtNodes = getNodeFromSel(\"segid \" + ligandSegID, dnad.nodesAtmSel, dnad.atomToNode)\n",
    "trgtNode = getNodeFromSel(\"segid \" + ligandSegID + \" and name P\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "print(trgtNodes, trgtNode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and Visualize Paths and Communities\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following cells will:\n",
    "\n",
    "- Open an NGL View window displaying the molecular system. If the molecular system is not loaded check for possible errors in loading the python modules.\n",
    "- Show the optimal path between two nodes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NGL View \n",
    "Here, we use the NGL View package to create **interactive** system visualizations. This helps us examine the results created in the previous step and understand the behavior of the system being studied. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.929252Z",
     "start_time": "2023-02-23T00:26:26.703660Z"
    }
   },
   "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.clear_representations()\n",
    "w.add_cartoon(\"protein\")\n",
    "w.add_hyperball(ligandSegID)\n",
    "w.add_hyperball(\"not protein\")\n",
    "\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NGL Selection Language\n",
    "The NGL environment has its own [selection language](https://proteinformatics.uni-leipzig.de/ngl/doc/index.html#User_manual/Usage/Selection_language), and to facilitate our usage of this package, DyNetAn provides the `getNGLSelFromNode` method to create an NGL-compatible selection string from a node selection.\n",
    "\n",
    "We also provide the `getPath` and `viewPath` methods to, respectively, extract a sequence of edges that connect two nodes and plot this selection in the NGL View window above for interactive visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.941101Z",
     "start_time": "2023-02-23T00:26:26.929928Z"
    }
   },
   "outputs": [],
   "source": [
    "# Select an arbitrary node in the enzyme\n",
    "enzNode = getNodeFromSel(\"segid ENZY and resname GLU and resid 115\", dnad.nodesAtmSel, dnad.atomToNode)\n",
    "\n",
    "# Create ball-and-stick representations for node\n",
    "NGL_sel = dna.toolkit.getNGLSelFromNode(enzNode[0], dnad.nodesAtmSel, atom=False)\n",
    "w.add_ball_and_stick(NGL_sel)\n",
    "\n",
    "# Get the optimal path connecting the node selected above to one of the target nondes.\n",
    "optpath = getPath( trgtNodes[1], enzNode[0], dnad.nodesAtmSel, dnad.preds)\n",
    "\n",
    "viewPath(w, optpath, dnad.distsAll, dnad.maxDirectDist, dnad.nodesAtmSel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notice how new residue and path representations were added to the NGL View window we openned above.\n",
    "You can add representations to the same window by using the same view object `w`. If you want to add more paths, for example, simply call the `viewPath` method with other initial and end nodes.\n",
    "\n",
    "After you are done examining the system, simply close the NGL View window with the following `close` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:26.943397Z",
     "start_time": "2023-02-23T00:26:26.941732Z"
    }
   },
   "outputs": [],
   "source": [
    "w.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following cells will:\n",
    "- Open an NGL View window displaying the system.\n",
    "- Show the contacts with the highest betweenness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.136106Z",
     "start_time": "2023-02-23T00:26:26.944171Z"
    }
   },
   "outputs": [],
   "source": [
    "# http://proteinformatics.charite.de/ngl/doc/#User_manual/Usage/Selection_language\n",
    "\n",
    "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.clear_representations()\n",
    "w.add_cartoon(\"protein\")\n",
    "w.add_hyperball(ligandSegID)\n",
    "w.add_hyperball(\"not protein\")\n",
    "\n",
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.160433Z",
     "start_time": "2023-02-23T00:26:27.136940Z"
    }
   },
   "outputs": [],
   "source": [
    "####\n",
    "# Plot the 10 highest betweennesses (edges most used by optimal paths)\n",
    "####\n",
    "\n",
    "w.clear_representations()\n",
    "w.add_cartoon(\"protein\")\n",
    "w.add_hyperball(ligandSegID)\n",
    "w.add_hyperball(\"not protein\")\n",
    "\n",
    "window = 0\n",
    "for k,v in islice(dnad.btws[window].items(),10):\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": "2023-02-23T00:26:27.162911Z",
     "start_time": "2023-02-23T00:26:27.161204Z"
    }
   },
   "outputs": [],
   "source": [
    "####\n",
    "# Save the current visualization in a figue.\n",
    "####\n",
    "\n",
    "#w.download_image(filename='nglview_frame.png', factor=15, trim=True, antialias=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.165416Z",
     "start_time": "2023-02-23T00:26:27.163705Z"
    }
   },
   "outputs": [],
   "source": [
    "####\n",
    "# Close the current visualization.\n",
    "####\n",
    "\n",
    "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 cannot guarantee that the same community will be assigned the same ID in different windows of the same simulation.\n",
    "\n",
    "We finally rank the communities by modularity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.172129Z",
     "start_time": "2023-02-23T00:26:27.166077Z"
    }
   },
   "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": "2023-02-23T00:26:27.184278Z",
     "start_time": "2023-02-23T00:26:27.172843Z"
    }
   },
   "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": "2023-02-23T00:26:27.190901Z",
     "start_time": "2023-02-23T00:26:27.185954Z"
    }
   },
   "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": "2023-02-23T00:26:27.205862Z",
     "start_time": "2023-02-23T00:26:27.192026Z"
    }
   },
   "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": "2023-02-23T00:26:27.210814Z",
     "start_time": "2023-02-23T00:26:27.206880Z"
    }
   },
   "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": "2023-02-23T00:26:27.217644Z",
     "start_time": "2023-02-23T00:26:27.211630Z"
    }
   },
   "outputs": [],
   "source": [
    "# Checks that target nodes are classified in ALL windows\n",
    "nodeCommDF.loc[ nodeCommDF[\"Node\"].isin(trgtNodes) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.223632Z",
     "start_time": "2023-02-23T00:26:27.218441Z"
    }
   },
   "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": "2023-02-23T00:26:27.229401Z",
     "start_time": "2023-02-23T00:26:27.224497Z"
    }
   },
   "outputs": [],
   "source": [
    "nodeCommDFmelt.loc[nodeCommDFmelt[\"Node\"].isin(trgtNodes)].groupby(\"Node\")[\"Cluster\"].apply(np.unique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.232878Z",
     "start_time": "2023-02-23T00:26:27.230139Z"
    }
   },
   "outputs": [],
   "source": [
    "trgtClusters = np.unique( nodeCommDFmelt.loc[nodeCommDFmelt[\"Node\"].isin(trgtNodes), \"Cluster\"].values )\n",
    "trgtClusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.424926Z",
     "start_time": "2023-02-23T00:26:27.233631Z"
    }
   },
   "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": "2023-02-23T00:26:27.431321Z",
     "start_time": "2023-02-23T00:26:27.427003Z"
    }
   },
   "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": "2023-02-23T00:26:27.440982Z",
     "start_time": "2023-02-23T00:26:27.432049Z"
    }
   },
   "outputs": [],
   "source": [
    "# Get all nodes that make contact with target nodes in any window\n",
    "contactNodes = np.unique( np.where( dnad.corrMatAll[:,trgtNodes,:] > 0 )[2] )\n",
    "contactNodesTrgts = list(trgtNodes)\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(trgtNodes, columns=[\"trgtNodes\"]).to_csv(\n",
    "    os.path.join(workDir, \"trgtNodes.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": "2023-02-23T00:26:27.443834Z",
     "start_time": "2023-02-23T00:26:27.441835Z"
    }
   },
   "outputs": [],
   "source": [
    "# Prepares variable names for multi-system comparisons. In this tutorial, we only have one system.\n",
    "system1 = \"OMP\"\n",
    "cDF = nodeCommDFmelt\n",
    "cDF[\"system\"] = system1\n",
    "refWindow1 = refWindow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:27.447158Z",
     "start_time": "2023-02-23T00:26:27.444647Z"
    }
   },
   "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": "2023-02-23T00:26:27.479826Z",
     "start_time": "2023-02-23T00:26:27.447856Z"
    }
   },
   "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": "2023-02-23T00:26:27.483326Z",
     "start_time": "2023-02-23T00:26:27.480743Z"
    }
   },
   "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": "2023-02-23T00:26:27.685525Z",
     "start_time": "2023-02-23T00:26:27.484072Z"
    }
   },
   "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": "2023-02-23T00:26:27.693466Z",
     "start_time": "2023-02-23T00:26:27.686305Z"
    }
   },
   "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(\"protein\")\n",
    "    w.add_hyperball(ligandSegID)\n",
    "    showCommunityByTarget(w, nodeCommDF, trgtNodes, 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": "2023-02-23T00:26:27.696180Z",
     "start_time": "2023-02-23T00:26:27.694285Z"
    }
   },
   "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": "2023-02-23T00:26:27.979976Z",
     "start_time": "2023-02-23T00:26:27.696932Z"
    }
   },
   "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": "2023-02-23T00:26:28.174256Z",
     "start_time": "2023-02-23T00:26:27.980855Z"
    }
   },
   "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": "2023-02-23T00:26:28.860001Z",
     "start_time": "2023-02-23T00:26:28.175198Z"
    }
   },
   "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": "2023-02-23T00:26:44.369019Z",
     "start_time": "2023-02-23T00:26:28.861073Z"
    }
   },
   "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": "2023-02-23T00:26:44.372031Z",
     "start_time": "2023-02-23T00:26:44.370013Z"
    }
   },
   "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": "2023-02-23T00:26:44.374523Z",
     "start_time": "2023-02-23T00:26:44.372787Z"
    }
   },
   "outputs": [],
   "source": [
    "plotFilePrefix = \"OMP_\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:44.886862Z",
     "start_time": "2023-02-23T00:26:44.375252Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,plotFilePrefix -w 800 -h 450\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": "2023-02-23T00:26:45.420322Z",
     "start_time": "2023-02-23T00:26:44.887832Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,trgtNodes,plotFilePrefix -w 800 -h 450\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% trgtNodes))  )\n",
    "\n",
    "breaksNodePlot = dt[Window == refWindow,][ Node %in% trgtNodes, NodePlot ]\n",
    "\n",
    "labelsNode = dt[Window == refWindow,][ Node %in% trgtNodes, 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",
    "\n",
    "ggsave(plotPath, p, device=\"png\")\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T00:26:45.975055Z",
     "start_time": "2023-02-23T00:26:45.421255Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,trgtNodes,contactNodesTrgts,plotFilePrefix -w 800 -h 450\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% trgtNodes))  )\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": "2023-02-23T00:26:46.375010Z",
     "start_time": "2023-02-23T00:26:45.976044Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i workDir,refWindow,trgtNode,trgtNodes,plotFilePrefix,ligandSegID -w 600 -h 450\n",
    "\n",
    "trgtNode = as.integer(trgtNode)\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% trgtNodes, ]\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(trgtNodes), label=dt[ Node %in% as.factor(trgtNodes), unique(resid)]) +\n",
    "    scale_y_discrete(limits=dt[Node == trgtNode]$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": "2023-02-23T00:26:46.475531Z",
     "start_time": "2023-02-23T00:26:46.376085Z"
    }
   },
   "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": "2023-02-23T02:59:44.099239Z",
     "start_time": "2023-02-23T02:59:43.858582Z"
    }
   },
   "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 trgtNodes  for i in range(dnad.numNodes)]),\n",
    "                 columns=[\"node\", \"targets\"])\n",
    "    \n",
    "    df[\"distances\"]  = df.apply( lambda row: dnad.distsAll[window, \n",
    "                                                           int(row[\"node\"]), \n",
    "                                                           int(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(int(row[\"node\"]), \n",
    "                                                         int(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(int(row[\"node\"]), \n",
    "                                                       int(row[\"targets\"]), dnad.numNodes,\n",
    "                                                       dnad.nodeDists, distType), axis=1)\n",
    "    df[\"path\"]       = df.apply( lambda row: list(getPath(int(row[\"node\"]), \n",
    "                                                          int(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": "2023-02-23T02:59:50.127052Z",
     "start_time": "2023-02-23T02:59:50.119502Z"
    }
   },
   "outputs": [],
   "source": [
    "# df.head()\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T02:59:52.093646Z",
     "start_time": "2023-02-23T02:59:52.080088Z"
    },
    "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": "2023-02-23T03:00:00.405892Z",
     "start_time": "2023-02-23T03:00:00.399142Z"
    }
   },
   "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": "2023-02-23T03:00:00.836823Z",
     "start_time": "2023-02-23T03:00:00.823635Z"
    },
    "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": "2023-02-23T03:00:04.046649Z",
     "start_time": "2023-02-23T03:00:02.815319Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "carCorMat[\"iRes\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"i\"]) & \n",
    "                               (df[\"targets\"] == trgtNode[0])][\"resid\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"jRes\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"j\"]) & \n",
    "                               (df[\"targets\"] == trgtNode[0])][\"resid\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"iType\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"i\"]) & \n",
    "                               (df[\"targets\"] == trgtNode[0])][\"type\"].iloc[0] , axis=1)\n",
    "\n",
    "carCorMat[\"jType\"] = carCorMat.apply(lambda row: df[ (df[\"node\"] == row[\"j\"]) & \n",
    "                               (df[\"targets\"] == trgtNode[0])][\"type\"].iloc[0] , axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:00:04.614064Z",
     "start_time": "2023-02-23T03:00:04.607584Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "carCorMat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:00:05.330548Z",
     "start_time": "2023-02-23T03:00:05.319840Z"
    }
   },
   "outputs": [],
   "source": [
    "carCorMat.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:00:05.642870Z",
     "start_time": "2023-02-23T03:00:05.639480Z"
    }
   },
   "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": "2023-02-23T03:00:08.186056Z",
     "start_time": "2023-02-23T03:00:06.983455Z"
    }
   },
   "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": "2023-02-23T03:00:11.131535Z",
     "start_time": "2023-02-23T03:00:11.062934Z"
    }
   },
   "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": "2023-02-23T03:00:27.523117Z",
     "start_time": "2023-02-23T03:00:27.232989Z"
    },
    "scrolled": false
   },
   "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 stite 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": "2023-02-23T03:00:32.741393Z",
     "start_time": "2023-02-23T03:00:32.263244Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i carCorMatInterface,workDir,ligandSegID -w 800 -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 != \"TIP3\"][Corr > corrCutoff][,iRes]\n",
    "dt <- dt[ iRes %in% iResList, ]\n",
    "dt <- dt[jType != \"TIP3\", .SD[which.max(Corr)] , by=.(iRes,jRes)]\n",
    "\n",
    "dt <- dt[, iRes := sapply(strsplit(iRes, \"_\"), '[', 1) ]\n",
    "\n",
    "colours <- c(\"Aminoacid\" = \"red\", \n",
    "          \"TIP3\" = \"blue\", \n",
    "          \"Nucleotide\" = \"darkorange\")\n",
    "colours <- c(colours, setNames(\"purple\",ligandSegID))\n",
    "\n",
    "shapes <- c(\"Aminoacid\" = 19,\n",
    "          \"TIP3\" = 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": [
    "### Vizualize "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:00:40.990157Z",
     "start_time": "2023-02-23T03:00:40.759194Z"
    }
   },
   "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": "2023-02-23T03:01:15.872111Z",
     "start_time": "2023-02-23T03:01:15.672866Z"
    }
   },
   "outputs": [],
   "source": [
    "# Paint all connections to the ligand\n",
    "window = 1\n",
    "\n",
    "w.clear_representations()\n",
    "w.add_cartoon(\"protein\")\n",
    "w.add_hyperball(ligandSegID)\n",
    "\n",
    "tmpDF = carCorMat[(carCorMat[\"i\"].isin(trgtNodes)) | (carCorMat[\"j\"].isin(trgtNodes))]\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": "2023-02-23T03:01:24.194302Z",
     "start_time": "2023-02-23T03:01:24.191636Z"
    }
   },
   "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": "2023-02-23T03:01:27.450845Z",
     "start_time": "2023-02-23T03:01:26.615555Z"
    }
   },
   "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": "2023-02-23T03:01:30.691051Z",
     "start_time": "2023-02-23T03:01:30.687080Z"
    }
   },
   "outputs": [],
   "source": [
    "srcNode = getNodeFromSel(\"resname VAL and resid 11 and segid ENZY\",dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "print(\"Source node:\", srcNode)\n",
    "\n",
    "trgNode = getNodeFromSel(\"resname ALA and resid 69 and segid ENZY\",dnad.nodesAtmSel, dnad.atomToNode)[0]\n",
    "print(\"Target node:\", trgNode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:01:31.515466Z",
     "start_time": "2023-02-23T03:01:31.511825Z"
    }
   },
   "outputs": [],
   "source": [
    "getSelFromNode(58,dnad.nodesAtmSel)"
   ]
  },
  {
   "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": "2023-02-23T03:01:32.633861Z",
     "start_time": "2023-02-23T03:01:32.631525Z"
    }
   },
   "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",
    "\n",
    "nodesForPaths = [ [0,58], [0,60] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:01:53.613552Z",
     "start_time": "2023-02-23T03:01:53.520946Z"
    }
   },
   "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": {
    "ExecuteTime": {
     "end_time": "2023-02-23T03:01:55.740613Z",
     "start_time": "2023-02-23T03:01:55.737058Z"
    }
   },
   "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": "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 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}