Scientifica_CaImAn_Workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo notebook for analysis of two-photon calcium imaging data\n",
    "\n",
    "The notebook is based on the CaImAn [demo pipeline notebook](https://github.com/flatironinstitute/CaImAn/blob/master/demos/notebooks/demo_pipeline.ipynb) and has been adapted for imaging data acquired with the Scientifica microscope in the Helmchen lab. The notebook shows how to do the following steps:\n",
    "1. Import data\n",
    "2. Specify parameters\n",
    "3. Run motion correction\n",
    "4. Evaluate motion correction\n",
    "5. Run CNMF source extraction\n",
    "6. Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bokeh.plotting as bpl\n",
    "import cv2\n",
    "import glob\n",
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os, platform, re\n",
    "from tifffile import TiffFile\n",
    "\n",
    "try:\n",
    "    cv2.setNumThreads(0)\n",
    "except():\n",
    "    pass\n",
    "\n",
    "try:\n",
    "    if __IPYTHON__:\n",
    "        # this is used for debugging purposes only. allows to reload classes\n",
    "        # when changed\n",
    "        get_ipython().magic('load_ext autoreload')\n",
    "        get_ipython().magic('autoreload 2')\n",
    "except NameError:\n",
    "    pass\n",
    "\n",
    "import caiman as cm\n",
    "from caiman.motion_correction import MotionCorrect\n",
    "from caiman.source_extraction.cnmf import cnmf as cnmf\n",
    "from caiman.source_extraction.cnmf import params as params\n",
    "from caiman.utils.utils import download_demo\n",
    "from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour\n",
    "bpl.output_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on Linux we have to add the caiman folder to Pythonpath\n",
    "if platform.system() == 'Linux':\n",
    "    sys.path.append(os.path.expanduser('~/caiman'))\n",
    "# environment variables for parallel processing\n",
    "os.environ['MKL_NUM_THREADS']='1'\n",
    "os.environ['OPENBLAS_NUM_THREADS']='1'\n",
    "os.environ['VECLIB_MAXIMUM_THREADS']='1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up logger (optional)\n",
    "You can log to a file using the filename parameter, or make the output more or less verbose by setting level to `logging.DEBUG`, `logging.INFO`, `logging.WARNING`, or `logging.ERROR`. A filename argument can also be passed to store the log file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Select data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_folder = '/Volumes/RI_DATA/Locomotion_Pia'\n",
    "day_folder = '20181022'\n",
    "animal_folder = 'mouse4Day2Reg'\n",
    "tif_basename = 'mouse4Day2Reg_' # this should match the start of the filename for all TIF files to be processed in this folder\n",
    "\n",
    "display_movie = False # display the first few movies \n",
    "\n",
    "data_folder = os.path.join(data_folder, day_folder, animal_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create list of files to be processed\n",
    "tif_files = [os.path.join(data_folder, x) for x in os.listdir(data_folder) if x.startswith(tif_basename) and x.endswith('.tif')]\n",
    "\n",
    "print('Processing %1.0f files:' % (len(tif_files)))\n",
    "print(*tif_files[:10], sep='\\n')\n",
    "if len(tif_files) > 10:\n",
    "    print('...')\n",
    "    print(*tif_files[-5:], sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get frame rate from TIFF header (ScanImage tags)\n",
    "tif_hdr = dict()\n",
    "with TiffFile(tif_files[0]) as tif:\n",
    "    for page in tif.pages:\n",
    "        for tag in page.tags.values():\n",
    "            tif_hdr[tag.name] = tag.value\n",
    "        break\n",
    "si_tags = tif_hdr['Software']\n",
    "p = re.compile('scanFrameRate = (\\d{1,3}.\\d{1,6})')\n",
    "frame_rate = float(p.search(si_tags).group(1))\n",
    "\n",
    "print('Frame rate: %1.4f Hz (based on %s header)' % (frame_rate, os.path.basename(tif_files[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Play some movies (optional)\n",
    "This will require loading the movie in memory which in general is not needed by the pipeline. Useful to check if data is imported correctly.\n",
    "\n",
    "Displaying the movie uses the OpenCV library. Press `q` to close the video panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if display_movie:\n",
    "    fnames = tif_files[:2]\n",
    "    m_orig = cm.load_movie_chain(fnames)\n",
    "    ds_ratio = 0.2\n",
    "    m_orig.resize(1, 1, ds_ratio)\n",
    "    print('Prepared concatenated movie for files:')\n",
    "    print(*tif_files[:2], sep='\\n')\n",
    "    print('Run the next cell to play the movie (may take a while)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Todo: debug\n",
    "if display_movie:\n",
    "    m_orig.play(q_max=99.5, fr=frame_rate, magnification=2, backend='notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Specify parameters\n",
    "We set some parameters that are relevant for motion correction, processing with CNMF and component quality evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dataset dependent parameters\n",
    "fr = frame_rate                     # imaging rate in frames per second\n",
    "decay_time = 0.4                    # length of a typical transient in seconds\n",
    "\n",
    "# motion correction parameters\n",
    "strides = (96, 96)          # start a new patch for pw-rigid motion correction every x pixels\n",
    "overlaps = (48, 48)         # overlap between pathes (size of patch strides+overlaps)\n",
    "max_shifts = (24,24)          # maximum allowed rigid shifts (in pixels)\n",
    "max_deviation_rigid = 6     # maximum shifts deviation allowed for patch with respect to rigid shifts\n",
    "pw_rigid = True             # flag for performing non-rigid motion correction\n",
    "splits_rig = 20\n",
    "splits_els = 20\n",
    "\n",
    "# parameters for source extraction and deconvolution\n",
    "p = 1                       # order of the autoregressive system\n",
    "gnb = 2                     # number of global background components\n",
    "merge_thresh = 0.8          # merging threshold, max correlation allowed\n",
    "rf = 50                     # half-size of the patches in pixels. e.g., if rf=25, patches are 50x50\n",
    "stride_cnmf = 24             # amount of overlap between the patches in pixels\n",
    "K = 4                       # number of components per patch\n",
    "gSig = [12, 12]               # expected half size of neurons in pixels\n",
    "method_init = 'greedy_roi'  # initialization method (if analyzing dendritic data using 'sparse_nmf')\n",
    "ssub = 1                    # spatial subsampling during initialization\n",
    "tsub = 1                    # temporal subsampling during intialization\n",
    "\n",
    "# parameters for component evaluation\n",
    "min_SNR = 2.0               # signal to noise ratio for accepting a component\n",
    "rval_thr = 0.85              # space correlation threshold for accepting a component\n",
    "cnn_thr = 0.99              # threshold for CNN based classifier\n",
    "cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create a parameters object\n",
    "You can create a parameters object by passing all the parameters as a single dictionary. Parameters not defined in the dictionary will assume their default values. The resulting `params` object is a collection of subdictionaries pertaining to the dataset to be analyzed `(params.data)`, motion correction `(params.motion)`, data pre-processing `(params.preprocess)`, initialization `(params.init)`, patch processing `(params.patch)`, spatial and temporal component `(params.spatial), (params.temporal)`, quality evaluation `(params.quality)` and online processing `(params.online)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames = tif_files[:2] # select files\n",
    "opts_dict = {'fnames': fnames,\n",
    "            'fr': fr,\n",
    "            'decay_time': decay_time,\n",
    "            'strides': strides,\n",
    "            'overlaps': overlaps,\n",
    "            'max_shifts': max_shifts,\n",
    "            'max_deviation_rigid': max_deviation_rigid,\n",
    "            'splits_rig': splits_rig,\n",
    "            'splits_els': splits_els,\n",
    "            'pw_rigid': pw_rigid,\n",
    "            'p': 1,\n",
    "            'nb': gnb,\n",
    "            'rf': rf,\n",
    "            'K': K, \n",
    "            'stride': stride_cnmf,\n",
    "            'method_init': method_init,\n",
    "            'rolling_sum': True,\n",
    "            'only_init': True,\n",
    "            'ssub': ssub,\n",
    "            'tsub': tsub,\n",
    "            'min_SNR': min_SNR,\n",
    "            'rval_thr': rval_thr,\n",
    "            'use_cnn': True,\n",
    "            'min_cnn_thr': cnn_thr,\n",
    "            'cnn_lowest': cnn_lowest}\n",
    "\n",
    "opts = params.CNMFParams(params_dict=opts_dict)\n",
    "opts.set('motion', {'splits_els': splits_els, 'splits_rig': splits_rig})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opts.get_group('motion')"
   ]
  },
  {
   "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 the cell below. The variable `backend` determines the type of cluster used. The default value `'local'` 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": [
    "n_processes = 8 # None to select automatically\n",
    "if 'dview' in locals():\n",
    "    cm.stop_server(dview=dview)\n",
    "c, dview, n_processes = cm.cluster.setup_cluster(\n",
    "    backend='local', n_processes=n_processes, single_thread=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Run Motion Correction\n",
    "First we create a motion correction object with the parameters specified. Note that the file is not loaded in memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now perform motion correction. From the movie above we see that the dateset exhibits non-uniform motion. We will perform piecewise rigid motion correction using the NoRMCorre algorithm. This has already been selected by setting `pw_rigid=True` when defining the parameters object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "mc.motion_correct(save_movie=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_els = cm.load(mc.fname_tot_els)\n",
    "border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0 \n",
    "    # maximum shift to be used for trimming against NaNs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the results by comparing the original movie. A more detailed presentation of the motion correction method can be found in the [demo motion correction](./demo_motion_correction.ipynb) notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if display_movie:\n",
    "    m_orig = cm.load_movie_chain(fnames)\n",
    "    ds_ratio = 0.2\n",
    "    cm.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie,\n",
    "                    m_els.resize(1, 1, ds_ratio)], \n",
    "                   axis=2).play(fr=frame_rate, gain=15, magnification=2, offset=0, backend='notebook')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc.mmap_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Memory mapping \n",
    "\n",
    "The cell below memory maps the file in order `'C'` and then loads the new memory mapped file. The saved files from motion correction are memory mapped files stored in `'F'` order. Their paths are stored in `mc.mmap_file`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# memory map the file in order 'C'\n",
    "fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C',\n",
    "                           border_to_0=border_to_0) # exclude borders\n",
    "\n",
    "# now load the file\n",
    "Yr, dims, T = cm.load_memmap(fname_new)\n",
    "images = np.reshape(Yr.T, [T] + list(dims), order='F') \n",
    "    #load frames in python format (T x X x Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now restart the cluster to clean up the memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm.stop_server(dview=dview)\n",
    "c, dview, n_processes = cm.cluster.setup_cluster(\n",
    "    backend='local', n_processes=n_processes, single_thread=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run CNMF on patches in parallel\n",
    "- The FOV is split is different overlapping patches that are subsequently processed in parallel by the CNMF algorithm.\n",
    "- The results from all the patches are merged with special attention to idendtified components on the border.\n",
    "- The results are then refined by additional CNMF iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "# First extract spatial and temporal components on patches and combine them\n",
    "# for this step deconvolution is turned off (p=0)\n",
    "opts.set('temporal', {'p': 0})\n",
    "cnm = cnmf.CNMF(n_processes, params=opts, dview=dview)\n",
    "cnm = cnm.fit(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspecting the results\n",
    "Briefly inspect the results by plotting contours of identified components against correlation image.\n",
    "The results of the algorithm are stored in the object `cnm.estimates`. More information can be found in the definition of the `estimates` object and in the [CaImAn Wiki](https://github.com/flatironinstitute/CaImAn/wiki/Interpreting-Results)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Cn = cm.local_correlations(images.transpose(1,2,0))\n",
    "Cn[np.isnan(Cn)] = 0\n",
    "# cnm.estimates.plot_contours_nb(img=Cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Re-run (seeded) CNMF  on the full Field of View  \n",
    "You can re-run the CNMF algorithm seeded on just the selected components from the previous step. Be careful, because components rejected on the previous step will not be recovered here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "cnm.params.set('temporal', {'p': p})\n",
    "cnm2 = cnm.refit(images, dview=dview)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnm2"
   ]
  },
  {
   "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\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% COMPONENT EVALUATION\n",
    "# the components are evaluated in three 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",
    "#   c) each shape passes a CNN based classifier\n",
    "\n",
    "cnm.estimates.evaluate_components(images, cnm.params, dview=dview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot contours of selected and rejected components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% PLOT COMPONENTS\n",
    "cnm.estimates.plot_contours_nb(img=Cn, 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.nb_view_components(img=Cn, idx=cnm.estimates.idx_components)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rejected components\n",
    "if len(cnm2.estimates.idx_components_bad) > 0:\n",
    "    cnm2.estimates.nb_view_components(img=Cn, idx=cnm2.estimates.idx_components_bad)\n",
    "else:\n",
    "    print(\"No components were rejected.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract DF/F values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Extract DF/F values\n",
    "cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select only high quality components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnm2.estimates.select_components(use_object=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display final results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cnm2.estimates.nb_view_components(img=Cn, denoised_color='red')\n",
    "print('you may need to change the data rate to generate this one: use jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10 before opening jupyter notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Closing, saving, and creating denoised version\n",
    "### You can save an hdf5 file with all the fields of the cnmf object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_results = False\n",
    "if save_results:\n",
    "    cnm2.save('analysis_results.hdf5')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stop cluster and clean up LOG files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% STOP CLUSTER and clean up log files\n",
    "cm.stop_server(dview=dview)\n",
    "log_files = glob.glob('*_LOG_*')\n",
    "for log_file in log_files:\n",
    "    os.remove(log_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View movie with the results\n",
    "We can inspect the denoised results by reconstructing the movie and playing alongside the original data and the resulting (amplified) residual movie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnm2.estimates.play_movie(images, q_max=99.9, gain_res=2,\n",
    "                                  magnification=2,\n",
    "                                  bpx=border_to_0,\n",
    "                                  include_bck=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The denoised movie can also be explicitly constructed using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% reconstruct denoised movie\n",
    "denoised = cm.movie(cnm2.estimates.A.dot(cnm2.estimates.C) + \\\n",
    "                    cnm2.estimates.b.dot(cnm2.estimates.f)).reshape(dims + (-1,), order='F').transpose([2, 0, 1])"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "caiman",
   "language": "python",
   "name": "caiman"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}