in2scienceUK'20_Task.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Task: Neural circuits that make us turn right or left.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMnzgZ/B/guVmXGC8KoVjo9",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/zanemit/materials4in2scienceUK/blob/master/in2scienceUK'20_Task.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DIuPm-Us7ELQ",
        "colab_type": "text"
      },
      "source": [
        "# Investigating the neural mechanisms that make us turn right or left"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YYMVwAhk6wiB",
        "colab_type": "text"
      },
      "source": [
        "##Introduction to the tutorial\n",
        "\n",
        "This tutorial is designed to supplement the 'Neuroscience of Locomotion' webinar by examining the neural circuits involved in left/right turning behaviour. The notebook will guide you through analysis of data modelled on a publication by the laboratory of Ole Kiehn (University of Copenhagen) earlier this year. \n",
        "\n",
        "**!!! Please run the tutorial in Google Colab by clicking on 'OPEN IN COLAB' on the top of this page**\n",
        "\n",
        "The tutorial is written in Python and is thoroughly annotated, so prior coding experience is not required. However, you are encouraged to read through the code and try to understand the logic behind it. You can run the code by clicking on the 'play' button in the top left corner of each cell. \n",
        "\n",
        "As part of the tutorial, you will generate several plots and be asked to interpret them. Being able to make sense of data is an essential skill in science. Reading the background information that is provided throughout the tutorial may also be helpful. \n",
        "\n",
        "There are three tasks, each with associated questions in **bold**. Tasks 1 and 2 can be carried out by simply running the code. In Task 3, you will be asked to fill in a few blanks in the code. That should be possible by examining analogous pieces of code elsewhere in this notebook. \n",
        "\n",
        "The goals of this tutorial are:\n",
        "* to give you an opportunity to think about, and engage in, science on a deep level\n",
        "* to give you insight into the process of data analysis that neuroscientists face on an almost daily basis. \n",
        "\n",
        "The questions in this tutorial are NOT EASY! Some of them might even pose a challenge to senior undergraduate students, so you should not worry if you do not know the answers to some (or, indeed, any) of them! In fact, feeling puzzled and not having answers is also part of being a scientist!\n",
        "\n",
        "Finally, please jot down any questions you have as you go through the tutorial, as you will be encouraged to ask them during the follow-up session."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tx0p5jmya8M2",
        "colab_type": "text"
      },
      "source": [
        "## Setup and data retrieval\n",
        "Run the cell below to:\n",
        "* import the Python packages necessary for this tutorial: `numpy`, `pandas`, and `matplotlib` are some of the most widely used packages in science (and beyond)\n",
        "* create a folder called `in2scienceTaskOutput` in our current working directory (`/content` in Colab by default, you can run `os.getcwd()` to check) for storing data\n",
        "* set several `matplotlib` parameters to standardise the appearance of our plots\n",
        "* download the data "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MziWq3Mebt_C",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# import packages\n",
        "import os\n",
        "import numpy as np\n",
        "import scipy.stats \n",
        "import pandas as pd\n",
        "import seaborn as sns \n",
        "from matplotlib import pyplot as plt \n",
        "from statsmodels.sandbox.stats.multicomp import multipletests\n",
        "from statsmodels.formula.api import ols\n",
        "from statsmodels.stats.anova import anova_lm\n",
        "%matplotlib inline\n",
        "\n",
        "# create a folder \n",
        "SPATH = \"in2scienceTaskOutput/\" #define folder name\n",
        "!mkdir -p \"{SPATH}\" #makes the folder in our current working directory\n",
        "\n",
        "# set matplotlib parameters\n",
        "plt.rcParams['figure.dpi'] = 140\n",
        "plt.rcParams['savefig.dpi'] = 300\n",
        "plt.rcParams['savefig.facecolor'] = (1,1,1,0)\n",
        "plt.rcParams['savefig.bbox'] = \"tight\"\n",
        "plt.rcParams['font.size'] = 10\n",
        "plt.rcParams['pdf.fonttype'] = 42\n",
        "plt.rcParams['xtick.labelsize'] = 10\n",
        "plt.rcParams['ytick.labelsize'] = 10\n",
        "plt.rcParams['axes.labelsize'] = 12\n",
        "\n",
        "# download and extract the data\n",
        "!wget -nc -O \"{SPATH}tutorial_data.zip\" \"https://ndownloader.figshare.com/files/23106284\"\n",
        "!unzip -d \"{SPATH}\" -n \"{SPATH}tutorial_data.zip\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k_MxHmydaoCv",
        "colab_type": "text"
      },
      "source": [
        "##Research context\n",
        "In the webinar, we learned about midbrain neurons that initiate full-bodied locomotion in a straight line. Stimulating those neurons on one side of the brain activates spinal cord circuits bilaterally. In other words, it is a symmetric start command. \n",
        "\n",
        "A related question concerns our ability to direct our movement to the left or right side. One might imagine how such directional control could be achieved through asymmetric activation of the central pattern generators in the spinal cord. Indeed, experiments on reduced preparations, such as an extruded spinal cord in a dish, have supported the existence of such asymmetrically controllable modules. Furthermore, anatomical data has suggested a specific class of neurons in the gigantocellular nucleus (Gi) of the reticular formation as plausible generators of the aforementioned asymmetry. \n",
        "\n",
        "This provides a good starting point for experiments in awake, locomoting animals that could ultimately lead to a thorough understand of turning behaviours."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UvIEe2k02_vn",
        "colab_type": "text"
      },
      "source": [
        "##Experiment 1: chemogenetic activation of specific Gi neurons\n",
        "\n",
        "To investigate the role of genetically defined Gi neurons in turning behaviour, one can artificially activate those neurons using _chemogenetics_ as we saw in the webinar. Thus, the mouse was made to express a specific receptor (hM3Dq) selectively in the Gi neurons of interest on _one side_ (right side) of the brain. When activated by clozapine-N-oxide (CNO), this receptor excites the respective neurons, making them fire action potentials. \n",
        "\n",
        "CNO was delivered into the membrane that lines the abdominal cavity (an _intraperitoneal_ route of administration). From there it can enter the bloodstream and reach the brain. In some sessions, the mouse was injected with saline instead of CNO.\n",
        "\n",
        "In this experiment, a mouse was placed in an empty [cylindrical arena](https://raw.githubusercontent.com/zanemit/images4in2scienceUK/master/cylinder.jpg) and left to explore the arena for 60 minutes. The mouse was filmed from above to monitor its turning behaviour.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-8-c9Mb0q3j_",
        "colab_type": "text"
      },
      "source": [
        "#### A note on data preprocessing\n",
        "Video files were processed with a deep learning algorithm called _DeepLabCut_ that is routinely used in neuroscience to efficiently extract the xy coordinates of specific body parts. Here, we care about the xy coordinates of two body parts: nose `(Nx,Ny)` and tail base `(Tx,Ty)`.\n",
        "\n",
        "Once the xy coordinates of those body parts were obtained, the position of body centre `(Bx,By)` was computed as the midpoint of the line segment between nose and tail base. \n",
        "\n",
        "This allows quantification of body orientation by considering the vector from the body centre to tail base `BT` in each camera frame. $$BT = (Bx-Tx,By-Ty)$$  Angular changes in body orientation can be calculated from the dot product of vectors `BT(t)` and `BT(t+1)` where `t` refers to the time point (camera frame) in the experiment. $$cos\\theta = \\frac{BT(t)\\cdot BT(t+1)}{\\lvert BT(t)\\rvert \\cdot \\lvert BT(t+1) \\rvert}$$\n",
        "\n",
        "A visual illustration of this process can be found [here](https://github.com/zanemit/images4in2scienceUK/blob/master/turning_angle.jpg?raw=true).\n",
        "\n",
        "Now we have a list of angular changes in body orientation corresponding to each camera frame. Remember that our goal is to identify neural circuits that have a strong effect on turning behaviour. Therefore, it is convenient to discretise the data by setting a threshold for the angular change in body orientation that is to be considered a proper turn. Here we define this threshold as 50 degrees, and end up with three possible outcomes at every time point: \n",
        "* right turn (ipsilateral), $50° < \\theta < 180°$, designated as `1` in the data\n",
        "* left turn (contralateral), $-180° < \\theta < -50°$, designated as `-1` in the data\n",
        "* no turn $-50° < \\theta < 50°$, designated as `0` in the data\n",
        "\n",
        "**!!! Note that ipsi- and contra-lateral sides are defined relative to the side of the brain that we perturb with chemogenetics.**\n",
        "\n",
        "The data we use in this tutorial has already been preprocessed and discretised as described above. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y85Jr-qvoNkd",
        "colab_type": "text"
      },
      "source": [
        "### Task 1: examine data from a single session\n",
        "Here we load the data from one session with saline, and one session with CNO, in the same mouse (id: `167.3`). \n",
        "\n",
        "We examine the general structure of the data table by calling `head()`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L7k8yMAwhTIi",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "path = SPATH + \"Gi_hM3Dq_167.3_3.csv\" #define full path to file\n",
        "df1 = pd.read_csv(path, sep = ';') #load data\n",
        "df1.head() #print the first five rows"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "47brOkORbs0v",
        "colab_type": "text"
      },
      "source": [
        "We have the timestamp of the camera frame (in minutes) and the turning outcome: right, left, or no turn (see 'Note on preprocessing' above) in saline and CNO sessions.\n",
        "\n",
        "Let us examine the relative prevalence of ipsilateral (right) or contralateral (left) turns in the two session types. \n",
        "\n",
        "We will do this over 5 minute long rolling windows. We start by defining a custom function `moving_average`. It can seem complicated if you are just starting to learn Python. Therefore, you are welcome to have a look at it, but a detail understanding of it is required to answer the questions in this tutorial."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dUfZTU_ohTE4",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def moving_average(array, window): \n",
        "  '''\n",
        "  A custom function to compute the relative prevalence of ipsi-\n",
        "  and contra-lateral turns over rolling windows.\n",
        "  :param, array - 1D data array\n",
        "  :param, window - rolling window size\n",
        "  ''' \n",
        "  #define empty lists to store the computed values\n",
        "  ipsi_fracs = []\n",
        "  contra_fracs = []\n",
        "\n",
        "  #use a loop to compute moving averages\n",
        "  for t in range(len(array)-win): \n",
        "    #get a subset of data correspond to our window\n",
        "    data = array[t:(t+win)] \n",
        "\n",
        "    #calculate the number of each type of turns\n",
        "    ipsi = list(data).count(1) \n",
        "    contra = list(data).count(-1)\n",
        "\n",
        "    #express each type of turns as fractions of total turns\n",
        "    ipsi_fraction = ipsi / (contra+ipsi)\n",
        "    contra_fraction = 1 - ipsi_fraction\n",
        "\n",
        "    #append these fractions to our (initially empty) lists\n",
        "    ipsi_fracs.append(ipsi_fraction)\n",
        "    contra_fracs.append(contra_fraction)\n",
        "  \n",
        "  #return the lists of turn fractions\n",
        "  return (ipsi_fracs, contra_fracs)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AXRnBExavveq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fps = 1 #frames per second\n",
        "win = fps * 60 * 5 #rolling window size\n",
        "\n",
        "cno_ipsi, cno_contra = moving_average(array=df1['turn_cno'], window=win) #run our custom-defined function on the CNO data\n",
        "saline_ipsi, saline_contra = moving_average(array=df1['turn_sal'], window=win) #run our custom-defined function on the saline data\n",
        "\n",
        "time = df1['time'][int(win/2):int(len(df1['time'])-(win/2))] #get timestamps corresponding to the midpoint of each moving window\n",
        "\n",
        "fig = plt.figure(figsize=(8,4)) #generate a figure template\n",
        "\n",
        "#plot the CNO data\n",
        "axes0 = fig.add_axes([0,.05,.4,.8]) #add a subplot, specifying the left & bottom margins, and subplot width & height respectively\n",
        "axes0.plot(time, cno_ipsi, label = 'ipsi') #plot the fraction of ipsilateral turns against time\n",
        "axes0.plot(time, cno_contra, label = 'contra') #plot the fraction of contralateral turns against time\n",
        "axes0.set_title('CNO') #set title of the subplot\n",
        "axes0.set_ylabel('fraction of turns') #set a label for the y-axis\n",
        "axes0.set_xlabel('time (min)') #set a label for the x-axis\n",
        "axes0.set_ylim(-0.05,1.05) #set the range of the y-axis\n",
        "axes0.legend() #add a legend to the plot\n",
        "\n",
        "#plot the saline data\n",
        "axes1 = fig.add_axes([0.5,.05,.4,.8]) #as before, but notice the difference in the parameter specifying the left margin\n",
        "axes1.plot(time, saline_ipsi, label = 'ipsi')\n",
        "axes1.plot(time, saline_contra, label = 'contra')\n",
        "axes1.set_title('Saline')\n",
        "axes1.set_ylabel('fraction of turns')\n",
        "axes1.set_xlabel('time (min)')\n",
        "axes1.set_ylim(-0.05,1.05)\n",
        "axes1.legend()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KcGGvvvff9aQ",
        "colab_type": "text"
      },
      "source": [
        "**Q-1.1. What is the reason for doing the experiment with saline?**\n",
        "\n",
        "**Q-1.2. Look at the CNO plot. Why is there a gradual change in the fraction of turns at the beginning of the session?** (hint: think about the route of CNO administration)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xytriOD238RC",
        "colab_type": "text"
      },
      "source": [
        "### Task 2: examine data from all sessions with neuronal activation\n",
        "The same experiment was performed with twelve mice. Each mouse performed three sessions with CNO and three sessions with saline. We will analyse only the final 10 minutes of each session to ensure that the effect of our experimental treatment has reached a plateau (see the plot generated in Task 1 for why this is a good idea).\n",
        "\n",
        "Here we load data from all mice and all sessions.\n",
        "\n",
        "As before, we examine the data table by calling `head()`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "21S6tFOQ4tWM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "path = SPATH + \"Gi_hM3Dq_summarised.csv\"\n",
        "df2 = pd.read_csv(path, sep = ';')\n",
        "df2.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QqrUE-jh6hgo",
        "colab_type": "text"
      },
      "source": [
        "This dataset looks different. Column `turn_count` contains the total number of turns observed in the time window (the last 10 minutes) from each session that we wish to analyse.\n",
        "\n",
        "Note that we are only seeing the first five rows of our data table. To understand our dataset better, let us have a look at the range of values in each column."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KLGfMZav7zTO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print('Types of experiment:', np.unique(df2['exp_type']))\n",
        "print('Mice:', np.unique(df2['mouse_id']))\n",
        "print('Sessions:', np.unique(df2['session']))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p4ROiq0x9TIO",
        "colab_type": "text"
      },
      "source": [
        "This verifies that we are still dealing with the same two types of experiment (CNO and saline), and the same two possible outcomes (ipsilateral and contralateral turns). We also see the IDs of all twelve mice used in the experiment, and that each mouse completed up to three sessions.\n",
        "\n",
        "Next we will summarise the data using [boxplots](https://seaborn.pydata.org/_images/seaborn-boxplot-2.png). The advantage of boxplots (over [barplots](https://seaborn.pydata.org/_images/seaborn-barplot-1.png), for example) is that they provide an indication of the spread of the data and point out any outliers.\n",
        "\n",
        "In our boxplot, we want each data point to signify the number of turns per mouse averaged across the three sessions. Therefore, we define another custom function `get_average_per_mouse` that uses `for` loops to iterate over mouse IDs and experiment types to calculate the mean value for each `mouse_id/exp_type` pair."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xqjH2DXz-Gli",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def get_average_per_mouse(array):\n",
        "  summarised_data = [] #create an empty list to store our computed values in\n",
        "  labels = [] #create an empty list to store the correspond exp_type labels in\n",
        "\n",
        "  # run a loop to iterate over mouse IDs\n",
        "  for mouse in np.unique(array['mouse_id']):\n",
        "    data = array[array['mouse_id'] == mouse] #get the subset of data from a given mouse\n",
        "    \n",
        "    #iterate over experiment types and compute the average value per type\n",
        "    aves_per_mouse = [np.mean(data[data['exp_type'] == exp]['turn_count']) for exp in np.unique(data['exp_type'])]\n",
        "    \n",
        "    summarised_data.append(aves_per_mouse) #store the computed values\n",
        "    labels.append(list(np.unique(data['exp_type']))) #store the corresponding 'exp_type' labels\n",
        "\n",
        "  return (np.array(summarised_data), np.array(labels)) #return the average values and their corresponing 'exp_type' labels"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4OSvWtwHUpFi",
        "colab_type": "text"
      },
      "source": [
        "Now, let us run this function on our data, and inspect its output. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EsWnf4J0Left",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "df2summary, labels = get_average_per_mouse(df2)\n",
        "print('data table:\\n',df2summary)\n",
        "print('shape:', np.shape(df2summary))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_3J6sM0yUz1O",
        "colab_type": "text"
      },
      "source": [
        "How do we make sense of this? We see that the data array contains 12 rows and 4 columns, corresponding to the number of mice and experiment types respectively. The order of mouse IDs and experiment types is the same across all rows and columns, and is determined by the `np.unique()` function we used inside our `get_average_per_mouse` function to extract unique values from our data for iteration purposes.\n",
        "\n",
        "Notice that our function returns another array `labels` that contains the `exp_type` label corresponding to each data point. Below, we combine these two arrays in a DataFrame because it makes generation of boxplots very straightforward."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4F5jYyr3NwKJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#combine the arrays into a DataFrame\n",
        "df2boxplot = pd.DataFrame(zip(df2summary.flatten(),labels.flatten()), columns = ['number of turns','experiment type'])\n",
        "\n",
        "fig = plt.figure(figsize=(4,4))\n",
        "axes = fig.add_axes([0,.05,.8,.8])\n",
        "\n",
        "#generate the boxplot\n",
        "axes = sns.boxplot(x = 'experiment type', y = 'number of turns', data = df2boxplot, showfliers = True, fliersize = 1)\n",
        "\n",
        "#add individual data points to the boxplot \n",
        "axes = sns.swarmplot(x = 'experiment type', y = \"number of turns\", data = df2boxplot, color = \".25\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BLsa0fPVuUeA",
        "colab_type": "text"
      },
      "source": [
        "**Q-2.1. Why might it be better to plot the average number of counts per mouse (one point per mouse) rather than the number of counts in each individual session (three points per mouse)?**\n",
        "\n",
        "**Q-2.2. Based on the plot, what does artificial activation of Gi neurons on one side of the brain tell us about turning behaviour?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RMvC4y0HvAOu",
        "colab_type": "text"
      },
      "source": [
        "Although one category clearly stands out, we cannot draw any conclusions yet. In science, we must always apply statistical analysis to determine if any two sets of data are _significantly different_. Although the details of statistical tests are beyond the scope of this tutorial, a few general ideas should be noted:\n",
        "* statistical tests address a _null hypothesis_ which states that there is no significant difference between the data distributions of interest and that any observed differences are due to insufficient sampling or experimental error\n",
        "* statistical tests produce a _test statistic_ which compares the data with what we would expect if the null hypothesis was true\n",
        "* the likelihood that all our data comes from a single distribution, _assuming that the null hypothesis is true_, is known as the _p-value_; you can think of it as the probability that, were you to repeat the experiment again and again, your data distributions would be at least as different from each other as currently observed _if the null hypothesis was true_\n",
        "* traditionally, distributions are considered statistically different if $p < 0.05$. This means that, if the null hypothesis was true and you repeated the experiment many times, the data distributions would be at least as different from each other as they are in your current experiment less than 5% of the time. Therefore, one may argue that the premise of the above statement _that the null hypothesis is true_ might not be true after all, and that the data distributions are statistically different\n",
        "* the $p < 0.05$ threshold is arbitrary and widely considered to be too lenient, so a good rule of thumb is to confidently claim statistical significance only if the p-value is orders of magnitude below 5%\n",
        "\n",
        "Whew! That is probably a lot to take in! Statistical analysis is by no means straightforward, and even well-seasoned scientists sometimes struggle. \n",
        "\n",
        "Let us see how this works in practice. \n",
        "\n",
        "To decide what statistical test to use, we must think about our data a bit more. For convenience of generating a boxplot, we have been talking about four experiment types. However, a more meaningful way to categorise our data would be to consider two types of treatment (CNO vs. saline) and two types of turn (ipsilateral vs. contralateral). The most appropriate statistical test for this type of comparison is two-way Analysis of Variance (ANOVA). It will generate p-values for the following three questions:\n",
        "* the significance of treatment\n",
        "* the significane of turning direction\n",
        "* the significance of treatment-turning direction interaction (_this will be explained below_)\n",
        "\n",
        "First we need to reorganise the data a little bit. You do not have to worry about the details in the code here."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qlIddsgCji9V",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#reorganise the data along with the appropriate treatment and turn labels into a new DataFrame\n",
        "data = df2summary.T.flatten()\n",
        "treatment = np.hstack((['cno']*int(len(data)/2), ['sal']*int(len(data)/2)))\n",
        "turn = np.hstack((['contra']*int(len(data)/4), ['ipsi']*int(len(data)/4),['contra']*int(len(data)/4), ['ipsi']*int(len(data)/4)))\n",
        "df2stats = pd.DataFrame(zip(data,treatment,turn), columns = ['data','treatment','turn'])\n",
        "\n",
        "#run two-way ANOVA\n",
        "formula = 'data ~ C(treatment) + C(turn) + C(treatment):C(turn)'\n",
        "model = ols(formula, df2stats).fit()\n",
        "aov_table = anova_lm(model, typ=2)\n",
        "\n",
        "aov_table"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "76NIAck5omi2",
        "colab_type": "text"
      },
      "source": [
        "The three p-values are printed in the column labeled `PR(>F)`. As you can see, all three are much lower than our $p < 0.05$ threshold. Therefore:\n",
        "* treatment (CNO vs saline) has a significant effect on the number of turns\n",
        "* turn direction (ipsi vs contra) has a significant effect on the number of turns\n",
        "* treatment-turn direction has a significant effect on the number of turns\n",
        "\n",
        "The third statement is the most important one. It means that _the effect of turn direction depends on treatment_. If this interaction term is significant, the other two terms become irrelevant. In other words, interdependence of treatment and turn direction means that it no longer makes sense to consider the effects of these two factors separately.\n",
        "\n",
        "The statistical analysis does not end here, however. We would also like to know exactly which categories are different from one another. To find that out, we will run statistical tests – Student's t-tests in this case – on all pairs of categories addressed by the interaction effect in the two-way ANOVA. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LHG4hYyVn-im",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#split the data into separate variables for convenience\n",
        "cno_contra = df2boxplot[df2boxplot['experiment type'] == 'cno_contra']['number of turns']\n",
        "cno_ipsi = df2boxplot[df2boxplot['experiment type'] == 'cno_ipsi']['number of turns']\n",
        "sal_contra = df2boxplot[df2boxplot['experiment type'] == 'sal_contra']['number of turns']\n",
        "sal_ipsi = df2boxplot[df2boxplot['experiment type'] == 'sal_ipsi']['number of turns']\n",
        "\n",
        "ttests = {} #create another empty dictionary to store the p values in\n",
        "\n",
        "#run t-tests to test how the effect of turn direction depends on treatment\n",
        "ttests['cno_contra_and_cno_ipsi'] = scipy.stats.ttest_ind(cno_contra, cno_ipsi)[1]\n",
        "ttests['sal_contra_and_sal_ipsi'] = scipy.stats.ttest_ind(sal_contra, sal_ipsi)[1]\n",
        "\n",
        "#run t-test to test how the effect of treatment depends on turn direction\n",
        "ttests['cno_contra_and_sal_contra'] = scipy.stats.ttest_ind(cno_contra, sal_contra)[1]\n",
        "ttests['cno_ipsi_and_sal_ipsi'] = scipy.stats.ttest_ind(cno_ipsi, sal_ipsi)[1]\n",
        "\n",
        "ttests"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zjWb0lm5vlCM",
        "colab_type": "text"
      },
      "source": [
        "These are some pretty small p-values. Unfortunately though, this is still not enough. \n",
        "\n",
        "Our $p < 0.05$ significance threshold means that, in any given statistical test, there is a 5% probability of getting a \"significant\" result simply by chance. Therefore, if we ran statistical tests on a hundred pairs of categories, it is likely that five of our tests would turn out \"significant\" just by chance too. To avoid such false positives, we must adjust our signficance threshold. In other words, we must require a lower p-value to consider a result significant. \n",
        "\n",
        "One approach is called a Bonferroni multiple comparisons test, and we apply that below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5in99MYDwkUa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#compute the Bonferroni-adjusted p-values\n",
        "adjusted = multipletests(pvals = list(ttests.values()), alpha = 0.05, method = 'bonferroni')\n",
        "\n",
        "#arrange these p-values in a table for convenient inspection\n",
        "ttests_corrected = pd.DataFrame(zip(ttests.keys(), adjusted[0], adjusted[1]), columns = ['comparison','significance','p-value'])\n",
        "\n",
        "ttests_corrected"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lLtRvwKse1ht",
        "colab_type": "text"
      },
      "source": [
        "**Q-2.3. What do these four p-values tell us about our results? How does it compare with your expectations?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DoARCDDOgX4k",
        "colab_type": "text"
      },
      "source": [
        "##Experiment 2: chemogenetic inhibition of specific Gi neurons\n",
        "\n",
        "A complementary way for investigating the role of genetically defined Gi neurons in turning behaviour is to artificially inhibit those neurons during the same cylinder assay as before. There are many ways to achieve neuronal inhibition, differing by the timescale of their action. \n",
        "\n",
        "In the webinar, we saw the use of diphtheria toxin to permanently 'inhibit' neurons. Here, mice were made to express a specific receptor (hM4Di) selectively in the Gi neurons of interest on _one side_ (right side) of the brain. As you might guess based on the name of this receptor, it is quite similar to the one used to activate neurons in Experiment 1. Indeed, this receptor is also sensitive to CNO. Yet, subtle differences in its biochemical structure mean that, upon activation by CNO, this receptor inhibits the neurons it is expressed in, making them less likely to fire action potentials.\n",
        "\n",
        "As before, the mouse was injected with saline instead of CNO in some sessions. Video data was preprocessed as before."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xESDNsTDkZyg",
        "colab_type": "text"
      },
      "source": [
        "### Task 3: examine data from all sessions with neuronal inhibition\n",
        "Twelve other mice individually performed six sessions in the cylinder: three sessions with CNO and three sessions with saline. We will analyse only the final 10 minutes of each session as before.\n",
        "\n",
        "Here we load data from all mice and all sessions like in Task 2. \n",
        "\n",
        "We analyse data in the exact same manner. You will therefore be asked to type in some lines of code yourself. Do not fear though! All you need to do is find the analogous part of code in Task 2 and change the names of a few variables.\n",
        "\n",
        "**Before you procede, make sure you have run the `get_average_per_mouse` function that we defined in Task 2!**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UADSw4Obk8Fz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "path = SPATH + \"Gi_hM4Di_summarised.csv\"\n",
        "df3 = pd.read_csv(path, sep = ';')\n",
        "df3.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "boLdH5EZuadC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#average the data per mouse across sessions\n",
        "df3summary, labels = get_average_per_mouse(df3)\n",
        "\n",
        "#combine the arrays into a DataFrame\n",
        "df3boxplot = pd.DataFrame(zip(df3summary.flatten(),labels.flatten()), columns = ['number of turns','experiment type'])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oBbPMAm5uRqe",
        "colab_type": "text"
      },
      "source": [
        "**Complete two lines of code below to generate a boxplot and add data points to the plot as we did in Task 2!**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BQBegd2P3e8t",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fig = plt.figure(figsize=(4,4))\n",
        "axes = fig.add_axes([0,.05,.8,.8])\n",
        "\n",
        "#generate the boxplot\n",
        "axes = \n",
        "\n",
        "#add individual data points to the boxplot \n",
        "axes = "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l1ibDWByuvER",
        "colab_type": "text"
      },
      "source": [
        "**Q-3.1. How is this plot different from that in Task 2? Can you explain why?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pB5UjJWBvJuS",
        "colab_type": "text"
      },
      "source": [
        "Let us continue by running the same statistical tests as above. You do not need to write any code in this section."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vK3sfKJKvmJQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#reorganise the data along with the appropriate treatment and turn labels into a new DataFrame\n",
        "data = df3summary.T.flatten()\n",
        "treatment = np.hstack((['cno']*int(len(data)/2), ['sal']*int(len(data)/2)))\n",
        "turn = np.hstack((['contra']*int(len(data)/4), ['ipsi']*int(len(data)/4),['contra']*int(len(data)/4), ['ipsi']*int(len(data)/4)))\n",
        "df3stats = pd.DataFrame(zip(data,treatment,turn), columns = ['data','treatment','turn'])\n",
        "\n",
        "#run two-way ANOVA\n",
        "formula = 'data ~ C(treatment) + C(turn) + C(treatment):C(turn)'\n",
        "model = ols(formula, df3stats).fit()\n",
        "aov_table = anova_lm(model, typ=2)\n",
        "\n",
        "aov_table"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "batrkD6OwL7S",
        "colab_type": "text"
      },
      "source": [
        "**Q-3.2. What can you learn from the p-values produced by this two-way ANOVA test?**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lgWYzhDQw1TJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#split the data into separate variables for convenience\n",
        "cno_contra = df3boxplot[df3boxplot['experiment type'] == 'cno_contra']['number of turns']\n",
        "cno_ipsi = df3boxplot[df3boxplot['experiment type'] == 'cno_ipsi']['number of turns']\n",
        "sal_contra = df3boxplot[df3boxplot['experiment type'] == 'sal_contra']['number of turns']\n",
        "sal_ipsi = df3boxplot[df3boxplot['experiment type'] == 'sal_ipsi']['number of turns']\n",
        "\n",
        "ttests = {} #create another empty dictionary to store the p values in\n",
        "\n",
        "#run t-tests to test how the effect of turn direction depends on treatment\n",
        "ttests['cno_contra_and_cno_ipsi'] = scipy.stats.ttest_ind(cno_contra, cno_ipsi)[1]\n",
        "ttests['sal_contra_and_sal_ipsi'] = scipy.stats.ttest_ind(sal_contra, sal_ipsi)[1]\n",
        "\n",
        "#run t-test to test how the effect of treatment depends on turn direction\n",
        "ttests['cno_contra_and_sal_contra'] = scipy.stats.ttest_ind(cno_contra, sal_contra)[1]\n",
        "ttests['cno_ipsi_and_sal_ipsi'] = scipy.stats.ttest_ind(cno_ipsi, sal_ipsi)[1]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tsQa9Sbqw-Tj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#compute the Bonferroni-adjusted p-values\n",
        "adjusted = multipletests(pvals = list(ttests.values()), alpha = 0.05, method = 'bonferroni')\n",
        "\n",
        "#arrange these p-values in a table for convenient inspection\n",
        "ttests_corrected = pd.DataFrame(zip(ttests.keys(), adjusted[0], adjusted[1]), columns = ['comparison','significance','p-value'])\n",
        "\n",
        "ttests_corrected"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ug_dIvnxa1K",
        "colab_type": "text"
      },
      "source": [
        "**Q-3.3. Based on your plot and these p-values, what does artificial inhibition of Gi neurons on one side of the brain tell us about turning behaviour?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9u2VJQmkyClL",
        "colab_type": "text"
      },
      "source": [
        "## Conclusions\n",
        "Well done on completing this tutorial! In many cases, analysing data takes up the majority of neuroscientists' daily lives, and now you have had a (hopefully enjoyable) glimpse into those!\n",
        "\n",
        "At the end of analysis, it is always important to bring all experiments together and consider the bigger picture.\n",
        "\n",
        "**Q-4.1. What can you conclude about the involvement of Gi neurons in turning behaviour based on both experiments in this tutorial?**\n",
        "\n",
        "**Q-4.2. Can you suggest any other experiments that should be done to further validate your current conclusions?**\n",
        "\n",
        "**Q-4.3. As an aspiring systems neuroscientist, what else might you want to learn about the neural circuits mediating turning behaviour to be able to claim a full understanding of how the central nervous system makes us turn right or left?**\n",
        "\n",
        "I hope you had fun doing this tutorial. See you in the follow-up session!"
      ]
    }
  ]
}