Pipeline
Here we use the provided demo script demo_large_data_1p.m as an example to show the whole pipeline of the analysis. Don't freak out by seeing such a long demo script. 50% of them are just parameter specification. As for the rest of them, they are just single line codes for multiple modules of CNMF-E. As long as you are familiar with the algorithm (Its not necessary, but I highly recommend you to understand the whole big picture), this is pretty easy to understand.
CNMF-E doesn't do motion correction, so please make sure that all your data has been motioned corrected before running CNMF-E. There are two ways of specifying data in CNMF-E.
- set the value of the variable nam with the file path
- set nam-[] and CNMF-E will open a graphical window to select data
%% choose data
neuron = Sources2D();
nam = './data_endoscope.tif'; % this demo data is very small, here we just use it as an example
nam = neuron.select_data(nam); %if nam is [], then select data interactively Here CNMF-E creates a Sources2D class object for saving all results and it can call multiple function directly. See illustrations below.
CNMF-E saves all parameters into a struct variable neuron.options before running the whole analysis. You can modify those values to adapt your data. The detailed explanations to all options were described in https://github.com/zhoupc/CNMF_E/wiki/Parameters.
Please don't be scared by these options. There are less than 10 parameters that are frequently used. In most cases, the default values should be working well. More parameters just give you the flexibility to handling complex data. There is no need to worry about them.
%% parameters
% ------------------------- COMPUTATION ------------------------- %
pars_envs = struct('memory_size_to_use', 8, ... % GB, memory space you allow to use in MATLAB
'memory_size_per_patch', 0.3, ... % GB, space for loading data within one patch
'patch_dims', [64, 64]); %GB, patch size
% ------------------------- SPATIAL ------------------------- %
gSig = 3; % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
gSiz = 11; % pixel, neuron diameter
ssub = 1; % spatial downsampling factor
with_dendrites = false; % with dendrites or not
if with_dendrites
% determine the search locations by dilating the current neuron shapes
updateA_search_method = 'dilate'; %#ok<UNRCH>
updateA_bSiz = 20;
updateA_dist = neuron.options.dist;
else
% determine the search locations by selecting a round area
updateA_search_method = 'ellipse'; %#ok<UNRCH>
updateA_dist = 5;
updateA_bSiz = neuron.options.dist;
end
spatial_constraints = struct('connected', true, 'circular', true); % you can include following constraints: 'circular'
% ------------------------- TEMPORAL ------------------------- %
Fs = 6; % frame rate
tsub = 1; % temporal downsampling factor
deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
'smin', -3, ... % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
'optimize_pars', true, ... % optimize AR coefficients
'optimize_b', true);% optimize the baseline);
nk = 5; % detrending the slow fluctuation. usually 1 is fine (no detrending)
% when changed, try some integers smaller than total_frame/(Fs*30)
detrend_method = 'spline'; % compute the local minimum as an estimation of trend.
% ------------------------- BACKGROUND ------------------------- %
bg_model = 'ring'; % model of the background {'ring', 'svd'(default), 'nmf'}
nb = 1; % number of background sources for each patch (only be used in SVD and NMF model)
bg_neuron_factor = 1.5;
ring_radius = round(bg_neuron_factor * gSiz); % when the ring model used, it is the radius of the ring used in the background model.
%otherwise, it's just the width of the overlapping area
% ------------------------- MERGING ------------------------- %
show_merge = false; % if true, manually verify the merging step
merge_thr = 0.85; % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
method_dist = 'max'; % method for computing neuron distances {'mean', 'max'}
dmin = 5; % minimum distances between two neurons. it is used together with merge_thr
dmin_only = 1; % merge neurons if their distances are smaller than dmin_only.
% ------------------------- INITIALIZATION ------------------------- %
K = []; % maximum number of neurons per patch. when K=[], take as many as possible.
min_corr = 0.8; % minimum local correlation for a seeding pixel
min_pnr = 10; % minimum peak-to-noise ratio for a seeding pixel
min_pixel = 4^2; % minimum number of nonzero pixels for each neuron
bd = 0; % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
frame_range = []; % when [], uses all frames
save_initialization = false; % save the initialization procedure as a video.
use_parallel = true; % use parallel computation for parallel computing
show_init = true; % show initialization results
choose_params = false; % manually choose parameters
center_psf = true; % set the value as true when the background fluctuation is large (usually 1p data)
% set the value as false when the background fluctuation is small (2p)
% ---------------------- WITH MANUAL INTERVENTION -------------------- %
with_manual_intervention = true;
% ------------------------- FINAL RESULTS ------------------------- %
save_demixed = true; % save the demixed file or not
kt = 1; % frame intervals
% ------------------------- UPDATE ALL ------------------------- %
neuron.updateParams('gSig', gSig, ... % -------- spatial --------
'gSiz', gSiz, ...
'ring_radius', ring_radius, ...
'ssub', ssub, ...
'search_method', updateA_search_method, ...
'bSiz', updateA_bSiz, ...
'dist', updateA_bSiz, ...
'spatial_constraints', spatial_constraints, ...
'tsub', tsub, ... % -------- temporal --------
'deconv_options', deconv_options, ... '
'nk', nk, ...
'detrend_method', detrend_method, ...
'background_model', bg_model, ... % -------- background --------
'nb', nb, ...
'ring_radius', ring_radius, ...
'merge_thr', merge_thr, ... % -------- merging ---------
'dmin', dmin, ...
'method_dist', method_dist, ...
'min_corr', min_corr, ... % ----- initialization -----
'min_pnr', min_pnr, ...
'min_pixel', min_pixel, ...
'bd', bd, ...
'center_psf', center_psf);
neuron.Fs = Fs;
This step loads the motion corrected video in different formats (*.tif, *.hdf5, *.avi) and saved it as a mat file. It divides the whole field of view into small patches for efficient loading in the following analysis.
%% distribute data and be ready to run source extraction
neuron.getReady(pars_envs); After doing so, following things happens:
- a folder named as filename_sources_extraction is created
- a mat file was saved into the folder and the file name is saved into neuron.P.mat_data
%% initialize neurons from the video data within a selected temporal range
if choose_params
% change parameters for optimized initialization
[gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
end
[center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
if show_init
h_init=figure();
imagesc(Cn, [0, 1]); colormap gray;
hold on;
plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
endIt will create a new folder within the filename_sources_extraction folder and named as LOGS_current_time. All the results related to this analysis will be put into this folder. If there are some previous analysis, you have the option to load previous results.
The log file within the folder records all steps in the analysis.
%% estimation of background components
neuron.update_background_parallel(use_parallel);
neuron_init = neuron.copy(); [center_res, Cn_res, PNR_res] =neuron.initComponents_residual_parallel([], save_initialization, use_parallel);
if show_init
figure(h_init);
plot(center_res(:, 2), center_res(:, 1), '.g', 'markersize', 10);
end
neuron_init_res = neuron.copy();
neuron.update_temporal_parallel(use_parallel);
%% estimate the spatial components
neuron.update_spatial_parallel(use_parallel);
%% deal with bad neurons
tags = neuron.tag_neurons_parallel(); % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
neuron.delete(tags>0);
%% merge neurons
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_close_neighbors(show_merge, dmin_only); neuron.orderROIs('decay_time'); % order neurons in different ways {'snr', 'decay_time', 'mean', 'circularity'}
if with_manual_intervention
neuron.viewNeurons([], neuron.C_raw);
end
neuron.update_background_parallel(use_parallel);
neuron.update_temporal_parallel(use_parallel);
neuron.update_spatial_parallel(use_parallel);
K = size(neuron.A,2);
tags = neuron.tag_neurons_parallel(); % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
neuron.delete(tags>0);
neuron.merge_neurons_dist_corr(show_merge);
neuron.merge_close_neighbors(show_merge, dmin_only);
if K~=size(neuron.A,2)
neuron.update_temporal_parallel(use_parallel);
neuron.update_spatial_parallel(use_parallel);
neuron.update_background_parallel(use_parallel);
endThis step is very important. It's good for you to save the results and continue the analysis in the fulter
%% save the workspace for future analysis
neuron.save_workspace(); %% show neuron contours
Coor = neuron.show_contours(0.8);
%% create a video for displaying the
neuron.show_demixed_video(save_demixed, kt);
%% save neurons shapes
neuron.save_neurons();