Connectome-analysis-JGU-course-exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a61b7712",
   "metadata": {},
   "source": [
    "Disclaimer\n",
    "-- The following script requieres internet connection to access the data set --"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39c351c1-1cac-42ab-85f1-83e20bf9491c",
   "metadata": {},
   "source": [
    "# Analysis of connectomics data in *Drosophila Melanogaster*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "855d051b-2cdc-4f94-889c-27c9cbde4ff3",
   "metadata": {},
   "source": [
    "In the following cells of code, we will access, analyse and plot electrone microscopy data from the Drosophila brain. For that, we will use a published \n",
    "<a href=\"https://neuprint-examples.janelia.org/?dataset=medulla7column&qt=findneurons&qr%5B0%5D%5Bcode%5D=fn&qr%5B0%5D%5Bds%5D=medulla7column&qr%5B0%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B0%5D%5Bpm%5D%5Binput_ROIs%5D%5B0%5D=proximal&qr%5B0%5D%5Bpm%5D%5Boutput_ROIs%5D%5B0%5D=proximal&qr%5B0%5D%5BvisProps%5D%5BrowsPerPage%5D=25&qr%5B1%5D%5Bcode%5D=ng&qr%5B1%5D%5Bds%5D=medulla7column&qr%5B1%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B1%5D%5Bpm%5D%5Bskip%5D=true&qr%5B1%5D%5Bpm%5D%5BbodyIds%5D=72142&qr%5B2%5D%5Bcode%5D=sk&qr%5B2%5D%5Bds%5D=medulla7column&qr%5B2%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B2%5D%5Bpm%5D%5Bskip%5D=true&qr%5B2%5D%5Bpm%5D%5BbodyIds%5D=72142&qr%5B2%5D%5Bpm%5D%5Bcoordinates%5D=-1857.405450505682%2C8945.456303462368%2C5389.841777784641%2C4390.418714324823%2C3113.443324386624%2C4544.509250182286&qr%5B3%5D%5Bcode%5D=sc&qr%5B3%5D%5Bds%5D=medulla7column&qr%5B3%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B3%5D%5Bpm%5D%5Bfind_inputs%5D=true&qr%5B3%5D%5Bpm%5D%5Bneuron_id%5D=72142&qr%5B3%5D%5BvisProps%5D%5BpaginateExpansion%5D=true&tab=0&ftab=&q=1\" target=\"_blank\">data set</a>\n",
    "containing information of seven medulla columns in the fly optic lobe, which is maintained by Janelia Research Campus. \n",
    "\n",
    "\n",
    "But before that, a bit of background:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "413331ee-6c77-488d-a870-b94018f62e91",
   "metadata": {},
   "source": [
    "## Introduction to the data set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ff6b5e9-b5ea-4806-ae88-56bd7ae5463b",
   "metadata": {},
   "source": [
    "The data set was first published by Takemura et al, 2013 and later used in publications in Takemura et al, 2015, and 2017. The following images highlights the brain area of the connectome: **the medulla** in the optic lobe (left) and the reconstruction volume (right) [Takemura et al., 2013]\n",
    "\n",
    "<img src=\"./images/Takemura1.png\"\n",
    "     alt=\"Optic lobe\"\n",
    "     width=\"500\" height=\"300\" style=\"display: block; margin: 0 auto\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f80725f2-e975-41f4-980a-36f3eacbf81e",
   "metadata": {},
   "source": [
    "The data set consist of thousands of micrometric brain slices (**a**) which has been imaged using electrone microscopy (EM). In each 2D slice, EM experts recognized and annotated many features, for example, synapses between neurons (**c**), and merged all segments that belong together in 2D (**b**) and 3D (**d**) to obtain all the neuronal reconstruaction.\n",
    "<img src=\"./images/Takemura2.png\"\n",
    "     alt=\"EM reconstruction\"\n",
    "     width=\"350\" height=\"350\" style=\"display: block; margin: 0 auto\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13f88ed2-c2b3-440e-8e73-b8dd48fe5847",
   "metadata": {},
   "source": [
    "All synaptic connections between different neuronal types can be summarized in connectivity tables (**a**) and the strength of a connection, based on the total or relative number of synapses, can be used to visualize the neuronal connectivity in different ways, for example, as connectivity clusters or connection fractions (**b**)\n",
    "\n",
    "<img src=\"./images/Takemura3.png\"\n",
    "     alt=\"Connectivity matrix\"\n",
    "     width=\"400\" height=\"400\" style=\"float: left ; width: 45%\"/>\n",
    "<img src=\"./images/Takemura4.png\"\n",
    "     alt=\"CLuster analysis\"\n",
    "     width=\"300\" height=\"400\" style=\"float: rigth ; width: 45%\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31dcab99-be42-4f01-9d76-274d76f2dbea",
   "metadata": {},
   "source": [
    "## Loading the data set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db334d0f-0a78-45a3-b9c2-863561fea581",
   "metadata": {},
   "source": [
    "The dataset **FIB-25** contains the raw EM data for the 7 medulla column FIB-SEM reconstruction, as reported in the paper *Synaptic circuits and their variations within different columns in the visual system of Drosophila* by Takemura.\n",
    "\n",
    "For accessing the dataset, you need a google account and a personal TOKEN.\n",
    "\n",
    "To obtain your authorization token, follow these steps:\n",
    "\n",
    "1. Navigate your web browser to the neuprint server\n",
    "<a href=\"https://neuprint-examples.janelia.org/?dataset=medulla7column&qt=findneurons&qr%5B0%5D%5Bcode%5D=fn&qr%5B0%5D%5Bds%5D=medulla7column&qr%5B0%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B0%5D%5Bpm%5D%5Binput_ROIs%5D%5B0%5D=proximal&qr%5B0%5D%5Bpm%5D%5Boutput_ROIs%5D%5B0%5D=proximal&qr%5B0%5D%5BvisProps%5D%5BrowsPerPage%5D=25&qr%5B1%5D%5Bcode%5D=ng&qr%5B1%5D%5Bds%5D=medulla7column&qr%5B1%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B1%5D%5Bpm%5D%5Bskip%5D=true&qr%5B1%5D%5Bpm%5D%5BbodyIds%5D=72142&qr%5B2%5D%5Bcode%5D=sk&qr%5B2%5D%5Bds%5D=medulla7column&qr%5B2%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B2%5D%5Bpm%5D%5Bskip%5D=true&qr%5B2%5D%5Bpm%5D%5BbodyIds%5D=72142&qr%5B2%5D%5Bpm%5D%5Bcoordinates%5D=-1857.405450505682%2C8945.456303462368%2C5389.841777784641%2C4390.418714324823%2C3113.443324386624%2C4544.509250182286&qr%5B3%5D%5Bcode%5D=sc&qr%5B3%5D%5Bds%5D=medulla7column&qr%5B3%5D%5Bpm%5D%5Bdataset%5D=medulla7column&qr%5B3%5D%5Bpm%5D%5Bfind_inputs%5D=true&qr%5B3%5D%5Bpm%5D%5Bneuron_id%5D=72142&qr%5B3%5D%5BvisProps%5D%5BpaginateExpansion%5D=true&tab=0&ftab=&q=1\" target=\"_blank\">address</a>.\n",
    "2. Log in.\n",
    "3. Using the account menu in the upper right-hand corner, select “Account” as shown in the screenshot below.\n",
    "4. Copy the entire auth token."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ecfc024-af4d-4708-a944-b825ceac11e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuprint import Client\n",
    "\n",
    "TOKEN = 'eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJlbWFpbCI6InNlYmFzdGlhbi5tb2xpbmEub2JhbmRvQGdtYWlsLmNvbSIsImxldmVsIjoibm9hdXRoIiwiaW1hZ2UtdXJsIjoiaHR0cHM6Ly9saDMuZ29vZ2xldXNlcmNvbnRlbnQuY29tL2EtL0FPaDE0R2hWZjMxR2RHeURzYmtfUW5qdW00b1U4SVZ5QTBEcXNVaXdNZ1ZrTEE9czk2LWM_c3o9NTA_c3o9NTAiLCJleHAiOjE4MzA5NTQ5MTB9.YUJV-C5VoOZ-huLvc73EhWa6KWnejmemqzl9V-OrBKs'\n",
    "c = Client('neuprint-examples.janelia.org', dataset='medulla7column', token= TOKEN)\n",
    "c.fetch_version()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc29b55f",
   "metadata": {},
   "source": [
    "<h1 style=\"text-align: center;\">NeuPrint model </h1>\n",
    "<img src=\"./images/Nodes_in_NeuPrint.png\"\n",
    "     alt=\"NeuPrint model\"\n",
    "     width=\"300\" height=\"400\" style=\"display: block; margin: 0 auto\" />\n",
    "\n",
    "<figcaption align = \"center\">The data is embedded in a graph model from which relationships between elements can be retrieved \n",
    "</figcaption>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace05ef3",
   "metadata": {},
   "source": [
    "Relationships between nodes are fetched using clauses of the \n",
    "<a href=\"https://neo4j.com/docs/cypher-manual/current/clauses/\" target=\"_blank\">Cypher query language</a> developed by Neo4j. \n",
    "\n",
    "**Neo4j** is a graph database that Neuprint is built on, and stores massive amounts of data as linked lists\n",
    "and index tables. Tables are dataframes and can be modified or analyzed using <a href= \"https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.html\" target=\"_blanck\"> pandas.</a>\n",
    "\n",
    "More information about Neo4j and NeuPrint and examples of cypher clauses can be found <a href=\"https://neuprint.janelia.org/public/neuprintuserguide.pdf\" target=\"_blank\"> here. </a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8c0bad1",
   "metadata": {},
   "source": [
    "## Fetching Neurons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ca5d9a3-7f3e-4deb-8933-23ac5fa48a54",
   "metadata": {},
   "source": [
    "In **NeuPrint**, each neuron is annotated with the list of regions (ROIs) it intersects, along with the synapse counts in each.\n",
    "The ROIs comprise a hierarchy, with smaller ROIs nested within larger ROIs. \n",
    "Furthermore, primary ROIs are guaranteed not to overlap, and they roughly tile the entire brain (with some gaps)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f1eeb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuprint import fetch_roi_hierarchy\n",
    "\n",
    "# Show the ROI hierarchy, with primary ROIs marked with '*'\n",
    "print(fetch_roi_hierarchy(False, mark_primary=True, format='text'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3954f7f2",
   "metadata": {},
   "source": [
    "Lets use now **Cypher clauses** to access an specific neuron as node in the data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8de3d38",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Fetching a specific type of neurons with cypher clauses ##\n",
    "\n",
    "#This query (q) will return all neurons in the dataset\n",
    "#that of type = 'Tm9'\n",
    "#Results are ordered by total synaptic sites (pre+post).\n",
    "q = \"\"\"\\\n",
    "    MATCH (n :Neuron)\n",
    "    WHERE n.type = 'Tm9'\n",
    "    RETURN n.bodyId AS bodyId, n.type as type, n.instance AS instance, n.pre AS numpre, n.post AS numpost\n",
    "    ORDER BY n.pre + n.post DESC\n",
    "\"\"\"\n",
    "results = c.fetch_custom(q)\n",
    "display(results)\n",
    "display(f'Data type: {type(results)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a17ce2c",
   "metadata": {},
   "source": [
    "<h3 style=\"text-align: center;\">The connectome for direction selectivity (DS)</h3>\n",
    "\n",
    "<img src=\"./images/Takemura5.png\"\n",
    "     alt=\"The connectome for direction selectivity\"\n",
    "     width=\"500\" height=\"500\" style=\"display: block; margin: 0 auto\" />\n",
    "\n",
    "<figcaption align = \"center\"><b>The connectome for direction selectivity (DS).</b> \n",
    "<br>A) Hassenstein and Reichert model for DS. B) Balow-Lewick model for DS. C) Neurons involved in the direction-selective detection of moving signals. D) Seven columns in the medulla connectome. Around the centered \"home\" columns there are six more named by different letters.</figcaption>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a379078b",
   "metadata": {},
   "source": [
    "<h3 style=\"text-align: center;\">Synaptic site:</h3>\n",
    "\n",
    "<img src=\"./images/Synaptic_site.png\"\n",
    "     alt=\"Synaptic site\"\n",
    "     width=\"300\" height=\"300\" style=\"display: block; margin: 0 auto\" />\n",
    "\n",
    "<figcaption align = \"center\"><b>Example of a synaptic site in the Drosophila.</b> \n",
    "<br>The arrow points to the pre-synaptic site called a T-bar. \n",
    "<br>Multiple outputs, called PSDs,exist, indicated by the triangles.</figcaption>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "991b74ed",
   "metadata": {},
   "source": [
    "<span style='color:Blue'> **EXERCISE 1:**\n",
    "<br> - How many unique Tm9s do exist in the dataset?\n",
    "<br> - Which is the mean number of postynaptic sites of a Tm9 neuron?\n",
    "<br> - Is there variability in the number of total postsynaptic sites across columns?\n",
    "<br><br>TIP: use pandas dataframe functions </span> \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab69f710",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unique Tm9s\n",
    "un_Tm9 = results[\"instance\"].unique()\n",
    "print(f\"There is a total of {len(un_Tm9)} unique Tm9\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de65d412",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean postsynaptic sites\n",
    "mean_post = results[\"numpost\"].mean()\n",
    "print(f\"The mean #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {mean_post}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8ad9714",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is there variability across columns?\n",
    "std_post = results[\"numpost\"].std()\n",
    "print(f\"The std #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {std_post}\")\n",
    "\n",
    "sem_post = results[\"numpost\"].sem()\n",
    "print(f\"The sem #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {sem_post}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d14293d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Just the main seven columns\n",
    "q = \"\"\"\\\n",
    "    \n",
    "    MATCH (n :Neuron)\n",
    "    WHERE n.bodyId IN [19642, 137381, 48678, 17016, 74030, 74625,11783]\n",
    "    RETURN n.bodyId AS bodyId, n.type as type, n.instance AS instance, n.pre AS numpre, n.post AS numpost\n",
    "    ORDER BY n.pre + n.post DESC\n",
    "\"\"\"\n",
    "results = c.fetch_custom(q)\n",
    "display(results)\n",
    "display(f'Data type: {type(results)}')\n",
    "\n",
    "un_Tm9 = results[\"instance\"].unique()\n",
    "\n",
    "# Mean postsynaptic sites\n",
    "mean_post = results[\"numpost\"].mean()\n",
    "print(f\"The mean #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {mean_post}\")\n",
    "# Is there variability across columns?\n",
    "std_post = results[\"numpost\"].std()\n",
    "print(f\"The std #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {std_post}\")\n",
    "\n",
    "sem_post = results[\"numpost\"].sem()\n",
    "print(f\"The sem #of postsynaptic sites across {len(un_Tm9)} Tm9s is:  {sem_post}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff662f03",
   "metadata": {},
   "source": [
    "When specifying the content of a cypher clause, one can use \n",
    "<a href=\"https://docs.python.org/3/library/re.html\" target=\"_blanck\"> regular expressions (regex). </a>\n",
    "A regular expression is a string of characters defining a search pattern. Specific characters are reserved for special use.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dc8805b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## This query will return all neurons in the dataset\n",
    "## that of type containing '.*XXX.*' in the string (Using regular expressions)\n",
    "## Results are ordered by total synaptic sites (pre+post).\n",
    "\n",
    "q = \"\"\"\\\n",
    "    MATCH (n :Neuron)\n",
    "    WHERE n.instance =~ '.* home'\n",
    "    RETURN n.bodyId AS bodyId, n.type as type, n.instance AS instance, n.pre AS numpre, n.post AS numpost\n",
    "    ORDER BY n.pre + n.post DESC\n",
    "\"\"\"\n",
    "results = c.fetch_custom(q)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0e5c066",
   "metadata": {},
   "source": [
    "<span style='color:Blue'> **EXCERCISE 2:**\n",
    "<br>- Analyse the variability of the total number of neurons present in the main seven columns of the dataset.\n",
    "<br>- Which is the columns with more unique reconstructed neurons?\n",
    "<br><br>TIP: write a for loop to query the different columns using regular expressions</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f257755",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Excercise\n",
    "\n",
    "q_all = \"\"\"\\\n",
    "    MATCH (n :Neuron)\n",
    "    WHERE n.instance =~ '.*{}'\n",
    "    RETURN n.bodyId AS bodyId, n.type as type, n.instance AS instance, n.pre AS numpre, n.post AS numpost\n",
    "    ORDER BY n.pre + n.post DESC\n",
    "\"\"\"\n",
    "l = []\n",
    "\n",
    "for neuron_columns in [\"home\",\"-A\", \"-B\",\"-C\",\"-D\",\"-E\",\"-F\"]:\n",
    "    results = c.fetch_custom(q_all.format(neuron_columns))\n",
    "    l.append(len(results))\n",
    "    \n",
    "print(l)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3ef0f81",
   "metadata": {},
   "source": [
    "## Fetching connections"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "676416a0",
   "metadata": {},
   "source": [
    "<img src=\"./images/Connects_to_NeuPrint.png\"\n",
    "     alt=\"Connection between neurons\"\n",
    "     width=\"500\" height=\"150\" style=\"display: block; margin: 0 auto\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f847996",
   "metadata": {},
   "source": [
    "The usefullness of a connectome is not only to map the different neurons (nodes in the graph) that exist at specific locations but the relationships (connections, edges in the graph) that exist between them, meaning how neurons influence each other in the circuit. One of the ways for a neuron to influence another (e.g., to excite or inhibite) is via a synaptic connection. \n",
    "\n",
    "Pairs of neurons build one or more synaptic connections between each other. It is assumed that the higher the number of synapses between two neurons the stronger their influence of a presynaptic neuron to its postsynaptic partner. This strenght of connection is usually refered as a connnection *weitgh*. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "042b189d",
   "metadata": {},
   "source": [
    "Let´s query from the dataset the synaptic weigths between a specific neuron and all its **presynaptic** partners: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f872142",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "q = \"\"\"\\\n",
    "MATCH (a:Neuron)-[w:ConnectsTo]->(b:Neuron)\n",
    "WHERE b.bodyId = 19642\n",
    "RETURN a.instance as instance, b.bodyId, b.instance, w.weight AS W\n",
    "ORDER BY W DESC\n",
    "\"\"\"\n",
    "results = c.fetch_custom(q)\n",
    "display(results.head())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa7dd9bc",
   "metadata": {},
   "source": [
    "<span style='color:Blue'> **EXCERCISE 3:**\n",
    "<br> - Analyse the variability of the main **seven** input weitghs to Tm9 across the **seven** columns. \n",
    "<br>- Plot both the absolut and the relative weitgh of these inputs in a bar plot with error bars or in a box plot.\n",
    "<br><br>TIP: Use a loop to first save the values of the different columns and store these on a dataframe, list or dictionary.\n",
    "<br>TIP: For plotting, you could use the <a href=\"https://matplotlib.org/\" target=\"_blanck\">  matplotlib</a> library or the <a href=\"https://seaborn.pydata.org//\" target=\"_blanck\">  seaborn</a> library.</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb4d0fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "#Excercise\n",
    "\n",
    "id_list = [19642, 137381, 48678, 17016, 74030, 74625,11783]\n",
    "\n",
    "q = \"\"\"\\\n",
    "MATCH (a:Neuron)-[w:ConnectsTo]->(b:Neuron)\n",
    "WHERE b.bodyId = {}\n",
    "RETURN a.instance as instance, b.bodyId, b.instance, w.weight AS W\n",
    "ORDER BY W DESC\n",
    "\"\"\"\n",
    "\n",
    "l_l = []\n",
    "i_l = []\n",
    "\n",
    "for ID in id_list:\n",
    "    results = c.fetch_custom(q.format(ID))\n",
    "    in_list = results[\"instance\"][0:7].tolist()\n",
    "    w_list = results[\"W\"][0:7].tolist()\n",
    "    l_l.append(w_list)\n",
    "    i_l.append(in_list)\n",
    "    \n",
    "\n",
    "#print(l_l)\n",
    "#print(i_l)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd6567d",
   "metadata": {},
   "source": [
    "## Data visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eb5725a",
   "metadata": {},
   "source": [
    "### The DS ON-pathway"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae6d1937",
   "metadata": {},
   "source": [
    "There are many ways to visualize connectomics data, from absolut or relative weitghs between neuronal pairs to representation of whole circuit connections. Next, we will explore some ways of visualizing data in order to uncover main input or output neurons of a the **direction selective (DS) circuit for ON signals**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a37d0b8c",
   "metadata": {},
   "source": [
    "The DS ON-pathway is defined as the collecion of neurons that are involved in the detection of moving light increments (ON). It is a contrast selective pathway with the first ON-selective neurons located in the optic lobe medulla. In a nutshell, photoreceptors in the retina depolarize to ON and hyperpolarize to OFF, transmitting both information of ON and OFF signals to lamina (L) neurons. The first order interneurons L1,L2 and L3 in the medulla invert the sign of the photoreceptor input and all depolarize to OFF and hyperpolarize to ON. Then, these L neurons pass their signals to medulla neurons which become selective and depolarize either to ON or to OFF.\n",
    "\n",
    "Medulla neurons in the DS ON-pathway, including the ON-depolarizing Mi1, Tm3 or Mi4, and the OFF-depolarizing Mi9 have in common that they respond to a contrast change, are not direction selective, and are main input of T4, the first ON direction selective neuron. \n",
    "\n",
    "<img src=\"./images/ON-pathway-DS-Strother-2017.png\"\n",
    "     alt=\"DS ON-Pathway\"\n",
    "     width=\"800\" height=\"400\" style=\"display: block; margin: 0 auto\"/>\n",
    "\n",
    "<figcaption align = \"center\"><b>T4 Neurons Receive an Array of Input Signals.</b> \n",
    "<br>(A) A schematic model for the fly visual system showing the primary columnar neuron types that contribute to the computation of directionally selective signals. \n",
    "<br>(B) Two classical models for the computation of directional selectivity from non-selective input signals.\n",
    "<br> Figure from Strother et al., 2017</figcaption>\n",
    "\n",
    "T4 cells project their dendrites in the proximal medulla and their axon terminals in the lobula plate, where they end up in four different layers that respond to four different directions of motions for a given z-depth, having in totall 6 different cell populations covering 6 different directions of motion. These differential axonal projection into four different layers has allowed to subdivide T4 cells in four subtypes: T4a, T4b, T4c and T4d.\n",
    "\n",
    "\n",
    "<img src=\"./images/T4-T5_six_subtypes_henning_etal_2022.png\"\n",
    "     alt=\"T4-T5 six subtypes\"\n",
    "     width=\"500\" height=\"250\" style=\"float: left ; width: 49%\"/>\n",
    "<img src=\"./images/T4-T5_four_subtypes_per_layer_henning_etal_2022.png\"\n",
    "     alt=\"T4-T5 four subtypes per layer\"\n",
    "     width=\"500\" height=\"250\" style=\"float: right ; width: 49%\" />\n",
    "\n",
    "<figcaption align = \"center\"><b>T4- T5 subtypes</b> \n",
    "<br>Left: Circular histograms of tuning preference. Black vectors depict average tuning per layer.\n",
    "<br>Right: One layer of the lobula plate with ROIs color-coded according to their directional tuning. Note that this is a different color code as compared to the one highlighting subtypes used in previous figures.\n",
    "<br> Figure  panels from Henning et al., 2022</figcaption>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29e07b47",
   "metadata": {},
   "source": [
    "## Data visualization: the OFF-pathway"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac0e58b3",
   "metadata": {},
   "source": [
    "### The DS OFF-pathway"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1317fe4",
   "metadata": {},
   "source": [
    "There are many ways to visualize connectomics data, from absolut or relative weitghs between neuronal pairs to representation of whole circuit connections. Next, we will explore some ways of visualizing data in order to uncover main connectivity of those neurons belonging to the direction selective circuit for light decrements (OFF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "884d4ffc",
   "metadata": {},
   "source": [
    "The DS OFF-pathway is defined as the collecion of neurons that are involved in the detection of moving light decrements (OFF). It is a contrast selective pathway with the first OFF-selective neurons located in the optic lobe medulla. In a nutshell, photoreceptors in the retina depolarize to ON and hyperpolarize to OFF, transmitting both information of ON and OFF signals to lamina (L) neurons. The first order interneurons L1,L2 and L3 in the medulla invert the sign of the photoreceptor input and all depolarize to OFF and hyperpolarize to ON. Then, these L neurons pass their signals to medulla neurons which become selective and depolarize either to ON or to OFF.\n",
    "\n",
    "Medulla neurons in the DS OFF-pathway, including the OFF-depolarizing Tm1, Tm2, Tm4, Tm9 and CT1, have in common that they respond to a contrast change, are not direction selective, and are main inputs of T5 cells, the first OFF direction selective neurons. \n",
    "\n",
    "Since T5 cells are located in the lobula and lobula plate neuropils, they are missing in the medulla connectome from the FIB25 data set. Nonetheless, one can explore the main input synapses of OFF medulla neurons and try to predict their role in circuitry.\n",
    "\n",
    "The **FIB19** data set,  reported in the paper \"Comparisons between the ON- and OFF-edge motion pathways in the Drosophila brain\" by Shinomiya, et al. 2019, contains detailed connectivity information for inputs to both the T4 and T5 pathways. This data set can be acced in the same way as the FIB25/7-medulla columns using NeuPrint:\n",
    "\n",
    "```python\n",
    "c = Client('neuprint.janelia.org', dataset='fib19', token= TOKEN)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef029324",
   "metadata": {},
   "source": [
    "<img src=\"./images/EM_reconstruction_shinomiya_2019.jpg\"\n",
    "     alt=\"FIB19 data set\"\n",
    "     width=\"600\" height=\"300\" style=\"display: block; margin: 0 auto\"/>\n",
    "\n",
    "<figcaption align = \"center\"><b>Horizontal section of the four optic lobe neuropils and their reconstructed neurons.</b> \n",
    "<br>(A) In the ON-edge pathway, the T4 cells receive inputs at the tenth layer of the medulla (M10) from upstream cells, whereas in the OFF-edge pathway, the T5 cells receive inputs at the first layer of the lobula (Lo1), before projecting to the lobula plate. Here, in a horizontal re-sliced section of the dataset, Mi1 and Tm1/Tm2 cells are shown as representative input neurons to the T4 and T5 cells. LA, lamina; ME, medulla; LO, lobula; LOP, lobula plate. \n",
    "<br> (B) Axon terminals of Tm cells and T5 cells in Lo1. The area roughly corresponds to the boxed inset in (A). T5 cell dendrites (in yellow) demarcate the Lo1 layer. Terminals of Tm1, Tm2, Tm9, and CT1 neurons form packed columnar units and have extensive synaptic contacts with T5 cells. Tm4 cells form no part of these columnar units but run between the columns, and also have synapses with T5 cells in Lo1. \n",
    "<br> (C) 3D reconstruction of the neurons. One representative T5 cell (yellow) is shown along with Tm cells and a CT1 terminal. The proximal tip of the axon terminal of Tm2 (arrowhead) extends slightly beyond the boundary between Lo2 and Lo1. Tm4 extends its axon further and terminates in Lo4 (Fischbach and Dittrich, 1989).</figcaption>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c301555a",
   "metadata": {},
   "source": [
    "### Connectivity matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c54befd9",
   "metadata": {},
   "source": [
    "For a given circuit (e.g, all  neurons on a given medulla column), a usefull way to fisrt visualize groups of interconnected neurons is to plot a connectivity matrix. This is is no more as a heat-map depicting the absolute or relative weitgh between pairs of neurons in a neuronal 2D matrix containing all neuronal types.\n",
    "\n",
    "<img src=\"./images/Takemura3.png\"\n",
    "     alt=\"Connectivity matrix\"\n",
    "     width=\"400\" height=\"400\" style=\"display: block; margin: 0 auto\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a315773",
   "metadata": {},
   "source": [
    "Next, we will generate a connectivity map for the home-column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f05f74aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "## This query will return all neurons in the dataset\n",
    "## that of type containing '.*XXX.*' in the string (Using regular expressions)\n",
    "## Results are ordered by total synaptic sites (pre+post).\n",
    "\n",
    "_column = 'home_column'\n",
    "\n",
    "q = \"\"\"\\\n",
    "    MATCH (n :Neuron)\n",
    "    WHERE n.instance =~ '.* home'\n",
    "    RETURN n.bodyId AS bodyId, n.type as type, n.instance AS instance, n.pre AS numpre, n.post AS numpost\n",
    "    ORDER BY n.pre + n.post DESC\n",
    "\"\"\"\n",
    "results = c.fetch_custom(q)\n",
    "display(results.head())\n",
    "# Keeping a list of all neurons of interest to use it later\n",
    "neuron_ids = results['bodyId'].tolist()\n",
    "neuron_type = results['type'].tolist()\n",
    "neuron_instance = results['instance'].tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47a41f9b",
   "metadata": {},
   "source": [
    "Neuprint has functions to get all connectivity information of a group of neurons via their IDs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94af7b66",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Getting the connectivity data\n",
    "\n",
    "from neuprint import fetch_adjacencies,  merge_neuron_properties\n",
    "\n",
    "# Example: Fetch all upstream and dowsntream connections of a set of neurons # AKA INPUTS AND OUTPUTS\n",
    "neuron_df, conn_df = fetch_adjacencies(neuron_ids, neuron_ids)\n",
    "\n",
    "conn_df = merge_neuron_properties(neuron_df, conn_df, ['type', 'instance']);\n",
    "display(conn_df.head())\n",
    "\n",
    "# Ordering data alphabetically\n",
    "conn_df.sort_values('type_pre', ascending=True, inplace=True)\n",
    "display(conn_df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c5931c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating the connectivity matrix\n",
    "\n",
    "from neuprint.utils import connection_table_to_matrix\n",
    "\n",
    "matrix = connection_table_to_matrix(conn_df, ('type_pre', 'type_post'),weight_col='weight')\n",
    "#print(matrix.shape)\n",
    "matrix.sort_index(axis=1, inplace=True)\n",
    "display(matrix.iloc[:10, :10])\n",
    "\n",
    "matrix.index = matrix.index.astype(str)\n",
    "matrix.columns = matrix.columns.astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3be6b2c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the connectibity matrix\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set(rc = {'figure.figsize':(15,8)})\n",
    "sns.heatmap(matrix, annot=False)\n",
    "plt.title(_column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc099184",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saving data in a csv file for later use\n",
    "import os\n",
    "\n",
    "# Changing some column names\n",
    "conn_df_save = conn_df.rename(columns = {'type_pre':'PreSynapticNeuron', 'type_post':'PostSynapticNeuron', 'weight':'N'}, inplace = False)\n",
    "\n",
    "data_path = 'YOUR\\PATH\\fib25'\n",
    "file_name = _column + '_fib25.csv'\n",
    "savePath = os.path.join(data_path,file_name)\n",
    "conn_df_save.to_csv(savePath)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "446019a7",
   "metadata": {},
   "source": [
    "<span style='color:Blue'>**EXCERCISE 4 FOR THE ON-PATHWAY:** \n",
    "<br>- Fetch the connectivity and plot the connectivity matrix between all Mi1s and T4 subtype a.\n",
    "<br>- Answer the following question: on average, how many different Mi1 contact a single T4a with more than 15 synapses?\n",
    "<br><br>TIP: fetch the connections for the different neurons (Mi1 and T4a) using the function fetch_neurons from neuprint.\n",
    "<br>TIP: to answer the question, use pandas dataframe grammar to filter the fetched connectivity as desired. </span> "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b30942ac",
   "metadata": {},
   "source": [
    "<span style='color:Blue'> **EXCERCISE 4 FOR THE OFF-PATHWAY:**\n",
    "<br>- Make Tm9 inputs become visible in a connectivity matrix. You can do this by changing either hetmap colors, diplaying the synaptic count, changing the scale of the heatmat (e.g. normalized values, log2, etc)\n",
    "<br>- Make the same plot for two more columns. Do the results differ?\n",
    "<br><br>TIP: explore options in the seaborn.heatmap function\n",
    "<br>TIP: you could modified the numbers in the pandas dataframe before generating the connectivity matrix. </span> \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "628526dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scaling the weitgh column\n",
    "# Normalizing to the max value:\n",
    "max_w = conn_df[\"weight\"].max() #\n",
    "norm_ls = conn_df[\"weight\"]/max_w # Normalization step\n",
    "conn_df[\"weight_norm\"] = norm_ls\n",
    "\n",
    "# Weigths to log2 \n",
    "import numpy as np\n",
    "log2_ls = np.log2(conn_df[\"weight\"])\n",
    "\n",
    "display(conn_df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9cf6d39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating the connectivity matrix\n",
    "\n",
    "from neuprint.utils import connection_table_to_matrix\n",
    "\n",
    "matrix = connection_table_to_matrix(conn_df, ('type_pre', 'type_post'),weight_col='weight_norm')\n",
    "#print(matrix.shape)\n",
    "matrix.sort_index(axis=1, inplace=True)\n",
    "matrix.index = matrix.index.astype(str)\n",
    "matrix.columns = matrix.columns.astype(str)\n",
    "\n",
    "# Plotting the connectibity matrix\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set(rc = {'figure.figsize':(15,8)})\n",
    "sns.heatmap(matrix, annot=False)\n",
    "plt.title(_column)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa9fc449",
   "metadata": {},
   "source": [
    "<span style='color:Blue'>**BONUS EXCERCISE 5 FOR THE OFF-PATHWAY:**\n",
    "<br>- Make an average connectivity matrix with data from all columns. Only include pair of connected neurons that are found across all seven columns. Discard the ones that only found in some. </span> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cf1ce81",
   "metadata": {},
   "source": [
    "### Connectivity graph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe85c09",
   "metadata": {},
   "source": [
    "Connections between pair of neurons locally inform us about which inputs or ouputs are important from the prespective of a single neuron. However, it gives little information about global interections that involved many more connections and neuronal types.\n",
    "\n",
    "For that, we can perform more complex analysis using graph theory."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75423f72",
   "metadata": {},
   "source": [
    "### Graph theory in a nutshell"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53c1f92c",
   "metadata": {},
   "source": [
    "A graph in graph theory is no more than a group of nodes that are connected to each other via edges. The nodes can represent many things in nature, from personal profiles in social networks to  modems in internet traffic. In our cases, it is useful to consider neuronal types as single nodes in a graph. The edges in a graph represent the strenght of connections between nodes via distance. This means that the closer two nodes are, the shorter the edge distance and the stronger their reationship.\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./images/graph_theory1.png\"\n",
    "     alt=\"Graph image\"\n",
    "     width=\"400\" height=\"200\" style=\"float: left; width: 50%\" />\n",
    "<img src=\"./images/graph_theory2.png\"\n",
    "     alt=\"Grapg function\"\n",
    "     width=\"400\" height=\"200\" style=\"float: right; width: 50%\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ac5cf5",
   "metadata": {},
   "source": [
    "### Creating a graph: wiring diagram"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b003f8d7",
   "metadata": {},
   "source": [
    "The current code to generate a graph based on connectomics and further analysis data needs to define some parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfc2ec99",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Users parameters\n",
    "user_parameters = {}\n",
    "user_parameters['dataPath']= r'D:\\fib25' # path for the folder with the stored data , r'YOUR\\PATH\\fib25'\n",
    "user_parameters['file']='home_column_fib25.csv' # data file\n",
    "user_parameters['column']='home column'# naming\n",
    "user_parameters['graph']= 'Fab25_data_7medulla_columns' # naming\n",
    "user_parameters['node_of_interest']= 'Tm9'  # for input-output plotting\n",
    "user_parameters['synaptic_stength_filter'] = 0 # to filter counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27138e39",
   "metadata": {},
   "source": [
    "Let's load the **previously saved file** that we used to generate the connectivity matrix and keep working with it. We use the package *pandas* to create a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7434e163",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os \n",
    "\n",
    "# Loading a data set\n",
    "filePath = os.path.join(user_parameters['dataPath'], user_parameters['file'])\n",
    "Table = pd.read_csv(filePath)\n",
    "Table = Table[Table['N'] >user_parameters['synaptic_stength_filter']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02d436c6",
   "metadata": {},
   "source": [
    "With the information on the table, containing all connection counts for each pair, we know can create a graph, which basically should contain a list of nodes and a list of edges connecting nodes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22ed10ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Creating list of nodes\n",
    "presynaptic_neurons = Table.PreSynapticNeuron.unique()\n",
    "postsynaptic_neurons = Table.PostSynapticNeuron.unique()\n",
    "neurons = list(np.unique(np.concatenate((presynaptic_neurons,postsynaptic_neurons), axis = 0)))\n",
    "for i,n in enumerate(neurons):\n",
    "    if n[-4:] == 'home':\n",
    "        neurons[i] = n[0:-4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d9d0328",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initialyzing and filling the graph with nodes and edges\n",
    "from network_analysis.core_functions_and_classes_network_analysis import Graph\n",
    "\n",
    "\n",
    "customGraph = Graph() # empty graph\n",
    "neuron_list = neurons\n",
    "\n",
    "# Adding nodes in the graph\n",
    "for n,neuron in enumerate(neuron_list):\n",
    "        customGraph.addNode(neuron)\n",
    "        customGraph.addNode_num(n)\n",
    "        \n",
    "# Adding edges in the graph        \n",
    "for index, row in Table.iterrows():\n",
    "    # accessing data using column names\n",
    "    if row['PreSynapticNeuron'][-4:] == 'home':\n",
    "        _pre = row['PreSynapticNeuron'][0:-4]\n",
    "        _post = row['PostSynapticNeuron'][0:-4]\n",
    "    else:\n",
    "        _pre = row['PreSynapticNeuron']\n",
    "        _post = row['PostSynapticNeuron']\n",
    "\n",
    "    pair = _pre ,_post \n",
    "        \n",
    "    if pair in customGraph.distances: # if the edge has been added, update the N.\n",
    "        temp_N = customGraph.distances[pair] + row['N']\n",
    "        customGraph.addEdge(pair[0],pair[1],temp_N)\n",
    "    else:\n",
    "        customGraph.addEdge(pair[0],pair[1],row['N']) # adding edge for the first time\n",
    "        # print('%f + %f = %f' % (customGraph.distances[pair],row['N'],temp_N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd14cb38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graph with networkx package \n",
    "\n",
    "# documentation: https://networkx.org/documentation/stable/reference/introduction.html\n",
    "# examples: https://algorithmx-python.readthedocs.io/en/latest/examples/networkx-examples.html\n",
    "\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Defing what synaptic weight means for us\n",
    "Weights = customGraph.distances # here, synaptic weight = absolut edge distance between two nodes\n",
    "for key, value in Weights.items():\n",
    "    Weights[key]= round(value) \n",
    "\n",
    "# We can change how is are the weigths (here same as edges) in different ways based on a scientific assumption\n",
    "edges_plot = [(k[0], k[1], {'weight': v}) for k, v in Weights.items()]# absolut edge distance between two nodes\n",
    "edges = [(k[0], k[1], {'weight': 1/v}) for k, v in Weights.items()] # 1/absolut edge distance between two nodes\n",
    "\n",
    "# Defining a graoh with the package networkx\n",
    "G_plot = nx.DiGraph()\n",
    "G = nx.DiGraph()\n",
    "\n",
    "# each edge here is already a tuple of the form (node1, node2, {'weight': weight})\n",
    "G_plot.add_edges_from(edges_plot)\n",
    "G.add_edges_from(edges)\n",
    "\n",
    "# Chossing a way to display the position of each node\n",
    "pos = nx.circular_layout(G_plot) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f1e64b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualizing the graph\n",
    "cm = 1/2.54  # centimeters in inches\n",
    "graph_fig,axes= plt.subplots(figsize=(20*cm, 20*cm))\n",
    "graph_fig.suptitle(user_parameters['file'])\n",
    "\n",
    "## nodes\n",
    "nx.draw_networkx_nodes(G_plot,pos,node_size=600)\n",
    "## labels\n",
    "nx.draw_networkx_labels(G_plot,pos,font_size=12,font_family='sans-serif')\n",
    "\n",
    "## edges\n",
    "#nx.draw_networkx_edges(G_plot,pos,edgelist=edges, width=2,connectionstyle=\"arc3,rad=-0.2\")\n",
    "nx.draw_networkx_edges(G_plot,pos,edgelist=edges, width=2)\n",
    "\n",
    "## weights\n",
    "#labels = nx.get_edge_attributes(G_plot,'weight')\n",
    "# for key, value in labels .items():\n",
    "#     labels [key]= round(1/value)\n",
    "# nx.draw_networkx_edge_labels(G_plot,pos,edge_labels=labels,label_pos = 0.3,font_size=10)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6591bee3",
   "metadata": {},
   "source": [
    "<span style='color:Blue'>**EXERCISE 6 - Circuit graph:**\n",
    "<br>- Try out different layouts from the networkx package \n",
    "<br>- Is there any to visualize more connected group of neurons? \n",
    "<br><br>TIP: consult <a href=\"https://networkx.org/documentation/stable/reference/drawing.html?highlight=layout#module-networkx.drawing.layout\" target=\"_blanck\">networkx</a> documentation. </span> "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5fa08c3",
   "metadata": {},
   "source": [
    "### The shortest path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc563680",
   "metadata": {},
   "source": [
    "The shortest-path algorithm calculates the shortest path from a start node to each node of a connected graph. Developed in 1956 by Edsger W. Dijsktra, it is the basis for all the apps that show you a shortest route from one place to another. In your, we can use it to show how information flows in the circtuitry and assume that the shortest path between two neurons could be the main way of information transfer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "181b98b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choosing start and end nodes\n",
    "# Nodes of interest to check all single paths downstream of a node and between two nodes\n",
    "user_parameters['start_node']= 'L3' \n",
    "user_parameters['last_node']= 'Tm9'\n",
    "\n",
    "# with NetwrokX\n",
    "length, path = nx.single_source_dijkstra(G, user_parameters['start_node'])\n",
    "\n",
    "# Paths downstream of a node\n",
    "path_str = ''.join('%s: %s ; ' % (k,v) for k,v in path.items())\n",
    "message_str = '\\n >>>>>>> All paths: \\n\\n%s can reach a neuron via [shortest path]: \\n\\n%s ' % (user_parameters['start_node'], path_str)\n",
    "print(message_str)\n",
    "\n",
    "# Paths between two nodes\n",
    "path_str = ''.join('%s: %s ; ' % (k,v) for k,v in path.items() if k == user_parameters['last_node'])\n",
    "message_str = '\\n >>>>>>> Single path: \\n\\n%s can reach %s via [shortest path]: \\n\\n%s ' % (user_parameters['start_node'],user_parameters['last_node'], path_str)\n",
    "print(message_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f68c3ab",
   "metadata": {},
   "source": [
    "<span style='color:Blue'>**EXERCISE 7 - Paths:**\n",
    "<br>- Print all the paths between L3 and Tm9.\n",
    "<br>- Which is the path with the highest number of synapses after normalized to the total number of edges involved?\n",
    "<br><br>TIP: use the function nx.all_simple_paths </span> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cca3362",
   "metadata": {},
   "source": [
    "### Centrality measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a69fda",
   "metadata": {},
   "source": [
    "In graph theory and network analysis, indicators of centrality assign numbers or rankings to nodes within a graph corresponding to their network position. Applications include identifying the most influential person(s) in a social network, key infrastructure nodes in the Internet or urban networks, super-spreaders of disease, and brain networks.\n",
    "\n",
    "Basic measurements of centrality include:\n",
    "1. **Degree centrality**: is defined as the number of links incident upon a node (i.e., the number of ties that a node has). The degree can be interpreted in terms of the immediate risk of a node for catching whatever is flowing through the network (such as a virus, or some information). In the case of a directed network (where ties have direction), we usually define two separate measures of degree centrality, namely indegree and outdegree.\n",
    "\n",
    "2. **Eigenvector centrality**: Also known as eigencentrality. It is a measure of the influence of a node in a network. It assigns relative scores to all nodes in the network based on the concept that connections to high-scoring nodes contribute more to the score of the node in question than equal connections to low-scoring nodes.\n",
    "\n",
    "3. **Closeness centrality**: the normalized closeness centrality (or closeness) of a node is the average length of the shortest path between the node and all other nodes in the graph. Thus the more central a node is, the closer it is to all other nodes. Taking distances from or to all other nodes is irrelevant in undirected graphs, whereas it can produce totally different results in directed graphs (e.g. a website can have a high closeness centrality from outgoing link, but low closeness centrality from incoming links).\n",
    "\n",
    "4. **Betweeness centrality**: Betweenness centrality quantifies the number of times a node acts as a bridge along the shortest path between two other nodes. It was introduced as a measure for quantifying the control of a human on the communication between other humans in a social network by Linton Freeman.[26] In his conception, vertices that have a high probability to occur on a randomly chosen shortest path between two randomly chosen vertices have a high betweenness.\n",
    "\n",
    "Node degree and eignencentrality assume that just the existance of a connection betweem nodees (bynary rule), is what matters, not considering the actual distance in the connections. For closeness and betweenss centrality, the distances matter, since they involve shortest path calculations.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61cd9dbe",
   "metadata": {},
   "source": [
    "<h4 style=\"float: left ; width: 25%\"> Degree centrality </h4> \n",
    "<h4 style=\"float: left ; width: 25%\"> Eignevector centrality </h4>\n",
    "<h4 style=\"float: left ; width: 25%\"> Closeness centrality </h4>\n",
    "<h4 style=\"float: left ; width: 25%\"> Betweenness centrality </h4>\n",
    "<img src=\"./images/degree_centrality.png\"\n",
    "     alt=\"Degree centrality\"\n",
    "     width=\"100\" height=\"100\" style=\"float: left ; width: 25%\"/>\n",
    "<img src=\"./images/eigencentrality.png\"\n",
    "     alt=\"Eignevector centrality\"\n",
    "     width=\"100\" height=\"100\" style=\"float: left ; width: 25%\"/>\n",
    "<img src=\"./images/closeness_centrality.png\"\n",
    "     alt=\"Closeness centrality\"\n",
    "     width=\"100\" height=\"100\" style=\"float: left ; width: 25%\"/>\n",
    "<img src=\"./images/betweenness_centrality.png\"\n",
    "     alt=\"Betweenness centrality\"\n",
    "     width=\"100\" height=\"100\" style=\"float: left ; width: 25%\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23193143",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Centrality measure\n",
    "from networkx.algorithms.centrality import closeness_centrality\n",
    "from networkx.algorithms.centrality import betweenness_centrality\n",
    "from networkx.algorithms.centrality import degree_centrality\n",
    "from networkx.algorithms.centrality import eigenvector_centrality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "475ea031",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating a common dataframe with all measures\n",
    "\n",
    "degree_dict = degree_centrality(G)\n",
    "eigenvector_dict = eigenvector_centrality(G)\n",
    "pagerank_dict = nx.pagerank(G, alpha = 0.8)\n",
    "betweenness_dict = betweenness_centrality(G, k=None, normalized=True, weight=\"weight\") # Use distances as weight?\n",
    "closeness_dict = closeness_centrality(G, u=None, distance=\"weight\", wf_improved=True) # Inputs distances?\n",
    "\n",
    "degree_df =pd.DataFrame(degree_dict, index=[0]) \n",
    "eigenvector_df =pd.DataFrame(eigenvector_dict, index=[0])\n",
    "pagerank_df =pd.DataFrame(pagerank_dict, index=[0])  \n",
    "betweenness_df =pd.DataFrame(betweenness_dict, index=[0]) \n",
    "closeness_df =pd.DataFrame(closeness_dict, index=[0]) \n",
    "\n",
    "degree_df = degree_df.T\n",
    "eigenvector_df = eigenvector_df.T\n",
    "pagerank_df = pagerank_df.T\n",
    "betweenness_df = betweenness_df.T\n",
    "closeness_df = closeness_df.T\n",
    "\n",
    "centrality_df = pd.concat([degree_df , eigenvector_df, pagerank_df,betweenness_df,closeness_df ], axis=1)\n",
    "centrality_df.columns = ['Degree', 'Eigenvector', 'Pagerank eigenvector', 'Betweenness', 'Closeness']\n",
    "centrality_df.reset_index(level=0, inplace=True)\n",
    "centrality_df.rename(columns = {'index':'Neuron'}, inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c7318d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plotting centrality measures\n",
    "import seaborn as sns\n",
    "\n",
    "# Barplots\n",
    "bar_fig_c, bar_axes_c = plt.subplots(nrows= 2,ncols=2,figsize=(20*cm, 20*cm)) # All together\n",
    "_fontsize = 6\n",
    "\n",
    "cur_df = centrality_df\n",
    "\n",
    "_color = 'tomato'\n",
    "sns.barplot(ax = bar_axes_c[0,0],x = 'Neuron', y = 'Degree', data = cur_df, order=cur_df.sort_values('Degree',ascending = False).Neuron,\n",
    "            color=_color)\n",
    "bar_axes_c[0,0].set_xticklabels(cur_df.sort_values('Degree',ascending = False).Neuron.tolist(), rotation = 90, fontsize = _fontsize)\n",
    "\n",
    "\n",
    "_color = 'purple'\n",
    "sns.barplot(ax = bar_axes_c[0,1],x = 'Neuron', y = 'Pagerank eigenvector', data = cur_df, order=cur_df.sort_values('Pagerank eigenvector',ascending = False).Neuron,\n",
    "            color=_color)\n",
    "bar_axes_c[0,1].set_xticklabels(cur_df.sort_values('Pagerank eigenvector',ascending = False).Neuron.tolist(), rotation = 90, fontsize = _fontsize)\n",
    "\n",
    "\n",
    "_color = 'teal'\n",
    "sns.barplot(ax = bar_axes_c[1,0],x = 'Neuron', y = 'Betweenness', data = cur_df, order=cur_df.sort_values('Betweenness',ascending = False).Neuron,\n",
    "            color=_color)\n",
    "bar_axes_c[1,0].set_xticklabels(cur_df.sort_values('Betweenness',ascending = False).Neuron.tolist(), rotation = 90, fontsize = _fontsize)\n",
    "\n",
    "\n",
    "_color = 'gold'\n",
    "sns.barplot(ax = bar_axes_c[1,1],x = 'Neuron', y = 'Closeness', data = cur_df, order=cur_df.sort_values('Closeness',ascending = False).Neuron,\n",
    "            color=_color)\n",
    "bar_axes_c[1,1].set_xticklabels(cur_df.sort_values('Closeness',ascending = False).Neuron.tolist(), rotation = 90, fontsize = _fontsize)\n",
    "\n",
    "\n",
    "bar_axes_c[0,0].spines['right'].set_visible(False)\n",
    "bar_axes_c[0,1].spines['right'].set_visible(False)\n",
    "bar_axes_c[1,0].spines['right'].set_visible(False)\n",
    "bar_axes_c[1,1].spines['right'].set_visible(False)\n",
    "bar_axes_c[0,0].spines['top'].set_visible(False)\n",
    "bar_axes_c[0,1].spines['top'].set_visible(False)\n",
    "bar_axes_c[1,0].spines['top'].set_visible(False)\n",
    "bar_axes_c[1,1].spines['top'].set_visible(False)\n",
    "\n",
    "\n",
    "_title = user_parameters['graph'] + ': ' + user_parameters['column'] + ' Centrality measures'\n",
    "bar_fig_c.suptitle(_title, fontsize = 12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6de7089",
   "metadata": {},
   "source": [
    "<span style='color:Blue'>**EXERCISE 8 - Centrality measures:**\n",
    "<br>- Plot the centrality measures only for medulla neurons, dicarting lamina inputs (L neurons) and output neurons (T4 cells)\n",
    "<br>- Plot centrality measures only for neurons belonging to the OFF-pathway.\n",
    "<br><br>TIP: use the DataFrame.Series.isin(values) function to include or discard neurons based on their name (here a Series) </span> "
   ]
  }
 ],
 "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.7.4"
  },
  "vscode": {
   "interpreter": {
    "hash": "6d46b44854ce120c1d45df5260ad03a52822abdd2850770c79221850bab6cd24"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}