https://doi.org/10.5281/zenodo.3597474
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.\n",
"When using your own data, be sure to check that you've specified the correct dtype (default is int16) and that the data is in row-major (or 'C') order, the default for NumPy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import urllib.request\n",
"import zipfile\n",
"from pathlib import Path\n",
"from tqdm import tqdm\n",
"\n",
"from kilosort.utils import DOWNLOADS_DIR\n",
"\n",
"\n",
"# NOTE: Be sure to update this filepath if you want the data downloaded to\n",
"# a specific location.\n",
"SAVE_PATH = DOWNLOADS_DIR / '.test_data' / 'ZFM-02370_mini.imec0.ap.short.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",
"\n",
"def download_url(url, output_path):\n",
" # Download zip-compressed data file.\n",
" zip_file = Path(output_path).with_suffix('.zip')\n",
" with DownloadProgressBar(unit='B', unit_scale=True,\n",
" miniters=1, desc=url.split('/')[-1]) as t:\n",
" urllib.request.urlretrieve(url, filename=zip_file, reporthook=t.update_to)\n",
" # Unzip to specified `output_path`.\n",
" with zipfile.ZipFile(zip_file, \"r\") as zip_ref:\n",
" zip_ref.extractall(output_path.parent)\n",
" # Remove zip archive after unzipping.\n",
" zip_file.unlink()\n",
"\n",
"\n",
"## CROPPED DATASET\n",
"URL = 'https://osf.io/download/67effd64f74150d8738b7f34/'\n",
"download_url(URL, SAVE_PATH)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"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 = {'filename': SAVE_PATH, 'n_chan_bin': 385}\n",
"\n",
"ops, st, clu, tF, Wall, similar_templates, is_ref, est_contam_rate, kept_spikes = \\\n",
" run_kilosort(\n",
" settings=settings, probe_name='NeuroPix1_default.mat',\n",
" # save_preprocessed_copy=True\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`run_kilosort` will also accept a list of paths for the `filename` setting, in which case the files will be concatenated in time in the order provided. The built-in `glob` library provides methods for easily collecting lists of file names, but make sure you check the ordering."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you want to save a pre-processed copy of the data (including whitening, high-pass filtering, and drift correction), you can set `save_preprocessed_copy = True` in the arguments for `run_kilosort`. Alternatively, `kilosort.io.save_prepocessing` can be used as a standalone utility to generate the same copy from saved sorting results, but this will not update options for Phy. By default, results are saved in the same directory as the binary data in the `kilosort4` subdirectory."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from kilosort.io import save_preprocessing, load_ops\n",
"\n",
"# NOTE: This will only create the .dat file, it will *NOT* update options for Phy.\n",
"# If you want to use this with Phy, you will need to modify `params.py`\n",
"# in the results directory to point to this new file. Additionally,\n",
"# you must set `hp_filtered=True` and `dtype='int16'` in `params.py``.\n",
"ops_path = SAVE_PATH.parent / 'kilosort4' / 'ops.npy'\n",
"ops = load_ops(ops_path)\n",
"save_preprocessing(SAVE_PATH.parent / 'temp_wh.dat', ops, bfile_path=SAVE_PATH)"
]
},
{
"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": 5,
"metadata": {},
"outputs": [],
"source": [
"from pathlib import Path\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from kilosort.io import load_ops\n",
"\n",
"\n",
"# outputs saved to results_dir\n",
"results_dir = SAVE_PATH.parent / 'kilosort4'\n",
"ops = load_ops(results_dir / 'ops.npy')\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()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.21"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "c450ffd893003221542b409439a3c6dd3d0fae98f595b9e4724fd711cbdc16b9"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
