New files

CellsortApplyFilter.m

@@ -0,0 +1,77 @@
+function cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+% cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+%
+%CellsortApplyFilter
+% Read in movie data and output signals corresponding to specified spatial
+% filters
+%
+% Inputs:
+%     fn - file name of TIFF movie file
+%     ica_segments - nIC x X matrix of ICA spatial filters
+%     flims - optional two-element vector of frame limits to be read
+%     movm - mean fluorescence image
+%     subtractmean - boolean specifying whether or not to subtract the mean
+%     fluorescence of each time frame
+%
+% Outputs:
+%     cell_sig - cellular signals
+%
+% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+% Email: eran@post.harvard.edu, mschnitz@stanford.edu
+%
+
+if (nargin<3)||isempty(flims)
+    nt = tiff_frames(fn);
+    flims = [1,nt];
+else
+    nt = diff(flims)+1;
+end
+if nargin<5
+    subtractmean = 1;
+end
+
+[pixw,pixh] = size(imread(fn,1));
+if (nargin<4)||isempty(movm)
+    movm = ones(pixw,pixh);
+else
+    movm = double(movm);
+end
+movm(movm==0) = 1; % Just in case there are black areas in the average image
+k=0;
+
+cell_sig = zeros(size(ica_segments,1), nt);
+ica_segments = reshape(ica_segments, [], pixw*pixh);
+
+fprintf('Loading %5g frames for %d ROIs.\n', nt, size(ica_segments,1))
+while k<nt
+    ntcurr = min(500, nt-k);
+    mov = zeros(pixw, pixh, ntcurr);
+    for j=1:ntcurr
+        movcurr = imread(fn, j+k+flims(1)-1);
+        mov(:,:,j) = movcurr;
+    end
+    mov = (mov ./ repmat(movm, [1,1,ntcurr])) - 1; % Normalize by background and subtract mean
+
+    if subtractmean
+        % Subtract the mean of each frame
+        movtm = mean(mean(mov,1),2);
+        mov = mov - repmat(movtm,[pixw,pixh,1]);
+    end
+
+    mov = reshape(mov, pixw*pixh, ntcurr);
+    cell_sig(:, k+[1:ntcurr]) = ica_segments*mov;
+
+    k=k+ntcurr;
+    fprintf('Loaded %3.0f frames; ', k)
+    toc
+end
+
+function j = tiff_frames(fn)
+%
+% n = tiff_frames(filename)
+%
+% Returns the number of slices in a TIFF stack.
+%
+% Modified April 9, 2013 for compatibility with MATLAB 2012b
+
+j = length(imfinfo(fn));

CellsortChoosePCs.m

@@ -0,0 +1,113 @@
+function [PCuse] = CellsortChoosePCs(fn, mixedfilters)
+% [PCuse] = CellsortChoosePCs(fn, mixedfilters)
+%
+% Allows the user to select which principal components will be kept
+% following dimensional reduction.
+%
+% Inputs:
+%   fn - movie file name. Must be in TIFF format.
+%   mixedfilters - N x X matrix of N spatial signal mixtures sampled at X
+%   spatial points.
+%
+% Outputs:
+%   PCuse - vector of indices of the PCs to be kept for dimensional
+%   reduction
+%
+% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+% Email: eran@post.harvard.edu, mschnitz@stanford.edu
+%
+
+fprintf('-------------- CellsortChoosePCs %s -------------- \n', date)
+
+[pixw,pixh] = size(imread(fn,1));
+
+npcs = 20; % Number of PCs to display concurrently
+currpcs = [1:npcs];
+PCf = [];
+while isempty(PCf)
+    showpcs(currpcs, mixedfilters, pixw, pixh)
+    yl = ylim;
+    xl = xlim;
+    set(gca,'Units','pixels')
+    title(['Choose first PC; showing PCs [',num2str(currpcs(1)),':',num2str(currpcs(end)),']'])
+    PCf = input('Number of first PC to retain, Klow (''b/f'' to scroll backwards/forwards)): ','s');
+    if PCf=='b'
+        currpcs = currpcs - min(npcs,currpcs(1)-1);
+        PCf = [];
+    elseif (PCf=='f')
+        currpcs = currpcs+npcs;
+        if nnz(currpcs>size(mixedfilters,2))
+            currpcs = [-npcs+1:0]+size(mixedfilters,2);
+            fprintf('Reached end of stored PCs.\n')
+        end
+        PCf = [];
+    else
+        PCf = str2num(PCf);
+    end
+end
+PCl=[];
+currpcs = [PCf:PCf+npcs-1];
+while isempty(PCl)
+    showpcs(currpcs, mixedfilters, pixw, pixh)
+    title(['Choose last PC; showing PCs [',num2str(currpcs(1)),':',num2str(currpcs(end)),']'])
+    PCl = input('Number of last PC to retain, Khigh (''b/f'' to scroll backwards/forwards): ','s');
+    if PCl=='b'
+        currpcs = currpcs - min(npcs,currpcs(1)-1);
+        PCl = [];
+    elseif (PCl=='f')
+        currpcs = currpcs+npcs;
+        if nnz(currpcs>size(mixedfilters,2))
+            currpcs = [-npcs+1:0]+size(mixedfilters,2);
+            fprintf('Reached end of stored PCs.\n')
+        end
+        PCl = [];
+    else
+        PCl = str2num(PCl);
+    end
+end
+currpcs = [PCf:PCl];
+PCbad=[];
+showpcs(currpcs, mixedfilters, pixw, pixh)
+
+PCuse = setdiff(currpcs, PCbad);
+showpcs(PCuse, mixedfilters, pixw, pixh)
+
+fprintf('  Retaining PCs in the range [Klow - Khigh] = [%d - %d].\n', PCf,PCl)
+
+function showpcs(usepcs, Efull, pixw, pixh)
+
+if nargin<3
+    fprintf('Assuming movie frames are square.\n')
+    pixw = sqrt(size(Efull,1));
+    pixh = sqrt(size(Efull,1));
+end
+if isempty(usepcs)
+    usepcs = [1:size(Efull,2)];
+end
+
+if ndims(Efull)>=3
+    Efull = reshape(Efull, pixw*pixh, []);
+end
+for j=usepcs
+    Efull(:,j) = zscore(Efull(:,j));
+end
+pcs = reshape(Efull(:,usepcs), pixw, pixh, []);
+pcs = permute(pcs, [1, 2, 4, 3]);
+montage(pcs)
+colormap(hot)
+axis on
+xl = xlim;
+yl = ylim;
+nw = ceil(xl(2)/pixh)-1;
+nh = ceil(yl(2)/pixw)-1;
+set(gca,'YTick',[pixw:pixw:yl(2)],'YTickLabel',  num2str(usepcs(min([0:nh]*nw+1, length(usepcs)))'), ...
+    'XTick',[pixh:pixh:xl(2)], ...
+    'XTickLabel',num2str(usepcs([(nh-1)*nw+1:length(usepcs)])'), 'XAxisLocation','bottom','LineWidth',2)
+grid on
+formataxes
+caxis([-1,1]*7)
+
+function formataxes
+
+set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out')
+set(gcf,'Color','w','PaperPositionMode','auto')

CellsortFindspikes.m

@@ -0,0 +1,54 @@
+function [spmat, spt, spc] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+% [spmat, spt, spc, zsig] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+%
+% CELLSORT
+%  Deconvolve signal and find spikes using a threshold
+%
+% Inputs:
+%   ica_sig - nIC x T matrix of ICA temporal signals
+%   thresh - threshold for spike detection
+%   dt - time step
+%   deconvtau - time constant for temporal deconvolution (Butterworth
+%   filter); if deconvtau=0 or [], no deconvolution is performed
+%   normalization - type of normalization to apply to ica_sig; 0 - no
+%   normalization; 1 - divide by standard deviation and subtract mean
+%
+% Outputs:
+%   spmat - nIC x T sparse binary matrix, containing 1 at the time frame of each
+%   spike
+%   spt - list of all spike times
+%   spc - list of the indices of cells for each spike
+%
+% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+% Email: eran@post.harvard.edu, mschnitz@stanford.edu
+%
+
+if size(ica_sig,2)==1
+    ica_sig = ica_sig';
+end
+
+if (nargin>=3)&&(deconvtau>0)
+    dsig = diff(ica_sig,1,2);
+    sig = ica_sig/deconvtau + [dsig(:,1),dsig]/dt;
+else
+    sig = ica_sig;
+end
+
+if (nargin<2)
+    thresh=3;
+    fprintf('Using threshold = 3 s.d. \n')
+end
+switch normalization
+    case 0 % Absolute units
+        zsig = sig';
+    case 1 % Standard-deviation
+        zsig = zscore(sig');
+end
+pp1=[zsig(1,:);zsig(1:end-1,:)];
+pp2=[zsig(2:end,:);zsig(end,:)];
+spmat = sparse((zsig>=thresh)&(zsig-pp1>=0)&(zsig-pp2>=0));
+
+if nargout>1
+    [spt,spc] = find(spmat);
+    spt = spt*dt;
+end