Skip to content

Commit

Permalink
save oasis_matlab as a folder directly, not submodule
Browse files Browse the repository at this point in the history
  • Loading branch information
zhoupc committed Nov 25, 2018
1 parent 561c14c commit 5450956
Show file tree
Hide file tree
Showing 100 changed files with 10,814 additions and 5 deletions.
3 changes: 0 additions & 3 deletions .gitmodules
View file
Original file line number Diff line number Diff line change
@@ -1,3 +0,0 @@
[submodule "OASIS_matlab"]
path = OASIS_matlab
url = https://github.com/zhoupc/OASIS_matlab.git
1 change: 0 additions & 1 deletion OASIS_matlab
Submodule OASIS_matlab deleted from 0a1251
674 changes: 674 additions & 0 deletions OASIS_matlab/subfolder_tmp/LICENSE
View file

Large diffs are not rendered by default.

255 changes: 255 additions & 0 deletions OASIS_matlab/subfolder_tmp/README.html
View file

Large diffs are not rendered by default.

33 changes: 33 additions & 0 deletions OASIS_matlab/subfolder_tmp/README.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
# OASIS: Fast online deconvolution of calcium imaging data

Code accompanying the paper "[Fast online deconvolution of calcium imaging data](http://journals.plos.org/ploscompbiol/article?rev=2&id=10.1371/journal.pcbi.1005423)".

Python: https://github.com/j-friedrich/OASIS

MATLAB: https://github.com/zhoupc/OASIS_matlab

[A brief summary of the FOOPSI approach for denoising & deconvolving calcium imaging](https://github.com/zhoupc/OASIS_matlab/blob/master/document/FOOPSI_.md)
## Get Started
### Installation
add OASIS function to the search path of MATLAB

`>> oasis_setup`

### Process yoru data
There is also a high-level function **deconvolveCa.m** for the ease of calling different methods. You only need to specify parameters and then denoise & deconvolve your raw trace. For example
```matlab
[c, s, options] = deconvolveCa(y, 'foopsi', 'ar1', 'smin', -3, ...'optimize_pars', true, 'optimize_b', true)
```
In this example, we deconvolve the raw trace $y$ using FOOPSI model and constrain the spike size to be $3\times $ noise levels. The AR coefficients and the baseline were updated automatically.

For more options, check the **examples/** folder and see the comments in *deconvolveCa.m*.

### Reproduce figures in the paper
You can reproduce the figures in the paper [3] using the following command (replace * with figure index).

`>> run examples/paper/fig*.m `

## Copyright

Pengcheng Zhou @ Colubmia University
2018
349 changes: 349 additions & 0 deletions OASIS_matlab/subfolder_tmp/deconvolveCa.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,349 @@
function [c, s, options] = deconvolveCa(y, varargin)
%% infer the most likely discretized spike train underlying an fluorescence trace
%% Solves mutliple formulation of the problem
% 1) FOOPSI,
% mininize_{c,s} 1/2 * norm(y-c,2)^2 + lambda * norm(s,1)
% subject to c>=0, s>=0, s=Gc
% 2) constrained FOOPSI
% minimize_{c,s} norm(s, q)
% subject to norm(y-c,2) <= sn*sqrt(T), c>=0, s>=0, s=Gc
% where q is either 1 or 0, rendering the problem convex or non-convex.
% 3) hard threshrinkage
% minimize_{c,s} 1/2 * norm(y-c, 2)^2
% subjec to c>=0, s=Gc, s=0 or s>=smin
% 4) Nonnegative least square problem (NNLS)
% min_{s} norm(y - s*h, 2)^2 + lambda * norm(s,1)
% subject to s>=0

%% inputs:
% y: T x 1 vector, fluorescence trace
% varargin: variable input arguments
% type: string, defines the model of the deconvolution kernel. possible
% options are:
% 'ar1': auto-regressive model with order p=1
% 'ar2': auto-regressive model with order p=2
% 'exp2': the convolution kernel is modeled as the difference of two
% exponential functions -
% h(t) = (exp(-t/tau_d) - exp(-t/tau_r)) / (tau_d-tau_r)
% 'kernel': a vector of the convolution kernel
% pars: parameters for the specified convolution kernel. it has
% different shapes for differrent types of the convolution model:
% 'ar1': scalar
% 'ar2': 2 x 1 vector, [r_1, r_2]
% 'exp2': 2 x 1 vector, [tau_r, tau_d]
% 'kernel': maxISI x 1 vector, the kernel.
% sn: scalar, standard deviation of the noise distribution. If no
% values is give, then sn is estimated from the data based on power
% spectual density method.
% b: fluorescence baseline vlaues. default is 0
% optimize_pars: estimate the parameters of the convolution kernel. default: 0
% optimize_b: estimate the baseline. default: 0
% lambda: penalty parameter
% method: methods for running deconvolution. {'foopsi',
% 'constrained_foopsi' (default), 'thresholded'},

%% outputs:
% c: T x 1 vector, denoised trace
% s: T x 1 vector, deconvolved signal
% b: fluorescence baseline
% kernel: struct variable containing the parameters for the selected
% convolution model
% lambda: Optimal Lagrange multiplier for noise constraint under L1 penalty
% """olves the noise constrained sparse nonnegat

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% ported from the Python implementation from Johannes Friedrich

%% References
% Friedrich J et.al., NIPS 2016, Fast Active Set Method for Online Spike Inference from Calcium Imaging

%% input arguments
y = reshape(y, [], 1); % reshape the trace as a vector
options = parseinputs(varargin{:}); % parse input arguments
if isempty(y)
c = []; s = [];
return;
end
win = options.window; % length of the convolution kernel
% estimate the noise
if isempty(options.sn)
options.sn = GetSn(y);
end
% estimate time constant
if isempty(options.pars) || all(options.pars==0)
switch options.type
case 'ar1'
try
options.pars = estimate_time_constant(y, 1, options.sn);
catch
c = y*0;
s = c;
fprintf('fail to deconvolve the trace\n');
return;
end
if length(options.pars)~=1
c = zeros(size(y));
s = zeros(size(y));
options.pars = 0;
return;
end
case 'ar2'
options.pars = estimate_time_constant(y, 2, options.sn);
if length(options.pars)~=2
c = zeros(size(y));
s = zeros(size(y));
options.pars =[0,0];
return;
end
case 'exp2'
g = estimate_time_constant(y, 2, options.sn);
options.pars = ar2exp(g);
case 'kernel'
g = estimate_time_constant(y, 2, options.sn);
taus = ar2exp(g);
options.pars = exp2kernel(taus, options.win); % convolution kernel
end
end

%% run deconvolution
c = y;
s = y;
switch lower(options.method)
case 'foopsi' %% use FOOPSI
if strcmpi(options.type, 'ar1') % AR 1
if options.smin<0
options.smin = abs(options.smin)*options.sn;
end
gmax = exp(-1/options.max_tau);
[c, s, options.b, options.pars] = foopsi_oasisAR1(y-options.b, options.pars, options.lambda, ...
options.smin, options.optimize_b, options.optimize_pars, [], options.maxIter, ...
options.tau_range, gmax);
elseif strcmpi(options.type, 'ar2') % AR 2
if options.smin<0
options.smin = abs(options.smin)*options.sn/max_ht(options.pars);
end
[c, s, options.b, options.pars] = foopsi_oasisAR2(y-options.b, options.pars, options.lambda, ...
options.smin);
elseif strcmpi(options.type, 'exp2') % difference of two exponential functions
kernel = exp2kernel(options.pars, options.window);
[c, s] = onnls(y-options.b, kernel, options.lambda, ...
options.shift, options.window);
elseif strcmpi(options.type, 'kernel') % convolution kernel itself
[c, s] = onnls(y-options.b, options.pars, options.lambda, ...
options.shift, options.window);
else
disp('to be done');
end
case 'constrained'
if strcmpi(options.type, 'ar1') % AR1
[c, s, options.b, options.pars, options.lambda] = constrained_oasisAR1(y,...
options.pars, options.sn, options.optimize_b, options.optimize_pars, ...
[], options.maxIter, options.tau_range);
else
[cc, options.b, c1, options.pars, options.sn, s] = constrained_foopsi(y,[],[],options.pars,options.sn, ...
options.extra_params);
gd = max(roots([1,-options.pars'])); % decay time constant for initial concentration
gd_vec = gd.^((0:length(y)-1));
c = cc(:) + c1*gd_vec';
options.cin = c1;
end
case 'thresholded' %% Use hard-shrinkage method
if strcmpi(options.type, 'ar1')
[c, s, options.b, options.pars, options.smin] = thresholded_oasisAR1(y,...
options.pars, options.sn, options.optimize_b, options.optimize_pars, ...
[], options.maxIter, options.thresh_factor, options.p_noise, ...
options.tau_range);
% if and(options.smin==0, options.optimize_smin) % smin is given
% [c, s, options.b, options.pars, options.smin] = thresholded_oasisAR1(y,...
% options.pars, options.sn, options.optimize_b, options.optimize_pars, ...
% [], options.maxIter, options.thresh_factor);
% else
% [c, s] = oasisAR1(y-options.b, options.pars, options.lambda, ...
% options.smin);
% end
elseif strcmpi(options.type, 'ar2')
[c, s, options.b, options.pars, options.smin] = thresholded_oasisAR2(y,...
options.pars, options.sn, options.smin, options.optimize_b, options.optimize_pars, ...
[], options.maxIter, options.thresh_factor);
% if and(options.smin==0, options.optimize_smin) % smin is given
% [c, s, options.b, options.pars, options.smin] = thresholded_oasisAR2(y,...
% options.pars, options.sn, options.optimize_b, options.optimize_pars, ...
% [], options.maxIter, options.thresh_factor);
% else
% [c, s] = oasisAR2(y-options.b, options.pars, options.lambda, ...
% options.smin);
% end
elseif strcmpi(options.type, 'exp2') % difference of two exponential functions
d = options.pars(1);
r = options.pars(2);
options.pars = (exp(log(d)*(1:win)) - exp(log(r)*(1:win))) / (d-r); % convolution kernel
[c, s] = onnls(y-options.b, options.pars, options.lambda, ...
options.shift, options.window, [], [], [], options.smin);
elseif strcmpi(options.type, 'kernel') % convolution kernel itself
[c, s] = onnls(y-options.b, options.pars, options.lambda, ...
options.shift, options.window, [], [], [], options.smin);
else
disp('to be done');
end
case 'mcmc'
SAMP = cont_ca_sampler(y,options.extra_params);
options.extra_params = SAMP;
options.mcmc_results = SAMP;
plot_continuous_samples(SAMP,y);
end

% deal with large residual
if options.remove_large_residuals && strcmpi(options.method, 'foopsi')
ind = (abs(smooth(y-c, 3))>options.smin);
c(ind) = max(0, y(ind));
end


function options=parseinputs(varargin)
%% parse input variables

%% default options
options.type = 'ar1';
options.pars = [];
options.sn = [];
options.b = 0;
options.lambda = 0;
options.optimize_b = false;
options.optimize_pars = false;
options.optimize_smin = false;
options.method = 'constrained';
options.window = 200;
options.shift = 100;
options.smin = 0;
options.maxIter = 10;
options.thresh_factor = 1.0;
options.extra_params = [];
options.p_noise = 0.9999;
options.max_tau = 100;
options.tau_range = [];
options.remove_large_residuals = false;

if isempty(varargin)
return;
elseif isstruct(varargin{1}) && ~isempty(varargin{1})
tmp_options = varargin{1};
field_nams = fieldnames(tmp_options);
for m=1:length(field_nams)
eval(sprintf('options.%s=tmp_options.%s;', field_nams{m}, field_nams{m}));
end
k = 2;
else
k = 1;
end
%% parse all input arguments
while k<=nargin
if isempty(varargin{k})
k = k+1;
end
switch lower(varargin{k})
case {'ar1', 'ar2', 'exp2', 'kernel'}
% convolution kernel type
options.type = lower(varargin{k});
if (k<nargin) && (isnumeric(varargin{k+1}))
options.pars = varargin{k+1};
k = k + 1;
end
k = k + 1;
case 'pars'
% parameters for the kernel
options.pars = varargin{k+1};
k = k+2;
case 'sn'
% noise
options.sn = varargin{k+1};
k = k+2;
case 'b'
% baseline
options.b = varargin{k+1};
k = k+2;
case 'optimize_b'
% optimize the baseline
options.optimize_b = true;
if (k<nargin) && (islogical(varargin{k+1}))
options.optimize_b = varargin{k+1};
k = k + 1;
end
k = k+1;
case 'optimize_pars'
% optimize the parameters of the convolution kernel
options.optimize_pars = true;
if (k<nargin) && (islogical(varargin{k+1}))
options.optimize_pars = varargin{k+1};
k = k+1;
end
k = k + 1;

case 'optimize_smin'
% optimize the parameters of the convolution kernel
options.optimize_smin = true;
if (k<nargin) && (islogical(varargin{k+1}))
options.optimize_smin = varargin{k+1};
k = k+1;
end
k = k+1;
case 'lambda'
% penalty
options.lambda = varargin{k+1};
k = k+2;
case {'foopsi', 'constrained', 'thresholded', 'mcmc'}
% method for running the deconvolution
options.method = lower(varargin{k});
k = k+1;
if strcmpi(options.method, 'mcmc') && (k<=length(varargin)) && (~ischar(varargin{k}))
options.extra_params = varargin{k};
k = k+1;
end
case 'window'
% maximum length of the kernel
options.window = varargin{k+1};
k = k+2;
case 'shift'
% number of frames by which to shift window from on run of NNLS
% to the next
options.shift = varargin{k+1};
k = k+2;
case 'smin'
% number of frames by which to shift window from on run of NNLS
% to the next
options.smin = varargin{k+1};
k = k+2;
case 'maxiter'
% number of frames by which to shift window from on run of NNLS
% to the next
options.maxIter = varargin{k+1};
k = k+2;
case 'thresh_factor'
% number of frames by which to shift window from on run of NNLS
% to the next
options.thresh_factor = varargin{k+1};
k = k+2;
case 'p_noise'
% number of frames by which to shift window from on run of NNLS
% to the next
options.p_noise = varargin{k+1};
k = k+2;
case 'tau_range'
options.tau_range = varargin{k+1};
k = k+2;
case 'remove_large_residuals'
% remove large residuals by setting c(t) = y(t)
options.remove_large_residuals = true;
if (k<nargin) && (islogical(varargin{k+1}))
options.remove_large_residuals = varargin{k+1};
k = k+1;
end
k = k+1;
otherwise
k = k+1;
end
end

%% correct some wrong inputs
if strcmpi(options.type, 'kernel')
options.window = numel(options.pars);
end
Loading

0 comments on commit 5450956

Please sign in to comment.