ET_manuscript_Figure_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5849609",
   "metadata": {},
   "source": [
    "### FIGURE 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bb3baf58",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import nglui\n",
    "from statistics import mean\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import cloudvolume\n",
    "from caveclient import CAVEclient\n",
    "#import pcg_skel\n",
    "import tqdm\n",
    "#from meshparty import meshwork\n",
    "import datetime\n",
    "\n",
    "client = CAVEclient('minnie65_phase3_v1')\n",
    "client.info.get_datastack_info()\n",
    "\n",
    "\n",
    "pd.options.display.max_rows = None\n",
    "pd.options.display.max_columns = None\n",
    "\n",
    "#Load dataframe\n",
    "synapse_table = pd.read_feather(\"ET_extended_synapse_table_revision.feather\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aef8efe3-3475-4c40-bfc7-ad9be9876abe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duplicate entries:\n",
      "[183966021]\n",
      "original synapse table size\n",
      "(5530, 32)\n",
      " synapse table size without duplicates\n",
      "(5525, 32)\n"
     ]
    }
   ],
   "source": [
    "# check duplicate entries \n",
    "duplicates = synapse_table[synapse_table['synapse_id'].duplicated(keep=False)]\n",
    "\n",
    "print(\"Duplicate entries:\")\n",
    "print(duplicates.synapse_id.unique())\n",
    "\n",
    "# Remove duplicate rows \n",
    "\n",
    "print(\"original synapse table size\")\n",
    "print(synapse_table.shape)\n",
    "\n",
    "synapse_table = synapse_table.drop_duplicates(subset='synapse_id')\n",
    "\n",
    "print(\" synapse table size without duplicates\")\n",
    "print(synapse_table.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e364aa2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# MAIN\n",
    "#'y' location of cortical surface\n",
    "surface_y_column =[84534, 85689, 86053, 87800, 89421, 90105, 82884, 81677, 86242]\n",
    "average_surface_location = mean(surface_y_column)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e8568a35",
   "metadata": {},
   "outputs": [],
   "source": [
    "#root IDs of ET neurons\n",
    "V1_neurons_rootID = synapse_table[synapse_table['pre_soma_area'] == 'v1']['pre_pt_root_id'].unique()\n",
    "HVA_neurons_rootID = synapse_table[synapse_table['pre_soma_area'] == 'hva']['pre_pt_root_id'].unique()\n",
    "All_neurons_rootID = list(V1_neurons_rootID) + list(HVA_neurons_rootID)\n",
    "\n",
    "#nucleus IDs of ET neurons\n",
    "V1_neurons_nucleusID = synapse_table[synapse_table['pre_soma_area'] == 'v1']['pre_nucleus_id'].unique()\n",
    "HVA_neurons_nucleusID = synapse_table[synapse_table['pre_soma_area'] == 'hva']['pre_nucleus_id'].unique()\n",
    "All_neurons_nucleusID = list(V1_neurons_nucleusID) + list(HVA_neurons_nucleusID)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "71b26dfd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                    ID  total_inputs\n",
      "0   864691136974876956         14717\n",
      "1   864691135273655569         11823\n",
      "2   864691136578820884         12659\n",
      "3   864691135397023777          9122\n",
      "4   864691136578814228         13407\n",
      "5   864691135517376211         11047\n",
      "6   864691136991202453         15326\n",
      "7   864691135489514810         26797\n",
      "8   864691136579127060         11562\n",
      "9   864691135503182685         11128\n",
      "10  864691135339044582         10355\n",
      "11  864691135464714565         10147\n",
      "mean:  13174\n",
      "stdev:  4468.219067916085\n"
     ]
    }
   ],
   "source": [
    "#Statistics of ET Neurons\n",
    "\n",
    "\n",
    "# synaptic input to ET neurons\n",
    "\n",
    "import statistics\n",
    "import numpy\n",
    "\n",
    "et_table = 'bodor_pt_cells'\n",
    "et_df = client.materialize.query_table(et_table)\n",
    "\n",
    "pt_rootIDs = list(et_df.pt_root_id.unique())\n",
    "\n",
    "\n",
    "ET_total_inputs = []\n",
    "\n",
    "for ii,pt_root_id in enumerate(pt_rootIDs): \n",
    "\n",
    "    ET_individuals = {\n",
    "        'ID': pt_root_id,\n",
    "        'total_inputs': len(client.materialize.synapse_query(post_ids=pt_root_id))\n",
    "    }\n",
    "    ET_total_inputs.append(ET_individuals)\n",
    "\n",
    "ET_total_inputs_df = pd.DataFrame(ET_total_inputs)\n",
    "\n",
    "print(ET_total_inputs_df)\n",
    "print('mean: ', mean(ET_total_inputs_df['total_inputs'].values))\n",
    "print('stdev: ', numpy.std(ET_total_inputs_df['total_inputs'].values))\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1e6fa4b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "0it [00:00, ?it/s]\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'meshwork' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[6], line 21\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m _, row \u001b[38;5;129;01min\u001b[39;00m tqdm\u001b[38;5;241m.\u001b[39mtqdm(et_df\u001b[38;5;241m.\u001b[39miterrows()):\n\u001b[1;32m     17\u001b[0m \n\u001b[1;32m     18\u001b[0m \n\u001b[1;32m     19\u001b[0m     \u001b[38;5;66;03m#print(row)\u001b[39;00m\n\u001b[1;32m     20\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mexists(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mskeldir\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m/\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mrow[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpt_root_id\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m.h5\u001b[39m\u001b[38;5;124m\"\u001b[39m):\n\u001b[0;32m---> 21\u001b[0m         nrns[row[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpt_root_id\u001b[39m\u001b[38;5;124m\"\u001b[39m]] \u001b[38;5;241m=\u001b[39m \u001b[43mmeshwork\u001b[49m\u001b[38;5;241m.\u001b[39mload_meshwork(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mskeldir\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m/\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mrow[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpt_root_id\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m.h5\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     23\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m     25\u001b[0m         nrns[row[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpt_root_id\u001b[39m\u001b[38;5;124m\"\u001b[39m]] \u001b[38;5;241m=\u001b[39m pcg_skel\u001b[38;5;241m.\u001b[39mcoord_space_meshwork(\n\u001b[1;32m     26\u001b[0m             row[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpt_root_id\u001b[39m\u001b[38;5;124m\"\u001b[39m],\n\u001b[1;32m     27\u001b[0m             client\u001b[38;5;241m=\u001b[39mclient,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     33\u001b[0m             timestamp \u001b[38;5;241m=\u001b[39m now,\n\u001b[1;32m     34\u001b[0m         )\n",
      "\u001b[0;31mNameError\u001b[0m: name 'meshwork' is not defined"
     ]
    }
   ],
   "source": [
    "#Statistics of ET Neurons\n",
    "\n",
    "# Axon length of ET neurons\n",
    "\n",
    "skeldir = 'skeletons'\n",
    "now = client.materialize.get_timestamp()\n",
    "\n",
    "et_table = 'bodor_pt_cells'\n",
    "et_df = client.materialize.query_table(et_table)\n",
    "\n",
    "ET_total_axon = []\n",
    "\n",
    "# Build the skeletons from dataframe\n",
    "nrns = {}\n",
    "\n",
    "for _, row in tqdm.tqdm(et_df.iterrows()):\n",
    "\n",
    "\n",
    "    #print(row)\n",
    "    if os.path.exists(f\"{skeldir}/{row['pt_root_id']}.h5\"):\n",
    "        nrns[row[\"pt_root_id\"]] = meshwork.load_meshwork(f\"{skeldir}/{row['pt_root_id']}.h5\")\n",
    "   \n",
    "    else:\n",
    "    \n",
    "        nrns[row[\"pt_root_id\"]] = pcg_skel.coord_space_meshwork(\n",
    "            row[\"pt_root_id\"],\n",
    "            client=client,\n",
    "            root_point=row[\"pt_position\"],\n",
    "            root_point_resolution=[4, 4, 40],\n",
    "            collapse_soma=True,\n",
    "            synapses=\"all\",\n",
    "            synapse_table=client.info.get_datastack_info().get(\"synapse_table\"),\n",
    "            timestamp = now,\n",
    "        )\n",
    "\n",
    "        nrns[row[\"pt_root_id\"]].save_meshwork(f\"{skeldir}/{row['pt_root_id']}.h5\")           \n",
    "# Get the axons and calculate length\n",
    "\n",
    "\n",
    "for rid, nrn in nrns.items():\n",
    "\n",
    "    print(f\"Processing neuron with ID: {rid}\")\n",
    "    \n",
    "    is_axon = meshwork.algorithms.split_axon_by_annotation(\n",
    "        nrn,\n",
    "        'pre_syn',\n",
    "        'post_syn',\n",
    "        return_quality=False\n",
    "    )\n",
    "    nrn.anno.add_annotations('is_axon', is_axon, mask=True)\n",
    "    \n",
    "    with nrn.mask_context( nrn.anno.is_axon.mesh_mask ):\n",
    "        print(nrn.path_length() / 1_000)\n",
    "        ET_individuals = {\n",
    "        'ID': {rid},\n",
    "        'total_axon': nrn.path_length() / 1_000\n",
    "        }\n",
    "    \n",
    "    ET_total_axon.append(ET_individuals)\n",
    "\n",
    "ET_total_axon_df = pd.DataFrame(ET_total_axon)\n",
    "\n",
    "#print(ET_total_axon_df)\n",
    "#print('mean: ', mean(ET_total_axon_df['total_axon'].values))\n",
    "#print('stdev: ', numpy.std(ET_total_axon_df['total_axon'].values))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2700c255",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ANALYSIS - GENERAL\n",
    "\n",
    "#NUMBER OF SYNAPSES\n",
    "\n",
    "synapse_numbers = { 'total_synapses' : len(synapse_table),\n",
    "                   \n",
    "                    'Synapses_with_cells_with_soma_in_the_volume' : len(synapse_table[(synapse_table['num_soma'] == 1)]),\n",
    "                    \n",
    "                    'synapses_with_class_labels' : len(synapse_table[pd.notna(synapse_table['consensus_class']) &\n",
    "                                                                            (synapse_table['manual_class'] != 'unclear')]),\n",
    "\n",
    "                    'synapses_with_class_labels_manual' : len(synapse_table[(pd.notna(synapse_table['manual_class'])) &\n",
    "                                                                            (synapse_table['manual_class'] != 'unclear')]),\n",
    "                   \n",
    "                    'synapses_without_class_labels' : len(synapse_table[(pd.isna(synapse_table['consensus_class']))]),\n",
    "                   \n",
    "                    'synapses_with_subclass_labels' : len(synapse_table[pd.notna(synapse_table['consensus_subclass'])]),\n",
    "                    \n",
    "                    'synapses with_subclass_labels_manual' : len(synapse_table[(pd.notna(synapse_table['manual_subclass'])) &\n",
    "                                                                            (synapse_table['manual_subclass'] != 'excitatory') & \n",
    "                                                                            (synapse_table['manual_subclass'] != 'inhibitory')]),\n",
    "                    \n",
    "                    'synapses_without_subclass_labels' : len(synapse_table[pd.isna(synapse_table['consensus_subclass'])]),\n",
    "                   \n",
    "                    'synapses_with_orphans' : len(synapse_table[synapse_table['num_soma'] == 0]),\n",
    "                                                 \n",
    "                    'synapses_with_orphans_with_class_labels' : len(synapse_table[(synapse_table['num_soma'] == 0) &\n",
    "                                                                    pd.notna(synapse_table['consensus_class'])]),\n",
    "                   \n",
    "                    'synapses_with_orphans_without_class_labels' : len(synapse_table[(synapse_table['num_soma'] == 0) &\n",
    "                                                                    pd.isna(synapse_table['consensus_class'])]),\n",
    "                   \n",
    "                    'synapses_with_orphan_spines': len(synapse_table[synapse_table['orphan'] == 'single_spine']),\n",
    "                   \n",
    "                    'synapses_with_orphan_dendrites_without_class_labels': len(synapse_table[(synapse_table['num_soma'] == 0) &\n",
    "                                                                    (pd.isna(synapse_table['consensus_class'])) &\n",
    "                                                                    (synapse_table['orphan'] == 'dendrite')]),\n",
    "                   \n",
    "                    'synapses_with_orphans_with_subclass_labels' : len(synapse_table[(synapse_table['num_soma'] == 0) &\n",
    "                                                                    pd.notna(synapse_table['consensus_subclass'])]),\n",
    "                    'synapses_with_orphans_without_subclass_labels' : len(synapse_table[(synapse_table['num_soma'] == 0) &\n",
    "                                                                    pd.isna(synapse_table['consensus_subclass'])]),\n",
    "\n",
    "                                   \n",
    "                    'synapses_with_multisoma' : len(synapse_table[synapse_table['num_soma'] > 1]),\n",
    "                    'synapses_with_multisoma_with_class_labels' : len(synapse_table[(synapse_table['num_soma'] > 1) &\n",
    "                                                                    pd.notna(synapse_table['consensus_class'])]),\n",
    "                    'synapses_with_multisoma_without_class_labels' : len(synapse_table[(synapse_table['num_soma'] > 1) &\n",
    "                                                                    pd.isna(synapse_table['consensus_class'])]),\n",
    "}\n",
    "\n",
    "\n",
    "#NUMBER OF SOMATA\n",
    "\n",
    "connection_numbers = {\n",
    "                    'total_connections' : len(synapse_table[synapse_table['num_soma'] == 1].post_pt_root_id.unique()),\n",
    "                   \n",
    "                    'connections_with_cells_with_soma' : len(synapse_table[(synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "                    \n",
    "                    'connections_with_class_labels' : len(synapse_table[pd.notna(synapse_table['consensus_class']) &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "    \n",
    "                    'connections_with_class_labels_manual' : len(synapse_table[(pd.notna(synapse_table['manual_class'])) &\n",
    "                                                             (synapse_table['manual_class'] != 'unclear') &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "                    \n",
    "                    'connections_without_class_labels' : len(synapse_table[(pd.isna(synapse_table['consensus_class'])) &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "                \n",
    "                    'with_subclass_labels' : len(synapse_table[(pd.notna(synapse_table['consensus_subclass'])) &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "                    'with_subclass_labels_manual' : len(synapse_table[(pd.notna(synapse_table['manual_subclass'])) &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "                    'without_subclass_labels' : len(synapse_table[(pd.isna(synapse_table['consensus_subclass'])) &\n",
    "                                                             (synapse_table['num_soma'] == 1)].post_pt_root_id.unique()),\n",
    "   \n",
    "}\n",
    "\n",
    "connection_numbers, synapse_numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d13e97c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ANALYSIS - FIGURE 3.e and f\n",
    "#MAKE STATS DATAFRAME - All synapses and connections\n",
    "\n",
    "np.seterr(divide='ignore', invalid='ignore')\n",
    "\n",
    "percentage = []\n",
    "values = []\n",
    "\n",
    "#IDs of presynaptic PT neurons\n",
    "pre_soma_IDs = synapse_table['pre_nucleus_id'].unique()\n",
    "\n",
    "\n",
    "for ii,pre_soma_ID in enumerate(pre_soma_IDs):        \n",
    "        \n",
    "   \n",
    "    stat_values={\n",
    "\n",
    "                    'ID': pre_soma_ID,\n",
    "                    \n",
    "        #SYNAPSES\n",
    "                    \n",
    "                    'all_syn_total': len(synapse_table[ (synapse_table['pre_nucleus_id'] == pre_soma_ID)]),\n",
    "                    'all_e_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'excitatory')]),\n",
    "                    'all_i_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'inhibitory')]),        \n",
    "                    'all_Undetermined_class_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_class'] == 'unclear') |\n",
    "                                      pd.isnull(synapse_table['consensus_class']))]),\n",
    "                    \n",
    "        \n",
    "                    'all_23P_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '23P')]),\n",
    "        \n",
    "                    'all_4P_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '4P')]),\n",
    "                    'all_5P-ET_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-ET')]),\n",
    "                    'all_5P-NP_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-NP')]), \n",
    "                    'all_5P-IT_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-IT')]), \n",
    "                \n",
    "                    'all_6P_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P')]),\n",
    "        \n",
    "                    'all_6P-CT_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P-CT')]),\n",
    "        \n",
    "                    'all_6P-IT_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P-IT')]),\n",
    "        \n",
    "                    'all_BC_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'BC')]),\n",
    "        \n",
    "                    'all_MC_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'MC')]),\n",
    "        \n",
    "                    'all_BPC_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'BPC')]),\n",
    "        \n",
    "                    'all_NGC_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'NGC')]),\n",
    "        \n",
    "                    'all_Undetermined_subclass_syn#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_subclass'] == 'unclear') |\n",
    "                                      pd.isnull(synapse_table['consensus_subclass']))]),\n",
    "        \n",
    "        \n",
    "        #CONNECTIONS\n",
    "                    'all_con_total': len(synapse_table[ (synapse_table['pre_nucleus_id'] == pre_soma_ID)]\n",
    "                                         ['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_e_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'excitatory')]\n",
    "                                      ['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_i_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'inhibitory')]\n",
    "                                      ['post_pt_root_id'].unique()),\n",
    "        \n",
    "                    'all_Undetermined_class_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_class'] == 'unclear') |\n",
    "                                      pd.isnull(synapse_table['consensus_class']))]['post_pt_root_id'].unique()),                      \n",
    "        \n",
    "                    'all_23P_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '23P')]['post_pt_root_id'].unique()),\n",
    "        \n",
    "                    'all_4P_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '4P')]['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_5P-ET_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-ET')]['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_5P-NP_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-NP')]['post_pt_root_id'].unique()), \n",
    " \n",
    "                    'all_5P-IT_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '5P-IT')]['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_6P_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P')]['post_pt_root_id'].unique()),\n",
    "                \n",
    "                    'all_6P-CT_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P-CT')]['post_pt_root_id'].unique()),\n",
    "        \n",
    "                    'all_6P-IT_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == '6P-IT')]['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_BC_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'BC')]['post_pt_root_id'].unique()),\n",
    "                     \n",
    "                    'all_MC_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'MC')]['post_pt_root_id'].unique()),\n",
    "       \n",
    "                    'all_BPC_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'BPC')]['post_pt_root_id'].unique()),\n",
    "                    \n",
    "                    'all_NGC_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_subclass'] == 'NGC')]['post_pt_root_id'].unique()),\n",
    "        \n",
    "                    'all_Undetermined_subclass_con#': len(synapse_table[(synapse_table['pre_nucleus_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_subclass'] == 'unclear') |\n",
    "                                      pd.isnull(synapse_table['consensus_subclass']))]['post_pt_root_id'].unique()),\n",
    "                        \n",
    "        \n",
    "    }\n",
    "    values.append(stat_values)\n",
    "    \n",
    "    \n",
    "    stat_percentage={\n",
    "\n",
    "                    'ID': pre_soma_ID,\n",
    "                    'Soma_depth':  (client.materialize.query_table('nucleus_detection_v0', \n",
    "                               filter_equal_dict={'id':pre_soma_ID})['pt_position'].apply(lambda x: x[1]).iloc[0]\n",
    "                                                  - average_surface_location)*0.004,\n",
    "                   \n",
    "                      \n",
    "                    'i_syn': stat_values['all_i_syn#'] / stat_values['all_syn_total'],\n",
    "                    'e_syn': stat_values['all_e_syn#'] / stat_values['all_syn_total'],\n",
    "                    'Undetermined_class_syn': stat_values['all_Undetermined_class_syn#'] / stat_values['all_syn_total'],\n",
    "\n",
    "        \n",
    "                    'i_syn_determined': stat_values['all_i_syn#'] / (stat_values['all_i_syn#']+stat_values['all_e_syn#']),\n",
    "                    'e_syn_determined': stat_values['all_e_syn#'] / (stat_values['all_i_syn#']+stat_values['all_e_syn#']),\n",
    "        \n",
    "        \n",
    "                    '23P_syn': stat_values['all_23P_syn#'] / stat_values['all_syn_total'],\n",
    "                    '4P_syn': stat_values['all_4P_syn#'] / stat_values['all_syn_total'],\n",
    "                    '5P-ET_syn': stat_values['all_5P-ET_syn#'] / stat_values['all_syn_total'],\n",
    "                    '5P-IT_syn': stat_values['all_5P-IT_syn#'] / stat_values['all_syn_total'],\n",
    "                    '5P-NP_syn': stat_values['all_5P-NP_syn#'] / stat_values['all_syn_total'],\n",
    "                    '6P_syn': stat_values['all_6P_syn#'] / stat_values['all_syn_total'],\n",
    "                    '6P-IT_syn': stat_values['all_6P-IT_syn#'] / stat_values['all_syn_total'],\n",
    "                    '6P-CT_syn': stat_values['all_6P-CT_syn#'] / stat_values['all_syn_total'],\n",
    "                    'BC_syn': stat_values['all_BC_syn#'] / stat_values['all_syn_total'],\n",
    "                    'MC_syn': stat_values['all_MC_syn#'] / stat_values['all_syn_total'],\n",
    "                    'BPC_syn': stat_values['all_BPC_syn#'] / stat_values['all_syn_total'],\n",
    "                    'NGC_syn': stat_values['all_NGC_syn#'] / stat_values['all_syn_total'],\n",
    "                    'Undetermined_subclass_syn': stat_values['all_Undetermined_subclass_syn#'] / stat_values['all_syn_total'],\n",
    "        \n",
    "        \n",
    "                    'i_con': stat_values['all_i_con#'] / stat_values['all_con_total'],\n",
    "                    'e_con': stat_values['all_e_con#'] / stat_values['all_con_total'],\n",
    "                    'Undetermined_class_con': stat_values['all_Undetermined_class_con#'] / stat_values['all_con_total'],\n",
    "                   \n",
    "        \n",
    "                    '23P_con': stat_values['all_23P_con#'] / stat_values['all_con_total'],\n",
    "                    '4P_con': stat_values['all_4P_con#'] / stat_values['all_con_total'],\n",
    "                    '5P-ET_con': stat_values['all_5P-ET_con#'] / stat_values['all_con_total'],\n",
    "                    '5P-IT_con': stat_values['all_5P-IT_con#'] / stat_values['all_con_total'],\n",
    "                    '5P-NP_con': stat_values['all_5P-NP_con#'] / stat_values['all_con_total'],\n",
    "                    '6P_con': stat_values['all_6P_con#'] / stat_values['all_con_total'],\n",
    "                    '6P-IT_con': stat_values['all_6P-IT_con#'] / stat_values['all_con_total'],\n",
    "                    '6P-CT_con': stat_values['all_6P-CT_con#'] / stat_values['all_con_total'],\n",
    "                    'BC_con': stat_values['all_BC_con#'] / stat_values['all_con_total'],\n",
    "                    'MC_con': stat_values['all_MC_con#'] / stat_values['all_con_total'],\n",
    "                    'BPC_con': stat_values['all_BPC_con#'] / stat_values['all_syn_total'],\n",
    "                    'NGC_con': stat_values['all_NGC_con#'] / stat_values['all_con_total'],\n",
    "                    'Undetermined_subclass_con': stat_values['all_Undetermined_subclass_con#'] / stat_values['all_con_total'],\n",
    "                   \n",
    "    }\n",
    "    percentage.append(stat_percentage) \n",
    "    \n",
    "\n",
    "synapse_table_values = pd.DataFrame(values)\n",
    "synapse_table_percentage = pd.DataFrame(percentage)\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e55b2103",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "synapse_table_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb47cac9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Number of synapses in a connection\n",
    "\n",
    "syn_per_con_e = synapse_table.query('consensus_class ==\"excitatory\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_i = synapse_table.query('consensus_class ==\"inhibitory\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_23P = synapse_table.query('consensus_subclass ==\"23P\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_4P = synapse_table.query('consensus_subclass ==\"4P\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_5P_IT = synapse_table.query('consensus_subclass ==\"5P-IT\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_5P_ET = synapse_table.query('consensus_subclass ==\"5P-ET\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_5P_NP = synapse_table.query('consensus_subclass ==\"5P-NP\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_6P = synapse_table.query('consensus_subclass ==\"6P\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_BC = synapse_table.query('consensus_subclass ==\"BC\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "syn_per_con_MC = synapse_table.query('consensus_subclass ==\"MC\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].mean()\n",
    "\n",
    "\n",
    "print(\"Number of synapses onto excitatory cells with soma in the volume:\", all_syn_e_soma )\n",
    "print(\"Number of target excitatory cells with soma in the volume:\", all_con_e_soma )\n",
    "\n",
    "print(\"Number of synapses onto inhibitory  cells with soma in the volume:\", all_syn_i_soma )\n",
    "print(\"Number of target inhibitory cells with soma in the volume:\", all_con_i_soma )\n",
    "\n",
    "print(\"Number of synapses per connection with excitatry cell with soma in the volume:\", syn_per_con_e)\n",
    "print(\"Number of synapses per connection with inhibitory cell with soma in the volume:\", syn_per_con_i)\n",
    "\n",
    "print(\"Number of synapses per connection with 23P cell with soma in the volume:\", \n",
    "      syn_per_con_23P)\n",
    "\n",
    "print(\"Number of synapses per connection with 4P cell with soma in the volume:\", \n",
    "      syn_per_con_4P)\n",
    "\n",
    "print(\"Number of synapses per connection with 5P-IT cell with soma in the volume:\", \n",
    "      syn_per_con_5P_IT)\n",
    "\n",
    "print(\"Number of synapses per connection with 5P-ET cell with soma in the volume:\", \n",
    "      syn_per_con_5P_ET)\n",
    "\n",
    "print(\"Number of synapses per connection with 5P-NP cell with soma in the volume:\", \n",
    "      syn_per_con_5P_NP)\n",
    "\n",
    "print(\"Number of synapses per connection with 6P cell with soma in the volume:\", \n",
    "      syn_per_con_6P)\n",
    "\n",
    "print(\"Number of synapses per connection with BC cell with soma in the volume:\", \n",
    "      syn_per_con_BC)\n",
    "\n",
    "print(\"Number of synapses per connection with MC cell with soma in the volume:\", \n",
    "      syn_per_con_MC)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec41dfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "synapse_table.query('consensus_class ==\"inhibitory\" and num_soma==1').groupby(\n",
    "['pre_pt_root_id','post_pt_root_id']).count()['synapse_id'].std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42c83a7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#FIGURE 3.e\n",
    "#Swarm plot\n",
    "import seaborn as sns\n",
    "sns.set(style=\"whitegrid\")\n",
    "#tips = sns.load_dataset(\"tips\")\n",
    "#ax = sns.swarmplot(x=tips[\"total_bill\"])\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "\n",
    "ax = sns.swarmplot(y=\"i_syn_determined\",\n",
    "                   data=synapse_table_percentage,\n",
    "                  size=15)\n",
    "\n",
    "\n",
    "ax.set(ylim=(0, 1), yticks=np.arange(0, 1.1 , 0.1))\n",
    "\n",
    "\n",
    "\n",
    "plt.savefig('Fig3e_synapses_revision.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38a04fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#FIGURE 3.f\n",
    "#Swarm plot\n",
    "import seaborn as sns\n",
    "sns.set(style=\"whitegrid\")\n",
    "#tips = sns.load_dataset(\"tips\")\n",
    "#ax = sns.swarmplot(x=tips[\"total_bill\"])\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "\n",
    "ax = sns.swarmplot(y=\"i_con\",\n",
    "                   data=synapse_table_percentage,\n",
    "                  size=15)\n",
    "\n",
    "ax.set(ylim=(0, 1))\n",
    "ax.set(ylim=(0, 1), yticks=np.arange(0, 1.1 , 0.1))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "plt.savefig('Fig3f_synapses_revision.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f928a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ANALYSIS - FIGURE 3.g-h\n",
    "\n",
    "np.seterr(divide='ignore', invalid='ignore')\n",
    "\n",
    "percentage = []\n",
    "values = []\n",
    "\n",
    "#IDs of presynaptic PT neurons\n",
    "pre_soma_IDs = All_neurons_rootID\n",
    "\n",
    "\n",
    "for ii,pre_soma_ID in enumerate(pre_soma_IDs):        \n",
    "        \n",
    "   \n",
    "    stat_values={\n",
    "\n",
    "                    'ID': pre_soma_ID,\n",
    "                    'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                     pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "                  \n",
    "                    #SYNAPSES\n",
    "        \n",
    "                    #V1\n",
    "                    \n",
    "                    'all_total_syn#': len(synapse_table[ (synapse_table['pre_pt_root_id'] == pre_soma_ID)]),\n",
    "        \n",
    "                    'NO_AREA_total_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                                            pd.isnull(synapse_table['post_soma_area'])]),\n",
    "                     \n",
    "                    'V1_total_syn#': len(synapse_table[ (synapse_table['post_soma_area'] == 'v1') &\n",
    "                                                    (synapse_table['pre_pt_root_id'] == pre_soma_ID)]),                \n",
    "                    \n",
    "                    'V1_i_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'inhibitory') &\n",
    "                                                   (synapse_table['post_soma_area'] == 'v1')]),\n",
    "                   \n",
    "        \n",
    "                    'V1_e_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'excitatory') &\n",
    "                                                   (synapse_table['post_soma_area'] == 'v1')]),\n",
    "        \n",
    "                    \n",
    "                    'V1_Undetermined_class_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_class'] == 'Unsure') |\n",
    "                                      pd.isnull(synapse_table['consensus_class'])) &\n",
    "                                                                    (synapse_table['post_soma_area'] == 'v1')]),\n",
    "        \n",
    "\n",
    "        \n",
    "                    #HVA\n",
    "                    'HVA_total_syn#': len(synapse_table[ (synapse_table['post_soma_area'] == 'hva') &\n",
    "                                                    (synapse_table['pre_pt_root_id'] == pre_soma_ID)]),\n",
    "                    \n",
    "                    'HVA_i_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'inhibitory') &\n",
    "                                                   (synapse_table['post_soma_area'] == 'hva')]),\n",
    "        \n",
    "                    'HVA_e_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     (synapse_table['consensus_class'] == 'excitatory') &\n",
    "                                                   (synapse_table['post_soma_area'] == 'hva')]),\n",
    "        \n",
    "                    'HVA_Undetermined_class_syn#': len(synapse_table[(synapse_table['pre_pt_root_id'] == pre_soma_ID) &\n",
    "                                     ((synapse_table['consensus_class'] == 'Unsure') |\n",
    "                                      pd.isnull(synapse_table['consensus_class'])) &\n",
    "                                                                    (synapse_table['post_soma_area'] == 'hva')]),\n",
    "                    \n",
    "                   \n",
    "        \n",
    "        \n",
    "    }\n",
    "    values.append(stat_values)\n",
    "    \n",
    "    #print(stat_values['HVA_total_syn#'])\n",
    "    if stat_values['HVA_total_syn#'] > 10:\n",
    "        \n",
    "        local_area = synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0]\n",
    "        if local_area == ('v1'):\n",
    "            #print('aqui')\n",
    "\n",
    "            stat_percentage={\n",
    "\n",
    "                        'ID': pre_soma_ID,\n",
    "                        'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "\n",
    "                        #SYNAPSES\n",
    "\n",
    "                        #V1\n",
    "\n",
    "                        'location': 'local',\n",
    "                        'i_syn': stat_values['V1_i_syn#'] / stat_values['V1_total_syn#'],\n",
    "                        'e_syn': stat_values['V1_e_syn#'] / stat_values['V1_total_syn#'],\n",
    "                        'Undetermined_class_syn': stat_values['V1_Undetermined_class_syn#'] / stat_values['V1_total_syn#'],\n",
    "            }\n",
    "            percentage.append(stat_percentage)\n",
    "\n",
    "            if stat_values['HVA_total_syn#'] > 0:\n",
    "                stat_percentage={\n",
    "\n",
    "                        'ID': pre_soma_ID,\n",
    "                        'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "\n",
    "                        #SYNAPSES\n",
    "\n",
    "                        #HVA\n",
    "                        'location': 'inter-areal',\n",
    "                        'i_syn': stat_values['HVA_i_syn#'] / stat_values['HVA_total_syn#'],\n",
    "                        'e_syn': stat_values['HVA_e_syn#'] / stat_values['HVA_total_syn#'],\n",
    "                        'Undetermined_class_syn': stat_values['HVA_Undetermined_class_syn#'] / stat_values['HVA_total_syn#'],\n",
    "                }\n",
    "                percentage.append(stat_percentage)\n",
    "            else:\n",
    "                stat_percentage={\n",
    "\n",
    "                        'ID': pre_soma_ID,\n",
    "                        'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "\n",
    "                        #SYNAPSES\n",
    "\n",
    "                        #HVA\n",
    "                        'location': 'inter-areal',\n",
    "                        'i_syn': stat_values['HVA_i_syn#'],\n",
    "                        'e_syn': stat_values['HVA_e_syn#'],\n",
    "                        'Undetermined_class_syn': stat_values['HVA_Undetermined_class_syn#'],\n",
    "                }\n",
    "                percentage.append(stat_percentage)\n",
    "\n",
    "\n",
    "        else:       \n",
    "            stat_percentage={\n",
    "\n",
    "                        'ID': pre_soma_ID,\n",
    "                        'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "\n",
    "                        #SYNAPSES\n",
    "\n",
    "                        #V1\n",
    "\n",
    "                        'location': 'inter-areal',\n",
    "                        'i_syn': stat_values['V1_i_syn#'] / stat_values['V1_total_syn#'],\n",
    "                        'e_syn': stat_values['V1_e_syn#'] / stat_values['V1_total_syn#'],\n",
    "                        'Undetermined_class_syn': stat_values['V1_Undetermined_class_syn#'] / stat_values['V1_total_syn#'],\n",
    "            }\n",
    "            percentage.append(stat_percentage)\n",
    "\n",
    "            stat_percentage={\n",
    "\n",
    "                        'ID': pre_soma_ID,\n",
    "                        'local_area': synapse_table[synapse_table['pre_pt_root_id'] == \n",
    "                                                         pre_soma_ID]['pre_soma_area'].unique()[0],\n",
    "\n",
    "                        #SYNAPSES\n",
    "\n",
    "                        #HVA\n",
    "                        'location': 'local',\n",
    "                        'i_syn': stat_values['HVA_i_syn#'] / stat_values['HVA_total_syn#'],\n",
    "                        'e_syn': stat_values['HVA_e_syn#'] / stat_values['HVA_total_syn#'],\n",
    "                        'Undetermined_class_syn': stat_values['HVA_Undetermined_class_syn#'] / stat_values['HVA_total_syn#'],\n",
    "            }\n",
    "            percentage.append(stat_percentage)\n",
    "\n",
    "synapse_area_table_values = pd.DataFrame(values)\n",
    "synapse_area_table_percentage = pd.DataFrame(percentage)\n",
    "\n",
    "\n",
    "#stats_type.append(stat)   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca7bad35",
   "metadata": {},
   "outputs": [],
   "source": [
    "synapse_area_table_percentage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c6138d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#PLOT - Figure 3.g and h\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "#INDIVIDUAL CELLS\n",
    "\n",
    "# Create a pointplot using Seaborn\n",
    "sns.pointplot(x=\"location\", y=\"i_syn\", hue=\"ID\", data=synapse_area_table_percentage)\n",
    "\n",
    "# Set the plot title\n",
    "plt.title(\"Individual neurons\")\n",
    "\n",
    "#save figure\n",
    "plt.savefig('Fig3g_interareal_neurons_revision.eps')\n",
    "\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#CELLS GROUPED BY AREA OF THE SOMA\n",
    "\n",
    "# Create a pointplot using Seaborn\n",
    "sns.pointplot(x=\"location\", y=\"i_syn\", hue=\"local_area\", data=synapse_area_table_percentage)\n",
    "\n",
    "# Set the plot title\n",
    "plt.title(\"Area of soma location\")\n",
    "\n",
    "#save figure\n",
    "plt.savefig('Fig3g_interareal_areas_revision.eps')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#ALL CELLS\n",
    "\n",
    "# Create a pointplot using Seaborn\n",
    "sns.pointplot(x=\"location\", y=\"i_syn\", data=synapse_area_table_percentage)\n",
    "\n",
    "# Set the plot title\n",
    "plt.title(\"Area of soma location\")\n",
    "\n",
    "#save figure\n",
    "plt.savefig('Fig3g_interareal_all_revision.eps')\n",
    "\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26b9ac0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURE 3.i\n",
    "import seaborn as sns\n",
    "\n",
    "synapse_table_E_I = synapse_table[(synapse_table['consensus_class'] == 'excitatory') |\n",
    "                                  (synapse_table['consensus_class'] == 'inhibitory')\n",
    "                                 ]\n",
    "\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "kde_dict = {\n",
    "            'bw_adjust' : 0.5,\n",
    "            \n",
    "}\n",
    "\n",
    "#sns.histplot(synapse_table_NP_ET,\n",
    "#            x=\"size\" , hue = 'consensus_subclass', kde = True, kde_kws = kde_dict, fill=False, binwidth=50000, element=\"step\")\n",
    "\n",
    "#sns.histplot(synapse_table_E_I,\n",
    "#            x=\"dist_to_root\" , hue = 'consensus_class',kde = True, kde_kws = kde_dict, binwidth=20000 )\n",
    "\n",
    "# Define a custom color palette\n",
    "custom_palette = ['#c9badb','#d1e6c5']  # Converted RGB colors\n",
    "\n",
    "# Create the histplot with the specified custom colors\n",
    "sns.histplot(synapse_table_E_I,\n",
    "             x=\"dist_to_root\",\n",
    "             hue=\"consensus_class\",\n",
    "             kde=True,\n",
    "             kde_kws=kde_dict,\n",
    "             binwidth=20000,\n",
    "             palette=custom_palette)  # Apply custom colors\n",
    "\n",
    "\n",
    "ax.set(xlim=(0, 1500000))#, xticks=np.arange(0, 1600000 , 100000))\n",
    "ax.set_xticks(np.arange(0, 1500000 , 250000))\n",
    "#ax.set(xlim=(0, 50000))\n",
    "\n",
    "plt.savefig('Fig3i_path_length_revision.eps')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a0077e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURE 3.k\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "def euclidian_distance(row):\n",
    "    if row['num_soma'] == 1:\n",
    "          \n",
    "        #print(row['pre_soma_pt'])\n",
    "        x_pre = row['pre_soma_pt'][0] * 4\n",
    "        y_pre = row['pre_soma_pt'][1] * 4\n",
    "        z_pre = row['pre_soma_pt'][2] * 40\n",
    "        \n",
    "        \n",
    "        #post-synaptic neuron soma\n",
    "        #x_post = row['post_soma_pt'][0] * 4\n",
    "        #y_post = row['post_soma_pt'][1] * 4\n",
    "        #z_post = row['post_soma_pt'][2] * 40\n",
    "        \n",
    "        #synapse location\n",
    "        x_post = row['ctr_pt_position'][0] * 4\n",
    "        y_post = row['ctr_pt_position'][1] * 4\n",
    "        z_post = row['ctr_pt_position'][2] * 40\n",
    " \n",
    "        \n",
    "        \n",
    "        \n",
    "        #print(row['pre_pt_position']-row['post_pt_position'])\n",
    "        #print(row['pre_pt_position'])\n",
    "        return np.linalg.norm(np.array([x_pre, y_pre, z_pre])-np.array([x_post, y_post, z_post]))\n",
    "        #return np.linalg.norm(row['pre_pt_position']-row['post_pt_position'])\n",
    "        \n",
    "    else:\n",
    "        return None\n",
    "\n",
    "synapse_table['euclidian_dist'] = synapse_table.apply(euclidian_distance, axis=1)\n",
    "\n",
    "\n",
    "synapse_table_E_I = synapse_table[(synapse_table['consensus_class'] == 'excitatory') |\n",
    "                                  (synapse_table['consensus_class'] == 'inhibitory')\n",
    "                                 ]\n",
    "\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "kde_dict = {\n",
    "            'bw_adjust' : 0.5,\n",
    "            \n",
    "}\n",
    "\n",
    "# Define a custom color palette\n",
    "custom_palette = ['#c9badb','#d1e6c5']  # Converted RGB colors\n",
    "\n",
    "#sns.histplot(synapse_table_NP_ET,\n",
    "#            x=\"size\" , hue = 'consensus_subclass', kde = True, kde_kws = kde_dict, fill=False, binwidth=50000, element=\"step\")\n",
    "\n",
    "sns.histplot(synapse_table_E_I,\n",
    "            x=\"euclidian_dist\" , hue = 'consensus_class',kde = True, kde_kws = kde_dict, binwidth=10000,\n",
    "            palette=custom_palette)\n",
    "\n",
    "\n",
    "ax.set(xlim=(0, 800000))\n",
    "\n",
    "plt.savefig('Fig3k_euclidean_length_revision.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90b8c579",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j and 3.l\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "raw_df = pd.read_feather(\"ET_extended_synapse_table_revision.feather\")\n",
    "\n",
    "\n",
    "#pre-synaptic neurons and post synaptic cell subclasses\n",
    "\n",
    "studied_nuc_ids = {266839, 589294, 303216, 232635, 494888, 495010, 526436, 302951,\n",
    "       267033, 267029, 527784, 337966}\n",
    "\n",
    "synapse_table = raw_df[raw_df[\"pre_nucleus_id\"].isin(studied_nuc_ids)]\n",
    "\n",
    "#histogram metrics\n",
    "\n",
    "max_hist_dist = 2_000_000\n",
    "euc_dist_step = 100 * 1000\n",
    "path_dist_step = 100 * 1000\n",
    "\n",
    "#cache_dir = \"../data/\"\n",
    "#results_dir = os.path.join(cache_dir, \"results\", \"dendrite-clean-subclass\")\n",
    "\n",
    "#check data inclusion\n",
    "unplotted = synapse_table[synapse_table[\"dist_to_root\"] > max_hist_dist]\n",
    "if len(unplotted) > 0:\n",
    "    print(\"Some data is not captured by the histogram\")\n",
    "\n",
    "\n",
    "# define euc_dist_to_root column\n",
    "resolution = np.array([4., 4., 40.])\n",
    "\n",
    "def vx_to_nm(point):\n",
    "    return point * resolution\n",
    "\n",
    "ctr_locs = np.stack(synapse_table[\"ctr_pt_position\"].apply(vx_to_nm).values)\n",
    "#soma_locs = np.stack(synapse_table[\"pre_pt_position\"].apply(vx_to_nm).values)\n",
    "\n",
    "soma_locs = np.stack(synapse_table[\"pre_soma_pt\"].apply(vx_to_nm).values)\n",
    "synapse_table[\"euc_dist_to_root\"] = np.sqrt(np.sum((ctr_locs - soma_locs)**2, axis=1)) #* 1000  # in nm\n",
    "\n",
    "# rename column\n",
    "synapse_table[\"path_dist_to_root\"] = synapse_table[\"dist_to_root\"]\n",
    "synapse_table[[\"path_dist_to_root\", \"euc_dist_to_root\"]]\n",
    "\n",
    "\n",
    "# Create a dictionary called 'syn_dfs' with the key \"inh\" and a DataFrame containing inhibitory synapses.\n",
    "# It filters the 'synapse_table' DataFrame where the 'consensus_class' column is equal to \"inhibitory\".\n",
    "syn_dfs = {\"inh\": synapse_table[synapse_table[\"consensus_class\"] == \"inhibitory\"]}\n",
    "syn_dfs[\"exc\"] = synapse_table[synapse_table[\"consensus_class\"] == \"excitatory\"]\n",
    "\n",
    "\n",
    "# Determine whether to use counts (False) or density (True) for the histogram\n",
    "density = False\n",
    "\n",
    "# Define bin edges for path distance\n",
    "path_bin_edges = np.arange(0, max_hist_dist, path_dist_step)\n",
    "\n",
    "# Calculate bin centers for path distance\n",
    "path_bin_centers = (path_bin_edges[1:] + path_bin_edges[:-1]) / 2\n",
    "\n",
    "# Create a dictionary to store histograms of counts of each cell type by distance along the axon\n",
    "path_hists = dict()\n",
    "\n",
    "# Loop through each cell type in the 'syn_dfs' dictionary\n",
    "for cell_type in syn_dfs:\n",
    "    # Calculate the histogram of path distances for the current cell type\n",
    "    # and store only the counts (not the bin edges) in the 'path_hists' dictionary\n",
    "    path_hists[cell_type] = np.histogram(syn_dfs[cell_type][\"path_dist_to_root\"], bins=path_bin_edges, density=density)[0]\n",
    "\n",
    "# Define bin edges for Euclidean distance\n",
    "euc_bin_edges = np.arange(0, max_hist_dist, euc_dist_step)\n",
    "\n",
    "# Calculate bin centers for Euclidean distance\n",
    "euc_bin_centers = (path_bin_edges[1:] + path_bin_edges[:-1]) / 2\n",
    "\n",
    "# Create a dictionary to store histograms of counts of each cell type by distance along the axon\n",
    "euc_hists = dict()\n",
    "\n",
    "# Loop through each cell type in the 'syn_dfs' dictionary\n",
    "for cell_type in syn_dfs:\n",
    "    # Calculate the histogram of Euclidean distances for the current cell type\n",
    "    # and store only the counts (not the bin edges) in the 'euc_hists' dictionary\n",
    "    euc_hists[cell_type] = np.histogram(syn_dfs[cell_type][\"euc_dist_to_root\"], bins=euc_bin_edges, density=density)[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61f48eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j and 3.l\n",
    "\n",
    "# Continuously moving bucket to calculate synapse density of each type\n",
    "# This implementation may be slower, but it's functional\n",
    "# A more efficient approach might involve looping over synapses and incrementing counters directly\n",
    "\n",
    "# Number of steps for the moving bucket method\n",
    "nsteps = 1000\n",
    "\n",
    "# Define the bin locations and step size for Euclidean distance\n",
    "euc_bin_locs = np.linspace(euc_dist_step / 2, max_hist_dist - euc_dist_step / 2, nsteps)\n",
    "step_size = max_hist_dist / nsteps\n",
    "\n",
    "# Define the bin locations and step size for path distance\n",
    "path_bin_locs = np.linspace(path_dist_step / 2, max_hist_dist - path_dist_step / 2, nsteps)\n",
    "\n",
    "# Create dictionaries to store Euclidean and path densities for each cell type\n",
    "euc_densities = dict([(cell_type, np.zeros(nsteps)) for cell_type in syn_dfs])\n",
    "path_densities = dict([(cell_type, np.zeros(nsteps)) for cell_type in syn_dfs])\n",
    "\n",
    "# Loop through each step\n",
    "for i in range(nsteps):\n",
    "    # For Euclidean distance\n",
    "    euc_c = euc_bin_locs[i]  # Get the current bin center\n",
    "    euc_l = euc_c - euc_dist_step / 2  # Calculate the left bin edge\n",
    "    euc_r = euc_c + euc_dist_step / 2  # Calculate the right bin edge\n",
    "    for cell_type in syn_dfs:\n",
    "        # Calculate the density within the current bin for the current cell type\n",
    "        # by counting synapses that fall within the bin\n",
    "        euc_densities[cell_type][i] = sum((euc_l < syn_dfs[cell_type][\"euc_dist_to_root\"]).values & (syn_dfs[cell_type][\"euc_dist_to_root\"] < euc_r).values) / euc_dist_step\n",
    "\n",
    "    # For path distance\n",
    "    path_c = path_bin_locs[i]  # Get the current bin center\n",
    "    path_l = path_c - path_dist_step / 2  # Calculate the left bin edge\n",
    "    path_r = path_c + path_dist_step / 2  # Calculate the right bin edge\n",
    "    for cell_type in syn_dfs:\n",
    "        # Calculate the density within the current bin for the current cell type\n",
    "        # by counting synapses that fall within the bin\n",
    "        path_densities[cell_type][i] = sum((path_l < syn_dfs[cell_type][\"path_dist_to_root\"]).values & (syn_dfs[cell_type][\"path_dist_to_root\"] < path_r).values) / path_dist_step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec08b278",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j and 3.l\n",
    "\n",
    "# Create a dictionary to store the proportion of synapses for each cell type by path distances\n",
    "p_type_by_path_dists = dict()\n",
    "\n",
    "# Calculate the total density of synapses by summing densities for all cell types\n",
    "path_total_density = sum(path_densities[t] for t in path_densities)\n",
    "\n",
    "# Loop through each cell type's path densities to calculate proportions\n",
    "for cell_type in path_densities:\n",
    "    # Calculate the proportion of synapses for the current cell type by dividing its density\n",
    "    # by the total density across all cell types\n",
    "    p_type_by_path_dists[cell_type] = path_densities[cell_type] / path_total_density\n",
    "\n",
    "# Create a dictionary to store the proportion of synapses for each cell type by Euclidean distances\n",
    "p_type_by_euc_dists = dict()\n",
    "\n",
    "# Calculate the total density of synapses by summing densities for all cell types\n",
    "euc_total_density = sum(euc_densities[t] for t in euc_densities)\n",
    "\n",
    "# Loop through each cell type's Euclidean densities to calculate proportions\n",
    "for cell_type in euc_densities:\n",
    "    # Calculate the proportion of synapses for the current cell type by dividing its density\n",
    "    # by the total density across all cell types\n",
    "    p_type_by_euc_dists[cell_type] = euc_densities[cell_type] / euc_total_density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "487c2396",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j and 3.l\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Assuming you have already calculated the path_hists and euc_hists dictionaries\n",
    "\n",
    "# Plot histograms for path distance\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.title(\"Histograms of Path Distance\")\n",
    "for cell_type, hist in path_hists.items():\n",
    "    plt.hist(path_bin_centers, bins=path_bin_edges, weights=hist, label=cell_type, alpha=0.5)\n",
    "\n",
    "plt.xlabel(\"Path Distance\")\n",
    "plt.ylabel(\"Count or Density\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Plot histograms for Euclidean distance\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.title(\"Histograms of Euclidean Distance\")\n",
    "for cell_type, hist in euc_hists.items():\n",
    "    plt.hist(euc_bin_centers, bins=euc_bin_edges, weights=hist, label=cell_type, alpha=0.5)\n",
    "\n",
    "plt.xlabel(\"Euclidean Distance\")\n",
    "plt.ylabel(\"Count or Density\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c68715f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j and 3.l\n",
    "\n",
    "# path-length distance split by post-synaptice cell subclass\n",
    "plot_prefix = \"inh_together\"\n",
    "cell_types = [\"inh\"] + [\"exc\"]\n",
    "plot_colors = [\"indigo\"] + [\"lime\"]# * len(exc_cell_types)\n",
    "type_colors = dict(zip(cell_types, plot_colors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3106f882",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.j\n",
    "\n",
    "fig, axes = plt.subplots(1, len(cell_types), figsize=(3 * len(cell_types), 3), dpi=1000)\n",
    "\n",
    "for i, cell_type in enumerate(cell_types):\n",
    "    p = p_type_by_path_dists[cell_type]\n",
    "    axes[i].plot(path_bin_locs/1000, p, color=type_colors[cell_type], label=cell_type)\n",
    "    std_hat = np.sqrt(p * (1 - p) / (path_dist_step * path_total_density))\n",
    "    axes[i].fill_between(path_bin_locs/1000, p - 2 * std_hat, p + 2 * std_hat, color=type_colors[cell_type], alpha=0.2)\n",
    "    axes[i].legend()\n",
    "    axes[i].grid(True)\n",
    "    axes[i].set_ylim((0, 1.2))\n",
    "    #axes[i].set_xlim((0, max_hist_dist/1000))\n",
    "    axes[i].set_xlim((0, 1500))\n",
    "    axes[i].set_xticks(np.arange(0, 1500 , 250))\n",
    "    axes[i].set_ylabel(\"proportion of target cells\")\n",
    "    axes[i].set_xlabel(\"distance along axon ($\\mu m$)\")\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig('Fig3j_path_length_fraction_revision.svg')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3e7891f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIGURES 3.l\n",
    "\n",
    "# path-length distance split by post-synaptice cell subclass\n",
    "fig, axes = plt.subplots(1, len(cell_types), figsize=(3 * len(cell_types), 3), dpi=1000)\n",
    "\n",
    "for i, cell_type in enumerate(cell_types):\n",
    "    print(cell_type)\n",
    "    p = p_type_by_euc_dists[cell_type]\n",
    "    axes[i].plot(euc_bin_locs/1000, p, color=type_colors[cell_type], label=cell_type)\n",
    "    std_hat = np.sqrt(p * (1 - p) / (euc_dist_step * euc_total_density))\n",
    "    axes[i].fill_between(euc_bin_locs/1000, p - 2 * std_hat, p + 2 * std_hat, color=type_colors[cell_type], alpha=0.2)\n",
    "    axes[i].legend()\n",
    "    axes[i].set_ylim((0, 1.2))\n",
    "    axes[i].set_xlim((0, 800))\n",
    "    axes[i].set_xticks(np.arange(0, 800 , 100))\n",
    "    axes[i].set_ylabel(\"proportion of target cells\")\n",
    "    axes[i].set_xlabel(\"Euclidean distance ($\\mu m$)\")\n",
    "    axes[i].grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "print('ping')\n",
    "\n",
    "#plt.show()\n",
    "plt.savefig('Fig3l_euc_fraction_revision.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89ce408a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}