notebooks/rastermap_singleneurons.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/MouseLand/rastermap/blob/main/notebooks/rastermap_singleneurons.ipynb)\n",
    "\n",
    "# Rastermap sorting of 137 neurons\n",
    "\n",
    "We will use a processed version of the data from [Grosmark & Buzsaki, 2016](https://crcns.org/data-sets/hc/hc-11/about-hc-11). 137 neurons from rat hippocampal CA1 were recorded over several hours using eight bilateral silicon-probes. We selected the 33 minute period from the recording in which the rat traverses a linear track of length 1.6 meters. We binned the neural activity into time bins of length 200 ms. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will install the required packages, if not already installed. If on google colab, it will require you to click the \"RESTART RUNTIME\" button because we are updating numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install numpy>=1.24 # (required for google colab)\n",
    "!pip install rastermap\n",
    "!pip install matplotlib"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data and import libraries\n",
    "\n",
    "If not already downloaded, the following cell will automatically download the processed data stored [here](https://osf.io/szmw6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# importing rastermap\n",
    "# (this will be slow the first time since it is compiling the numba functions)\n",
    "from rastermap import Rastermap, utils\n",
    "from scipy.stats import zscore\n",
    "\n",
    "# download processed hippocampus recording from Grosmark & Buzsaki 2016\n",
    "filename = utils.download_data(data_type=\"hippocampus\")\n",
    "\n",
    "dat = np.load(filename)\n",
    "\n",
    "# spks is neurons by time\n",
    "# (each timepoint is 200 ms)\n",
    "spks = dat[\"spks\"]\n",
    "n_neurons, n_time = spks.shape\n",
    "print(f\"{n_neurons} neurons by {n_time} timepoints\")\n",
    "# zscore activity (each neuron activity trace is then mean 0 and standard-deviation 1)\n",
    "spks = zscore(spks, axis=1)\n",
    "\n",
    "# location of the rat and speed\n",
    "loc2d = dat[\"loc2d\"] # 2D location\n",
    "loc_signed = dat[\"loc_signed\"] # left runs are positive and right runs are negative\n",
    "speed = (np.diff(loc2d, axis=0)**2).sum(axis=1)**0.5\n",
    "speed = np.concatenate((np.zeros((1,)), speed), axis=0)\n",
    "\n",
    "# which neurons in the recording are pyramidal cells\n",
    "pyr_cells = dat[\"pyr_cells\"].astype(\"int\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run Rastermap\n",
    "\n",
    "Let's sort the single neurons with Rastermap, skipping clustering and upsampling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Rastermap(n_clusters=None, # None turns off clustering and sorts single neurons \n",
    "                  n_PCs=64, # use fewer PCs than neurons\n",
    "                  locality=0.1, # some locality in sorting (this is a value from 0-1)\n",
    "                  time_lag_window=15, # use future timepoints to compute correlation\n",
    "                  grid_upsample=0, # 0 turns off upsampling since we're using single neurons\n",
    "                ).fit(spks)\n",
    "y = model.embedding # neurons x 1\n",
    "isort = model.isort"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "compute tuning curves along linear corridor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pos = 15\n",
    "bins = np.arange(-1, 1+1./n_pos, 1./n_pos)\n",
    "ibin = np.digitize(loc_signed, bins) - 1\n",
    "n_bins = ibin.max()\n",
    "inan = np.isnan(loc_signed)\n",
    "ibin[inan] = -1\n",
    "tcurves = np.zeros((spks.shape[0], n_bins))\n",
    "for b in range(n_bins):\n",
    "    tcurves[:, b] = spks[:, ibin==b].mean(axis=1)\n",
    "tcurves -= tcurves.mean(axis=1, keepdims=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization\n",
    "\n",
    "Use the Rastermap sorting to visualize the neural activity and tuning curves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# timepoints to visualize\n",
    "xmin = 1000\n",
    "xmax = xmin + 2000\n",
    "\n",
    "# make figure with grid for easy plotting\n",
    "fig = plt.figure(figsize=(12,8), dpi=200)\n",
    "grid = plt.GridSpec(10, 24, figure=fig, wspace = 0.1, hspace = 0.4)\n",
    "\n",
    "# plot location\n",
    "ax = plt.subplot(grid[0, :-5])\n",
    "ax.plot(loc2d[xmin:xmax])#, color=kp_colors[0])\n",
    "ax.set_xlim([0, xmax-xmin])\n",
    "ax.axis(\"off\")\n",
    "ax.set_title(\"2D location\")\n",
    "\n",
    "# plot running speed\n",
    "ax = plt.subplot(grid[1, :-5])\n",
    "ax.plot(speed[xmin:xmax], color=0.5*np.ones(3))\n",
    "ax.set_xlim([0, xmax-xmin])\n",
    "ax.axis(\"off\")\n",
    "ax.set_title(\"running speed\")\n",
    "\n",
    "# plot sorted neural activity\n",
    "ax = plt.subplot(grid[2:, :-5])\n",
    "ax.imshow(spks[isort, xmin:xmax], cmap=\"gray_r\", vmin=0, vmax=1.2, aspect=\"auto\")\n",
    "ax.set_xlabel(\"time\")\n",
    "ax.set_ylabel(\"superneurons\")\n",
    "\n",
    "# excitatory cells in yellow, and inhibitory cells in dark blue\n",
    "# (could replace this with a colorbar or other property)\n",
    "ax = plt.subplot(grid[2:, -5])\n",
    "ax.imshow(pyr_cells[isort, np.newaxis],\n",
    "          cmap=\"viridis\", aspect=\"auto\")\n",
    "ax.axis(\"off\")\n",
    "\n",
    "# plot single-neuron tuning curves\n",
    "ax = plt.subplot(grid[2:, -4:])\n",
    "x = np.arange(0, n_pos)\n",
    "dy = 2\n",
    "xpad = n_pos/10\n",
    "nn = spks.shape[0]\n",
    "for t in range(len(tcurves)):\n",
    "    ax.plot(x, tcurves[isort[t], :n_pos]*dy - dy/2 + nn - t, \n",
    "            color=\"k\", lw=0.5)\n",
    "    ax.plot(x+n_pos+xpad, tcurves[isort[t], n_pos:]*dy - dy/2 + nn - t, \n",
    "            color=\"k\", lw=0.5)\n",
    "for j in range(2):\n",
    "    xstr = \"position\\n(left run)\" if j==0 else \"position\\n(right run)\"\n",
    "    ax.text(n_pos/2 + j*(n_pos+xpad), -14, xstr, ha=\"center\")\n",
    "    ax.text(j*(n_pos+xpad), -3, \"0\")\n",
    "    ax.text(n_pos + j*(n_pos+xpad), -3, \"1.6\", ha=\"right\")\n",
    "ax.set_ylim([0, nn])\n",
    "ax.axis(\"off\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Settings\n",
    "\n",
    "You can see all the rastermap settings with `Rastermap?`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rastermap?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outputs\n",
    "\n",
    "All the attributes assigned to the Rastermap `model` are listed with `Rastermap.fit?`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rastermap.fit?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}