demos/notebooks/demo_pipeline_cnmfE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pipeline for microendoscopic data processing in CaImAn using the CNMF-E algorithm\n",
    "This demo presents a complete pipeline for processing microendoscopic data using CaImAn. It includes:\n",
    "- Motion Correction using the NoRMCorre algorithm\n",
    "- Source extraction using the CNMF-E algorithm\n",
    "- Deconvolution using the OASIS algorithm\n",
    "\n",
    "Some basic visualization is also included. The demo illustrates how to `params`, `MoctionCorrection` and `cnmf` object for processing 1p microendoscopic data. For processing two-photon data consult the related `demo_pipeline.ipynb` demo. For more information see the companion CaImAn paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    get_ipython().run_line_magic('load_ext' , 'autoreload')\n",
    "    get_ipython().run_line_magic('autoreload', 2)\n",
    "except:\n",
    "    pass\n",
    "\n",
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "logging.basicConfig(format=\n",
    "                          \"%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s\",\n",
    "                    # filename=\"/tmp/caiman.log\",\n",
    "                    level=logging.WARNING)\n",
    "\n",
    "import caiman as cm\n",
    "from caiman.source_extraction import cnmf\n",
    "from caiman.utils.utils import download_demo\n",
    "from caiman.utils.visualization import inspect_correlation_pnr, nb_inspect_correlation_pnr\n",
    "from caiman.motion_correction import MotionCorrect\n",
    "from caiman.source_extraction.cnmf import params as params\n",
    "from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour\n",
    "import cv2\n",
    "\n",
    "try:\n",
    "    cv2.setNumThreads(0)\n",
    "except:\n",
    "    pass\n",
    "import bokeh.plotting as bpl\n",
    "import holoviews as hv\n",
    "bpl.output_notebook()\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select file(s) to be processed\n",
    "The `download_demo` function will download the specific file for you and return the complete path to the file which will be stored in your `caiman_data` directory. If you adapt this demo for your data make sure to pass the complete path to your file(s). Remember to pass the `fnames` variable as a list. Note that the memory requirement of the CNMF-E algorithm are much higher compared to the standard CNMF algorithm. Test the limits of your system before trying to process very large amounts of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames = ['data_endoscope.tif']  # filename to be processed\n",
    "fnames = [download_demo(fnames[0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup a cluster\n",
    "To enable parallel processing a (local) cluster needs to be set up. This is done with a cell below. The variable `backend` determines the type of cluster used. The default value `'multiprocessing'` uses the multiprocessing package. The `ipyparallel` option is also available. More information on these choices can be found [here](https://github.com/flatironinstitute/CaImAn/blob/master/CLUSTER.md). The resulting variable `dview` expresses the cluster option. If you use `dview=dview` in the downstream analysis then parallel processing will be used. If you use `dview=None` then no parallel processing will be employed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% start a cluster for parallel processing (if a cluster already exists it will be closed and a new session will be opened)\n",
    "if 'dview' in locals():\n",
    "    cm.stop_server(dview=dview)\n",
    "c, dview, n_processes = cm.cluster.setup_cluster(\n",
    "    backend='multiprocessing', n_processes=None, single_thread=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup some parameters\n",
    "We first set some parameters related to the data and motion correction and create a `params` object. We'll modify this object with additional settings later on. You can also set all the parameters at once as demonstrated in the `demo_pipeline.ipynb` notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# dataset dependent parameters\n",
    "frate = 10                       # movie frame rate\n",
    "decay_time = 0.4                 # length of a typical transient in seconds\n",
    "\n",
    "# motion correction parameters\n",
    "motion_correct = True    # flag for performing motion correction\n",
    "pw_rigid = False         # flag for performing piecewise-rigid motion correction (otherwise just rigid)\n",
    "gSig_filt = (3, 3)       # size of high pass spatial filtering, used in 1p data\n",
    "max_shifts = (5, 5)      # maximum allowed rigid shift\n",
    "strides = (48, 48)       # start a new patch for pw-rigid motion correction every x pixels\n",
    "overlaps = (24, 24)      # overlap between pathes (size of patch strides+overlaps)\n",
    "max_deviation_rigid = 3  # maximum deviation allowed for patch with respect to rigid shifts\n",
    "border_nan = 'copy'      # replicate values along the boundaries\n",
    "\n",
    "mc_dict = {\n",
    "    'fnames': fnames,\n",
    "    'fr': frate,\n",
    "    'decay_time': decay_time,\n",
    "    'pw_rigid': pw_rigid,\n",
    "    'max_shifts': max_shifts,\n",
    "    'gSig_filt': gSig_filt,\n",
    "    'strides': strides,\n",
    "    'overlaps': overlaps,\n",
    "    'max_deviation_rigid': max_deviation_rigid,\n",
    "    'border_nan': border_nan\n",
    "}\n",
    "\n",
    "opts = params.CNMFParams(params_dict=mc_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motion Correction\n",
    "The background signal in micro-endoscopic data is very strong and makes the motion correction challenging. \n",
    "As a first step the algorithm performs a high pass spatial filtering with a Gaussian kernel to remove the bulk of the background and enhance spatial landmarks. \n",
    "The size of the kernel is given from the parameter `gSig_filt`. If this is left to the default value of `None` then no spatial filtering is performed (default option, used in 2p data).\n",
    "After spatial filtering, the NoRMCorre algorithm is used to determine the motion in each frame. The inferred motion is then applied to the *original* data so no information is lost.\n",
    "\n",
    "The motion corrected files are saved in memory mapped format. If no motion correction is being performed, then the file gets directly memory mapped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if motion_correct:\n",
    "    # do motion correction rigid\n",
    "    mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion'))\n",
    "    mc.motion_correct(save_movie=True)\n",
    "    fname_mc = mc.fname_tot_els if pw_rigid else mc.fname_tot_rig\n",
    "    if pw_rigid:\n",
    "        bord_px = np.ceil(np.maximum(np.max(np.abs(mc.x_shifts_els)),\n",
    "                                     np.max(np.abs(mc.y_shifts_els)))).astype(int)\n",
    "    else:\n",
    "        bord_px = np.ceil(np.max(np.abs(mc.shifts_rig))).astype(int)\n",
    "        plt.subplot(1, 2, 1); plt.imshow(mc.total_template_rig)  # % plot template\n",
    "        plt.subplot(1, 2, 2); plt.plot(mc.shifts_rig)  # % plot rigid shifts\n",
    "        plt.legend(['x shifts', 'y shifts'])\n",
    "        plt.xlabel('frames')\n",
    "        plt.ylabel('pixels')\n",
    "\n",
    "    bord_px = 0 if border_nan == 'copy' else bord_px\n",
    "    fname_new = cm.save_memmap(fname_mc, base_name='memmap_', order='C',\n",
    "                               border_to_0=bord_px)\n",
    "else:  # if no motion correction just memory map the file\n",
    "    fname_new = cm.save_memmap(fnames, base_name='memmap_',\n",
    "                               order='C', border_to_0=0, dview=dview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load memory mapped file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load memory mappable file\n",
    "Yr, dims, T = cm.load_memmap(fname_new)\n",
    "images = Yr.T.reshape((T,) + dims, order='F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter setting for CNMF-E\n",
    "We now define some parameters for the source extraction step using the CNMF-E algorithm. \n",
    "We construct a new dictionary and use this to modify the *existing* `params` object,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for source extraction and deconvolution\n",
    "p = 1               # order of the autoregressive system\n",
    "K = None            # upper bound on number of components per patch, in general None\n",
    "gSig = (3, 3)       # gaussian width of a 2D gaussian kernel, which approximates a neuron\n",
    "gSiz = (13, 13)     # average diameter of a neuron, in general 4*gSig+1\n",
    "Ain = None          # possibility to seed with predetermined binary masks\n",
    "merge_thr = .7      # merging threshold, max correlation allowed\n",
    "rf = 40             # half-size of the patches in pixels. e.g., if rf=40, patches are 80x80\n",
    "stride_cnmf = 20    # amount of overlap between the patches in pixels\n",
    "#                     (keep it at least large as gSiz, i.e 4 times the neuron size gSig)\n",
    "tsub = 2            # downsampling factor in time for initialization,\n",
    "#                     increase if you have memory problems\n",
    "ssub = 1            # downsampling factor in space for initialization,\n",
    "#                     increase if you have memory problems\n",
    "#                     you can pass them here as boolean vectors\n",
    "low_rank_background = None  # None leaves background of each patch intact,\n",
    "#                     True performs global low-rank approximation if gnb>0\n",
    "gnb = 0             # number of background components (rank) if positive,\n",
    "#                     else exact ring model with following settings\n",
    "#                         gnb= 0: Return background as b and W\n",
    "#                         gnb=-1: Return full rank background B\n",
    "#                         gnb<-1: Don't return background\n",
    "nb_patch = 0        # number of background components (rank) per patch if gnb>0,\n",
    "#                     else it is set automatically\n",
    "min_corr = .8       # min peak value from correlation image\n",
    "min_pnr = 10        # min peak to noise ration from PNR image\n",
    "ssub_B = 2          # additional downsampling factor in space for background\n",
    "ring_size_factor = 1.4  # radius of ring is gSiz*ring_size_factor\n",
    "\n",
    "opts.change_params(params_dict={'method_init': 'corr_pnr',  # use this for 1 photon\n",
    "                                'K': K,\n",
    "                                'gSig': gSig,\n",
    "                                'gSiz': gSiz,\n",
    "                                'merge_thr': merge_thr,\n",
    "                                'p': p,\n",
    "                                'tsub': tsub,\n",
    "                                'ssub': ssub,\n",
    "                                'rf': rf,\n",
    "                                'stride': stride_cnmf,\n",
    "                                'only_init': True,    # set it to True to run CNMF-E\n",
    "                                'nb': gnb,\n",
    "                                'nb_patch': nb_patch,\n",
    "                                'method_deconvolution': 'oasis',       # could use 'cvxpy' alternatively\n",
    "                                'low_rank_background': low_rank_background,\n",
    "                                'update_background_components': True,  # sometimes setting to False improve the results\n",
    "                                'min_corr': min_corr,\n",
    "                                'min_pnr': min_pnr,\n",
    "                                'normalize_init': False,               # just leave as is\n",
    "                                'center_psf': True,                    # leave as is for 1 photon\n",
    "                                'ssub_B': ssub_B,\n",
    "                                'ring_size_factor': ring_size_factor,\n",
    "                                'del_duplicates': True,                # whether to remove duplicates from initialization\n",
    "                                'border_pix': bord_px})                # number of pixels to not consider in the borders)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect summary images and set parameters\n",
    "Check the optimal values of `min_corr` and `min_pnr` by moving slider in the figure that pops up. You can modify them in the `params` object. \n",
    "Note that computing the correlation pnr image can be computationally and memory demanding for large datasets. In this case you can compute\n",
    "only on a subset of the data (the results will not change). You can do that by changing `images[::1]` to `images[::5]` or something similar.\n",
    "This will compute the correlation pnr image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute some summary images (correlation and peak to noise)\n",
    "cn_filter, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile\n",
    "# inspect the summary images and set the parameters\n",
    "nb_inspect_correlation_pnr(cn_filter, pnr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can inspect the correlation and PNR images to select the threshold values for `min_corr` and `min_pnr`. The algorithm will look for components only in places where these value are above the specified thresholds. You can adjust the dynamic range in the plots shown above by choosing the selection tool (third button from the left) and selecting the desired region in the histogram plots on the right of each panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print parameters set above, modify them if necessary based on summary images\n",
    "print(min_corr) # min correlation of peak (from correlation image)\n",
    "print(min_pnr)  # min peak to noise ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the CNMF-E algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnm = cnmf.CNMF(n_processes=n_processes, dview=dview, Ain=Ain, params=opts)\n",
    "cnm.fit(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternate way to run the pipeline at once\n",
    "It is possible to run the combined steps of motion correction, memory mapping, and cnmf fitting in one step as shown below. The command is commented out since the analysis has already been performed. It is recommended that you familiriaze yourself with the various steps and the results of the various steps before using it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview)\n",
    "# cnm1.fit_file(motion_correct=motion_correct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Component Evaluation\n",
    "\n",
    "The processing in patches creates several spurious components. These are filtered out by evaluating each component using three different criteria:\n",
    "\n",
    "- the shape of each component must be correlated with the data at the corresponding location within the FOV\n",
    "- a minimum peak SNR is required over the length of a transient\n",
    "- each shape passes a CNN based classifier\n",
    "\n",
    "<img src=\"../../docs/img/evaluationcomponent.png\"/>\n",
    "After setting some parameters we again modify the existing `params` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% COMPONENT EVALUATION\n",
    "# the components are evaluated in two ways:\n",
    "#   a) the shape of each component must be correlated with the data\n",
    "#   b) a minimum peak SNR is required over the length of a transient\n",
    "# Note that here we do not use the CNN based classifier, because it was trained on 2p not 1p data\n",
    "\n",
    "min_SNR = 3            # adaptive way to set threshold on the transient size\n",
    "r_values_min = 0.85    # threshold on space consistency (if you lower more components\n",
    "#                        will be accepted, potentially with worst quality)\n",
    "cnm.params.set('quality', {'min_SNR': min_SNR,\n",
    "                           'rval_thr': r_values_min,\n",
    "                           'use_cnn': False})\n",
    "cnm.estimates.evaluate_components(images, cnm.params, dview=dview)\n",
    "\n",
    "print(' ***** ')\n",
    "print('Number of total components: ', len(cnm.estimates.C))\n",
    "print('Number of accepted components: ', len(cnm.estimates.idx_components))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do some plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#%% plot contour plots of accepted and rejected components\n",
    "cnm.estimates.plot_contours_nb(img=cn_filter, idx=cnm.estimates.idx_components)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View traces of accepted and rejected components. Note that if you get data rate error you can start Jupyter notebooks using:\n",
    "'jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# accepted components\n",
    "cnm.estimates.hv_view_components(img=cn_filter, idx=cnm.estimates.idx_components,\n",
    "                                denoised_color='red', cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rejected components\n",
    "cnm.estimates.hv_view_components(img=cn_filter, idx=cnm.estimates.idx_components_bad,\n",
    "                                denoised_color='red', cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stop cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm.stop_server(dview=dview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some instructive movies\n",
    "Play the reconstructed movie alongside the original movie and the (amplified) residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with background \n",
    "cnm.estimates.play_movie(images, q_max=99.5, magnification=2,\n",
    "                                 include_bck=True, gain_res=10, bpx=bord_px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# without background\n",
    "cnm.estimates.play_movie(images, q_max=99.9, magnification=2,\n",
    "                                 include_bck=False, gain_res=4, bpx=bord_px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}