tutorial_allosteric_pathways_including_cofactors.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Allosteric pathways with current flow analysis on protein-cofactor networks\n",
    "\n",
    "*This tutorial shows how to build and analyze networks that include protein residues and cofactors (e.g. lipids or small molecules).*\n",
    "\n",
    "***Note***: To build and analyze a residue interaction network of the isolated protein only, just skip the steps in Section 2b and 3a, and inputs named *interactor_atom_inds_file.npy* or *additional_interactor_**.\n",
    "\n",
    "## Citing this work\n",
    "The code and developments here are described in two papers. <br>\n",
    "\n",
    "\n",
    "**[1]** P.W. Kang, A.M. Westerlund, J. Shi, K. MacFarland White, A.K. Dou, A.H. Cui, J.R. Silva, L. Delemotte and J. Cui. <br>\n",
    "*Calmodulin acts as a state-dependent switch to control a cardiac potassium channel opening*. 2020<br><br>\n",
    "**[2]** A.M. Westerlund, O. Fleetwood, S. Perez-Conesa and L. Delemotte. <br>\n",
    "*Network analysis reveals how lipids and other cofactors influence membrane protein allostery*. 2020\n",
    "\n",
    "[1] is an applications-oriented paper describing how to analyze **residue interaction networks** of **isolated proteins**. <br>\n",
    "[2] is a methods-oriented paper of how to build and analyze **residue interaction networks** that include **proteins and cofactors**.\n",
    "\n",
    "\n",
    "## Short background\n",
    "A residue interaction network is typically obtained from the element-wise product of two matrices: <br>\n",
    "&emsp; 1) Contact map. <br>\n",
    "&emsp; 2) Correlation (of node fluctuations) map.\n",
    "\n",
    "For protein residue interaction networks, the node fluctuations correspond to protein residue fluctuations around an equilibrium position [1]. The method used to build contact and correlation maps which include cofactor nodes is described in details in [2]. \n",
    "\n",
    "### Contact map\n",
    "The contact map here is defined using a truncated Gaussian kernel $K$ to smooth the contacts. For a frame with given a distance $d$ between two nodes\n",
    "\n",
    "$$\n",
    "K(d) = \n",
    "\\begin{cases} \n",
    "1 & \\text{if } d \\le c \\\\\n",
    "\\exp (-\\frac{d^2}{2\\sigma^2}) / \\exp (-\\frac{c^2}{2\\sigma^2}) & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "By default, $c=0.45$ nm and $\\sigma=0.138$ nm. <br>\n",
    "The cutoff, $c=0.45$, ensures a contact if $d \\le 4.5$ Å. The standard deviation, $\\sigma=0.138$, is chosen such that $K(d=0.8 \\text{ nm}) = 10^{-5}$. <br><br>\n",
    "The final contact map is averaged over frames.\n",
    "\n",
    "### Correlation map\n",
    "The correlation of node (protein residues in the case of isolated proteins) fluctuations is calculated using mutual information.\n",
    "\n",
    "$$\n",
    "M_{ij} = H_i + H_j - H_{ij},\n",
    "$$\n",
    "where\n",
    "$$\n",
    "H_i = -\\int\\limits_X \\rho(x)\\ln \\rho(x).\n",
    "$$\n",
    "\n",
    "$\\rho_i(x)$ is the density of distances from the node equilibrium position. This is estimated with Gaussian mixture models and the Bayesian information criterion model selection. \n",
    "\n",
    "### Including cofactors in the network\n",
    "Cofactors, such as lipids and small molecules, are treated slighlty differently than protein residues. The details are described in [2]. Practically, cofactors are processesed and added to the network in separate steps than the protein residues. The network nodes that represent cofactors are called *interactor nodes*. The following is needed to add cofactors in the network:\n",
    "\n",
    "1. **Trajectory (and .pdb file) with protein and cofactors**: If the trajectory is centered on the protein, make sure that the other molecules are not split across simulation box boundaries. In gromacs, for example, this may be avoided in *gmx trjconv* by using the option *-pbc res*. <br>\n",
    "2. **Definition of interactors**: A cofactor may be described by one or several *interactors*. An interactor could e.g. be the lipid head group. We therefore have to specify which cofactor atoms form an interactor. More details are found in Section 2b. <br>\n",
    "3. **Contact map and fluctuations**: The practical details are outlined in Sections 2b and 3a.\n",
    "\n",
    "### Current flow analysis\n",
    "The networks are analyzed using a current flow analysis [3,4] framework. The code supports both current flow betweenness and current flow closeness analysis. In short, the current flow computes the net diffusion along edges between network nodes. The net throughput of a node is given by the sum over edges. \n",
    "\n",
    "Current flow betweenness is useful for identifying allosteric pathways [5,1]. Specifically, it shows how important each residue is for transmitting allosteric pathways from a source (allosteric site) to a sink (functional site). Current flow closeness centrality [3], instead indicates signaling efficiency within the network (using a \"distance\" measured in current flow).\n",
    "\n",
    "To perform current flow analysis, you need a contact map and a similarity map (e.g. mutual information or Pearson correlation). These are computed in Section 2-3. The practical details are described in Section 4.\n",
    "\n",
    "\n",
    "## Additional references\n",
    "[3] U. Brandes and D. Fleischer, Springer, Berlin, Heidelberg, 2005 <br>\n",
    "[4] M. E. J. Newman, Social Networks, 2005 <br>\n",
    "[5] W.M. Botello-Smith and Y. Luo, J. Chem. Theory Comput., 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import allopath\n",
    "import numpy as np\n",
    "\n",
    "# Set the trajectory that should be analyzed.\n",
    "structure=['input_data/my_system.pdb']\n",
    "trajs=['input_data/system_traj1.dcd','input_data/system_traj2.dcd']\n",
    "\n",
    "# Specify how many cores to run the calculations on.\n",
    "n_cores=4\n",
    "\n",
    "# Set the output directories (out_dir is where the main data will be saved, \n",
    "# while out_dir_MI will contain the MI matrix data, see below on how they are used).\n",
    "out_dir='Results_data/'\n",
    "out_dir_MI='Results_data/MI_data/'\n",
    "\n",
    "file_label='my_system'     # Simulation label which will be appended to filenames of all written files (optional)\n",
    "dt=1                        # Trajectory stride (default=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Semi-binary contact maps\n",
    "------------------------------------------\n",
    "### 2a. Protein residue contact map\n",
    "To compute the protein (only including protein residue-residue interactions) contact map we will use _ContactMap_.\n",
    "\n",
    "***allopath.ContactMap***(**self,** *topology_file*, \\**kwargs)\n",
    "\n",
    "where *kwargs* is a dictionary with the keyword arguments (https://docs.python.org/2/glossary.html). This means that to contruct a _ContactMap_ object we have to give at least the topology_file (_structure_) as input (but in principle we want the average over a trajectory):\n",
    "> CM = allopath.ContactMap(structure)\n",
    "\n",
    "We now create a dictionary, *kwargs*, to define the named/keyword arguments that should not assume default values, such as the trajectory, ie. you may include all named input arguments that you want to modify and remove those that you wish to keep at default value.\n",
    "\n",
    "List of input keyword parameters:\n",
    "* **trajectory_files**: Input trajectory files (.xtc, .dcd, etc)\n",
    "* **trajectory_file_directory**:Input directory with trajectory files (.xtc, .dcd, etc.). This will load all trajectory files in the specified directory (this is complementary to *trajectory_files*).\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files.\n",
    "* **out_directory**: The directory where data should be written.\n",
    "* **dt**: Trajectory stride.\n",
    "* **query**: Atom-selection used on the trajectory, e.g. \"protein and !(type H)\" or \"protein and name CA\". \n",
    "* **n_cores**: Number of jobs to run with joblib.\n",
    "* **cutoff**: Cutoff value, $c$, in the truncated Gaussian kernel. For distances < cutoff, the contact will be set to one (default $c=0.45$ nm, see \"Background: contact map\" for definition).\n",
    "* **std_dev**: Standard deviation value, $\\sigma$, in the truncated Gaussian kernel. (default $\\sigma=0.138$ nm => 1e-5 contact at 0.8 nm, see \"Background: contact map\" for definition)\n",
    "* **per_frame**: Whether or not to compute contact map per frame instead of averaging over the trajectory (default=False).\n",
    "* **start_frame**: Defines which frame to start calculations from. Used in combination with *per_frame*=True.\n",
    "* **end_frame**: Defines which frame to end calculations at. Used in combination with *per_frame*=True.\n",
    "* **ref_cmap_file**: File with reference cmap (e.g. average over all frames). Is used to make computations sparse/speed up calculation. Used in combination with *per_frame*=True. \n",
    "\n",
    "The default values are: <br>\n",
    "{'trajectory_files': '', <br>\n",
    "'trajectory_file_directory': '',  <br>\n",
    "'dt': 1,  <br>\n",
    "'n_cores': 4,  <br>\n",
    "'out_directory': '',  <br>\n",
    "'file_label': '',  <br>\n",
    "'cutoff': 0.45,  <br>\n",
    "'query': 'protein and !(type H)',  <br>\n",
    "'start_frame': 0,  <br>\n",
    "'end_frame': -1,  <br>\n",
    "'ref_cmap_file': '',  <br>\n",
    "'per_frame': False,  <br>\n",
    "'std_dev': 0.138} <br>\n",
    "\n",
    "Note that the trajectory files can either be given by explicitly naming them and inputting as *trajectory_files* (as we do with _trajs_, see below), or by simply inputting a directory containing all the '.xtc' or '.dcd' files that should be analyzed (*trajectory_file_directory*).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set inputs\n",
    "kwargs={\n",
    "    'trajectory_files': trajs,\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir,\n",
    "    'dt': dt,\n",
    "    'n_cores': n_cores\n",
    "}\n",
    "\n",
    "# Compute contact map and write to file\n",
    "CM = allopath.ContactMap(structure, **kwargs)\n",
    "CM.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Interactor node - protein residue contact map\n",
    "The contact map of interactor-interactor and interactor-protein residue node contacts wille be computed using *CofactorInteractors*.\n",
    "\n",
    "\n",
    "***allopath.CofactorInteractors***(**self,** *topology_file*, \\**kwargs)\n",
    "\n",
    "The *CofactorInteractors* is used to both compute the interactions that include cofactors, and the cofactor fluctuations. The cofactor fluctuations will be used as input to the MI calculations.\n",
    "\n",
    "List of input keyword parameters to create a contact map:\n",
    "* **trajectory_files**: Input trajectory files (.xtc, .dcd, etc)\n",
    "* **trajectory_file_directory**:Input directory with trajectory files (.xtc, .dcd, etc.). This will load all trajectory files in the specified directory.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files.\n",
    "* **out_directory**: The directory where data should be written.\n",
    "* **dt**: Trajectory stride.\n",
    "* **cofactor_domain_selection**: A file containing cofactor-interactor selections. Each row should list the atoms that make up an interactor. Example of a domain selection file content: <br><br>\n",
    "*resname POPC and name N C11 C12 C13 C14 C15 P O13 O14 O12 O11 C1 C2 <br>\n",
    "resname POPC and name O21 C21 C22 O22 C23 C24 C25 C26 C27 C28 C29 C210 C211 C212 C213 C214 C215 C216 C217 C218 <br>\n",
    "resname POPC and name C3 O31 C31 C32 O32 C33 C34 C35 C36 C37 C38 C39 C310 C311 C312 C313 C314 C315 C316* \n",
    "<br><br>\n",
    "* **cutoff**: Cutoff value, $c$, for binary residue-lipid contacts. For distances < cutoff, the contact will be set to one (default=0.45 nm).\n",
    "* **std_dev**: Standard deviation value, $\\sigma$, on the semi-binary Gaussian-kernel. (default=0.138 nm => 1e-5 contact at 0.8 nm)\n",
    "\n",
    "\n",
    "The default values are: <br>\n",
    "{'trajectory_files': '', <br>\n",
    "'trajectory_file_directory': '',  <br>\n",
    "'dt': 1,  <br>\n",
    "'out_directory': '',  <br>\n",
    "'file_label': '',  <br>\n",
    "'cofactor_domain_selection': '',  <br>\n",
    "'cofactor_interactor_inds': '', <br>\n",
    "'cofactor_interactor_coords':, '', <br>\n",
    "'compute_cofactor_interactor_fluctuations': False, <br>\n",
    "'cofactor_interactor_atom_inds': '', <br>\n",
    "'cutoff': 0.45,  <br>\n",
    "'std_dev': 0.138} <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set inputs\n",
    "cofactor_domain_selection_file='input_data/cofactor_domain_selection.txt'\n",
    "\n",
    "kwargs={\n",
    "    'trajectory_files': trajs,\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir,\n",
    "    'dt': dt,\n",
    "    'cofactor_domain_selection': cofactor_domain_selection_file\n",
    "}\n",
    "\n",
    "# Compute contact map and write to file\n",
    "CI = allopath.CofactorInteractors(structure, **kwargs)\n",
    "CI.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Mutual information\n",
    "-----------------------------------\n",
    "To compute mutual information (MI) between nodes we use *MutualInformation* and *CofactorInteractors*.\n",
    "\n",
    "The MI is done in **four** steps. <br>\n",
    "**(a)** Computing the interactor node fluctuations using *CofactorInteractors*. These will be given as input to *MutualInformation*.<br>\n",
    "**(b)** Computing the off-diagonal elements in the MI matrix using *MutualInformation*. Because this is computationally demanding, we can 1) use the contact map as input to ignore non-contacting residues and 2) split the matrix into blocks that can be processed in parallel (although we will do it in sequence in this tutorial). \n",
    "> We will divide the matrix into 4 blocks along the column and 4 blocks along the rows. As we include the diagonal blocks but use symmetry on off-diagonal blocks, we get *n_matrix_block_cols=4 and *n_blocks*= n_matrix_block_cols(n_matrix_block_cols-1)/2 + n_matrix_block_cols = 10 number of blocks. The input argument *i_block* should be between 1 and *n_blocks*, denoting which block should be constructed. <br>\n",
    "\n",
    "**(c)** Computing the diagonal elements in the MI matrix using *MutualInformation*. This requires *do_diagonal*=True as input. *Note: This is only needed if you normalize the mutual information in allopath.CurrentFlow.* (Section 4)<br>\n",
    "**(d)** Building the full off-diagonal matrix based on blocks.<br><br>\n",
    "*Note:* The calculations in **(b)** and **(c)** are time consuming, but they are structured so that they can be launched in parallel. **(d)** cannot be done until the calculations in **(b)** have finished."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. Computing interactor node fluctuations\n",
    "The interactor fluctuations will be computed using *CofactorInteractors*.\n",
    "\n",
    "\n",
    "***allopath.CofactorInteractors***(**self,** *topology_file*, \\**kwargs)\n",
    "\n",
    "As mentioned, *CofactorInteractors* is used to both compute the interactions that include cofactors, and the cofactor fluctuations. To compute interactor fluctuations, we need to set **compute_cofactor_interactor_fluctuations=True** in *kwargs*.\n",
    "\n",
    "List of input keyword parameters to compute interactor fluctuations:\n",
    "* **trajectory_files**: Input trajectory files (.xtc, .dcd, etc)\n",
    "* **trajectory_file_directory**:Input directory with trajectory files (.xtc, .dcd, etc.). This will load all trajectory files in the specified directory.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files.\n",
    "* **out_directory**: The directory where data should be written.\n",
    "* **dt**: Trajectory stride.\n",
    "* **cofactor_interactor_inds**: (generated when computing the interactor node contact map).\n",
    "* **cofactor_interactor_coords**:(generated when computing the interactor node contact map).\n",
    "* **compute_interactor_node_fluctuations**: Whether or not to compute the fluctuations. Default is False. Set to True.\n",
    "\n",
    "The default values are: <br>\n",
    "{'trajectory_files': '', <br>\n",
    "'trajectory_file_directory': '',  <br>\n",
    "'dt': 1,  <br>\n",
    "'out_directory': '',  <br>\n",
    "'file_label': '',  <br>\n",
    "'cofactor_domain_selection': '',  <br>\n",
    "'cofactor_interactor_inds': '', <br>\n",
    "'cofactor_interactor_coords':, '', <br>\n",
    "'compute_cofactor_interactor_fluctuations': False, <br>\n",
    "'cofactor_interactor_atom_inds': '', <br>\n",
    "'cutoff': 0.45,  <br>\n",
    "'std_dev': 0.138} <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set inputs\n",
    "cofactor_interactor_inds = out_dir+'cofactor_interactor_indices_'+file_label+'.npy'\n",
    "cofactor_interactor_coords = out_dir+'cofactor_interactor_coords_'+file_label+'.npy'\n",
    "\n",
    "kwargs={\n",
    "    'trajectory_files': trajs,\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir,\n",
    "    'dt': dt,\n",
    "    'cofactor_interactor_inds': cofactor_interactor_inds,\n",
    "    'cofactor_interactor_coords': cofactor_interactor_coords,\n",
    "    'compute_interactor_node_fluctuations': True\n",
    "}\n",
    "\n",
    "# Compute interactor node fluctuations and write to file\n",
    "CI = allopath.CofactorInteractors(structure, **kwargs)\n",
    "CI.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b. Computing off-diagonal elements\n",
    "The MI matrix is obtained with *MutualInformation*.\n",
    "\n",
    "***allopath.MutualInformation*** (**self,** *topology_file*, \\**kwargs)\n",
    "\n",
    "Similarly to *ContactMap* and *CofactorInteractors* it is in principle enough to input the structure.\n",
    "> MI = allopath.MutualInformation(structure)\n",
    "\n",
    "\n",
    "List of input keyword parameters:\n",
    "* **trajectory_files**: Input trajectory files (.xtc, .dcd, etc)\n",
    "* **trajectory_file_directory**:Input directory with trajectory files (.xtc, .dcd, etc.). This will load all trajectory files in the specified directory.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files.\n",
    "* **out_directory**: The directory where data should be written.\n",
    "* **dt**: Trajectory stride.\n",
    "* **n_cores**: Number of jobs to run with joblib.\n",
    "* **n_matrix_block_cols**: Number of blocks of the column of the MI matrix. Example: 4 blocks => 10 parts (upper triangle + diagonal). See part (a) above.\n",
    "* **i_block**: The matrix block for which MI should be calculated. See part (a) above.\n",
    "* **n_split_sets**: Number of sampled sets with the same size as the original data set to use for more accurate estimate of entropy. Can also be used to check unceratinty of the MI matrix.\n",
    "* **additional_interactor_protein_contacts**: The interactor contact map (computed in Section 2b).\n",
    "* **additional_interactor_fluctuations**: The interactor fluctuations (computed in Section 3a).\n",
    "* **n_components_range:** Array with the lower and upper limit of GMM components used to estimate densities. \n",
    "* **do_diagonal**: Whether or not to compute diagonal of residue-residue mutual information (default=False).\n",
    "\n",
    "\n",
    "The default values are: <br>\n",
    "{'trajectory_files': '', <br>\n",
    "'trajectory_file_directory': '',  <br>\n",
    "'dt': 1,  <br>\n",
    "'out_directory': '',  <br>\n",
    "'file_label': '',  <br>\n",
    "'n_cores': -1, <br>\n",
    "'contact_map_file': '',  <br>\n",
    "'i_block': 0, <br>\n",
    "'n_matrix_block_cols': 1 <br>\n",
    "'n_split_sets': 0, <br>\n",
    "'additional_interactor_protein_contacts': '', <br>\n",
    "'additional_interactor_fluctuations': '', <br>\n",
    "'n_components_range': [1,4], <br>\n",
    "'do_diagonal': False\n",
    "} <br>\n",
    "\n",
    "To compute the off-diagonal elements, we use the default *do_diagonal*=False and split the matrix into 10 blocks. We also do 10 bootstrap samplings to obtain a better entropy estimate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_blocks = 10\n",
    "n_cols = 4\n",
    "n_bootstraps = 10\n",
    "contact_map = out_dir+'distance_matrix_semi_bin_'+file_label+'.txt'\n",
    "additional_interactor_fluctuations = out_dir+'interactor_centroid_fluctuations_'+file_label+'.npy'\n",
    "additional_interactor_protein_contacts = out_dir+'cofactor_protein_residue_semi_binary_cmap_'+file_label+'.npy'\n",
    "\n",
    "n_components_range = [1,4]\n",
    "\n",
    "for i_block in range(1,n_blocks+1):\n",
    "    # Set inputs\n",
    "    kwargs={\n",
    "        'trajectory_files': trajs,\n",
    "        'dt': dt,\n",
    "        'contact_map_file': contact_map,\n",
    "        'additional_interactor_fluctuations': additional_interactor_fluctuations,\n",
    "        'additional_interactor_protein_contacts': additional_interactor_protein_contacts,\n",
    "        'i_block': i_block,\n",
    "        'n_matrix_block_cols': n_cols,\n",
    "        'n_split_sets': n_bootstraps,\n",
    "        'n_components_range': n_components_range,\n",
    "        'file_label': file_label,\n",
    "        'out_directory': out_dir_MI,\n",
    "        'n_cores': n_cores,\n",
    "    }\n",
    "\n",
    "    # Compute mutual information matrix\n",
    "    MI = allopath.MutualInformation(structure, **kwargs)\n",
    "    MI.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c. Computing diagonal elements \n",
    "To estimate the diagonal elements, we use the same inputs as above except setting *do_diagonal*=True. Moreover, the matrix is not divided into blocks since the diagonal is much faster to compute.\n",
    "\n",
    "***Note:*** *This step is only needed if you choose to normalize the mutual information in allopath.CurrentFlow (Section 4).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set inputs\n",
    "kwargs={\n",
    "    'trajectory_files': trajs,\n",
    "    'dt': dt,\n",
    "    'additional_interactor_fluctuations': additional_interactor_fluctuations,\n",
    "    'n_split_sets': n_bootstraps,\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir_MI,\n",
    "    'n_components_range': n_components_range,\n",
    "    'n_cores': n_cores,    \n",
    "    'do_diagonal': True\n",
    "}\n",
    "\n",
    "# Compute diagonal of the MI matrix\n",
    "MI = allopath.MutualInformation(structure, **kwargs)\n",
    "MI.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3d. Building matrix from blocks \n",
    "Next, the full MI matrix is built. \n",
    "\n",
    "***allopath.from_matrix.build_matrix*** (*base_file_name*, *n_blocks*, file_label='', out_directory='')\n",
    "\n",
    "We use the same parameters as above. \n",
    "\n",
    "List of input parameters:\n",
    "* **base_file_name**: the base name of each file to be processed. This is given by *base_file_name*=*path_to_data*+'res_res_MI_part_' .\n",
    "* **n_blocks**: Total number of generated matrix blocks.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files (default is '').\n",
    "* **out_directory**: The directory where data should be written to  (default is '').\n",
    "\n",
    "The input *base_file_name* is named after the files in \"Results_data/MI_data/\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_file_name=out_dir+'MI_data/res_res_MI_part_'\n",
    "\n",
    "# Set inputs\n",
    "kwargs={\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir+'MI_data/'\n",
    "}\n",
    "\n",
    "# Build matrix\n",
    "allopath.from_matrix_blocks.build_matrix(base_file_name, n_blocks, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Current flow analysis\n",
    "-----------------------------------\n",
    "Current flow analysis is done with *CurrentFlow*.\n",
    "\n",
    "***allopath.CurrentFlow*** (**self,** *similarity_map_filename*, *contact_map_filenames*, *sources_filename*, *sinks_filename*, \\**kwargs)\n",
    "\n",
    "To run current flow analysis in its simplest form, the files containing the similarity map (ie. our MI matrix), the contact map and the source and sink indices are needed.\n",
    "\n",
    "> allopath.CurrentFlow(similarity_map_filename, contact_map_filename, sources_filename, sinks_filename)\n",
    "\n",
    "Explanation of input (positional) parameters:\n",
    "* **similarity_map_filename**: File containing the similarity map (ie. the mutual information matrix).\n",
    "* **contact_map_filenames**: File containing the contact map(s). If multiple are given, one current flow profile per contact map will be computed (*Note: multiple network calculations are only supported for isolated-protein networks*).\n",
    "* **sources_filename**: File containing the residue indices of the sources.\n",
    "* **sinks_filenams**: File containing the residues indices of the sinks.\n",
    "\n",
    "Explanation of input keyword parameters:\n",
    "* **similarity_map_diagonal_filename**: File containing the diagonal elements of the mutual information matrix.\n",
    "* **additional_interactor_protein_contacts**: The interactor contact map (computed in Section 2b).\n",
    "* **out_directory**: The directory where data should be written.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced files.\n",
    "* **n_chains**: The number of (homomeric) chains/subunits in the main-protein (e.g. a tetrameric ion channel => n_chains = 4).\n",
    "* **n_cores**: Number of jobs to run with joblib.\n",
    "* **cheap_write**: If set to True, fewer files will be written.\n",
    "* **start_frame**: Used if multiple contact maps are supplied. *start_frame* is the index of the first frame to analyze.\n",
    "* **normalize_similarity_map**: Whether or not to normalize the similarity map with symmetric unertainty (*Note: applies to mutual information maps; Witten & Frank, 2005*)\n",
    "* **auxiliary_protein_indices**: Residue indices of auxiliary subunits. This is used when symmeterizing current flow over subunits (chains). The auxiliary subunits will also be averaged over chains, ie. one auxiliary subunit per chain is assumed. If there is no auxiliary subunit, just ignore this input to the current flow script.\n",
    "* **compute_current_flow_closeness**: Whether or not to compute current flow closeness instead of current flow betweenness.\n",
    "\n",
    "\n",
    "The default values are: <br>\n",
    "{'out_directory': '',  <br>\n",
    "'file_label': '',  <br>\n",
    "'similarity_map_diagonal_filename': '', <br>\n",
    "'n_chains': 1, <br>\n",
    "'n_cores': 1, <br>\n",
    "'cheap_write': False,  <br>\n",
    "'start_frame': 0, <br>\n",
    "'normalize_similarity_map': False, <br>\n",
    "'auxiliary_protein_indices': '', <br>\n",
    "'additional_interactor_protein_contacts': '',  <br>\n",
    "'compute_current_flow_closeness': False } <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "similarity_map = out_dir+'MI_data/res_res_MI_compressed_'+file_label+'.npy'\n",
    "similarity_map_diagonal = out_dir+'MI_data/diagonal_MI_'+file_label+'.npy'\n",
    "contact_maps = [out_dir+'distance_matrix_semi_bin_'+file_label+'.txt']\n",
    "additional_interactor_protein_contacts = out_dir+'cofactor_protein_residue_semi_binary_cmap_'+file_label+'.npy' \n",
    "\n",
    "n_chains=4\n",
    "\n",
    "source_inds='input_data/inds_sources.txt' \n",
    "sink_inds='input_data/inds_sinks.txt'   \n",
    "aux_inds='input_data/auxiliary_prot_inds.txt' \n",
    "\n",
    "compute_current_flow_closeness = False  # False (ie. default) => will compute current flow betweenness. \n",
    "                                        # Set this to True to compute current flow closeness centrality between each\n",
    "                                        # source and all sinks instead.\n",
    "\n",
    "kwargs={\n",
    "    'file_label': file_label,\n",
    "    'out_directory': out_dir,\n",
    "    'n_chains': n_chains,\n",
    "    'n_cores': n_cores,\n",
    "    'similarity_map_diagonal_filename': similarity_map_diagonal,\n",
    "    'normalize_similarity_map': False,\n",
    "    'auxiliary_protein_indices': aux_inds,\n",
    "    'additional_interactor_protein_contacts': additional_interactor_protein_contacts,\n",
    "    'compute_current_flow_closeness': compute_current_flow_closeness\n",
    "}\n",
    "\n",
    "CF = allopath.CurrentFlow(similarity_map, contact_maps, source_inds, sink_inds, **kwargs)\n",
    "CF.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Project current flow on structure \n",
    "----------------------------------------------------\n",
    "As a last step, we project the current flow onto the structure (PDB file) with *make_pdb*. The current flow of ech residue will be mapped to the beta-column in the PDB. This can be visualized in VMD by setting the \"Coloring method\" to \"Beta\" in \"Graphical Representations\".\n",
    "\n",
    "> ***allopath.make_pdb.project_current_flow***(*pdb_file*, *current_flow_file*, \\**kwargs)\n",
    "\n",
    "Explanation of input (positional arguments) parameters:\n",
    "* **pdb_file**: The .pdb file corresponding to the first trajectory frame. *Note: .gro does not work.*\n",
    "* **current_flow_file**: File containing the current flow. This is created by *CurrentFlow*, Section 4. **Note:** For homomultimers (using *n_chains > 1* in *CurrentFlow*), the file is  *out_dir+'average_current_flow_'+file_label+'.npy'*. For *n_chains = 1*, the file is *out_dir+'current_flow_betweenness_'+file_label+'.npy'*.\n",
    "\n",
    "Explanation of input keyword arguments:\n",
    "* **out_directory**: The directory where pdb should be written.\n",
    "* **file_label**: \"File end name\": label of the system that will be appended to the end of the produced pdb.\n",
    "* **max_min_normalize**: Whether or not to scale the current flow between 0 and 1.\n",
    "* **interactor_atom_inds_file**: The atom indices used to define the interactors (generated in Section 2b).\n",
    "\n",
    "The default values are: <br>\n",
    "{'out_directory': '', <br>\n",
    "'file_label': '', <br>\n",
    "'max_min_normalize': False,<br>\n",
    "'interactor_atom_inds_file': None }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_file = out_dir+'PDBs/current_flow_'+file_label+'.pdb'\n",
    "current_flow = out_dir+'average_current_flow_'+file_label+'.npy'\n",
    "\n",
    "interactor_atom_inds_file = out_dir+'cofactor_interactor_atom_indices_'+file_label+'.npy'\n",
    "\n",
    "kwargs={\n",
    "    'out_directory': out_dir+'PDBs/',\n",
    "    'file_label': file_label,\n",
    "    'interactor_atom_inds_file': interactor_atom_inds_file\n",
    "}\n",
    "\n",
    "# Create PDB with current flow values on the beta column\n",
    "allopath.make_pdb.project_current_flow(structure[0], current_flow, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}