14-Correlations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correlations\n",
    "\n",
    "In this notebook, we'll explore correlations in single cell two-photon imaging data as well as extracellular recordings.\n",
    "\n",
    "### By the end of this lesson, you'll be able to:\n",
    "- Pull signal correlation and distance information from the Brain Observatory dataset.\n",
    "- Use `imshow` and `scatter` to inspect and plot these data.\n",
    "- Grab data from the Allen Neuropixels dataset, and look for correlations between cells.\n",
    "\n",
    "### Table of Contents\n",
    "1. [Part One: Correlations in two-photon imaging data](#one)\n",
    "2. [Part Two: Correlations in extracellular recordings](#two)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Allen specific toolboxes\n",
    "import allensdk.brain_observatory.stimulus_info as stim_info\n",
    "from allensdk.core.brain_observatory_cache import BrainObservatoryCache\n",
    "\n",
    "# Additional 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": [
    "<a id=\"one\"></a>\n",
    "\n",
    "## Part I. Correlations in two-photon imaging data.\n",
    "\n",
    "### Step 1. Create a BrainObservatoryCache object and create a dataframe of experiments\n",
    "Just as we did in the Brain Observatory notebook, here we'll choose a visual_area and cre_line and look at all of the experiments that meet that criteria.\n",
    "\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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# We will create an instance of the Brain Observatory Cache as an object, \"boc.\"\n",
    "boc = BrainObservatoryCache(manifest_file='/datasets/allen-brain-observatory/visual-coding-2p/manifest.json')\n",
    "\n",
    "# Choose a visual area and cre-line\n",
    "visual_area = 'VISp'\n",
    "cre_line ='Cux2-CreERT2'\n",
    "\n",
    "# Get the experiment containers for that visual area & cre line\n",
    "exps = boc.get_experiment_containers(targeted_structures=[visual_area], cre_lines=[cre_line])\n",
    "\n",
    "# Make a dataframe and look at it to remind ourselves what's in here.\n",
    "exps = pd.DataFrame(exps)\n",
    "exps.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2. Choose a specific experiment\n",
    "For our exercise today, we'll take a look at correlations in a field of view from two-photon imaging data. Let's first look at data that was collected while mice were viewing natural scenes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the first experiment container id\n",
    "expt_container_id = exps['id'][0]\n",
    "\n",
    "# Get the first session id\n",
    "session_id = boc.get_ophys_experiments(\n",
    "    experiment_container_ids=[expt_container_id],\n",
    "    stimuli=['natural_scenes'])[0]['id']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3. Get the computed signal correlation \n",
    "To do so, we'll use `get_ophys_experiment_analysis()` to import an analysis object as `ns`. This object contains many computed metrics, such as the signal correlation between single neurons (the `signal_correlation` attribute of our `ns` object)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the natural scenes analysis\n",
    "ns = boc.get_ophys_experiment_analysis(ophys_experiment_id=session_id, stimulus_type='natural_scenes')\n",
    "\n",
    "# Get the signal correlation\n",
    "sc = ns.signal_correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Task</b>: Look at the shape of the <code>sc</code> object we created above. Why is it that shape? (Hint: look at the available methods and attributes from <code>ns</code> for additional information about this dataset.)\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Task</b>: Now that we have an idea of what this <code>sc</code> object is, visualize it using <code>plt.imshow()</code>. What should the x and y axes be labeled here?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4. Look at the raw data for this experiment\n",
    "So some of these cells are correlated more than other cells -- you might wonder whether how *far* the cells are from eachother actually predicts these correlations. To figure this out, we first need to calculate the distance between cells in our field of view. To remind us that we're looking down at a bunch of cells, let's look at the maximum projection of all of the imaging in this session."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the max projection of the data and show it\n",
    "dataset = boc.get_ophys_experiment_data(ophys_experiment_id=session_id)\n",
    "max_proj = dataset.get_max_projection()\n",
    "plt.imshow(max_proj, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5. Get the regions of interest\n",
    "In the block below, we'll get the regions of interest (ROIs) using the `get_roi_mask_array()` method on our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rois = dataset.get_roi_mask_array()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6. Calculate distance between ROIs\n",
    "Now we can calculate the distance between our ROIs, using simple geometric distance. This distance that we calculate will be in \"pixels\" of the image -- we'd need to know the size of our field of view in order to get this into actual units of distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_cells = rois.shape[0]\n",
    "loc_x = np.zeros((num_cells))\n",
    "loc_y = np.zeros((num_cells))\n",
    "\n",
    "for i in range(num_cells):\n",
    "    ind = np.where(rois[i])\n",
    "    loc_x[i] = np.mean(ind[1])\n",
    "    loc_y[i] = np.mean(ind[0])\n",
    "\n",
    "distance = np.zeros((num_cells, num_cells))\n",
    "for i in range(num_cells):\n",
    "    for j in range(num_cells):\n",
    "        distance[i, j] = np.sqrt( (loc_x[i]-loc_x[j])**2 + (loc_y[i]-loc_y[j])**2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 7. Plot the distance matrix between cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the distance matrix\n",
    "plt.imshow(distance)\n",
    "plt.xlabel(\"Cell #\")\n",
    "plt.ylabel(\"Cell #\")\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 8. Plot the distribution of distances.\n",
    "Below, we'll plot the distribution of distances. Importantly, we only count each distance once, so we'll use a numpy function called <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.triu_indices.html\">triu_indices</a> to get the x and y indices of the the values. Be sure to eliminate the diagonal itself)\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inds = np.triu_indices(num_cells, k=1)\n",
    "distance_vector = distance[inds[0], inds[1]]\n",
    "plt.hist(distance_vector, bins=100)\n",
    "plt.xlabel(\"Distance (pixels)\")\n",
    "plt.ylabel(\"Number of pairs\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 9. See if there is a relationship between distance and signal correlation!\n",
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Task</b>: Plot the distance versus the single cell correlation values, to see if they correlate. Add this where indicated.\n",
    "</div>\n",
    "\n",
    "After plotting your scatterplot, there is code below to calculate the linear least squares regression of these data, using `scipy.stats.linregress()` ([documentation here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html)).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "\n",
    "# YOUR SCATTER PLOT HERE\n",
    "\n",
    "# Calculate the linear regression\n",
    "slope, intercept, r_value, p_value, std_err = sp.stats.linregress(distance_vector, sc[inds[0], inds[1]])\n",
    "plt.plot(distance_vector,intercept+slope*distance[inds[0], inds[1]],\"k-\")\n",
    "print(\"r_value: %f    p_value: %f\" % (r_value, p_value))\n",
    "\n",
    "# Calculate the Pearson's correlation\n",
    "pearsons_r , pearsons_p = sp.stats.pearsonr(distance_vector, sc[inds[0], inds[1]])\n",
    "print(\"pearson's r: %f    pearson's p: %f\" % (pearsons_r, pearsons_p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"two\"></a>\n",
    "\n",
    "## Part II: Correlations in extracellular recordings \n",
    "\n",
    "We can also look at correlations between cells in extracellular recordings.\n",
    "\n",
    "### Step 1. Get the Neuropixels data.\n",
    "As always, we first need to get the data! To do so, we have two options:\n",
    "\n",
    "#### Option 1: From the data stored on the DataHub.\n",
    "This option requires that the DataHub is updated, and that you're running this from the DataHub."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# To do this with data stored on the datahub \n",
    "\n",
    "# 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Option 2: Directly from the Allen\n",
    "If we don't have the data on the DataHub, we can instead pull it directly from the Allen's data storage. This takes some time as well as some data storage space (just over 50 MB), so if possible, the above option should be chosen. See additional details in the documentation and Data Access notebook found [here](https://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Step 2. Look at what is available in this dataset, and pull relevant sessions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the sessions available in this dataset\n",
    "sessions = cache.get_session_table()\n",
    "sessions.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These data contain recordings across many visual areas (see column \"ecephys_structure_acroynyms\" above).\n",
    "\n",
    "For the sake of continuity, let's just look at experiments that contain the same area we looked at in Part I."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create a session list based on some criteria\n",
    "session_list = []\n",
    "\n",
    "for idx,structure_list in enumerate(sessions['ecephys_structure_acronyms']):\n",
    "    if visual_area in structure_list:\n",
    "        session_list.append(sessions.index[idx])   \n",
    "        \n",
    "print('There are '+str(len(session_list))+' sessions that meet this criteria:')\n",
    "print(session_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's just look at the first experiment session, and get that data using `get_session_data()`. We will then get the response to drifting gratings using the `get_stimulus_table()` method on our session object.\n",
    "\n",
    "<font color='red'>**Note**: <code>get_session_data()</code> will download a ~2 GB file. This will take a while, unfortunately.</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "session_id = session_list[0]\n",
    "session = cache.get_session_data(session_id)\n",
    "drift = session.get_stimulus_table('drifting_gratings')\n",
    "drift.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3. Get the spike counts for a drifting grating stimulus.\n",
    "We need to do a little bit of wrangling to get the spike counts for our drifting grating stimulus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the info about the drifting gratings session\n",
    "first_drift_id = drift.index.values[0]\n",
    "first_drift_duration = drift.loc[first_drift_id, \"stop_time\"] - drift.loc[first_drift_id, \"start_time\"]\n",
    "\n",
    "# construct the time domain at 10 ms resolution\n",
    "time_step = 1 / 100\n",
    "time_domain = np.arange(0.0, first_drift_duration + time_step, time_step)\n",
    "\n",
    "histograms_drift = session.presentationwise_spike_counts(\n",
    "    bin_edges=time_domain,\n",
    "    stimulus_presentation_ids=drift.index,\n",
    "    unit_ids=None\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spike_counts_drift = histograms_drift[0]\n",
    "max_len=spike_counts_drift.shape[0]\n",
    "\n",
    "# get two spike trains in drifting gratings activity\n",
    "spike_train_1_drift=spike_counts_drift[:max_len, 81]\n",
    "spike_train_2_drift=spike_counts_drift[:max_len, 83]\n",
    "\n",
    "# plot the spike trains for which the correlogram will be computed\n",
    "\n",
    "fig, ax = plt.subplots(1, 2,figsize=(10,4))\n",
    "\n",
    "ax[0].plot(spike_train_1_drift)\n",
    "ax[0].set_ylabel('Spike count 1')\n",
    "ax[0].set_xlabel('Time bins (10 ms)')\n",
    "\n",
    "ax[1].plot(spike_train_2_drift)\n",
    "ax[1].set_ylabel('Spike count 2')\n",
    "ax[1].set_xlabel('Time bins (10 ms)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the correlogram for spontaneous activity\n",
    "xcorr_drift=sp.signal.correlate(spike_train_1_drift,spike_train_2_drift)\n",
    "\n",
    "# time steps\n",
    "time_shift_drift=np.arange(-len(xcorr_drift)/2,len(xcorr_drift)/2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've computed the cross-correlations between our spike trains, we can plot the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the cross-correlations for drifting gratings\n",
    "\n",
    "plt.figure(figsize=(14,8))\n",
    "plt.plot(time_shift_drift,xcorr_drift)\n",
    "plt.ylabel('Signal correlation')\n",
    "plt.xlabel('Time steps (10 ms)')\n",
    "plt.title('Drifting gratings')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "<b>Task</b>:\n",
    "1. Find two cells that are clearly correlated in time.\n",
    "2. Change the brain area we are looking at -- are spikes more or less correlated there?\n",
    "3. Instead of looking for correlations in response to drifting gratings, look for responses during spontaneous sessions ('spontaneous' in <code>sessions.get_stimulus_table()</code>).\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About this notebook\n",
    "\n",
    "Exercises here inspired by [this Brain Observatory notebook](https://github.com/AllenInstitute/SWDB_2019/blob/master/DynamicBrain/solutions/Other/BrainObservatory_solutions.ipynb) and [this Neuropixels notebook](https://github.com/AllenInstitute/SWDB_2019/blob/master/DynamicBrain/solutions/Neuropixels_exercises_solutions.ipynb) by the Allen Institute."
   ]
  }
 ],
 "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
}