Step1_Process_Trajectory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Network Analysis Tutorial - Step 1\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/. \n",
    "\n",
    "In this tutorial, we will use the Dynamic Network Analysis python package to explore the interactions between the OMP decarboxylase enzyme and its substrate, identifying clusters of amino acid residues that form domains, and active site residues that are important for binding and enzymatic activity.\n",
    "\n",
    "The tutorial is divided in two jupyter notebooks. In this notebook, **Step 1**, we will analyze the MD trajectory for the OMP decarboxylase system and generate the network data used for analysis. \n",
    "\n",
    "The trajectory files have approximately 500MB in size, and must be downloaded [from this link](http://www.ks.uiuc.edu/~rcbernardi/NetworkAnalysis/DynamicNetworkAnalysis_MDdata.tar.gz) and placed in the *TutorialData* folder.\n",
    "\n",
    "In the accompanying notebook, **Step 2**, we will load the data generated in the Step 1, and create analysis plots and visualizations of the network. \n",
    "\n",
    "In **Step 3**, which is outside of the jupyter notebook environment, we will produce high-quality renderings of the network using the popular visualization software VMD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:26.970833Z",
     "start_time": "2020-04-14T20:02:20.399680Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load the python package\n",
    "import os\n",
    "import dynetan\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",
    "import argparse\n",
    "from operator import itemgetter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the object that processes MD trajectories.\n",
    "dnap = DNAproc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bPythonExport = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapResidueNames={'ALA':'A','CYS':'C','ASP':'D','GLU':'E','PHE':'F',\n",
    "                 'GLY':'G','HIS':'H','HSD':'H','HSE':'H','ILE':'I','LYS':'K','LEU':'L',\n",
    "                 'MET':'M','ASN':'N','PRO':'P','GLN':'Q','ARG':'R',\n",
    "                 'SER':'S','THR':'T','VAL':'V','TRP':'W','TYR':'Y',\n",
    "                 'MG':'Mg','POPC':'Popc',\n",
    "                 'ATP':'Atp','GTP':'Gtp',\n",
    "                 'NA':'Sod','SOD':'Sod','CLA':'Cl','CL':'Cl','POT':'Pot','K':'Pot',\n",
    "                 'SOL':'h2o','HOH':'h2o','WAT':'h2o','TIP':'h2o','H2O':'h2o',\n",
    "                }\n",
    "def name_node(dnap, node):\n",
    "    #i=dnap.nodesAtmSel[node].index\n",
    "    resname=dnap.nodesAtmSel[node].resname ; resid=dnap.nodesAtmSel[node].resid\n",
    "    return \"%s%s\" % (mapResidueNames[resname], resid)\n",
    "\n",
    "def clarify_duplicate_nodes(dictNames, dictSuffix):\n",
    "    \"\"\"\n",
    "    From two dicts with the same keys, add the respective suffix to all keys in the former that possess duplicate values.\n",
    "    \"\"\"\n",
    "    from itertools import chain\n",
    "    dictRev = {}\n",
    "    for k, v in dictNames.items():\n",
    "        dictRev.setdefault(v, set()).add(k)\n",
    "    setDuplicateKeys = set(chain.from_iterable( v for k, v in dictRev.items() if len(v) > 1))\n",
    "    \n",
    "    for k in setDuplicateKeys:\n",
    "        dictNames[k] = dictNames[k]+\"_\"+dictSuffix[k]\n",
    "    return dictNames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Node Groups\n",
    "\n",
    "In Network Analysis, each residue is represented by one or more \"nodes\", serving as proxies for groups of atoms form the original residue (Figure 1). This approach lowers computational cost and noise.\n",
    "\n",
    "For our purposes, we treat this as a coarse-graining procedure such that large molecules like ATP and POPC will gain multiple nodes to better distinguish between connections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = Default workflow control parameters.\n",
    "bHasInterfaceDefinitions = False\n",
    "bIncludeSolvent = False\n",
    "segIDs = [] ; h2oName = []\n",
    "numSampledFrames = 100\n",
    "strideDCD = 0\n",
    "cutoffDist = 4.5 ; contactPersistence = 0.75\n",
    "numCoresAvailable=1\n",
    "bInMem = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if bPythonExport:\n",
    "    parser = argparse.ArgumentParser(description='Process trajactories for Dynamic Network Analysis.'\n",
    "                                                 'Equivalent to Step 1 of the tutorial without using Jupyter.',\n",
    "                                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)\n",
    "    parser.add_argument('-o', '--output', type=str, dest='outputPrefix', default=None,\n",
    "                        help='Prefix used for all output files.')    \n",
    "    parser.add_argument('--top', type=str, dest='fileTopology', default='seg.psf',\n",
    "                        help='The name of a suitable file for MDAnalysis to generate a topology.')\n",
    "    parser.add_argument('fileTrajectories', type=str, metavar='N', nargs='+', default='sum.xtc',\n",
    "                        help='One or more trajectory files for MDAnalysis to read in. All frames are combined together.')\n",
    "    parser.add_argument('-n', '--num_windows', type=int, dest='numWinds', default=None,\n",
    "                        help='[Opt.] Split the total trajectory frames into N windows. By default equal to the number of files read.')\n",
    "    parser.add_argument('-f', '--num_frames', type=int, dest='numFrames', default=numSampledFrames,\n",
    "                        help='Number of frames to read for each window. Total number of frames will then be nWind*nFrames. '\n",
    "                       'This is computed by DNA as the total frames divided evenly.')\n",
    "    parser.add_argument('--dcd_stride', type=int, dest='strideDCD', default=strideDCD,\n",
    "                        help='The frame stride for the output dcd file, give to save disk space. Give 0 to save only the first and last frames.')\n",
    "#    = = = Very inconvenient to do within MDAnalysis. Outside design doc.\n",
    "#    parser.add_argument('--no_span', dest='bNoSpan', action='store_true',\n",
    "#                        help='When the number of frames is given, prevent windows from spanning across multiple trajectory files. '\n",
    "#                       'This will cause remainder frames to be discarded.')\n",
    "    parser.add_argument('--cpus', type=int, dest='nCPUs', default=numCoresAvailable,\n",
    "                            help='Set number of processes to enable partial parallelisation.')\n",
    "    #   =============== Analysis parameters\n",
    "    parser.add_argument('--segIDs', type=str, dest='segIDs', default=None,\n",
    "                        help='A list of comma-separated segment IDs of components that should be included in network analysis. '\n",
    "                             'By default, will try to predict and include segIDs that are not water or named after simply ions. '\n",
    "                             'Note: The water segment is defined in the next argument.')\n",
    "    parser.add_argument('--segIDs_solv', type=str, dest='segIDs_solv', default=None,\n",
    "                        help='Name of the water segment. A prediction will be attempted if not given '\n",
    "                             'by searching for possible names like SOL, HOH, TIP, WAT, H2O.')\n",
    "    parser.add_argument('--interfaceA', type=str, dest='interfaceA', default=None,\n",
    "                        help='Definition of side-A of an interface to be analysed, using MDAnalysis selection syntax.')\n",
    "    parser.add_argument('--interfaceB', type=str, dest='interfaceB', default=None,\n",
    "                        help='Definition of side-B of an interface to be analysed, using MDAnalysis selection syntax.')\n",
    "    parser.add_argument('--exclude_solvent', action='store_true', dest='bExcludeSolvent',\n",
    "                        help='Boolean to exclude potential solvent contacts in the analysis.')\n",
    "    parser.add_argument('--cutoff', type=float, dest='distCutoff', default=cutoffDist,\n",
    "                        help='Distance in Angstroms to consider two nodes as being in contact.'\n",
    "                             'Calculated as the nearest distance between all atoms of the node.')\n",
    "    parser.add_argument('--persistence', type=float, dest='ratioContactPersistence', default=contactPersistence,\n",
    "                        help='Minimum portion of frames for DNA to consider two nodes to be in sufficient contact,'\n",
    "                             'such that it will include this node pair as an edge in the output graph.')\n",
    "    parser.add_argument('--in_memory', action='store_true', dest='bInMem',\n",
    "                        help='Boolean to load the enitre MD trajectroy into RAM for increased performance.')\n",
    "\n",
    "    args = parser.parse_args()\n",
    "\n",
    "    numCoresAvailable=args.nCPUs\n",
    "\n",
    "    # = = = Analysis Configurations = = =\n",
    "    # Cutoff for contact map (In Angstroms)\n",
    "    cutoffDist = args.distCutoff\n",
    "    # Minimum contact persistance (In ratio of total trajectory frames)\n",
    "    contactPersistence = args.ratioContactPersistence\n",
    "    # Inclusion of solvent\n",
    "    bIncludeSolvent = not args.bExcludeSolvent\n",
    "\n",
    "    # = = = Optional interface definition. A dummy is other calculated.\n",
    "    if args.interfaceA is not None and args.interfaceB is not None:\n",
    "        bHasInterfaceDefinitions=True\n",
    "        seltextInterfaceA=args.interfaceA\n",
    "        seltextInterfaceB=args.interfaceB\n",
    "    \n",
    "    # = = Define selections of segments that will be included. = = \n",
    "    # = = Following the NAMD/VMD system, DynaNet uses a segment syntax instead of chain syntax to classify atoms.\n",
    "    #     Thus, segments will need to be defined in the topology file.\n",
    "    if args.segIDs is not None:\n",
    "        segIDs = args.segIDs.split(\",\")\n",
    "    if args.segIDs_solv is not None:\n",
    "        h2oName = args.segIDs_solv.split(\",\")\n",
    "    \n",
    "    topFile = args.fileTopology\n",
    "    trjFiles = args.fileTrajectories\n",
    "    bInMem = args.bInMem\n",
    "    \n",
    "    if args.outputPrefix is None:\n",
    "        pathToData = \"./results\" % (allele, temperature)\n",
    "        fileNameRoot = \"analysis\"\n",
    "    else:\n",
    "        pathToData, fileNameRoot = os.path.split(args.outputPrefix)\n",
    "        \n",
    "    fullPathRoot = os.path.join(pathToData, fileNameRoot)\n",
    "\n",
    "    # Number of windows created from full simulation.\n",
    "    if args.numWinds is not None:\n",
    "        numWinds = args.numWinds\n",
    "    else:\n",
    "        numWinds = len(args.fileTrajectories)\n",
    "        \n",
    "    # Sampled frames per window\n",
    "    numSampledFrames = args.numFrames\n",
    "    \n",
    "    # Output\n",
    "    strideDCD = args.strideDCD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = Change to the relevant work folder\n",
    "if not bPythonExport:\n",
    "    systemExample=\"caspase-1\"\n",
    "    if systemExample == \"UbqCHARMM\":\n",
    "        %cd /home/zharmad/host/projects/Ubq-md\n",
    "    elif systemExample == \"UbqAthi\":\n",
    "        %cd /home/zharmad/host/shared-colleague/Ubq-2017\n",
    "    elif systemExample == \"periplasmic\":\n",
    "        %cd /home/zharmad/projects/periplasmic/leucine-binding_protein\n",
    "    elif systemExample == \"caspase-1\":\n",
    "        %cd /home/zharmad/host/projects/caspase-1        \n",
    "    elif systemExample == \"CFTR\":\n",
    "        %cd /home/zharmad/projects/cftr/DyNetAn\n",
    "    else:\n",
    "        %cd .."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = Customised workflow control parameters for publication examples.\n",
    "if not bPythonExport:\n",
    "    bInMem=True\n",
    "    numCoresAvailable = 6\n",
    "    if systemExample == \"UbqCHARMM\":\n",
    "        # apo1 apo2 apo3\n",
    "        state=\"apo1\" ; workDir = \"./apo-md01\"\n",
    "        topFile = \"./tops/prot-segids.pdb\"\n",
    "        trjFiles= os.path.join(workDir, \"prot-masses.xtc\" )\n",
    "        numWinds = 5 ; numSampledFrames = 1000\n",
    "        pathToData = \"./dynetan\"\n",
    "        fileNameRoot = \"%s_%ix%i\" % (state, numWinds, numSampledFrames)\n",
    "        fullPathRoot = os.path.join(pathToData, fileNameRoot)\n",
    "        #segIDs = [\"UBQ\"] ; h2oName = [\"SOL\"]\n",
    "    elif systemExample == \"UbqAthi\":\n",
    "        # UbqI13V  UbqI23A  UbqI30A  UbqL43A  UbqL67A  UbqL69A  UbqV17A  UbqWT\n",
    "        state=\"UbqV17A\" ; workDir = \"./%s\" % state\n",
    "        topFile = os.path.join(workDir, \"reference.pdb\" )\n",
    "        trjFiles= os.path.join(workDir, \"traj-1ns.xtc\" )\n",
    "        numWinds = 5 ; numSampledFrames = 200\n",
    "        pathToData = \"./dynetan\"\n",
    "        fileNameRoot = \"%s_%ix%i\" % (state, numWinds, numSampledFrames)\n",
    "        fullPathRoot = os.path.join(pathToData, fileNameRoot)\n",
    "        #segIDs = [\"UBQ\"] ; h2oName = [\"SOL\"]        \n",
    "    elif systemExample == \"periplasmic\":\n",
    "        # apo holo-leu\n",
    "        state=\"apo\" ; workDir = \"./apo\"\n",
    "        topFile = \"./apo/tops/prot-segids.pdb\"\n",
    "        trjFiles= os.path.join(workDir, \"protcent.xtc\" )\n",
    "        numWinds = 5 ; numSampledFrames = 200\n",
    "        pathToData = \"./dynetan\"\n",
    "        fileNameRoot = \"%s_%ix%i\" % (state, numWinds, numSampledFrames)\n",
    "        fullPathRoot = os.path.join(pathToData, fileNameRoot)\n",
    "        #segIDs = [\"LBP\"] ; h2oName = [\"SOL\"]\n",
    "        if state != \"apo\":\n",
    "            bHasInterfaceDefinitions = True\n",
    "            seltextInterfaceA=\"resid 15 18 77 78 79 80 100 101 102 103 118 121 124 150 202 226 227 276\"\n",
    "            seltextInterfaceB=\"resid 347\"\n",
    "    elif systemExample == \"caspase-1\":\n",
    "        # off-state on-state\n",
    "        state=\"off-state\" ; workDir = \"./off-state\"\n",
    "        topFile = \"%s/apo-em.pdb\" % workDir\n",
    "        trjFiles= os.path.join(workDir, \"prod/centered.xtc\" )\n",
    "        numWinds = 2 ; numSampledFrames = 5\n",
    "        pathToData = \"./dynetan\"\n",
    "        fileNameRoot = \"%s_%ix%i\" % (state, numWinds, numSampledFrames)\n",
    "        fullPathRoot = os.path.join(pathToData, fileNameRoot)\n",
    "        #segIDs = [\"C1A\",\"C1B\",\"C2A\",\"C2B\"] ; h2oName = [\"SOL\"]\n",
    "        bIncludeSolvent=True\n",
    "    elif systemExample == \"CFTR\":\n",
    "        # Define mutant file IO locations. wt, P67L, E56K, R75Q, dF508 , S945L\n",
    "        allele=\"wt\" ; temperature=\"310K\" ; nRepl=6\n",
    "        workDir = \"./trajectories/%s/%s/\" % (allele, temperature)\n",
    "        topFile = os.path.join(workDir, \"1/clustered_d3.5_r0.50.pdb\")\n",
    "        trjFiles=[]\n",
    "        for i in range(1,nRepl+1):\n",
    "            trjFiles.append(os.path.join(workDir, \"%i/clustered_d3.5_r0.50.xtc\" %i))\n",
    "        pathToData = \"./results/%s/%s/\" % (allele, temperature)\n",
    "        fileNameRoot = \"1to%i\" % nRepl\n",
    "        fullPathRoot = os.path.join(pathToData, fileNameRoot)        \n",
    "        numWinds = 5 ; numSampledFrames = 200\n",
    "        #segIDs = [\"LAS\",\"TD1\",\"ND1\",\"RDO\",\"TD2\",\"ND2\",\"CTR\",\"LIP\",\"ATP\",\"POT\",\"CLA\"]\n",
    "        segIDs = [\"LAS\",\"TD1\",\"ND1\",\"RDO\",\"TD2\",\"ND2\",\"CTR\",\"CRY\"]\n",
    "        h2oName = [\"SOL\"]\n",
    "        bIncludeSolvent=True\n",
    "    else:\n",
    "        raise KeyboardInterrupt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if numCoresAvailable>1:\n",
    "    strBackend=\"openmp\"\n",
    "else:\n",
    "    strBackend=\"serial\"\n",
    "\n",
    "if not os.path.exists(pathToData):\n",
    "    os.makedirs(pathToData)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:26.989963Z",
     "start_time": "2020-04-14T20:02:26.973220Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Number of sampled frames for automatic selection of solvent and ions.\n",
    "# numAutoFrames = numSampledFrames*numWinds\n",
    "\n",
    "# Network Analysis will make one node per protein residue (in the alpha carbon)\n",
    "# For all other residues, the user must specify atom(s) that will represent a node.\n",
    "# ...\n",
    "# We also need to know the heavy atoms that compose each node group.\n",
    "\n",
    "customResNodes = {}\n",
    "usrNodeGroups = {}\n",
    "\n",
    "# Cater to HSD and HSE in the CHARMM forcefield.\n",
    "# Lip: Uses POPC as a reference. Node designations include:\n",
    "#      ...zwitterion centers, ester centers, and two additional carbon centers for each tail capturing ~7 atoms.\n",
    "#      ...heavy atom assignment based on charge groups in CHARMMM36\n",
    "#      ...a PSFGEN error means that the naming of carbons C216 -> 6C12.\n",
    "#customResNodes[\"HSD\"] = [\"CA\"]\n",
    "#usrNodeGroups[\"HSD\"] = {}\n",
    "#usrNodeGroups[\"HSD\"][\"CA\"] = set(\"N CA CB ND1 CG CE1 NE2 CD2 C O\".split())\n",
    "#customResNodes[\"HSE\"] = [\"CA\"]\n",
    "#usrNodeGroups[\"HSE\"] = {}\n",
    "#usrNodeGroups[\"HSE\"][\"CA\"] = set(\"N CA CB ND1 CG CE1 NE2 CD2 C O\".split())\n",
    "\n",
    "customResNodes[\"K\"] = [\"K\"]\n",
    "usrNodeGroups[\"K\"] = {}\n",
    "usrNodeGroups[\"K\"][\"K\"] = set(\"K\".split())\n",
    "\n",
    "customResNodes[\"CL\"] = [\"CL\"]\n",
    "usrNodeGroups[\"CL\"] = {}\n",
    "usrNodeGroups[\"CL\"][\"CL\"] = set(\"CL\".split())\n",
    "\n",
    "customResNodes[\"SOL\"] = [\"OW\"]\n",
    "usrNodeGroups[\"SOL\"] = {}\n",
    "usrNodeGroups[\"SOL\"][\"OW\"] = set(\"OW HW1 HW2\".split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#################################\n",
    "### Load info to object\n",
    "\n",
    "dnap.setNumWinds(numWinds)\n",
    "dnap.setNumSampledFrames(numSampledFrames)\n",
    "dnap.setCutoffDist(cutoffDist)\n",
    "dnap.setContactPersistence(contactPersistence)\n",
    "dnap.seth2oName(h2oName)\n",
    "dnap.setSegIDs(segIDs)\n",
    "\n",
    "dnap.setCustomResNodes(customResNodes)\n",
    "dnap.setUsrNodeGroups(usrNodeGroups)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the trajectory\n",
    "\n",
    "Our Generalized Network Analysis leverages the MDAnalysis package to create a *universe* that contains all the trajectory and system information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:28.159684Z",
     "start_time": "2020-04-14T20:02:26.993405Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.loadSystem(topFile,trjFiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = First check if requests are sane.\n",
    "if numWinds*numSampledFrames > dnap.workU.trajectory.n_frames:\n",
    "    print(\"= = ERROR: You have requested more simulation frames to be analysed than are available!\")\n",
    "    if bPythonExport:\n",
    "        import sys\n",
    "        sys.exit()\n",
    "    else:\n",
    "        raise KeyboardInterrupt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# = = = Attempt to predict segment IDs to be analysed, if not given.\n",
    "ionSegIDs   = ['ION','SOD','CLA','POT','HET','NA','K','CL']\n",
    "waterSegIDs = ['HOH','SOL','WAT','TIP','H2O']\n",
    "if segIDs == []:\n",
    "    # = = = grab every segment that is not named after likely simple ion or water segments.\n",
    "    soluteIDs = []\n",
    "    for s in dnap.workU.segments.segids:\n",
    "        if (not s in ionSegIDs) and (not s in waterSegIDs):\n",
    "            soluteIDs.append( s )\n",
    "    dnap.setSegIDs( soluteIDs )\n",
    "    print(\"= = Automated prediction has assigned %s solute segments to be included for analysis.\" % str(soluteIDs))\n",
    "if h2oName == []:\n",
    "    # = = = Search the segments for a name that sounds like water\n",
    "    for s in waterSegIDs:\n",
    "        if s in dnap.workU.segments.segids:\n",
    "            h2oName.append(s)\n",
    "            break\n",
    "    if h2oName != []:\n",
    "        dnap.seth2oName(h2oName)\n",
    "    else:\n",
    "        # = = = Is possibly a protein-only trajectory. Give a dummy argument.\n",
    "        dnap.seth2oName(['HOH'])\n",
    "    print(\"= = Automated prediction has assigned %s to be the water segment ID.\" % str(h2oName))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checks segments and residue names\n",
    "\n",
    "This is important to know if there are residues in the structure that we didn't know of, and need to be addresssed so that network analysis can create nodes in all selected residues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:28.485023Z",
     "start_time": "2020-04-14T20:02:28.187188Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.checkSystem()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatically identify crystallographic waters and ions.\n",
    "\n",
    "In this section we identify all residues which will be checked for connectivity with selected segments.\n",
    "\n",
    "- First, we define if all solvent molecules will be checked, or if just ions, ligands, lipids, and other molecules will be checked.\n",
    "\n",
    "- Second, we sample a small set of frames from the trajectory to select likely residues, then we check all trajectory to see if they are closer than the cutoff distance for at least x% of the simulation (where x is the \"contact persistence\" fraction of the trajectory).\n",
    "\n",
    "- Third, we load the trajectory of the relevant atoms to memory to speed up the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:35.528972Z",
     "start_time": "2020-04-14T20:02:28.494862Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.selectSystem(withSolvent=bIncludeSolvent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare network representation of the system\n",
    "\n",
    "Here we check that we know how to treat all types of residues in the final selection. Every residue will generate one or more nodes in the final network. Then we store the groups of atoms that define each node."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dnap.getU().residues[300]\n",
    "# dnap.getU().residues[300].resname\n",
    "# 'C13' in set.union(*dnap.resNodeGroups['POPC'].values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:36.062392Z",
     "start_time": "2020-04-14T20:02:35.535533Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.prepareNetwork()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Align the trajectory based on selected segments\n",
    "\n",
    "We align the trajectory to its first frame using heavy atoms (non-hydrogen) from the selected segments. In the process, we also transfer the trajectory to the computer memory, so that future analysis and manipulations are completed faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:02:37.698151Z",
     "start_time": "2020-04-14T20:02:36.082241Z"
    }
   },
   "outputs": [],
   "source": [
    "# If your system is too large, you can turn off the \"in memory\" option, at a cost for performance.\n",
    "dnap.alignTraj(inMemory=bInMem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select residues that are closer than 4.5A for more than 75% of simulation\n",
    "\n",
    "Creates an N-by-N matrix for all N nodes in the selected region, and automatically selected nodes (ions, solvent).\n",
    "\n",
    "The following cell defines efficient functions to run the analysis and create a contact matrix. We leverage both MDAnalysis parallel contact detection tools, as well as accelerated Numba and Cython function. After creating the contact matrix, we remove any automatically selected nodes that have insuficient persistance, and filter the contacts by (optionally) removing contacts between nodes in consecutive residues in a protein or nucleic chain.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Attention** For every time you start this Jupyter Notebook, the first time you execute this function may take significanlty longer (several seconds) to start. This is because we use *Cython* and *Numba* to compile functions \"on-demand\", and a new compilation may be necessary after the notebook is re-started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:25.577099Z",
     "start_time": "2020-04-14T20:02:37.700995Z"
    }
   },
   "outputs": [],
   "source": [
    "# To speed-up the contact matrix calculation, a larger stride can be selected, at a cost for precision.\n",
    "dnap.findContacts(stride=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing contacts between nodes in the same residue.\n",
    "\n",
    "The following function guarantees that there will be no \"self-contacts\" (contacts between a node and itself), and gives you the opportunity to remove contacts between nodes in consecutive residues (such as sequential amino acids in the same chain, removing back-bone interactions). \n",
    "\n",
    "The function also removes nodes that are never in contact with any other node in the system (such as the ends of flexible chains, or residues in flexible loops). This will automatically update the MDanalysis universe and related network informatio, such as number of nodes and atom-to-node mappings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:30.345737Z",
     "start_time": "2020-04-14T20:03:25.588031Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dnap.filterContacts(notSameRes=True, notConsecutiveRes=False, removeIsolatedNodes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Generalized Correlation with Python/Numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:43.328370Z",
     "start_time": "2020-04-14T20:03:30.348061Z"
    }
   },
   "outputs": [],
   "source": [
    "# We can calculate generalized correlaions in parallel using Python's multiprocessing package.\n",
    "dnap.calcCor(ncores=numCoresAvailable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate cartesian distances between all nodes in the selected system.\n",
    "\n",
    "Here, we will calculate the **shortest** distance between atoms in all pairs of nodes. It is similar to the contact matrix calculation, but we check all distances and keep the shortest one to use in our analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( dnap.atomToNode.shape )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:54.619509Z",
     "start_time": "2020-04-14T20:03:43.331634Z"
    }
   },
   "outputs": [],
   "source": [
    "# We can leverage MDanalysis parallelization options with backend=\"serial\" or backend=\"openmp\".\n",
    "# For very small systems, the serial can be faster!\n",
    "dnap.calcCartesian(backend=strBackend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Network Calculations\n",
    "Create a graph from our correlation matrix. Different properties are calculated:\n",
    "\n",
    "*Density* measures how connected the graph is compared to how connected it *could* be. It is the ratio between edges in the graph over all possible edges between all pairs of nodes.\n",
    "\n",
    "*Transitivity* maesures the triadic closure, comparing present triangles to possible triangles. In a triangle, if A is connected to B, and B connected to C, then A is connected to C.\n",
    "\n",
    "*Degree* measures the number of connections a node has.\n",
    "\n",
    "(Reference)[1]\n",
    "\n",
    "[1]:https://programminghistorian.org/en/lessons/exploring-and-analyzing-network-data-with-python#advanced-networkx-community-detection-with-modularity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:56.272811Z",
     "start_time": "2020-04-14T20:03:54.622084Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.calcGraphInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:56.283952Z",
     "start_time": "2020-04-14T20:03:56.276920Z"
    }
   },
   "outputs": [],
   "source": [
    "# Basic information of the network as interpreted as a graph.\n",
    "print( nx.info(dnap.nxGraphs[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:56.507667Z",
     "start_time": "2020-04-14T20:03:56.287113Z"
    }
   },
   "outputs": [],
   "source": [
    "# Both density and transitivity are scaled from 0 to 1\n",
    "for win in range(dnap.numWinds):\n",
    "    print(\"----- Window {} -----\".format(win))\n",
    "    print(\"Density:\", round( nx.density(dnap.nxGraphs[win]), 4) )\n",
    "    print(\"Transitivity:\", round( nx.transitivity(dnap.nxGraphs[win]), 4) )\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:03:56.530227Z",
     "start_time": "2020-04-14T20:03:56.511898Z"
    }
   },
   "outputs": [],
   "source": [
    "from operator import itemgetter\n",
    "\n",
    "# We can check the nodes that have the most connections in each window.\n",
    "for win in range(dnap.numWinds):\n",
    "    print(\"----- Window {} -----\".format(win))\n",
    "    \n",
    "    sorted_degree = sorted(dnap.getDegreeDict(win).items(), key=itemgetter(1), reverse=True)\n",
    "    \n",
    "    print(\"Top 10 nodes by degree: [node --> degree : selection]\")\n",
    "    for n,d in sorted_degree[:5]:\n",
    "        print(\"{0:>4} --> {1:>2} : {2}\".format(n, d, getSelFromNode(n, dnap.nodesAtmSel)))\n",
    "    \n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate optimal paths\n",
    "We choose the Floyd Warshall algorithm[1]. This uses the **correlations as weights** to calculate network distances and shortest distances.\n",
    "\n",
    "[1]:https://networkx.github.io/documentation/stable/reference/algorithms/generated/networkx.algorithms.shortest_paths.dense.floyd_warshall.html?highlight=warshall#networkx.algorithms.shortest_paths.dense.floyd_warshall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:03.692662Z",
     "start_time": "2020-04-14T20:03:56.533388Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.calcOptPaths(ncores=numCoresAvailable)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate betweenness\n",
    "\n",
    "We calculate both betweenness centrality[1] for edges and eigenvector centrality[2] for nodes.\n",
    "\n",
    "[1]:https://networkx.github.io/documentation/stable/reference/algorithms/generated/networkx.algorithms.centrality.edge_betweenness_centrality.html?highlight=betweenness#networkx.algorithms.centrality.edge_betweenness_centrality\n",
    "[2]:https://networkx.github.io/documentation/networkx-1.10/reference/generated/networkx.algorithms.centrality.eigenvector_centrality.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:04.371111Z",
     "start_time": "2020-04-14T20:04:03.698752Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dnap.calcBetween(ncores=numCoresAvailable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:04.383621Z",
     "start_time": "2020-04-14T20:04:04.374440Z"
    }
   },
   "outputs": [],
   "source": [
    "from itertools import islice\n",
    "\n",
    "# Pairs of nodes with highest Betweeness values, compared to their correlation values (in Window 0)\n",
    "for k,v in islice(dnap.btws[0].items(),5):\n",
    "    print(\"Nodes {} have betweenes {} and correlation {}.\".format(k, \n",
    "                                                                  round(v,3), \n",
    "                                                                  round(dnap.corrMatAll[0, k[0], k[1]], 3) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Turn to node centrality instead of edge centrality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:04.706385Z",
     "start_time": "2020-04-14T20:04:04.386209Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.calcEigenCentral()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate communities\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using **Louvain heuristices** is feasible. \n",
    "This method also maximizes the modularity of the network.\n",
    "\n",
    "http://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:04.985294Z",
     "start_time": "2020-04-14T20:04:04.712272Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.calcCommunities()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:04.997552Z",
     "start_time": "2020-04-14T20:04:04.988395Z"
    }
   },
   "outputs": [],
   "source": [
    "# Sort communities based on number of nodes\n",
    "for comIndx in dnap.nodesComm[0][\"commOrderSize\"]:\n",
    "    print(\"Modularity Class {0:>2}: {1:>3} nodes.\".format(comIndx, len(dnap.nodesComm[0][\"commNodes\"][comIndx])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:05.099322Z",
     "start_time": "2020-04-14T20:04:05.000447Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Sort communities based on the node with highest eigenvector centrality\n",
    "for comIndx in dnap.nodesComm[0][\"commOrderEigenCentr\"]:\n",
    "    print(\"Modularity Class {0} ({1} nodes) Sorted by Eigenvector Centrality:\".format(\n",
    "                                                                    comIndx, \n",
    "                                                                len(dnap.nodesComm[0][\"commNodes\"][comIndx])))\n",
    "    for node in dnap.nodesComm[0][\"commNodes\"][comIndx][:5]:\n",
    "        print(\"Name: {0:>4} | Degree: {1:>2} | Eigenvector Centrality: {2}\".format(\n",
    "            node, dnap.nxGraphs[win].nodes[node]['degree'], dnap.nxGraphs[win].nodes[node]['eigenvector']))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process Interface Residues\n",
    "\n",
    "We now find all nodes that are close to both selections chosen by the user. That may include amino acids in the interface, as well as ligands, waters and ions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:06.765269Z",
     "start_time": "2020-04-14T20:04:05.102432Z"
    }
   },
   "outputs": [],
   "source": [
    "if bHasInterfaceDefinitions:\n",
    "    dnap.interfaceAnalysis(selAstr=seltextInterfaceA, selBstr=seltextInterfaceB)\n",
    "else:\n",
    "    # Conduct analysis on dummy inputs since this is hard-coded.\n",
    "    dnap.interfaceAnalysis(selAstr=getSelFromNode(0,dnap.nodesAtmSel), selBstr=getSelFromNode(1,dnap.nodesAtmSel) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute additional graph properties here to save some time for the analysis in Step 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: btws is the same as the results derived from nx.edge_betweenness_centrality()\n",
    "for w in range(numWinds):\n",
    "    G = dnap.nxGraphs[w]\n",
    "\n",
    "    # Transfer edge-betweenness to graph.\n",
    "    for u,v in dnap.btws[w]:\n",
    "        G.edges[u,v]['btws']=dnap.btws[w][u,v]\n",
    "    # Note: edge-betweeness weighted clustering coefficient. \n",
    "    c = nx.clustering(G, weight='btws')\n",
    "    for x in range(dnap.numNodes):\n",
    "        G.nodes[x]['bwcc']=c[x]\n",
    "    # Node betweenness centrality as an alternative to eigenvector centrality.\n",
    "    c = nx.betweenness_centrality(dnap.nxGraphs[w])\n",
    "    for x in range(dnap.numNodes):\n",
    "        dnap.nxGraphs[w].nodes[x]['btws']=c[x]\n",
    "\n",
    "    # Set the name of the node for future display. Append atom names to residues that have multiple nodes.\n",
    "    # It's sometimes important to do things here in Step 1 when issues of duplicate nodes might arise.\n",
    "    nodeNames={} ; nodeSegIDs={} ; nodeAtomNames={}\n",
    "    for x in G.nodes():\n",
    "        #i=dnap.nodesAtmSel[x].index\n",
    "        nodeNames[x]  = name_node(dnap, x)\n",
    "        nodeSegIDs[x] = dnap.nodesAtmSel[x].segid\n",
    "        nodeAtomNames[x] = dnap.nodesAtmSel[x].name\n",
    "    nodeNames = clarify_duplicate_nodes( nodeNames, nodeAtomNames )\n",
    "    nx.set_node_attributes(G, nodeNames, \"name\")\n",
    "    nx.set_node_attributes(G, nodeSegIDs, \"segid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save data and reduced trajectory for visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dnap.saveData(fullPathRoot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:07.074514Z",
     "start_time": "2020-04-14T20:04:06.955460Z"
    }
   },
   "outputs": [],
   "source": [
    "# This function will save a reduced DCD trajectory with the heavy atoms used for network analysis\n",
    "# A smaller trajectory can be created by choosing a \"stride\" that sub-samples the original trajectory.\n",
    "# This function will also produce a PDB file so that information on atoms and residues can be loaded to\n",
    "#    visualization software such as VMD.\n",
    "if strideDCD == 0:\n",
    "    strideDCD = dnap.workU.trajectory.n_frames-1\n",
    "print(\"We will save {} heavy atoms and {} frames.\".format(dnap.workU.atoms.n_atoms, \n",
    "                                                          len(dnap.workU.trajectory[::strideDCD]) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-14T20:04:10.709069Z",
     "start_time": "2020-04-14T20:04:07.078899Z"
    }
   },
   "outputs": [],
   "source": [
    "dnap.saveReducedTraj(fullPathRoot, stride = strideDCD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MDAnalysis may print warnings regarding missing data fields, such as altLocs, icodes, occupancies, or tempfactor, which provide information commonly found in PDB files.\n",
    "The warnings are for your information and in the context of this tutorial they are expected and do not indicate a problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis\n",
    "\n",
    "The we have finished processing the trajectory and storing all related data. We can now move on to analysis of the network properties calculated here.\n",
    "\n",
    "**All analysis code was placed in a second tutorial notebook for clarity.**\n",
    "\n",
    "\n",
    "# ---- 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
}