demos/notebooks/demo_multisession_registration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multisession registration with CaImAn\n",
    "\n",
    "This notebook will help to demonstrate how to use CaImAn on movies recorded in multiple sessions. CaImAn has in-built functions that align movies from two or more sessions and try to recognize components that are imaged in some or all of these recordings.\n",
    "\n",
    "The basic function for this is `caiman.base.rois.register_ROIs()`. It takes two sets of spatial components and finds components present in both using an intersection over union metric and the Hungarian algorithm for optimal matching.\n",
    "`caiman.base.rois.register_multisession()` takes a list of spatial components, aligns sessions 1 and 2, keeps the union of the matched and unmatched components to register it with session 3 and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import get_ipython\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pickle\n",
    "\n",
    "from caiman.base.rois import register_multisession\n",
    "from caiman.utils import visualization\n",
    "from caiman.utils.utils import download_demo\n",
    "\n",
    "try:\n",
    "    if __IPYTHON__:\n",
    "        get_ipython().run_line_magic('load_ext', 'autoreload')\n",
    "        get_ipython().run_line_magic('autoreload', '2')\n",
    "except NameError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We provide an example file generated from data courtesy of Sue Ann Koay and David Tank (Princeton University). The file contains the spatial footprints derived from the CNMF analysis of the same FOV over six different days, as well as a template (correlation image) for each day. The `download_demo` command will automatically download the file and store it in your caiman_data folder the first time you run it. To use the demo in your own dataset you can set:\n",
    "\n",
    "```file_path = '/path/to/file'```\n",
    "\n",
    "or construct a list of spatial footprints and templates and use that to perform the registration as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load multisession data (spatial components and mean intensity templates) (should be replaced by actual data)\n",
    "file_path = download_demo('alignment.pickle')\n",
    "infile = open(file_path,'rb')\n",
    "data = pickle.load(infile)\n",
    "infile.close()\n",
    "\n",
    "spatial = data[0]\n",
    "templates = data[1]\n",
    "dims = templates[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `register_multisession()`\n",
    "\n",
    "The function `register_multisession()` requires 3 arguments:\n",
    "- `A`: A list of ndarrays or scipy.sparse.csc matrices with (# pixels X # component ROIs) for each session\n",
    "- `dims`: Dimensions of the FOV, needed to restore spatial components to a 2D image\n",
    "- `templates`: List of ndarray matrices of size `dims`, template image of each session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_union, assignments, matchings = register_multisession(A=spatial, dims=dims, templates=templates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns 3 variables for further analysis:\n",
    "- `spatial_union`: csc_matrix (# pixels X # total distinct components), the union of all ROIs across all sessions aligned to the FOV of the last session.\n",
    "- `assignments`: ndarray (# total distinct components X # sessions). `assignments[i,j]=k` means that component `k` from session `j` has been identified as component `i` from the union of all components, otherwise it takes a `NaN` value. Note that for each `i` there is at least one session index `j` where `assignments[i,j]!=NaN`.\n",
    "- `matchings`: list of (# sessions) lists. Saves `spatial_union` indices of individual components in each session. `matchings[j][k] = i` means that component `k` from session `j` is represented by component `i` in the union of all components `spatial_union`. In other words `assignments[matchings[j][k], j] = j`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post-alignment screening\n",
    "\n",
    "The three outputs can be used to filter components in various ways. For example we can find the components that were active in at least a given a number of sessions. For more examples, check [this script](https://github.com/flatironinstitute/CaImAn/blob/master/use_cases/eLife_scripts/figure_9/Figure_9_alignment.py) that reproduces the results of [Figure 9, as presented in our eLife paper](https://elifesciences.org/articles/38173#fig9)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter components by number of sessions the component could be found\n",
    "\n",
    "n_reg = 6  # minimal number of sessions that each component has to be registered in\n",
    "\n",
    "# Use number of non-NaNs in each row to filter out components that were not registered in enough sessions\n",
    "assignments_filtered = np.array(np.nan_to_num(assignments[np.sum(~np.isnan(assignments), axis=1) >= n_reg]), dtype=int);\n",
    "\n",
    "# Use filtered indices to select the corresponding spatial components\n",
    "spatial_filtered = spatial[0][:, assignments_filtered[:, 0]]\n",
    "\n",
    "# Plot spatial components of the selected components on the template of the last session\n",
    "visualization.plot_contours(spatial_filtered, templates[-1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining data of components over multiple sessions (optional)\n",
    "\n",
    "Now that all sessions are aligned and we have a list of re-registered neurons, we can use `assignments` and `matchings` to collect traces from neurons over different sessions.\n",
    "\n",
    "As an exercise, we can collect the traces of all neurons that were registered in all sessions. We already gathered the indices of these neurons in the previous cell in `assignments_filtered`. Assuming that traces of each session are saved in their own `CNMF` object collected in a list, we can iterate through `assignments_filtered` and use these indices to find the re-registered neurons in every session.\n",
    "\n",
    "Note: This notebook does not include the traces of the extracted neurons, only their spatial components. As such the loop below will produce an error if you uncomment it. However, it demonstrates how to use the results of the registration to in your own analysis to extract the traces of the same neurons across different sessions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# traces = np.zeros(assignments_filtered.shape, dtype=np.ndarray)\n",
    "# for i in range(traces.shape[0]):\n",
    "#     for j in range(traces.shape[1]):\n",
    "#         traces[i,j] = cnm_list[j].estimates.C[int(assignments_filtered[i,j])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the array `traces`, where element `traces[i,j] = k` is the temporal component of neuron `i` at session `j`. This can be performed with `F_dff` data or `S` spikes as well."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}