Skip to content

Commit

Permalink
initial commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@jovo
jovo committed Jul 19, 2011
0 parents commit 1d6153d
Show file tree
Hide file tree
Showing 9 changed files with 1,735 additions and 0 deletions.
Binary file added .DS_Store
Binary file not shown.
1 change: 1 addition & 0 deletions README
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
This is a repo containing the most current code for doing model-based spike train inference from calcium imaging. Manuscripts explaining the theory and providing some results on simulated and real data are available from the fast-oopsi, smc-oopsi, and pop-oopsi github repositories. Those repositories also contain code that you may run and data to download to play with things. Any question, ask them in the issues tab. Please let me know of any positive or negative experiences. Much Obliged, joshua
410 changes: 410 additions & 0 deletions fast_oopsi.m

Large diffs are not rendered by default.

237 changes: 237 additions & 0 deletions run_oopsi.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,237 @@
function varargout = run_oopsi(F,V,P)
% this function runs our various oopsi filters, saves the results, and
% plots the inferred spike trains. make sure that fast-oopsi and
% smc-oopsi repository are in your path if you intend to use them.
%
% to use the code, simply provide F, a vector of fluorescence observations,
% for each cell. the fast-oopsi code can handle a matrix input,
% corresponding to a set of pixels, for each time bin. smc-oopsi expects a
% 1D fluorescence trace.
%
% see documentation for fast-oopsi and smc-oopsi to determine how to set
% variables
%
% input
% F: fluorescence from a single neuron
% V: Variables necessary to define for code to function properly (optional)
% P: Parameters of the model (optional)
%
% possible outputs
% fast: fast-oopsi MAP estimate of spike train, argmax_{n\geq 0} P[n|F], (fast.n), parameter estimate (fast.P), and structure of variables for algorithm (fast.V)
% smc: smc-oopsi estimate of {P[X_t|F]}_{t<T}, where X={n,C} or {n,C,h}, (smc.E), parameter estimates (smc.P), and structure of variables for algorithm (fast.V)

%% set code Variables

if nargin < 2, V = struct; end % create structure for algorithmic variables, if none provided
if ~isfield(V,'fast_iter_max')
V.fast_iter_max = input('\nhow many iterations of fast-oopsi would you like to do [0,1,2,...]: ');
end
if ~isfield(V,'smc_iter_max')
V.smc_iter_max = input('\how many iterations of smc-oopsi would you like to do [0,1,2,...]: ');
end

if ~isfield(V,'dt'), % frame duration
fr = input('\nwhat was the frame rate for this movie (in Hz)?: ');
V.dt = 1/fr;
end

if ~isfield(V,'preprocess'),
V.preprocess = input('\ndo you want to high-pass filter [0=no, 1=yes]?: ');
end

if ~isfield(V,'n_max'), V.n_max = 2; end
if nargin < 3, P = struct; end % create structure for parameters, if none provided
if ~isfield(V,'plot'), V.plot = 1; end % whether to plot the fluorescence and spike trains
if ~isfield(V,'name'), % give data a unique, time-stamped name, if there is not one specified
lic = str2num(license); % jovo's license #
if lic == 273165, % if using jovo's computer, set data and fig folders
fdat = '~/Research/oopsi/meta-oopsi/data/jovo';
ffig = '~/Research/oopsi/meta-oopsi/figs/jovo';
else % else just use current dir
fdat = pwd;
ffig = pwd;
end
V.name = ['/oopsi_' datestr(clock,30)];
else
fdat = pwd;
ffig = pwd;
end

if ~isfield(V,'save'), V.save = 0; end % whether to save results and figs
if V.save == 1
if isfield(V,'dat_dir'), fdat=V.dat_dir; end
V.name_dat = [fdat V.name]; % filename for data
save(V.name_dat,'V')
end

%% preprocess - remove the lowest 10 frequencies

if V.preprocess==1
V.T = length(F);
f = detrend(F);
nfft = 2^nextpow2(V.T);
y = fft(f,nfft);
bw = 3;
y(1:bw) = 0; y(end-bw+1:end)=0;
iy = ifft(y,nfft);
F = z1(real(iy(1:V.T)));
end

%% infer spikes and estimate parameters

% infer spikes using fast-oopsi
if V.fast_iter_max > 0
fprintf('\nfast-oopsi\n')
[fast.n fast.P fast.V]= fast_oopsi(F,V,P);
if V.save, save(V.name_dat,'fast','-append'); end
end

stupid=1;
% infer spikes using smc-oopsi
if V.smc_iter_max > 0
fprintf('\nsmc-oopsi\n')
siz=size(F); if siz(1)>1, F=F'; end
if V.fast_iter_max > 0;
if ~isfield(P,'A'), P.A = 50; end % initialize jump in [Ca++] after spike
if ~isfield(P,'n'), P.n = 1; end % Hill coefficient
if ~isfield(P,'k_d'), P.k_d = 200; end % dissociation constant
if ~isfield(V,'T'), V.T = fast.V.T; end % number of time steps
if ~isfield(V,'dt'), V.dt = fast.V.dt; end % frame duration, aka, 1/(framte rate)

if ~exist('stupid')
bign1=find(fast.n>0.1);
bign0=bign1-1;
df=max((F(bign1)-F(bign0))./(F(bign0)));

P.C_init= P.C_0;
S0 = Hill_v1(P,P.C_0);
arg = S0 + df*(S0 + 1/13);
P.A = ((arg*P.k_d)./(1-arg)).^(1/P.n)-P.C_0;
end
P.C_0 = 0;
P.tau_c = fast.V.dt/(1-fast.P.gam); % time constant
nnorm = V.n_max*fast.n/max(fast.n); % normalize inferred spike train
C = filter(1,[1 -fast.P.gam],P.A*nnorm)'+P.C_0; % calcium concentration
C1 = [Hill_v1(P,C); ones(1,V.T)]; % for brevity
ab = C1'\F'; % estimate scalse and offset
P.alpha = ab(1); % fluorescence scale
P.beta = ab(2); % fluorescence offset
P.zeta = (mad(F-ab'*C1,1)*1.4785)^2;
P.gamma = P.zeta/5; % signal dependent noise
P.k = V.spikegen.EFGinv(0.01, P, V);

end
[smc.E smc.P smc.V] = smc_oopsi(F,V,P);
if V.save, save(V.name_dat,'smc','-append'); end
end

%% provide outputs for later analysis

if nargout == 1
if V.fast_iter_max > 0
varargout{1} = fast;
else
varargout{1} = smc;
end
elseif nargout == 2
if V.fast_iter_max>0 && V.smc_iter_max>0
varargout{1} = fast;
varargout{2} = smc;
else
if V.fast_iter_max>0
varargout{1} = fast;
varargout{2} = V;
else
varargout{1} = smc;
varargout{2} = V;
end
end
elseif nargout == 3
varargout{1} = fast;
varargout{2} = smc;
varargout{3} = V;
end

%% plot results

if V.plot
if isfield(V,'fig_dir'), ffig=V.fig_dir; end
V.name_fig = [ffig V.name]; % filename for figure
fig = figure(3);
clf,
V.T = length(F);
nrows = 1;
if V.fast_iter_max>0, nrows=nrows+1; end
if V.smc_iter_max>0, nrows=nrows+1; end
gray = [.75 .75 .75]; % define gray color
inter = 'tex'; % interpreter for axis labels
xlims = [1 V.T]; % xmin and xmax for current plot
fs = 12; % font size
ms = 5; % marker size for real spike
sw = 2; % spike width
lw = 1; % line width
xticks = 0:1/V.dt:V.T; % XTick positions
skip = round(length(xticks)/5);
xticks = xticks(1:skip:end);
tvec_o = xlims(1):xlims(2); % only plot stuff within xlims
if isfield(V,'true_n'), V.n=V.true_n; end
if isfield(V,'n'), spt=find(V.n); end

% plot fluorescence data
i=1; h(i)=subplot(nrows,1,i); hold on
plot(tvec_o,z1(F(tvec_o)),'-k','LineWidth',lw);
ylab=ylabel([{'Fluorescence'}],'Interpreter',inter,'FontSize',fs);
set(ylab,'Rotation',0,'HorizontalAlignment','right','verticalalignment','middle')
set(gca,'YTick',[],'YTickLabel',[])
set(gca,'XTick',xticks,'XTickLabel',[])
axis([xlims 0 1.1])

% plot fast-oopsi output
if V.fast_iter_max>0
i=i+1; h(i)=subplot(nrows,1,i); hold on,
n_fast=fast.n/max(fast.n);
spts=find(n_fast>1e-3);
stem(spts,n_fast(spts),'Marker','none','LineWidth',sw,'Color','k')
if isfield(V,'n'),
stem(spt,V.n(spt)/max(V.n(spt))+0.1,'Marker','v','MarkerSize',ms,'LineStyle','none','MarkerFaceColor',gray,'MarkerEdgeColor',gray)
end
axis([xlims 0 1.1])
hold off,
ylab=ylabel([{'fast'}; {'filter'}],'Interpreter',inter,'FontSize',fs);
set(ylab,'Rotation',0,'HorizontalAlignment','right','verticalalignment','middle')
set(gca,'YTick',0:2,'YTickLabel',[])
set(gca,'XTick',xticks,'XTickLabel',[])
box off
end

% plot smc-oopsi output
if V.smc_iter_max>0
i=i+1; h(i)=subplot(nrows,1,i); hold on,
spts=find(smc.E.nbar>1e-3);
stem(spts,smc.E.nbar(spts),'Marker','none','LineWidth',sw,'Color','k')
if isfield(V,'n'),
stem(spt,V.n(spt)+0.1,'Marker','v','MarkerSize',ms,'LineStyle','none','MarkerFaceColor',gray,'MarkerEdgeColor',gray)
end
axis([xlims 0 1.1])
hold off,
ylab=ylabel([{'smc'}; {'filter'}],'Interpreter',inter,'FontSize',fs);
set(ylab,'Rotation',0,'HorizontalAlignment','right','verticalalignment','middle')
set(gca,'YTick',0:2,'YTickLabel',[])
set(gca,'XTick',xticks,'XTickLabel',[])
box off
end

% label last subplot
set(gca,'XTick',xticks,'XTickLabel',round(xticks*V.dt*100)/100)
xlabel('Time (sec)','FontSize',fs)
linkaxes(h,'x')

% print fig
if V.save
wh=[7 3]; %width and height
set(gcf,'PaperSize',wh,'PaperPosition',[0 0 wh],'Color','w');
print('-depsc',V.name_fig)
print('-dpdf',V.name_fig)
saveas(fig,V.name_fig)
end
end
117 changes: 117 additions & 0 deletions smc_oopsi.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
function [M P V] = smc_oopsi(F,V,P)
% this function runs the SMC-EM on a fluorescence time-series, and outputs the inferred
% distributions and parameter estimates
%
% Inputs
% F: fluorescence time series
% V: structure of stuff necessary to run smc-em code
% P: structure of initial parameter estimates
%
% Outputs
% M: structure containing mean, variance, and percentiles of inferred distributions
% P: structure containing the final parameter estimates
% V: structure Variables for algorithm to run

if nargin < 2, V = struct; end
if ~isfield(V,'T'), V.T = length(F); end % # of observations
if ~isfield(V,'freq'), V.freq = 1; end % # time steps between observations
if ~isfield(V,'T_o'), V.T_o = V.T; end % # of observations
if ~isfield(V,'x'), V.x = ones(1,V.T); end % stimulus
if ~isfield(V,'scan'), V.scan = 0; end % epi or scan
if ~isfield(V,'name'), V.name ='oopsi'; end % name for output and figure
if ~isfield(V,'Nparticles'), V.Nparticles = 99; end % # particles
if ~isfield(V,'Nspikehist'), V.Nspikehist = 0; end % # of spike history terms
if ~isfield(V,'CaBaselineDrift'), V.CaBaselineDrift = 0; end %whether to include baseline drift in resting calcium concentration
if V.CaBaselineDrift && V.freq > 1
warning('baseline drift is not supported for intermittent sampling, using fixed baseline instead');
V.CaBaselineDrift = 0;
end
if ~isfield(V,'condsamp'), V.condsamp = 1; end % whether to use conditional sampler

if ~isfield(V,'ignorelik'), V.ignorelik = 1; end % epi or scan
if ~isfield(V,'true_n'), % if true spikes are not available
V.use_true_n = 0; % don't use them for anything
else
V.use_true_n = 1;
V.true_n = repmat(V.true_n',V.Nparticles,1);
end
if ~isfield(V,'smc_iter_max'), % max # of iterations before convergence
reply = str2double(input('\nhow many EM iterations would you like to perform \nto estimate parameters (0 means use default parameters): ', 's'));
V.smc_iter_max = reply;
end
if ~isfield(V,'dt'),
fr = input('what was the frame rate for this movie (in Hz)? ');
V.dt = 1/fr;
end

% set which parameters to estimate
if ~isfield(V,'est_c'), V.est_c = 1; end % tau_c, A, C_0
if ~isfield(V,'est_t'), V.est_t = 1; end % tau_c (useful when data is poor)
if ~isfield(V,'est_n'), V.est_n = 1; end % b,k
if ~isfield(V,'est_h'), V.est_h = 0; end % w
if ~isfield(V,'est_F'), V.est_F = 1; end % alpha, beta
if ~isfield(V,'smc_plot'), V.smc_plot = 1; end % plot results with each iteration

%% initialize model Parameters

if nargin < 3, P = struct; end
if ~isfield(P,'tau_c'), P.tau_c = 1; end % calcium decay time constant (sec)
if ~isfield(P,'A'), P.A = 50; end % change ins [Ca++] after a spike (\mu M)
if ~isfield(P,'C_0'), P.C_0 = 30; end % baseline [Ca++] (\mu M)
if ~isfield(P,'C_init'),P.C_init= 20; end % initial [Ca++] (\mu M)
if ~isfield(P,'sigma_c'),P.sigma_c= 10; end % standard deviation of noise (\mu M)
if ~isfield(P,'sigma_cr'),P.sigma_cr= 10; end % standard deviation of noise (\mu M)
if ~isfield(P,'n'), P.n = 1; end % hill equation exponent
if ~isfield(P,'k_d'), P.k_d = 200; end % hill coefficient

if ~isfield(V,'spikegen')
[sg,V] = get_spikegen_info('Bernoulli', P, V);
V.spikegen = sg;
end
if V.freq > 1 && ~strcmpi(V.spikegen.name, 'bernoulli')
warning('only binary spiking is supported with intermittent sampling, switching to bernoulli spike count distributions');
[sg,V] = get_spikegen_info('Bernoulli', P, V);
V.spikegen = sg;
end

if ~isfield(V, 'maxspikes'), V.maxspikes = 15; end
if ~isfield(P,'k'), % linear filter
k = str2double(input('approx. how many spikes underly this trace: ', 's'));
P.k = V.spikegen.EFGinv(k / V.T);
%P.k = log(-log(1-k/V.T)/V.dt);
end

if ~isfield(P,'alpha'), P.alpha = mean(F); end % scale of F
if ~isfield(P,'beta'), P.beta = min(F); end % offset of F
if ~isfield(P,'zeta'), P.zeta = P.alpha/5; end % constant variance
if ~isfield(P,'gamma'), P.gamma = P.zeta/5; end % scaled variance
if V.Nspikehist==1 % if there are spike history terms
if ~isfield(P,'omega'), P.omega = -1; end % weight
if ~isfield(P,'tau_h'), P.tau_h = 0.02; end % time constant
if ~isfield(P,'sigma_h'), P.sigma_h = 0; end % stan dev of noise
if ~isfield(P,'g'), P.g = V.dt/P.tau_h; end % for brevity
if ~isfield(P,'sig2_h'), P.sig2_h = P.sigma_h^2*V.dt; end % for brevity
end
if ~isfield(P,'a'), P.a = V.dt/P.tau_c; end % for brevity
if ~isfield(P,'sig2_c'),P.sig2_c= P.sigma_c^2*V.dt; end % for brevity
if ~isfield(P,'sig2_cr'),P.sig2_cr= P.sigma_cr^2*V.dt; end % for brevity

%% initialize other stuff
starttime = cputime;
P.lik = -inf; % we are trying to maximize the likelihood here
F = max(F,eps); % in case there are any zeros in the F time series
% P.k = [P.k; 1];
S = smc_oopsi_forward(F,V,P); % forward step
M = smc_oopsi_backward(S,V,P); % backward step
if V.smc_iter_max>1, P.conv=false; else P.conv=true; end

while P.conv==false;
P = smc_oopsi_m_step(V,S,M,P,F); % m step
S = smc_oopsi_forward(F,V,P); % forward step
M = smc_oopsi_backward(S,V,P); % backward step
end
fprintf('\n')

V.smc_iter_tot = length(P.lik);
V.smc_time = cputime-starttime;
V = orderfields(V);
Loading

0 comments on commit 1d6153d

Please sign in to comment.