15-NeuronsBehavior.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neurons & Behavior\n",
    "\n",
    "In this notebook, we'll look at behavior in the Allen Brain Observatory two-photon imaging dataset. Specifically, we'll see if we can use the activity of neurons to predict the animal's running speed.\n",
    "\n",
    "<font color='red'>**Note**: You need to run this notebook in the Docker container that includes allen-brain-observatory (often denoted by 'w/ allen'). Otherwise you will get a permissions error when you attempt to run the cell below.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Typical toolboxes\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy as sp\n",
    "\n",
    "# Plotting setup\n",
    "%config InlineBackend.figure_format = 'retina' # Improve the resolution of our plots!\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step I. Get raw fluorescence & running traces\n",
    "\n",
    "This dataset uses **2-photon calcium imaging** to measure neural activity in mice while they're running on a wheel and viewing visual stimuli. We don't have access to the animal's performance on the task, but we can still use the animal's movement on the wheel (a very rich measure of arousal, actually!) as an interesting behavioral output.\n",
    "\n",
    "First, let's get some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Allen specific toolboxes\n",
    "from allensdk.core.brain_observatory_cache import BrainObservatoryCache\n",
    "\n",
    "boc = BrainObservatoryCache(manifest_file='/datasets/allen-brain-observatory/visual-coding-2p/manifest.json')\n",
    "\n",
    "# Get the drifting gratings experiments for a specific cell ID\n",
    "cell_id = 662074297\n",
    "#cell_id  = 541513979 # An alternate cell ID you can try\n",
    "\n",
    "exps = boc.get_ophys_experiments(cell_specimen_ids=[cell_id], stimuli=['drifting_gratings'])\n",
    "\n",
    "# Let's just look at the first session\n",
    "session_id = exps[0]['id']\n",
    "data_set = boc.get_ophys_experiment_data(session_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than look at single cells, let's look at all of the cells in our field of view. We'll use the `get_dff_traces()` method to get the raw data traces.\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Task</b>: After getting the dff trace, calculate a mean, and plot the mean fluorescence versus time.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all of the fluorescence traces (ts = timestamps; dff = deltaF/F)\n",
    "ts, dff = data_set.get_dff_traces()\n",
    "\n",
    "# Get the average change in fluorescence\n",
    "\n",
    "# Plot the sum of all of the fluorescence activity\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's the imaged activity of all of the cells in our field of view, but does this activity correlate with the behavior of the animal? Specifically, let's look at the mouse's running speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the running speed from this dataset (dxcm = distance in cm; ts = timestamps)\n",
    "dxcm, ts = data_set.get_running_speed()\n",
    "\n",
    "# Plot it\n",
    "plt.figure(figsize=(14,5))\n",
    "plt.plot(ts,dxcm)\n",
    "plt.xlabel('Time (seconds)')\n",
    "plt.ylabel('Speed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2. Clean up the running data and visualize.\n",
    "\n",
    "Running activity is fairly coarse, so rather than try to predict the activity at 30Hz, let's bin the activity (summed across all neurons) into bins of 5 seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bin_size = 5 \n",
    "bin_stamps = np.arange(0,ts.max(),bin_size)\n",
    "num_bins = len(bin_stamps)\n",
    "\n",
    "run_bin = np.empty(num_bins)\n",
    "response_bin = np.empty(num_bins)\n",
    "\n",
    "for i in range(num_bins):\n",
    "    \n",
    "    # Get the values within our time bin and take a mean\n",
    "    run_bin[i] = dxcm[np.where((ts>i*bin_size)&(ts<((i+1)*bin_size)))].mean() \n",
    "    response_bin[i] = dff_mean[np.where((ts>i*bin_size)&(ts<((i+1)*bin_size)))].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at the `run_bin` we just created, we'll see there are some NaN values, since the running data is recorded for slightly less time than the 2P imaging. This might cause a problem later, so let's remove them. \n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Task</b>: Remove the NaNs from the running trace, and remove those same data points from the summed activity trace.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "response_bin = response_bin[~np.isnan(run_bin)]\n",
    "run_bin = run_bin[~np.isnan(run_bin)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot our binned running and response traces together to see if there is, visually, any correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale the response trace\n",
    "response_bin_scaled = response_bin*1000\n",
    "\n",
    "# Plot it\n",
    "plt.figure(figsize=(14,5))\n",
    "plt.plot(run_bin)\n",
    "plt.plot(response_bin_scaled)\n",
    "plt.xlabel('Time (seconds)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Predict fluorescence data from mouse's running speed\n",
    "\n",
    "It's tough to tell visually, but maybe if we run a linear regression we'll be able to see if there is a relationship. \n",
    "\n",
    "When we run this linear regression (using least squares method) on our data, we need to decide which data to *train* the model on, and which data to *test* it on. Let's train it on samples 200 onward, since that contains a lot of variability. We'll then see if we can use the model to predict the first 200 samples.\n",
    "\n",
    "**Side note**: We saw `linregress()` in the correlations notebook, which can work for these purposes. If you ever neet to run a regression with multivariate data, though, you'll need [lstsq](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html) in `scipy.linalg()`. See the commented lines of code below to see how this would work here. You'll notice they return similar, but slightly different, results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Subtract the mean from both response and running\n",
    "response_bin -= response_bin.mean()\n",
    "run_bin -= run_bin.mean()\n",
    "\n",
    "# Run regression on training data & determine the coefficient \n",
    "coef, intercept, r_value, p_value, std_err = sp.stats.linregress(run_bin[200:], response_bin[200:])\n",
    "\n",
    "\n",
    "# Alternative way to do this below, using a slightly different function\n",
    "# from scipy.linalg import lstsq\n",
    "# This ',None' business is because we have to append a column of to calculate intercept\n",
    "# solution, residues, rank, singular_values = lstsq(run_bin[200:,None], response_bin[200:,None])\n",
    "# coef = solution[0]\n",
    "\n",
    "print(coef)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the coefficient we generated to see if we can predict running activity\n",
    "predict = coef*run_bin[:200]\n",
    "\n",
    "# Plot the prediction\n",
    "plt.figure(figsize=(14,5))\n",
    "plt.plot(response_bin[:200])\n",
    "plt.plot(predict)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test how good this prediction is, we can use Pearson's *r*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pearsons_r , pearsons_p = sp.stats.pearsonr(response_bin[:200], predict)\n",
    "print(\"pearson's r: %f    pearson's p: %f\" % (pearsons_r, pearsons_p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4. Try this on other data!\n",
    "\n",
    "Not every cell will have a clear linear relationship between running and neural activity. However, there definitely *are* cells that have a clear relationship! We're more likely to find a clear linear relationship within certain cell types.\n",
    "\n",
    "Consider the following passage from [Yavorska & Weher (2016)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5040712/):\n",
    "\n",
    ">  In visual cortex, running depolarizes both PNs and inhibitory cells. The resulting increase in both excitation and inhibition in PNs reduces membrane potential variance, and leads to more stimulus-evoked spikes without any increase in spontaneous activity (Polack et al., 2013). Whereas cholinergic input affects membrane potential fluctuations during quiescent periods, the effect of locomotion on membrane potential variance is mostly dependent on noradrenergic input. Interestingly, SOM neurons do not show decreased membrane potential variability during running, suggesting a differential influence of norepinephrine on SOM neurons and PNs (Polack et al., 2013). Different classes of inhibitory neurons show marked differences in how they are modulated by locomotion in the visual cortex. VIP neurons are depolarized throughout the entire running period, while PV cells only respond transiently at the beginning. SOM neurons are typically suppressed during running, and fire mostly at the end of the running period (Fu et al., 2014). These results suggest that the effect of locomotion is mediated by a disinhibitory circuit, in which VIP cells inhibit SOM cells and thereby increase the activity of neighboring PNs. \n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>Task</b>: Check to see if there are either <b>single cells</b> in the same experiment above or <b>an experiment in a different cell type</b> have strong correlations with running by editing the code above. Remember to choose a reasonable set of the data to train your model on.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you're curious, you could also try this on the Neuropixels dataset. It will work fairly similarly! The code below will get you started, see [here](https://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html) for more documentation. Once you extract a session, you can use `session.running_speed()` or `session.get_pupil_data()`. See the [Neuropixels Cheat Sheet](https://brainmapportal-live-4cc80a57cd6e400d854-f7fdcae.divio-media.net/filer_public/0f/5d/0f5d22c9-f8f6-428c-9f7a-2983631e72b4/neuropixels_cheat_sheet_nov_2019.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Neuropixels Cache\n",
    "from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache\n",
    "\n",
    "# We have all of this data on the datahub! This is where it lives.\n",
    "manifest_path = '/datasets/allen-brain-observatory/visual-coding-neuropixels/ecephys-cache/manifest.json' \n",
    "\n",
    "# Create the EcephysProjectCache object\n",
    "cache = EcephysProjectCache.fixed(manifest=manifest_path)\n",
    "\n",
    "# Get the sessions available in this dataset\n",
    "sessions = cache.get_session_table()\n",
    "print('Total number of sessions: ' + str(len(sessions)))\n",
    "sessions.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To do this by pulling data from the Allen SDK\n",
    "\n",
    "# Import the Neuropixels Cache\n",
    "from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache\n",
    "\n",
    "manifest_path = 'neuropixels/manifest.json'\n",
    "cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>published_at</th>\n",
       "      <th>specimen_id</th>\n",
       "      <th>session_type</th>\n",
       "      <th>age_in_days</th>\n",
       "      <th>sex</th>\n",
       "      <th>full_genotype</th>\n",
       "      <th>unit_count</th>\n",
       "      <th>channel_count</th>\n",
       "      <th>probe_count</th>\n",
       "      <th>ecephys_structure_acronyms</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>715093703</th>\n",
       "      <td>2019-10-03T00:00:00Z</td>\n",
       "      <td>699733581</td>\n",
       "      <td>brain_observatory_1.1</td>\n",
       "      <td>118.0</td>\n",
       "      <td>M</td>\n",
       "      <td>Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt</td>\n",
       "      <td>884</td>\n",
       "      <td>2219</td>\n",
       "      <td>6</td>\n",
       "      <td>[CA1, VISrl, nan, PO, LP, LGd, CA3, DG, VISl, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>719161530</th>\n",
       "      <td>2019-10-03T00:00:00Z</td>\n",
       "      <td>703279284</td>\n",
       "      <td>brain_observatory_1.1</td>\n",
       "      <td>122.0</td>\n",
       "      <td>M</td>\n",
       "      <td>Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt</td>\n",
       "      <td>755</td>\n",
       "      <td>2214</td>\n",
       "      <td>6</td>\n",
       "      <td>[TH, Eth, APN, POL, LP, DG, CA1, VISpm, nan, N...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>721123822</th>\n",
       "      <td>2019-10-03T00:00:00Z</td>\n",
       "      <td>707296982</td>\n",
       "      <td>brain_observatory_1.1</td>\n",
       "      <td>125.0</td>\n",
       "      <td>M</td>\n",
       "      <td>Pvalb-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt</td>\n",
       "      <td>444</td>\n",
       "      <td>2229</td>\n",
       "      <td>6</td>\n",
       "      <td>[MB, SCig, PPT, NOT, DG, CA1, VISam, nan, LP, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>732592105</th>\n",
       "      <td>2019-10-03T00:00:00Z</td>\n",
       "      <td>717038288</td>\n",
       "      <td>brain_observatory_1.1</td>\n",
       "      <td>100.0</td>\n",
       "      <td>M</td>\n",
       "      <td>wt/wt</td>\n",
       "      <td>824</td>\n",
       "      <td>1847</td>\n",
       "      <td>5</td>\n",
       "      <td>[grey, VISpm, nan, VISp, VISl, VISal, VISrl]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>737581020</th>\n",
       "      <td>2019-10-03T00:00:00Z</td>\n",
       "      <td>718643567</td>\n",
       "      <td>brain_observatory_1.1</td>\n",
       "      <td>108.0</td>\n",
       "      <td>M</td>\n",
       "      <td>wt/wt</td>\n",
       "      <td>568</td>\n",
       "      <td>2218</td>\n",
       "      <td>6</td>\n",
       "      <td>[grey, VISmma, nan, VISpm, VISp, VISl, VISrl]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   published_at  specimen_id           session_type  \\\n",
       "id                                                                    \n",
       "715093703  2019-10-03T00:00:00Z    699733581  brain_observatory_1.1   \n",
       "719161530  2019-10-03T00:00:00Z    703279284  brain_observatory_1.1   \n",
       "721123822  2019-10-03T00:00:00Z    707296982  brain_observatory_1.1   \n",
       "732592105  2019-10-03T00:00:00Z    717038288  brain_observatory_1.1   \n",
       "737581020  2019-10-03T00:00:00Z    718643567  brain_observatory_1.1   \n",
       "\n",
       "           age_in_days sex                                    full_genotype  \\\n",
       "id                                                                            \n",
       "715093703        118.0   M    Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt   \n",
       "719161530        122.0   M    Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt   \n",
       "721123822        125.0   M  Pvalb-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt   \n",
       "732592105        100.0   M                                            wt/wt   \n",
       "737581020        108.0   M                                            wt/wt   \n",
       "\n",
       "           unit_count  channel_count  probe_count  \\\n",
       "id                                                  \n",
       "715093703         884           2219            6   \n",
       "719161530         755           2214            6   \n",
       "721123822         444           2229            6   \n",
       "732592105         824           1847            5   \n",
       "737581020         568           2218            6   \n",
       "\n",
       "                                  ecephys_structure_acronyms  \n",
       "id                                                            \n",
       "715093703  [CA1, VISrl, nan, PO, LP, LGd, CA3, DG, VISl, ...  \n",
       "719161530  [TH, Eth, APN, POL, LP, DG, CA1, VISpm, nan, N...  \n",
       "721123822  [MB, SCig, PPT, NOT, DG, CA1, VISam, nan, LP, ...  \n",
       "732592105       [grey, VISpm, nan, VISp, VISl, VISal, VISrl]  \n",
       "737581020      [grey, VISmma, nan, VISpm, VISp, VISl, VISrl]  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get the sessions available in this dataset\n",
    "sessions = cache.get_session_table()\n",
    "sessions.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's see if we can get data for the first experiment\n",
    "session = cache.get_session_data(715093703)\n",
    "drift = session.get_stimulus_table('drifting_gratings')\n",
    "drift.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}