docs/tutorials/basic_example.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example spike-sorting analysis with sample data\n",
    "\n",
    "This tutorial is also available as a [collab notebook](https://github.com/MouseLand/Kilosort/blob/main/docs/tutorials/kilosort4.ipynb)\n",
    "if you would like to try Kilosort4 without installing the code locally."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Download example data\n",
    "\n",
    "This is an example electrophysiological recording from the International Brain Laboratory, recorded using a Neuropixels 1.0 probe (all data [here](https://ibl.flatironinstitute.org/public/)). The full recording is over 4000 seconds long, and the cropped recording is 90 seconds long.\n",
    "\n",
    "Downloading the recording may take a few minutes. If it fails, please try running the cell again.\n",
    "\n",
    "You can alternatively use any .bin file. See the \"Loading other data formats\" tutorial for loading other file extensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import urllib.request\n",
    "from pathlib import Path\n",
    "from tqdm import tqdm\n",
    "\n",
    "# NOTE: Be sure to update this filepath if you want the data downloaded to\n",
    "#       a specific location.\n",
    "SAVE_PATH = Path('ZFM-02370_mini.imec0.ap.bin')\n",
    "\n",
    "class DownloadProgressBar(tqdm):\n",
    "    \"\"\" from https://stackoverflow.com/a/53877507 \"\"\"\n",
    "    def update_to(self, b=1, bsize=1, tsize=None):\n",
    "        if tsize is not None:\n",
    "            self.total = tsize\n",
    "        self.update(b * bsize - self.n)\n",
    "def download_url(url, output_path):\n",
    "    with DownloadProgressBar(unit='B', unit_scale=True,\n",
    "                             miniters=1, desc=url.split('/')[-1]) as t:\n",
    "        urllib.request.urlretrieve(url, filename=output_path, reporthook=t.update_to)\n",
    "\n",
    "## CROPPED DATASET\n",
    "URL = 'http://www.kilosort.org/downloads/ZFM-02370_mini.imec0.ap.bin'\n",
    "download_url(URL, SAVE_PATH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download channel maps for default probes\n",
    "from kilosort.utils import download_probes\n",
    "download_probes()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run kilosort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from kilosort import run_kilosort\n",
    "\n",
    "# NOTE: 'n_chan_bin' is a required setting, and should reflect the total number\n",
    "#       of channels in the binary file. For information on other available\n",
    "#       settings, see `kilosort.run_kilosort.default_settings`.\n",
    "settings = {'data_dir': SAVE_PATH.parent, 'n_chan_bin': 385}\n",
    "\n",
    "ops, st, clu, tF, Wall, similar_templates, is_ref, est_contam_rate = \\\n",
    "    run_kilosort(settings=settings, probe_name='neuropixPhase3B1_kilosortChanMap.mat')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the results\n",
    "\n",
    "Note: at this point, you can also load the results in `phy`.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "\n",
    "# outputs saved to results_dir\n",
    "results_dir = Path(settings['data_dir']).joinpath('kilosort4')\n",
    "ops = np.load(results_dir / 'ops.npy', allow_pickle=True).item()\n",
    "camps = pd.read_csv(results_dir / 'cluster_Amplitude.tsv', sep='\\t')['Amplitude'].values\n",
    "contam_pct = pd.read_csv(results_dir / 'cluster_ContamPct.tsv', sep='\\t')['ContamPct'].values\n",
    "chan_map =  np.load(results_dir / 'channel_map.npy')\n",
    "templates =  np.load(results_dir / 'templates.npy')\n",
    "chan_best = (templates**2).sum(axis=1).argmax(axis=-1)\n",
    "chan_best = chan_map[chan_best]\n",
    "amplitudes = np.load(results_dir / 'amplitudes.npy')\n",
    "st = np.load(results_dir / 'spike_times.npy')\n",
    "clu = np.load(results_dir / 'spike_clusters.npy')\n",
    "firing_rates = np.unique(clu, return_counts=True)[1] * 30000 / st.max()\n",
    "dshift = ops['dshift']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec, rcParams\n",
    "rcParams['axes.spines.top'] = False\n",
    "rcParams['axes.spines.right'] = False\n",
    "gray = .5 * np.ones(3)\n",
    "\n",
    "fig = plt.figure(figsize=(10,10), dpi=100)\n",
    "grid = gridspec.GridSpec(3, 3, figure=fig, hspace=0.5, wspace=0.5)\n",
    "\n",
    "ax = fig.add_subplot(grid[0,0])\n",
    "ax.plot(np.arange(0, ops['Nbatches'])*2, dshift);\n",
    "ax.set_xlabel('time (sec.)')\n",
    "ax.set_ylabel('drift (um)')\n",
    "\n",
    "ax = fig.add_subplot(grid[0,1:])\n",
    "t0 = 0 \n",
    "t1 = np.nonzero(st > ops['fs']*5)[0][0]\n",
    "ax.scatter(st[t0:t1]/30000., chan_best[clu[t0:t1]], s=0.5, color='k', alpha=0.25)\n",
    "ax.set_xlim([0, 5])\n",
    "ax.set_ylim([chan_map.max(), 0])\n",
    "ax.set_xlabel('time (sec.)')\n",
    "ax.set_ylabel('channel')\n",
    "ax.set_title('spikes from units')\n",
    "\n",
    "ax = fig.add_subplot(grid[1,0])\n",
    "nb=ax.hist(firing_rates, 20, color=gray)\n",
    "ax.set_xlabel('firing rate (Hz)')\n",
    "ax.set_ylabel('# of units')\n",
    "\n",
    "ax = fig.add_subplot(grid[1,1])\n",
    "nb=ax.hist(camps, 20, color=gray)\n",
    "ax.set_xlabel('amplitude')\n",
    "ax.set_ylabel('# of units')\n",
    "\n",
    "ax = fig.add_subplot(grid[1,2])\n",
    "nb=ax.hist(np.minimum(100, contam_pct), np.arange(0,105,5), color=gray)\n",
    "ax.plot([10, 10], [0, nb[0].max()], 'k--')\n",
    "ax.set_xlabel('% contamination')\n",
    "ax.set_ylabel('# of units')\n",
    "ax.set_title('< 10% = good units')\n",
    "\n",
    "for k in range(2):\n",
    "    ax = fig.add_subplot(grid[2,k])\n",
    "    is_ref = contam_pct<10.\n",
    "    ax.scatter(firing_rates[~is_ref], camps[~is_ref], s=3, color='r', label='mua', alpha=0.25)\n",
    "    ax.scatter(firing_rates[is_ref], camps[is_ref], s=3, color='b', label='good', alpha=0.25)\n",
    "    ax.set_ylabel('amplitude (a.u.)')\n",
    "    ax.set_xlabel('firing rate (Hz)')\n",
    "    ax.legend()\n",
    "    if k==1:\n",
    "        ax.set_xscale('log')\n",
    "        ax.set_yscale('log')\n",
    "        ax.set_title('loglog')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probe = ops['probe']\n",
    "# x and y position of probe sites\n",
    "xc, yc = probe['xc'], probe['yc']\n",
    "nc = 16 # number of channels to show\n",
    "good_units = np.nonzero(contam_pct <= 0.1)[0]\n",
    "mua_units = np.nonzero(contam_pct > 0.1)[0]\n",
    "\n",
    "\n",
    "gstr = ['good', 'mua']\n",
    "for j in range(2):\n",
    "    print(f'~~~~~~~~~~~~~~ {gstr[j]} units ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')\n",
    "    print('title = number of spikes from each unit')\n",
    "    units = good_units if j==0 else mua_units \n",
    "    fig = plt.figure(figsize=(12,3), dpi=150)\n",
    "    grid = gridspec.GridSpec(2,20, figure=fig, hspace=0.25, wspace=0.5)\n",
    "\n",
    "    for k in range(40):\n",
    "        wi = units[np.random.randint(len(units))]\n",
    "        wv = templates[wi].copy()  \n",
    "        cb = chan_best[wi]\n",
    "        nsp = (clu==wi).sum()\n",
    "        \n",
    "        ax = fig.add_subplot(grid[k//20, k%20])\n",
    "        n_chan = wv.shape[-1]\n",
    "        ic0 = max(0, cb-nc//2)\n",
    "        ic1 = min(n_chan, cb+nc//2)\n",
    "        wv = wv[:, ic0:ic1]\n",
    "        x0, y0 = xc[ic0:ic1], yc[ic0:ic1]\n",
    "\n",
    "        amp = 4\n",
    "        for ii, (xi,yi) in enumerate(zip(x0,y0)):\n",
    "            t = np.arange(-wv.shape[0]//2,wv.shape[0]//2,1,'float32')\n",
    "            t /= wv.shape[0] / 20\n",
    "            ax.plot(xi + t, yi + wv[:,ii]*amp, lw=0.5, color='k')\n",
    "\n",
    "        ax.set_title(f'{nsp}', fontsize='small')\n",
    "        ax.axis('off')\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.18 ('kilosort4')",
   "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.18"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c450ffd893003221542b409439a3c6dd3d0fae98f595b9e4724fd711cbdc16b9"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}