This notebook will help to demonstrate the process of CaImAn and how it uses different functions to denoise, deconvolve and demix neurons from a two-photon Calcium Imaging dataset. The demo shows how to construct the `params`, `MotionCorrect` and `cnmf` objects and call the relevant functions. You can also run a large part of the pipeline with a single method (`cnmf.fit_file`). See inside for details. Dataset couresy of Sue Ann Koay and David Tank (Princeton University) This demo pertains to two photon data. For a complete analysis pipeline for one photon microendoscopic data see demo_pipeline_cnmfE.ipynb
More information can be found in the companion paper.
In [4]:
import bokeh.plotting as bpl
import cv2
import glob
import logging
import matplotlib.pyplot as plt
import numpy as np
import os
try:
cv2.setNumThreads(0)
except():
pass
try:
if __IPYTHON__:
# this is used for debugging purposes only. allows to reload classes
# when changed
get_ipython().magic('load_ext autoreload')
get_ipython().magic('autoreload 2')
except NameError:
pass
import caiman as cm
from caiman.motion_correction import MotionCorrect
from caiman.source_extraction.cnmf import cnmf as cnmf
from caiman.source_extraction.cnmf import params as params
from caiman.utils.utils import download_demo
from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour
bpl.output_notebook()
from caiman.paths import caiman_datadir
In [5]:
logging.basicConfig(format=
"%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s",
# filename="/tmp/caiman.log",
level=logging.WARNING)
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 fname variable as a list.
In [6]:
fnames = [os.path.join(caiman_datadir(), 'example_movies/Sue_2x_3000_40_-46.nwb')] # filename to be created or processed
#fnames = ['Sue_2x_3000_40_-46.tif'] # filename to be processed
#if fnames[0] in ['Sue_2x_3000_40_-46.tif', 'demoMovie.tif']:
# fnames = [download_demo(fnames[0])]
In [7]:
fr = 15 # imaging rate in frames per second
decay_time = 0.4 # length of a typical transient in seconds
#%% load the file and save it in the NWB format (if it doesn't exist already)
if not os.path.exists(fnames[0]):
fnames_orig = 'Sue_2x_3000_40_-46.tif' # filename to be processed
if fnames_orig in ['Sue_2x_3000_40_-46.tif', 'demoMovie.tif']:
fnames_orig = [download_demo(fnames_orig)]
print(fnames_orig)
orig_movie = cm.load(fnames_orig, fr=fr)
# save file in NWB format with various additional info
orig_movie.save(fnames[0], sess_desc='test', identifier='demo 1',
exp_desc='demo movie', imaging_plane_description='single plane',
emission_lambda=520.0, indicator='GCAMP6f',
location='parietal cortex', starting_time=0.,
experimenter='Sue Ann Koay', lab_name='Tank Lab',
institution='Princeton U',
experiment_description='Experiment Description',
session_id='Session 1')
We set some parameters that are relevant to the file, and then parameters for motion correction, processing with CNMF and component quality evaluation. Note that the dataset Sue_2x_3000_40_-46.tif has been spatially downsampled by a factor of 2 and has a lower than usual spatial resolution (2um/pixel). As a result several parameters (gSig, strides, max_shifts, rf, stride_cnmf) have lower values (halved compared to a dataset with spatial resolution 1um/pixel).
In [8]:
# dataset dependent parameters
# these next to parameters have been defined in the previous cell
#fr = 15 # imaging rate in frames per second
#decay_time = 0.4 # length of a typical transient in seconds
# dataset dependent parameters
#fr = 30 # imaging rate in frames per second
#decay_time = 0.4 # length of a typical transient in seconds
# motion correction parameters
dxy = (2., 2.) # spatial resolution in x and y in (um per pixel)
# note the lower than usual spatial resolution here
max_shift_um = (12., 12.) # maximum shift in um
patch_motion_um = (100., 100.) # patch size for non-rigid correction in um
pw_rigid = True # flag for performing non-rigid motion correction
max_shifts = [int(a/b) for a, b in zip(max_shift_um, dxy)]
# start a new patch for pw-rigid motion correction every x pixels
strides = tuple([int(a/b) for a, b in zip(patch_motion_um, dxy)])
# overlap between pathes (size of patch in pixels: strides+overlaps)
overlaps = (24, 24)
# maximum deviation allowed for patch with respect to rigid shifts
max_deviation_rigid = 3
# parameters for source extraction and deconvolution
p = 1 # order of the autoregressive system
gnb = 2 # number of global background components
merge_thr = 0.85 # merging threshold, max correlation allowed
rf = 15 # half-size of the patches in pixels. e.g., if rf=25, patches are 50x50
stride_cnmf = 6 # amount of overlap between the patches in pixels
K = 4 # number of components per patch
gSig = [4, 4] # expected half size of neurons in pixels
method_init = 'greedy_roi' # initialization method (if analyzing dendritic data using 'sparse_nmf')
ssub = 2 # spatial subsampling during initialization
tsub = 2 # temporal subsampling during intialization
# parameters for component evaluation
min_SNR = 2.0 # signal to noise ratio for accepting a component
rval_thr = 0.85 # space correlation threshold for accepting a component
cnn_thr = 0.99 # threshold for CNN based classifier
cnn_lowest = 0.1 # neurons with cnn probability lower than this value are rejected
You can creating 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)
In [9]:
opts_dict = {'fnames': fnames,
'fr': fr,
'decay_time': decay_time,
'dxy': dxy,
'strides': strides,
'overlaps': overlaps,
'max_shifts': max_shifts,
'max_deviation_rigid': max_deviation_rigid,
'pw_rigid': pw_rigid,
'border_nan': 'copy',
'var_name_hdf5': 'acquisition/mov',
'p': p,
'nb': gnb,
'rf': rf,
'K': K,
'stride': stride_cnmf,
'method_init': method_init,
'rolling_sum': True,
'only_init': True,
'ssub': ssub,
'tsub': tsub,
'gSig': gSig,
'merge_thr': merge_thr,
'min_SNR': min_SNR,
'rval_thr': rval_thr,
'use_cnn': True,
'min_cnn_thr': cnn_thr,
'cnn_lowest': cnn_lowest}
opts = params.CNMFParams(params_dict=opts_dict)
In [10]:
display_movie = False
if display_movie:
m_orig = cm.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5'])
ds_ratio = 0.2
m_orig.resize(1, 1, ds_ratio).play(
q_max=99.5, fr=30, magnification=2)
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 'local' uses the multiprocessing package. The ipyparallel option is also available. More information on these choices can be found here. 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.
In [11]:
#%% start a cluster for parallel processing (if a cluster already exists it will be closed and a new session will be opened)
if 'dview' in locals():
cm.stop_server(dview=dview)
c, dview, n_processes = cm.cluster.setup_cluster(
backend='local', n_processes=None, single_thread=False)
In [12]:
# first we create a motion correction object with the parameters specified
mc = MotionCorrect(fnames, dview=dview, var_name_hdf5=opts.data['var_name_hdf5'], **opts.get_group('motion'))
# note that the file is not loaded in memory
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.
In [13]:
%%capture
#%% Run piecewise-rigid motion correction using NoRMCorre
mc.motion_correct(save_movie=True)
m_els = cm.load(mc.fname_tot_els)
border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0
# maximum shift to be used for trimming against NaNs
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 notebook.
In [14]:
#%% compare with original movie
display_movie = False
if display_movie:
m_orig = cm.load_movie_chain(fnames, var_name_hdf5=opts.data['var_name_hdf5'])
ds_ratio = 0.2
cm.concatenate([m_orig.resize(1, 1, ds_ratio) - mc.min_mov*mc.nonneg_movie,
m_els.resize(1, 1, ds_ratio)],
axis=2).play(fr=60, gain=15, magnification=2, offset=0) # press q to exit
In [15]:
#%% MEMORY MAPPING
# memory map the file in order 'C'
fname_new = cm.save_memmap(mc.mmap_file, base_name='memmap_', order='C',
border_to_0=border_to_0) # exclude borders
# now load the file
Yr, dims, T = cm.load_memmap(fname_new)
images = np.reshape(Yr.T, [T] + list(dims), order='F')
#load frames in python format (T x X x Y)
Now restart the cluster to clean up the memory
In [16]:
#%% restart cluster to clean up memory
cm.stop_server(dview=dview)
c, dview, n_processes = cm.cluster.setup_cluster(
backend='local', n_processes=None, single_thread=False)
In [17]:
%%capture
#%% RUN CNMF ON PATCHES
# First extract spatial and temporal components on patches and combine them
# for this step deconvolution is turned off (p=0)
opts.change_params({'p': 0,'n_processes': n_processes})
cnm = cnmf.CNMF(n_processes, params=opts, dview=dview)
cnm = cnm.fit(images)
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.
In [18]:
# cnm1 = cnmf.CNMF(n_processes, params=opts, dview=dview)
# cnm1.fit_file(motion_correct=True)
Briefly inspect the results by plotting contours of identified components against correlation image.
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 wiki.
In [19]:
#%% plot contours of found components
Cn = cm.local_correlations(images, swap_dim=False)
#Cn = cm.local_correlations(images.transpose(1,2,0))
Cn[np.isnan(Cn)] = 0
cnm.estimates.plot_contours_nb(img=Cn)
cnm.estimates.Cn = Cn
cnm.save(fname_new[:-4]+'hdf5')
In [20]:
%%capture
#%% RE-RUN seeded CNMF on accepted patches to refine and perform deconvolution
cnm.params.change_params({'p': p})
cnm2 = cnm.refit(images, dview=dview)
The processing in patches creates several spurious components. These are filtered out by evaluating each component using three different criteria:
In [21]:
#%% COMPONENT EVALUATION
# the components are evaluated in three ways:
# a) the shape of each component must be correlated with the data
# b) a minimum peak SNR is required over the length of a transient
# c) each shape passes a CNN based classifier
cnm2.estimates.evaluate_components(images, cnm2.params, dview=dview)
cnm2.estimates.Cn = Cn
cnm2.save(fname_new[:-4] + 'hdf5')
Plot contours of selected and rejected components
In [20]:
#%% PLOT COMPONENTS
#cnm2.estimates.plot_contours_nb(img=Cn, idx=cnm2.estimates.idx_components)
Out[20]:
View traces of accepted and rejected components. Note that if you get data rate error you can start Jupyter notebooks using: 'jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10'
In [21]:
# accepted components
#cnm2.estimates.view_components(images, img=Cn,idx=cnm2.estimates.idx_components)
#cnm2.estimates.nb_view_components(img=Cn, idx=cnm2.estimates.idx_components)
Out[21]:
In [22]:
# rejected components
#if len(cnm2.estimates.idx_components_bad) > 0:
# cnm2.estimates.view_components(images, img=Cn,idx=cnm2.estimates.idx_components_bad)
# cnm2.estimates.nb_view_components(img=Cn, idx=cnm2.estimates.idx_components_bad)
#else:
# print("No components were rejected.")
In [22]:
#%% Extract DF/F values
cnm2.estimates.detrend_df_f(quantileMin=8, frames_window=250)
Out[22]:
In [23]:
cnm2.estimates.select_components(use_object=True)
Out[23]:
In [25]:
#cnm2.estimates.view_components(img=Cn)
#cnm2.estimates.nb_view_components(img=Cn, denoised_color='red')
#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')
In [25]:
save_results = True
if save_results:
#cnm2.save('analysis_results.hdf5')
cnm2.estimates.save_NWB(fnames[0], imaging_rate=fr)
In [27]:
#%% STOP CLUSTER and clean up log files
cm.stop_server(dview=dview)
log_files = glob.glob('*_LOG_*')
for log_file in log_files:
os.remove(log_file)
In [28]:
cnm2.estimates.play_movie(images, q_max=99.9, gain_res=2,
magnification=2,
bpx=border_to_0,
include_bck=False)
The denoised movie can also be explicitly constructed using:
In [ ]:
#%% reconstruct denoised movie
denoised = cm.movie(cnm2.estimates.A.dot(cnm2.estimates.C) + \
cnm2.estimates.b.dot(cnm2.estimates.f)).reshape(dims + (-1,), order='F').transpose([2, 0, 1])