diff --git a/CellsortApplyFilter.m b/CellsortApplyFilter.m new file mode 100644 index 0000000..5e44e1b --- /dev/null +++ b/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 ksize(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') diff --git a/CellsortFindspikes.m b/CellsortFindspikes.m new file mode 100644 index 0000000..78e88c3 --- /dev/null +++ b/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 diff --git a/CellsortICA.m b/CellsortICA.m new file mode 100644 index 0000000..872aad6 --- /dev/null +++ b/CellsortICA.m @@ -0,0 +1,155 @@ +function [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, ... + mixedfilters, CovEvals, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds) +% [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, CovEvals, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds) +% +%CELLSORT +% Perform ICA with a standard set of parameters, including skewness as the +% objective function +% +% Inputs: +% mixedsig - N x T matrix of N temporal signal mixtures sampled at T +% points. +% mixedfilters - N x X x Y array of N spatial signal mixtures sampled at +% X x Y spatial points. +% CovEvals - eigenvalues of the covariance matrix +% PCuse - vector of indices of the components to be included. If empty, +% use all the components +% mu - parameter (between 0 and 1) specifying weight of temporal +% information in spatio-temporal ICA +% nIC - number of ICs to derive +% termtol - termination tolerance; fractional change in output at which +% to end iteration of the fixed point algorithm. +% maxrounds - maximum number of rounds of iterations +% +% Outputs: +% ica_sig - nIC x T matrix of ICA temporal signals +% ica_filters - nIC x X x Y array of ICA spatial filters +% ica_A - nIC x N orthogonal unmixing matrix to convert the input to output signals +% numiter - number of rounds of iteration before termination +% +% Routine is based on the fastICA package (Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo +% Hyvärinen, http://www.cis.hut.fi/projects/ica/fastica) +% +% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 +% Email: eran@post.harvard.edu, mschnitz@stanford.edu +% + +fprintf('-------------- CellsortICA %s -------------- \n', date) + +if (nargin<4) || isempty(PCuse) + PCuse = [1:size(mixedsig,1)]; +end +if (nargin<6) || isempty(nIC) + nIC = length(PCuse); +end +if (nargin<7) || isempty(ica_A_guess) + ica_A_guess = randn(length(PCuse), nIC); +end +if (nargin<8) || isempty(termtol) + termtol = 1e-6; +end +if (nargin<9) || isempty(maxrounds) + maxrounds = 100; +end +if isempty(mu)||(mu>1)||(mu<0) + error('Spatio-temporal parameter, mu, must be between 0 and 1.') +end + +% Check that ica_A_guess is the right size +if size(ica_A_guess,1)~= length(PCuse) || size(ica_A_guess,2)~=nIC + error('Initial guess for ica_A is the wrong size.') +end +if nIC>length(PCuse) + error('Cannot estimate more ICs than the number of PCs.') +end + +[pixw,pixh] = size(mixedfilters(:,:,1)); +npix = pixw*pixh; + +% Select PCs +if mu > 0 || ~isempty(mixedsig) + mixedsig = mixedsig(PCuse,:); +end +if mu < 1 || ~isempty(mixedfilters) + mixedfilters = reshape(mixedfilters(:,:,PCuse),npix,length(PCuse)); +end +CovEvals = CovEvals(PCuse); + +% Center the data by removing the mean of each PC +mixedmean = mean(mixedsig,2); +mixedsig = mixedsig - mixedmean * ones(1, size(mixedsig,2)); + +% Create concatenated data for spatio-temporal ICA +nx = size(mixedfilters,1); +nt = size(mixedsig,2); +if mu == 1 + % Pure temporal ICA + sig_use = mixedsig; +elseif mu == 0 + % Pure spatial ICA + sig_use = mixedfilters'; +else + % Spatial-temporal ICA + sig_use = [(1-mu)*mixedfilters', mu*mixedsig]; + sig_use = sig_use / sqrt(1-2*mu+2*mu^2); % This normalization ensures that, if both mixedfilters and mixedsig have unit covariance, then so will sig_use +end + +% Perform ICA +[ica_A, numiter] = fpica_standardica(sig_use, nIC, ica_A_guess, termtol, maxrounds); + +% Sort ICs according to skewness of the temporal component +ica_W = ica_A'; + +ica_sig = ica_W * mixedsig; +ica_filters = reshape((mixedfilters*diag(CovEvals.^(-1/2))*ica_A)', nIC, nx); % This is the matrix of the generators of the ICs +ica_filters = ica_filters / npix^2; + +icskew = skewness(ica_sig'); +[icskew, ICord] = sort(icskew, 'descend'); +ica_A = ica_A(:,ICord); +ica_sig = ica_sig(ICord,:); +ica_filters = ica_filters(ICord,:); +ica_filters = reshape(ica_filters, nIC, pixw, pixh); + +% Note that with these definitions of ica_filters and ica_sig, we can decompose +% the sphered and original movie data matrices as: +% mov_sphere ~ mixedfilters * mixedsig = ica_filters * ica_sig = (mixedfilters*ica_A') * (ica_A*mixedsig), +% mov ~ mixedfilters * pca_D * mixedsig. +% This gives: +% ica_filters = mixedfilters * ica_A' = mov * mixedsig' * inv(diag(pca_D.^(1/2)) * ica_A' +% ica_sig = ica_A * mixedsig = ica_A * inv(diag(pca_D.^(1/2))) * mixedfilters' * mov + + function [B, iternum] = fpica_standardica(X, nIC, ica_A_guess, termtol, maxrounds) + + numSamples = size(X,2); + + B = ica_A_guess; + BOld = zeros(size(B)); + + iternum = 0; + minAbsCos = 0; + + errvec = zeros(maxrounds,1); + while (iternum < maxrounds) && ((1 - minAbsCos)>termtol) + iternum = iternum + 1; + + % Symmetric orthogonalization. + B = (X * ((X' * B) .^ 2)) / numSamples; + B = B * real(inv(B' * B)^(1/2)); + + % Test for termination condition. + minAbsCos = min(abs(diag(B' * BOld))); + + BOld = B; + errvec(iternum) = (1 - minAbsCos); + end + + if iternum1)||(mu<0) + error('Spatio-temporal parameter, mu, must be between 0 and 1.') +end + +% Check that ica_A_guess is the right size +if size(ica_A_guess,1)~= length(PCuse) || size(ica_A_guess,2)~=nIC + error('Initial guess for ica_A is the wrong size.') +end +if nIC>length(PCuse) + error('Cannot estimate more ICs than the number of PCs.') +end + +[pixw,pixh] = size(mixedfilters(:,:,1)); +npix = pixw*pixh; + +% Select PCs +if mu > 0 || ~isempty(mixedsig) + mixedsig = mixedsig(PCuse,:); +end +if mu < 1 || ~isempty(mixedfilters) + mixedfilters = reshape(mixedfilters(:,:,PCuse),npix,length(PCuse)); +end +CovEvals = CovEvals(PCuse); + +% Center the data by removing the mean of each PC +mixedmean = mean(mixedsig,2); +mixedsig = mixedsig - mixedmean * ones(1, size(mixedsig,2)); + +% Create concatenated data for spatio-temporal ICA +nx = size(mixedfilters,1); +nt = size(mixedsig,2); +if mu == 1 + % Pure temporal ICA + sig_use = mixedsig; +elseif mu == 0 + % Pure spatial ICA + sig_use = mixedfilters'; +else + % Spatial-temporal ICA + sig_use = [(1-mu)*mixedfilters', mu*mixedsig]; + sig_use = sig_use / sqrt(1-2*mu+2*mu^2); % This normalization ensures that, if both mixedfilters and mixedsig have unit covariance, then so will sig_use +end + +% Calculate a 3-tensor containing all the moments of the data we will need +SkewMoms = zeros(length(PCuse),length(PCuse),length(PCuse)); +for i1=1:length(PCuse) + for i2=i1:length(PCuse) + for i3=i2:length(PCuse) + SkewMoms(i1,i2,i3) = sum(sig_use(i1,:).*sig_use(i2,:).*sig_use(i3,:)); + SkewMoms(i3,i1,i2) = SkewMoms(i1,i2,i3); + SkewMoms(i2,i3,i1) = SkewMoms(i1,i2,i3); + SkewMoms(i1,i3,i2) = SkewMoms(i1,i2,i3); + SkewMoms(i2,i1,i3) = SkewMoms(i1,i2,i3); + SkewMoms(i3,i2,i1) = SkewMoms(i1,i2,i3); + end + end +end +clear sig_use + +% Perform ICA +[ica_A, numiter] = fpica_standardica(SkewMoms, nIC, ica_A_guess, termtol, maxrounds); + +% Sort ICs according to skewness of the temporal component +ica_W = ica_A'; + +ica_sig = ica_W * mixedsig; +ica_filters = reshape((mixedfilters*diag(CovEvals.^(-1/2))*ica_A)', nIC, nx); % This is the matrix of the generators of the ICs +ica_filters = ica_filters / npix^2; + +icskew = skewness(ica_sig'); +[icskew, ICord] = sort(icskew, 'descend'); +ica_A = ica_A(:,ICord); +ica_sig = ica_sig(ICord,:); +ica_filters = ica_filters(ICord,:); +ica_filters = reshape(ica_filters, nIC, pixw, pixh); + +% Note that with these definitions of ica_filters and ica_sig, we can decompose +% the sphered and original movie data matrices as: +% mov_sphere ~ mixedfilters * mixedsig = ica_filters * ica_sig = (mixedfilters*ica_A') * (ica_A*mixedsig), +% mov ~ mixedfilters * pca_D * mixedsig. +% This gives: +% ica_filters = mixedfilters * ica_A' = mov * mixedsig' * inv(diag(pca_D.^(1/2)) * ica_A' +% ica_sig = ica_A * mixedsig = ica_A * inv(diag(pca_D.^(1/2))) * mixedfilters' * mov + + function [B, iternum] = fpica_standardica(SkewMoms, nIC, ica_A_guess, termtol, maxrounds) + + B = ica_A_guess; + BOld = zeros(size(B)); + nPC = size(B, 1); + + iternum = 0; + minAbsCos = 0; + + errvec = zeros(maxrounds,1); + while (iternum < maxrounds) && ((1 - minAbsCos)>termtol) + iternum = iternum + 1; + + % Symmetric orthogonalization. + % B = (X * ((X' * B) .^ 2)) / numSamples; + Bnew = B; + for i=1:nPC + Bnew(i,:) = trace(B' * SkewMoms(:,:,i) * B); + end + B = B * real(inv(B' * B)^(1/2)); + + % Test for termination condition. + minAbsCos = min(abs(diag(B' * BOld))); + + BOld = B; + errvec(iternum) = (1 - minAbsCos); + end + + if iternum 1) +% spc - vector of spike indices (which cell) (needed for plottype > 1) +% +% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 +% Email: eran@post.harvard.edu, mschnitz@stanford.edu +% + +colord=[ 0 0 1.0000 + 0 0.4000 0 + 1.0000 0 0 + 0 0.7500 0.7500 + 0.7500 0 0.7500 + 0.8, 0.5, 0 + 0 0 0.5 + 0 0.85 0]; + +% Check input arguments +nIC = size(ica_sig,1); +if nargin<5 || isempty(tlims) + tlims = [0, size(ica_sig,2)*dt]; % seconds +end +if nargin<8 || isempty(plottype) + plottype = 1; +end +if nargin<9 || isempty(ICuse) + ICuse = [1:nIC]; +end +if size(ICuse,2)==1 + ICuse = ICuse'; +end + +% Reshape the filters +[pixw,pixh] = size(f0); +if size(ica_filters,1)==nIC + ica_filters = reshape(ica_filters,nIC,pixw,pixh); +elseif size(ica_filters,2)==nIC + ica_filters = reshape(ica_filters,nIC,pixw,pixh); +end + +switch mode + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'series'} + colmax = 20; % Maximum # of ICs in one column + ncols = ceil(length(ICuse)/colmax); + if plottype==4 + ncols=1; + end + nrows = ceil(length(ICuse)/ncols); + + if size(ica_filters(:,:,1))==size(f0(:,:,1)) + ica_filters = permute(ica_filters,[3,1,2]); + end + + subplot(1,3*ncols,[2:3]) + tlims(2) = min(tlims(2),size(ica_sig,2)*dt); + tlims(1) = max(tlims(1),0); + + clf + f_pos = get(gcf,'Position'); + f_pos(4) = max([f_pos(4),500,50*nrows]); + f_pos(3) = max(400*ncols,0.9*f_pos(4)); + + colormap(hot) + colord=get(gca,'ColorOrder'); + ll=0; + filtax = []; + if ~isempty(ica_filters) + for k=0:ncols-1 + jj=3*k; + nrows_curr = min(nrows,length(ICuse)-k*nrows); + for j=1:nrows_curr + filtax= [filtax,subplot(nrows_curr + (plottype==4),3*ncols, jj+1)]; + jj=jj+3*ncols; + ll=ll+1; + imagesc(squeeze(ica_filters(ICuse(ll),:,:))) + axis image tight off + end + end + end + + ax = []; + for j=0:ncols-1 + ax(j+1)=subplot(1,3*ncols,3*j+[2:3]); + ICuseuse = ICuse([1+j*nrows:min(length(ICuse),(j+1)*nrows)]); + if plottype<=2 + complot(ica_sig, ICuseuse, dt) + end + formataxes + if plottype>=2 + [spcICuse,spcord] = ismember(spc,ICuseuse); + spc_curr = spcord(spcICuse&(spt<=tlims(2))&(spt>=tlims(1))); + spt_curr = spt(spcICuse&(spt<=tlims(2))&(spt>=tlims(1))); + alpha = diff(ylim)/length(ICuseuse); + switch plottype + case 2 + hold on + scatter(spt_curr,-alpha*(spc_curr-1)+2,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled') + hold off + case 3 + cla + scatter(spt_curr,-alpha*(spc_curr-1)+0.4*alpha,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled') + yticks=sort(-alpha*([1:length(ICuseuse)]-1)+0.4*alpha); + yl = sort(-alpha*[length(ICuseuse)-1,-1]); + set(gca,'YTick',yticks) + set(gca,'YTicklabel',num2str(fliplr(ICuseuse)'),'YLim',yl) + set(gca,'YLim',sort(-alpha*[length(ICuseuse)-1,-0.6])) + + case 4 + ax4=[]; + tvec = [ratebin/2:ratebin:size(ica_sig,2)*dt]; + binsize = ratebin*ones(size(tvec)); + binsize(end) = size(ica_sig,2)*dt-tvec(end-1)-ratebin/2; + fprintf('IC mean rate\n') + sprm = []; sprs = []; + for jj=1:length(ICuseuse) + sptj = spt_curr(spc_curr==jj); + spn = hist(sptj,tvec); + spr = spn./binsize; + sperr = sqrt(spn)./binsize; + + ax4(jj) = subplot(length(ICuseuse)+1,3,[3*jj-1,3*jj]); + errorbar(tvec,spr,sperr,'Color',colord(mod(jj-1,size(colord,1))+1,:),'LineWidth',1) + hold on + bar(tvec,spr) + hold off + axis tight + formataxes + ylabel({['IC',int2str(ICuseuse(jj))],[num2str(mean(spr),2),'']}) + sprm(jj) = sum(spn)/sum(binsize); + fprintf('%3.3f\n', sprm(jj)) + box off + + end + histax = subplot(length(ICuseuse)+1,3, 3*length(ICuseuse)+1); + hist(sprm,[0:0.1:1]) + formataxes + xlabel('Spike rate (Hz)') + ylabel('# Cells') + xlim([0,1]) + + spn = hist(spt_curr,tvec); + spr = spn./(binsize*length(ICuseuse)); + sperr = sqrt(spn)./(binsize*length(ICuseuse)); %Standard error for measurement of mean from Poisson process + + ax4(jj+1) = subplot(length(ICuseuse)+1,3,3*length(ICuseuse)+[2,3]); + errorbar(tvec,spr,sperr,'k','LineWidth',3) + formataxes + set(ax4,'XLim',tlims) + set(ax4(1:jj-1),'XTick',[]) + ylabel({'Pop. mean',[num2str(mean(spr),2),' Hz']},'FontAngle','i') + fprintf('Pop. ave\tPop. SD\n') + fprintf('%3.3f\t\t%3.3f\n', mean(sprm), std(sprm)) + xlabel('Time (s)','FontAngle','i') + box off + ax(j+1) = ax4(1); + case 5 + % Scatter plot of spikes, showing synchrony + cla + [nsp, nsp_t] = hist(spt_curr, unique(spt_curr)); + nsp_t = nsp_t(nsp>1); + nsp = nsp(nsp>1); + for jj = 1:length(nsp_t) + plot( nsp_t(jj)*ones(nsp(jj),1), spc_curr(spt_curr==nsp_t(jj)), 'o-', ... + 'Color', colord(mod(nsp(jj)-2,size(colord,1))+1,:), 'LineWidth', 2, 'MarkerSize',10) + hold on + end + scatter(spt_curr,-spc_curr,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled') + hold off + yticks=sort([1:max(spc_curr)]); + yl = sort([max(spc_curr)+1,min(spc_curr)-1]); + set(gca,'YTick',yticks) + set(gca,'YTicklabel',num2str(fliplr(ICuse)'),'YLim',yl) + end + end + formataxes + xlabel('Time (s)') + xlim(tlims) + yl = ylim; + drawnow + end + set(gcf,'Color','w','PaperPositionMode','auto') + + %%%% + % Resize plots to appropriate size + if (plottype<4)&(length(ICuse)>=3) + bigpos = get(ax(1),'Position'); + aheight = 0.9*bigpos(4)/nrows; + for k=1:length(filtax) + axpos = get(filtax(k),'Position'); + axpos(3) = aheight; + axpos(4) = aheight; + set(filtax(k),'Position',axpos) + end + + set(gcf,'Units','normalized') + fpos = get(gcf,'Position'); + for j=1:ncols + axpos = get(ax(j),'OuterPosition'); + filtpos = get(filtax(1+(j-1)*nrows),'Position'); + axpos(1) = filtpos(1) + filtpos(3)*1.1; + set(ax(j),'OuterPosition',axpos,'ActivePositionProperty','outerposition') + axpos = get(ax(j),'Position'); + end + set(gcf,'Resize','on','Units','characters') + end + + for j=1:ncols + ax = [ax,axes('Position',get(ax(j),'Position'),'XAxisLocation','top','Color','none')]; + if plottype==4 + xt = get(ax4(end),'XTick'); + else + xt = get(ax(j),'XTick'); + end + xlim(tlims) + formataxes + set(gca,'YTick',[],'XTick',xt,'XTickLabel',num2str(xt'/dt, '%15.0f')) + xlabel('Frame number') + axes(ax(j)) + box on + end + linkaxes(ax,'xy') + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'contour'} + f_pos = get(gcf,'Position'); + if f_pos(4)>f_pos(3) + f_pos(4) = 0.5*f_pos(3); + set(gcf,'Position',f_pos); + end + set(gcf,'Renderer','zbuffer','RendererMode','manual') + + subplot(1,2,2) + + clf + colormap(gray) + + set(gcf,'DefaultAxesColorOrder',colord) + subplot(1,2,1) + + crange = [min(min(min(ica_filters(ICuse,:,:)))),max(max(max(ica_filters(ICuse,:,:))))]; + contourlevel = crange(2) - diff(crange)*[1,1]*0.8; + + cla + if ndims(f0)==2 + imagesc(f0) + else + image(f0) + end + cax = caxis; + sigsm = 1; + shading flat + hold on + for j=1:length(ICuse) + ica_filtersuse = gaussblur(squeeze(ica_filters(ICuse(j),:,:)), sigsm); + contour(ica_filtersuse, [1,1]*(mean(ica_filtersuse(:))+4*std(ica_filtersuse(:))), ... + 'Color',colord(mod(j-1,size(colord,1))+1,:),'LineWidth',2) + end + for j=1:length(ICuse) + ica_filtersuse = gaussblur(squeeze(ica_filters(ICuse(j),:,:)), sigsm); + + % Write the number at the cell center + [ypeak, xpeak] = find(ica_filtersuse == max(max(ica_filtersuse)),1); + text(xpeak,ypeak,num2str(j), 'horizontalalignment','c','verticalalignment','m','color','y') + end + hold off + caxis(cax) + formataxes + axis image tight off + title('Avg of movie, with contours of ICs') + + ax = subplot(1,2,2); + if plottype<=2 + complot(ica_sig, ICuse, dt) + end + set(gca,'ColorOrder',colord) + if plottype>=2 + [spcICuse,spcord] = ismember(spc,ICuse); + spc = spcord(spcICuse&(spt<=tlims(2))&(spt>=tlims(1))); + spt = spt(spcICuse&(spt<=tlims(2))&(spt>=tlims(1))); + alph = diff(ylim)/length(ICuse); + switch plottype + case 2 + % Traces with spikes + hold on + scatter(spt,-alph*(spc-1)+0.4*alph,20,colord(mod(spc-1,size(colord,1))+1,:),'o','filled') + hold off + case 3 + % Spike raster only + cla + scatter(spt,-alph*(spc-1)+0.4*alph,20,colord(mod(spc-1,size(colord,1))+1,:),'o','filled') + yticks=sort(-alph*([1:length(ICuse)]-1)+0.4*alph); + yl = sort(-alph*[length(ICuse)-1,-1]); + set(gca,'YTick',yticks) + set(gca,'YTicklabel',num2str(fliplr(ICuse)'),'YLim',yl) + set(gca,'YLim',sort(-alph*[length(ICuse)-1,-0.6])) + case 4 + % Binned spike rate + tvec = [ratebin/2:ratebin:size(ica_sig,2)*dt]; + binsize = ratebin*ones(size(tvec)); + binsize(end) = size(ica_sig,2)*dt-tvec(end-1)-ratebin/2; + for j=1:length(ICuse) + sptj = spt(spc==j); + spn = hist(sptj,tvec); + spr(j,:) = spn./binsize; + sperr(j,:) = sqrt(spn)./binsize; + end + spn = hist(spt,tvec); + mspr = spn./(binsize*length(ICuse)); + msperr = sqrt(spn)./(binsize*length(ICuse)); + + hold on + plot(tvec,mspr,'k','LineWidth',3); + patch([tvec(end:-1:1),tvec],[mspr(end:-1:1)+msperr(end:-1:1),mspr-msperr],0.5*[1,1,1],'FaceAlpha',1) + hold off + axis tight + formataxes + end + end + formataxes + xlim(tlims) + xlabel('Time (s)','FontAngle','i') + if plottype<=3 + ylabel('IC #','FontAngle','i') + else + ylabel('Spike rate (Hz)','FontAngle','i') + end + set(gcf,'Color','w','PaperPositionMode','auto') + set(gca,'yticklabel',num2str(fliplr([1:length(ICuse)])')) + + axes('Position',get(ax,'Position'),'XAxisLocation','top','Color','none') + xt = get(ax,'XTick'); + xlim(tlims) + formataxes + set(gca,'YTick',[], ... + 'XTick',xt,'XTickLabel',num2str(xt'/dt, '%15.0f')) + xlabel('Frame number') + axes(ax) + box on + +end + +%%%%%%%%%%%%%%%%%%%%% +function complot(sig, ICuse, dt) + +for i = 1:length(ICuse) + zsig(i, :) = zscore(sig(ICuse(i),:)); +end + +alpha = mean(max(zsig')-min(zsig')); +if islogical(zsig) + alpha = 1.5*alpha; +end + +zsig2 = zsig; +for i = 1:size(ICuse,2) + zsig2(i,:) = zsig(i,:) - alpha*(i-1)*ones(size(zsig(1,:))); +end + +tvec = [1:size(zsig,2)]*dt; +if islogical(zsig) + plot(tvec, zsig2','LineWidth',1) +else + plot(tvec, zsig2','LineWidth',1) +end +axis tight + +set(gca,'YTick',(-size(zsig,1)+1)*alpha:alpha:0); +set(gca,'YTicklabel',fliplr(ICuse)); + + +function formataxes + +set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out') +set(gcf,'Color','w','PaperPositionMode','auto') + +function fout = gaussblur(fin, smpix) +% +% Blur an image with a Gaussian kernel of s.d. smpix +% + +if ndims(fin)==2 + [x,y] = meshgrid([-ceil(3*smpix):ceil(3*smpix)]); + smfilt = exp(-(x.^2+y.^2)/(2*smpix^2)); + smfilt = smfilt/sum(smfilt(:)); + + fout = imfilter(fin, smfilt, 'replicate', 'same'); +else + [x,y] = meshgrid([-ceil(smpix):ceil(smpix)]); + smfilt = exp(-(x.^2+y.^2)/(2*smpix^2)); + smfilt = smfilt/sum(smfilt(:)); + + fout = imfilter(fin, smfilt, 'replicate', 'same'); +end diff --git a/CellsortPCA.m b/CellsortPCA.m new file mode 100644 index 0000000..7017444 --- /dev/null +++ b/CellsortPCA.m @@ -0,0 +1,327 @@ +function [mixedsig, mixedfilters, CovEvals, covtrace, movm, ... + movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes) +% [mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes) +% +% CELLSORT +% Read TIFF movie data and perform singular-value decomposition (SVD) +% dimensional reduction. +% +% Inputs: +% fn - movie file name. Must be in TIFF format. +% flims - 2-element vector specifying the endpoints of the range of +% frames to be analyzed. If empty, default is to analyze all movie +% frames. +% nPCs - number of principal components to be returned +% dsamp - optional downsampling factor. If scalar, specifies temporal +% downsampling factor. If two-element vector, entries specify temporal +% and spatial downsampling, respectively. +% outputdir - directory in which to store output .mat files +% badframes - optional list of indices of movie frames to be excluded +% from analysis +% +% Outputs: +% mixedsig - N x T matrix of N temporal signal mixtures sampled at T +% points. +% mixedfilters - N x X x Y array of N spatial signal mixtures sampled at +% X x Y spatial points. +% CovEvals - largest eigenvalues of the covariance matrix +% covtrace - trace of covariance matrix, corresponding to the sum of all +% eigenvalues (not just the largest few) +% movm - average of all movie time frames at each pixel +% movtm - average of all movie pixels at each time frame, after +% normalizing each pixel deltaF/F +% +% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 +% Email: eran@post.harvard.edu, mschnitz@stanford.edu +% + +% +% June 2016 (Cellsort 1.4): Fixed an inconsistency between create_tcov and +% creat_xcov. In the current version, both methods agree with the previous +% version of create_tcov. +% + +tic +fprintf('-------------- CellsortPCA %s: %s -------------- \n', date, fn) + +%----------------------- +% Check inputs +if isempty(dir(fn)) + error('Invalid input file name.') +end +if (nargin<2)||(isempty(flims)) + nt_full = tiff_frames(fn); + flims = [1,nt_full]; +end + +useframes = setdiff((flims(1):flims(2)), badframes); +nt = length(useframes); + +if nargin<3 || isempty(nPCs) + nPCs = min(150, nt); +end +if nargin<4 || isempty(dsamp) + dsamp = [1,1]; +end +if nargin<5 || isempty(outputdir) + outputdir = [pwd,'/cellsort_preprocessed_data/']; +end +if nargin<6 + badframes = []; +end +if isempty(dir(outputdir)) + mkdir(pwd, '/cellsort_preprocessed_data/') +end +if outputdir(end)~='/'; + outputdir = [outputdir, '/']; +end + +% Downsampling +if length(dsamp)==1 + dsamp_time = dsamp(1); + dsamp_space = 1; +else + dsamp_time = dsamp(1); + dsamp_space = dsamp(2); % Spatial downsample +end + +[fpath, fname] = fileparts(fn); +if isempty(badframes) + fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)), '_', date,'.mat']; +else + fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)),'_selframes_', date,'.mat']; +end +if ~isempty(dir(fnmat)) + fprintf('CELLSORT: Movie %s already processed;', ... + fn) + forceload = input(' Re-load data? [0-no/1-yes] '); + if isempty(forceload) || forceload==0 + load(fnmat) + return + end +end + +fncovmat = [outputdir, fname, '_cov_', num2str(flims(1)), ',', num2str(flims(2)), '_', date,'.mat']; + +[pixw,pixh] = size(imread(fn,1)); +npix = pixw*pixh; + +fprintf(' %d pixels x %d time frames;', npix, nt) +if nt1 + firstframe = imresize(firstframe, size(mov(:,:,1)),'bilinear'); +end + +%------------ +% Save the output data +save(fnmat,'mixedfilters','CovEvals','mixedsig', ... + 'movm','movtm','covtrace') +fprintf(' CellsortPCA: saving data and exiting; ') +toc + + function [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp) + %----------------------- + % Load movie data to compute the spatial covariance matrix + + npix = pixw*pixh; + + % Downsampling + if length(dsamp)==1 + dsamp_time = dsamp(1); + dsamp_space = 1; + else + dsamp_time = dsamp(1); + dsamp_space = dsamp(2); % Spatial downsample + end + + if (dsamp_space==1) + mov = zeros(pixw, pixh, nt); + for jjind=1:length(useframes) + jj = useframes(jjind); + mov(:,:,jjind) = imread(fn,jj); + if mod(jjind,500)==1 + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + end + end + else + [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' )); + mov = zeros(pixw_dsamp, pixh_dsamp, nt); + for jjind=1:length(useframes) + jj = useframes(jjind); + mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' ); + if mod(jjind,500)==1 + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + end + end + end + + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + mov = reshape(mov, npix, nt); + + % DFoF normalization of each pixel + movm = mean(mov,2); % Average over time + movmzero = (movm==0); + movm(movmzero) = 1; + mov = mov ./ (movm * ones(1,nt)) - 1; % Compute Delta F/F + mov(movmzero, :) = 0; + + if dsamp_time>1 + mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2); + mov = downsample(mov', dsamp)'; + end + + movtm = mean(mov,1); % Average over space + mov = mov - ones(size(mov,1),1)*movtm; + + covmat = (mov*mov')/size(mov,2); + covmat = covmat * size(mov,2)/size(mov,1); % Rescale to gree with create_tcov + toc + end + + function [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp) + %----------------------- + % Load movie data to compute the temporal covariance matrix + npix = pixw*pixh; + + % Downsampling + if length(dsamp)==1 + dsamp_time = dsamp(1); + dsamp_space = 1; + else + dsamp_time = dsamp(1); + dsamp_space = dsamp(2); % Spatial downsample + end + + if (dsamp_space==1) + mov = zeros(pixw, pixh, nt); + for jjind=1:length(useframes) + jj = useframes(jjind); + mov(:,:,jjind) = imread(fn,jj); + if mod(jjind,500)==1 + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + end + end + else + [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' )); + mov = zeros(pixw_dsamp, pixh_dsamp, nt); + for jjind=1:length(useframes) + jj = useframes(jjind); + mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' ); + if mod(jjind,500)==1 + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + end + end + end + + fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt) + toc + mov = reshape(mov, npix, nt); + + % DFoF normalization of each pixel + movm = mean(mov,2); % Average over time + movmzero = (movm==0); % Avoid dividing by zero + movm(movmzero) = 1; + mov = mov ./ (movm * ones(1,nt)) - 1; + mov(movmzero, :) = 0; + + if dsamp_time>1 + mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2); + mov = downsample(mov', dsamp)'; + end + + c1 = (mov'*mov)/npix; + movtm = mean(mov,1); % Average over space + covmat = c1 - movtm'*movtm; + clear c1 + end + + function [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix) + %----------------------- + % Perform SVD + + covtrace = trace(covmat) / npix; + + opts.disp = 0; + opts.issym = 'true'; + if nPCs0); + CovEvals = CovEvals(CovEvals>0); + end + + mixedsig = mixedsig' * nt; + CovEvals = CovEvals / npix; + + percentvar = 100*sum(CovEvals)/covtrace; + fprintf([' First ',num2str(nPCs),' PCs contain ',num2str(percentvar,3),'%% of the variance.\n']) + end + + function [mixedfilters] = reload_moviedata(npix, mov, mixedsig, CovEvals) + %----------------------- + % Re-load movie data + nPCs = size(mixedsig,1); + + Sinv = inv(diag(CovEvals.^(1/2))); + + movtm = mean(mov,1); % Average over space + movuse = mov - ones(npix,1) * movtm; + mixedfilters = reshape(movuse * mixedsig' * Sinv, npix, nPCs); + 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)); + end +end \ No newline at end of file diff --git a/CellsortPlotPCspectrum.m b/CellsortPlotPCspectrum.m new file mode 100644 index 0000000..a7e362a --- /dev/null +++ b/CellsortPlotPCspectrum.m @@ -0,0 +1,81 @@ +function CellsortPlotPCspectrum(fn, CovEvals, PCuse) +% CellsortPlotPCspectrum(fn, CovEvals, PCuse) +% +% Plot the principal component (PC) spectrum and compare with the +% corresponding random-matrix noise floor +% +% Inputs: +% fn - movie file name. Must be in TIFF format. +% CovEvals - eigenvalues of the covariance matrix +% PCuse - [optional] - indices of PCs included in dimensionally reduced +% data set +% +% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 +% Email: eran@post.harvard.edu, mschnitz@stanford.edu +% + +if nargin<3 + PCuse = []; +end + +[pixw,pixh] = size(imread(fn,1)); +npix = pixw*pixh; +nt = tiff_frames(fn); + +% Random matrix prediction (Sengupta & Mitra) +p1 = npix; % Number of pixels +q1 = nt; % Number of time frames +q = max(p1,q1); +p = min(p1,q1); +sigma = 1; +lmax = sigma*sqrt(p+q + 2*sqrt(p*q)); +lmin = sigma*sqrt(p+q - 2*sqrt(p*q)); +lambda = [lmin: (lmax-lmin)/100.0123423421: lmax]; +rho = (1./(pi*lambda*(sigma^2))).*sqrt((lmax^2-lambda.^2).*(lambda.^2-lmin^2)); +rho(isnan(rho)) = 0; +rhocdf = cumsum(rho)/sum(rho); +noiseigs = interp1(rhocdf, lambda, [p:-1:1]'/p, 'linear', 'extrap').^2 ; + +% Normalize the PC spectrum +normrank = min(nt-1,length(CovEvals)); +pca_norm = CovEvals*noiseigs(normrank) / (CovEvals(normrank)*noiseigs(1)); + +clf +plot(pca_norm, 'o-', 'Color', [1,1,1]*0.3, 'MarkerFaceColor', [1,1,1]*0.3, 'LineWidth',2) +hold on +plot(noiseigs / noiseigs(1), 'b-', 'LineWidth',2) +plot(2*noiseigs / noiseigs(1), 'b--', 'LineWidth',2) +if ~isempty(PCuse) + plot(PCuse, pca_norm(PCuse), 'rs', 'LineWidth',2) +end +hold off +formataxes +set(gca,'XScale','log','YScale','log', 'Color','none') +xlabel('PC rank') +ylabel('Normalized variance') +axis tight +if isempty(PCuse) + legend('Data variance','Noise floor','2 x Noise floor') +else + legend('Data variance','Noise floor','2 x Noise floor','Retained PCs') +end + +fntitle = fn; +fntitle(fn=='_') = ' '; +title(fntitle) + +function formataxes + +set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out') +set(gcf,'Color','w','PaperPositionMode','auto') + + +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)); diff --git a/CellsortSegmentation.m b/CellsortSegmentation.m new file mode 100644 index 0000000..18bd7fa --- /dev/null +++ b/CellsortSegmentation.m @@ -0,0 +1,152 @@ +function [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting) +% [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting) +% +%CellsortSegmentation +% Segment spatial filters derived by ICA +% +% Inputs: +% ica_filters - X x Y x nIC matrix of ICA spatial filters +% smwidth - standard deviation of Gaussian smoothing kernel (pixels) +% thresh - threshold for spatial filters (standard deviations) +% arealims - 2-element vector specifying the minimum and maximum area +% (in pixels) of segments to be retained; if only one element is +% specified, use this as the minimum area +% plotting - [0,1] whether or not to show filters +% +% Outputs: +% ica_segments - segmented spatial filters +% segmentabel - indices of the ICA filters from which each segment was derived +% segcentroid - X,Y centroid, in pixels, of each segment +% +% Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 +% Email: eran@post.harvard.edu, mschnitz@stanford.edu +% + +tic +fprintf('-------------- CellsortSegmentation %s -------------- \n', date) + +if (nargin<3)||isempty(thresh) + thresh = 2; +end +if (nargin<4)||isempty(arealims) + arealims = 200; +end +if (nargin<5)||isempty(plotting) + plotting = 0; +end +[nic,pixw,pixh] = size(ica_filters); + +ica_filtersorig = ica_filters / abs(std(ica_filters(:))); +ica_filters = (ica_filters - mean(ica_filters(:)))/abs(std(ica_filters(:))); +if smwidth>0 + % Smooth mixing filter with a Gaussian of s.d. smwidth pixels + smrange = max(5,3*smwidth); + [x,y] = meshgrid([-smrange:smrange]); + + smy = 1; smx = 1; + ica_filtersfilt = exp(-((x/smx).^2 + (y/smy).^2)/(2*smwidth^2)); + + ica_filtersfilt = ica_filtersfilt/sum(ica_filtersfilt(:)); + ica_filtersbw = false(pixw,pixh,nic); + tic + for j = 1:size(ica_filters,1) + ica_filtersuse = ica_filters(j,:,:); + ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:))); + ica_filtersbw(:,:,j) = (imfilter(ica_filtersuse, ica_filtersfilt, 'replicate', 'same') > thresh); + end +else + ica_filtersbw = (permute(ica_filters,[2,3,1]) > thresh); + ica_filtersfilt = 1; +end + +tic +if plotting + clf + set(gcf,'Color','w') + colormap(gray) + + subplot(223) + imagesc(squeeze(sum(ica_filters,1))) + axis image off + hold on +end +ica_filterslabel = []; +ica_segments = []; +k=0; +L=[]; +segmentlabel = []; +segcentroid = []; +[x,y] = meshgrid([1:pixh], [1:pixw]); +for j = 1:nic + % Label contiguous components + L = bwlabel(ica_filtersbw(:,:,j), 4); + Lu = 1:max(L(:)); + + % Delete small components + Larea = struct2array(regionprops(L, 'area')); + Lcent = regionprops(L, 'Centroid'); + + if length(arealims)==2 + Lbig = Lu( (Larea >= arealims(1))&(Larea <= arealims(2))); + Lsmall = Lu((Larea < arealims(1))|(Larea > arealims(2))); + else + Lbig = Lu(Larea >= arealims(1)); + Lsmall = Lu(Larea < arealims(1)); + end + + L(ismember(L,Lsmall)) = 0; + + for jj = 1:length(Lbig) + segcentroid(jj+k,:) = Lcent(Lbig(jj)).Centroid; + end + + ica_filtersuse = squeeze(ica_filtersorig(j,:,:)); + for jj = 1:length(Lbig) + ica_segments(jj+k,:,:) = ica_filtersuse .* ( 0*(L==0) + (L==Lbig(jj)) ); % Exclude background + end + + if plotting && ~isempty(Lbig) + if smwidth>0 + subplot(2,2,2) + ica_filtersuse = squeeze(ica_filters(j,:,:)); + ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:))); + imagesc(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'),[-1,4]) + hold on + contour(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'), [1,1]*thresh, 'k') + hold off + hc = colorbar('Position',[0.9189 0.6331 0.0331 0.2253]); + ylabel(hc,'Std. dev.') + title(['IC ',num2str(j),' smoothed']) + axis image off + + subplot(2,2,1) + else + subplot(211) + end + imagesc(squeeze(ica_filters(j,:,:))) + title(['IC ',num2str(j),' original']) + axis image off + + colord = lines(k+length(Lbig)); + for jj = 1:length(Lbig) + subplot(223) + contour(ica_filtersbw(:,:,j), [1,1]*0.5, 'color',colord(jj+k,:),'linewidth',2) + hold on + text(segcentroid(jj+k,1), segcentroid(jj+k,2), num2str(jj+k), 'horizontalalignment','c', 'verticalalignment','m') + set(gca, 'ydir','reverse','tickdir','out') + axis image + xlim([0,pixw]); ylim([0,pixh]) + + subplot(224) + imagesc(squeeze(ica_segments(jj+k,:,:))) + hold on + plot(segcentroid(jj+k,1), segcentroid(jj+k,2), 'bo') + hold off + axis image off + title(['Segment ',num2str(jj+k)]) + drawnow + end + end + k = size(ica_segments,1); +end +toc diff --git a/README.md b/README.md new file mode 100644 index 0000000..18aac11 --- /dev/null +++ b/README.md @@ -0,0 +1,26 @@ +CellSort 1.1 + +Independent Component Analysis for Optical Imaging Data Sets + +Copyright Eran Mukamel, 2009 email: eran@post.harvard.edu + +INSTALLATION ---------------------- +Add the folder "CellSort 1.1" to the MATLAB path by +selecting "Set Path" in the file menu. + +INSTRUCTIONS ---------------------- D + +ocumentation about each of the routines in this +package is available in the "doc" directory. For more information, see "Automated analysis +of cellular signals from large-scale calcium imaging data" by Eran Mukamel, Axel Nimmerjahn +and Mark Schnitzer, NEURON (2009). + + +UPDATE 1.3 - APRIL 2013 ------------------------------- Fixed a problem with the function +tiff_frames that altered its behavior in MATLAB 2012. + +UPDATE 1.4 - JUNE 2016 ------------------------------- +Fixed an inconsistency in CellsortPCA between create_tcov and +creat_xcov. In the current version, both methods agree with the previous +version of create_tcov. + diff --git a/README.rtf b/README.rtf new file mode 100644 index 0000000..47da4be --- /dev/null +++ b/README.rtf @@ -0,0 +1,25 @@ +{\rtf1\ansi\ansicpg1252\cocoartf1038\cocoasubrtf360 +{\fonttbl\f0\fswiss\fcharset0 Helvetica;} +{\colortbl;\red255\green255\blue255;} +\margl1440\margr1440\vieww9000\viewh8400\viewkind0 +\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\ql\qnatural\pardirnatural + +\f0\fs24 \cf0 CellSort 1.1\ +\ +Independent Component Analysis for Optical Imaging Data Sets\ +\ +Copyright Eran Mukamel, 2009\ +email: eran@post.harvard.edu\ +\ +INSTALLATION\ +----------------------\ +Add the folder "CellSort 1.1" to the MATLAB path by selecting "Set Path" in the file menu.\ +\ +INSTRUCTIONS\ +----------------------\ +Documentation about each of the routines in this package is available in the "doc" directory. For more information, see "Automated analysis of cellular signals from large-scale calcium imaging data" by Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, NEURON (2009).\ +\ +\ +UPDATE 1.3 - APRIL 2013\ +-------------------------------\ +Fixed a problem with the function tiff_frames that altered its behavior in MATLAB 2012.} \ No newline at end of file diff --git a/doc/CellsortApplyFilter.html b/doc/CellsortApplyFilter.html new file mode 100644 index 0000000..bd8b196 --- /dev/null +++ b/doc/CellsortApplyFilter.html @@ -0,0 +1,144 @@ + + + + Description of CellsortApplyFilter + + + + + + + + + +
CellSort 1.0 > CellsortApplyFilter.m
+ + + +

CellsortApplyFilter +

+ +

PURPOSE ^

+
cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+ +

SYNOPSIS ^

+
function cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+ +

DESCRIPTION ^

+
 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
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+0002 % cell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
+0003 %
+0004 %CellsortApplyFilter
+0005 % Read in movie data and output signals corresponding to specified spatial
+0006 % filters
+0007 %
+0008 % Inputs:
+0009 %     fn - file name of TIFF movie file
+0010 %     ica_segments - nIC x X matrix of ICA spatial filters
+0011 %     flims - optional two-element vector of frame limits to be read
+0012 %     movm - mean fluorescence image
+0013 %     subtractmean - boolean specifying whether or not to subtract the mean
+0014 %     fluorescence of each time frame
+0015 %
+0016 % Outputs:
+0017 %     cell_sig - cellular signals
+0018 %
+0019 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0020 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0021 
+0022 if (nargin<3)||isempty(flims)
+0023     nt = tiff_frames(fn);
+0024     flims = [1,nt];
+0025 else
+0026     nt = diff(flims)+1;
+0027 end
+0028 if nargin<5
+0029     subtractmean = 1;
+0030 end
+0031 
+0032 [pixw,pixh] = size(imread(fn,1));
+0033 if (nargin<4)||isempty(movm)
+0034     movm = ones(pixw,pixh);
+0035 else
+0036     movm = double(movm);
+0037 end
+0038 movm(movm==0) = 1; % Just in case there are black areas in the average image
+0039 k=0;
+0040 
+0041 cell_sig = zeros(size(ica_segments,1), nt);
+0042 ica_segments = reshape(ica_segments, [], pixw*pixh);
+0043 
+0044 fprintf('Loading %5g frames for %d ROIs.\n', nt, size(ica_segments,1))
+0045 while k<nt
+0046     ntcurr = min(500, nt-k);
+0047     mov = zeros(pixw, pixh, ntcurr);
+0048     for j=1:ntcurr
+0049         movcurr = imread(fn, j+k+flims(1)-1);
+0050         mov(:,:,j) = movcurr;
+0051     end
+0052     mov = (mov ./ repmat(movm, [1,1,ntcurr])) - 1; % Normalize by background and subtract mean
+0053 
+0054     if subtractmean
+0055         % Subtract the mean of each frame
+0056         movtm = mean(mean(mov,1),2);
+0057         mov = mov - repmat(movtm,[pixw,pixh,1]);
+0058     end
+0059 
+0060     mov = reshape(mov, pixw*pixh, ntcurr);
+0061     cell_sig(:, k+[1:ntcurr]) = ica_segments*mov;
+0062 
+0063     k=k+ntcurr;
+0064     fprintf('Loaded %3.0f frames; ', k)
+0065     toc
+0066 end
+0067 
+0068 
+0069     function j = tiff_frames(fn)
+0070         %
+0071         % n = tiff_frames(filename)
+0072         %
+0073         % Returns the number of slices in a TIFF stack.
+0074         %
+0075         % Modified April 9, 2013 for compatibility with MATLAB 2012b
+0076 
+0077         j = length(imfinfo(fn));
+0078 end
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortChoosePCs.html b/doc/CellsortChoosePCs.html new file mode 100644 index 0000000..40a1591 --- /dev/null +++ b/doc/CellsortChoosePCs.html @@ -0,0 +1,175 @@ + + + + Description of CellsortChoosePCs + + + + + + + + + +
CellSort 1.0 > CellsortChoosePCs.m
+ + + +

CellsortChoosePCs +

+ +

PURPOSE ^

+
[PCuse] = CellsortChoosePCs(fn, mixedfilters)
+ +

SYNOPSIS ^

+
function [PCuse] = CellsortChoosePCs(fn, mixedfilters)
+ +

DESCRIPTION ^

+
 [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
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function [PCuse] = CellsortChoosePCs(fn, mixedfilters)
+0002 % [PCuse] = CellsortChoosePCs(fn, mixedfilters)
+0003 %
+0004 % Allows the user to select which principal components will be kept
+0005 % following dimensional reduction.
+0006 %
+0007 % Inputs:
+0008 %   fn - movie file name. Must be in TIFF format.
+0009 %   mixedfilters - N x X matrix of N spatial signal mixtures sampled at X
+0010 %   spatial points.
+0011 %
+0012 % Outputs:
+0013 %   PCuse - vector of indices of the PCs to be kept for dimensional
+0014 %   reduction
+0015 %
+0016 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0017 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0018 
+0019 fprintf('-------------- CellsortICA %s -------------- \n', date)
+0020 
+0021 [pixw,pixh] = size(imread(fn,1));
+0022 
+0023 npcs = 20; % Number of PCs to display concurrently
+0024 currpcs = [1:npcs];
+0025 PCf = [];
+0026 while isempty(PCf)
+0027     showpcs(currpcs, mixedfilters, pixw, pixh)
+0028     yl = ylim;
+0029     xl = xlim;
+0030     set(gca,'Units','pixels')
+0031     title(['Choose first PC; showing PCs [',num2str(currpcs(1)),':',num2str(currpcs(end)),']'])
+0032     PCf = input('Number of first PC to retain, Klow (''b/f'' to scroll backwards/forwards)): ','s');
+0033     if PCf=='b'
+0034         currpcs = currpcs - min(npcs,currpcs(1)-1);
+0035         PCf = [];
+0036     elseif (PCf=='f')
+0037         currpcs = currpcs+npcs;
+0038         if nnz(currpcs>size(mixedfilters,2))
+0039             currpcs = [-npcs+1:0]+size(mixedfilters,2);
+0040             fprintf('Reached end of stored PCs.\n')
+0041         end
+0042         PCf = [];
+0043     else
+0044         PCf = str2num(PCf);
+0045     end
+0046 end
+0047 PCl=[];
+0048 currpcs = [PCf:PCf+npcs-1];
+0049 while isempty(PCl)
+0050     showpcs(currpcs, mixedfilters, pixw, pixh)
+0051     title(['Choose last PC; showing PCs [',num2str(currpcs(1)),':',num2str(currpcs(end)),']'])
+0052     PCl = input('Number of last PC to retain, Khigh (''b/f'' to scroll backwards/forwards): ','s');
+0053     if PCl=='b'
+0054         currpcs = currpcs - min(npcs,currpcs(1)-1);
+0055         PCl = [];
+0056     elseif (PCl=='f')
+0057         currpcs = currpcs+npcs;
+0058         if nnz(currpcs>size(mixedfilters,2))
+0059             currpcs = [-npcs+1:0]+size(mixedfilters,2);
+0060             fprintf('Reached end of stored PCs.\n')
+0061         end
+0062         PCl = [];
+0063     else
+0064         PCl = str2num(PCl);
+0065     end
+0066 end
+0067 currpcs = [PCf:PCl];
+0068 PCbad=[];
+0069 showpcs(currpcs, mixedfilters, pixw, pixh)
+0070 
+0071 PCuse = setdiff(currpcs, PCbad);
+0072 showpcs(PCuse, mixedfilters, pixw, pixh)
+0073 
+0074 fprintf('  Retaining PCs in the range [Klow - Khigh] = [%d - %d].\n', PCf,PCl)
+0075 
+0076 function showpcs(usepcs, Efull, pixw, pixh)
+0077 
+0078 if nargin<3
+0079     fprintf('Assuming movie frames are square.\n')
+0080     pixw = sqrt(size(Efull,1));
+0081     pixh = sqrt(size(Efull,1));
+0082 end
+0083 if isempty(usepcs)
+0084     usepcs = [1:size(Efull,2)];
+0085 end
+0086 
+0087 if ndims(Efull)>=3
+0088     Efull = reshape(Efull, pixw*pixh, []);
+0089 end
+0090 for j=usepcs
+0091     Efull(:,j) = zscore(Efull(:,j));
+0092 end
+0093 pcs = reshape(Efull(:,usepcs), pixw, pixh, []);
+0094 pcs = permute(pcs, [1, 2, 4, 3]);
+0095 montage(pcs)
+0096 colormap(hot)
+0097 axis on
+0098 xl = xlim;
+0099 yl = ylim;
+0100 nw = ceil(xl(2)/pixh)-1;
+0101 nh = ceil(yl(2)/pixw)-1;
+0102 set(gca,'YTick',[pixw:pixw:yl(2)],'YTickLabel',  num2str(usepcs(min([0:nh]*nw+1, length(usepcs)))'), ...
+0103     'XTick',[pixh:pixh:xl(2)], ...
+0104     'XTickLabel',num2str(usepcs([(nh-1)*nw+1:length(usepcs)])'), 'XAxisLocation','bottom','LineWidth',2)
+0105 grid on
+0106 formataxes
+0107 caxis([-1,1]*7)
+0108 
+0109 function formataxes
+0110 
+0111 set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out')
+0112 set(gcf,'Color','w','PaperPositionMode','auto')
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortFindspikes.html b/doc/CellsortFindspikes.html new file mode 100644 index 0000000..31dd734 --- /dev/null +++ b/doc/CellsortFindspikes.html @@ -0,0 +1,120 @@ + + + + Description of CellsortFindspikes + + + + + + + + + +
CellSort 1.0 > CellsortFindspikes.m
+ + + +

CellsortFindspikes +

+ +

PURPOSE ^

+
[spmat, spt, spc, zsig] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+ +

SYNOPSIS ^

+
function [spmat, spt, spc] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+ +

DESCRIPTION ^

+
 [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
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + + +

SOURCE CODE ^

+
0001 function [spmat, spt, spc] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+0002 % [spmat, spt, spc, zsig] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+0003 %
+0004 % CELLSORT
+0005 %  Deconvolve signal and find spikes using a threshold
+0006 %
+0007 % Inputs:
+0008 %   ica_sig - nIC x T matrix of ICA temporal signals
+0009 %   thresh - threshold for spike detection
+0010 %   dt - time step
+0011 %   deconvtau - time constant for temporal deconvolution (Butterworth
+0012 %   filter); if deconvtau=0 or [], no deconvolution is performed
+0013 %   normalization - type of normalization to apply to ica_sig; 0 - no
+0014 %   normalization; 1 - divide by standard deviation and subtract mean
+0015 %
+0016 % Outputs:
+0017 %   spmat - nIC x T sparse binary matrix, containing 1 at the time frame of each
+0018 %   spike
+0019 %   spt - list of all spike times
+0020 %   spc - list of the indices of cells for each spike
+0021 %
+0022 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0023 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0024 
+0025 if size(ica_sig,2)==1
+0026     ica_sig = ica_sig';
+0027 end
+0028 
+0029 if (nargin>=3)&&(deconvtau>0)
+0030     dsig = diff(ica_sig,1,2);
+0031     sig = ica_sig/deconvtau + [dsig(:,1),dsig]/dt;
+0032 else
+0033     sig = ica_sig;
+0034 end
+0035 
+0036 if (nargin<2)
+0037     thresh=3;
+0038     fprintf('Using threshold = 3 s.d. \n')
+0039 end
+0040 switch normalization
+0041     case 0 % Absolute units
+0042         zsig = sig';
+0043     case 1 % Standard-deviation
+0044         zsig = zscore(sig');
+0045 end
+0046 pp1=[zsig(1,:);zsig(1:end-1,:)];
+0047 pp2=[zsig(2:end,:);zsig(end,:)];
+0048 spmat = sparse((zsig>=thresh)&(zsig-pp1>=0)&(zsig-pp2>=0));
+0049 
+0050 if nargout>1
+0051     [spt,spc] = find(spmat);
+0052     spt = spt*dt;
+0053 end
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortICA.html b/doc/CellsortICA.html new file mode 100644 index 0000000..9a396e9 --- /dev/null +++ b/doc/CellsortICA.html @@ -0,0 +1,233 @@ + + + + Description of CellsortICA + + + + + + + + + +
CellSort 1.0 > CellsortICA.m
+ + + +

CellsortICA +

+ +

PURPOSE ^

+
[ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
+ +

SYNOPSIS ^

+
function [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig,mixedfilters, CovEvals, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
+ +

DESCRIPTION ^

+
 [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
+
+CELLSORT
+ Perform ICA with a standard set of parameters, including skewness as the
+ objective function
+
+ Inputs:
+   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
+   points.
+   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
+   X x Y spatial points.
+   CovEvals - eigenvalues of the covariance matrix
+   PCuse - vector of indices of the components to be included. If empty,
+   use all the components
+   mu - parameter (between 0 and 1) specifying weight of temporal
+   information in spatio-temporal ICA
+   nIC - number of ICs to derive
+   termtol - termination tolerance; fractional change in output at which
+   to end iteration of the fixed point algorithm.
+   maxrounds - maximum number of rounds of iterations
+
+ Outputs:
+     ica_sig - nIC x T matrix of ICA temporal signals
+     ica_filters - nIC x X x Y array of ICA spatial filters
+     ica_A - nIC x N orthogonal unmixing matrix to convert the input to output signals
+     numiter - number of rounds of iteration before termination
+
+ Routine is based on the fastICA package (Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo
+ Hyvärinen, http://www.cis.hut.fi/projects/ica/fastica)
+
+ Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+ Email: eran@post.harvard.edu, mschnitz@stanford.edu
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, ...
+0002     mixedfilters, CovEvals, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
+0003 % [ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
+0004 %
+0005 %CELLSORT
+0006 % Perform ICA with a standard set of parameters, including skewness as the
+0007 % objective function
+0008 %
+0009 % Inputs:
+0010 %   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
+0011 %   points.
+0012 %   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
+0013 %   X x Y spatial points.
+0014 %   CovEvals - eigenvalues of the covariance matrix
+0015 %   PCuse - vector of indices of the components to be included. If empty,
+0016 %   use all the components
+0017 %   mu - parameter (between 0 and 1) specifying weight of temporal
+0018 %   information in spatio-temporal ICA
+0019 %   nIC - number of ICs to derive
+0020 %   termtol - termination tolerance; fractional change in output at which
+0021 %   to end iteration of the fixed point algorithm.
+0022 %   maxrounds - maximum number of rounds of iterations
+0023 %
+0024 % Outputs:
+0025 %     ica_sig - nIC x T matrix of ICA temporal signals
+0026 %     ica_filters - nIC x X x Y array of ICA spatial filters
+0027 %     ica_A - nIC x N orthogonal unmixing matrix to convert the input to output signals
+0028 %     numiter - number of rounds of iteration before termination
+0029 %
+0030 % Routine is based on the fastICA package (Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo
+0031 % Hyvärinen, http://www.cis.hut.fi/projects/ica/fastica)
+0032 %
+0033 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0034 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0035 
+0036 fprintf('-------------- CellsortICA %s -------------- \n', date)
+0037 
+0038 if (nargin<4) || isempty(PCuse)
+0039     PCuse = [1:size(mixedsig,1)];
+0040 end
+0041 if (nargin<6) || isempty(nIC)
+0042     nIC = length(PCuse);
+0043 end
+0044 if (nargin<7) || isempty(ica_A_guess)
+0045     ica_A_guess = randn(length(PCuse), nIC);
+0046 end
+0047 if (nargin<8) || isempty(termtol)
+0048     termtol = 1e-6;
+0049 end
+0050 if (nargin<9) || isempty(maxrounds)
+0051     maxrounds = 100;
+0052 end
+0053 if isempty(mu)||(mu>1)||(mu<0)
+0054     error('Spatio-temporal parameter, mu, must be between 0 and 1.')
+0055 end
+0056 
+0057 % Check that ica_A_guess is the right size
+0058 if size(ica_A_guess,1)~= length(PCuse) || size(ica_A_guess,2)~=nIC
+0059     error('Initial guess for ica_A is the wrong size.')
+0060 end
+0061 if nIC>length(PCuse)
+0062     error('Cannot estimate more ICs than the number of PCs.')
+0063 end
+0064     
+0065 [pixw,pixh] = size(mixedfilters(:,:,1));
+0066 npix = pixw*pixh;
+0067 
+0068 % Select PCs
+0069 if mu > 0 || ~isempty(mixedsig)
+0070     mixedsig = mixedsig(PCuse,:);
+0071 end
+0072 if mu < 1 || ~isempty(mixedfilters)
+0073     mixedfilters = reshape(mixedfilters(:,:,PCuse),npix,length(PCuse));
+0074 end
+0075 CovEvals = CovEvals(PCuse);
+0076 
+0077 % Center the data by removing the mean of each PC
+0078 mixedmean = mean(mixedsig,2);
+0079 mixedsig = mixedsig - mixedmean * ones(1, size(mixedsig,2));
+0080 
+0081 % Create concatenated data for spatio-temporal ICA
+0082 nx = size(mixedfilters,1);
+0083 nt = size(mixedsig,2);
+0084 if mu == 1
+0085     % Pure temporal ICA
+0086     sig_use = mixedsig;
+0087 elseif mu == 0
+0088     % Pure spatial ICA
+0089     sig_use = mixedfilters';
+0090 else
+0091     % Spatial-temporal ICA
+0092     sig_use = [(1-mu)*mixedfilters', mu*mixedsig];
+0093     sig_use = sig_use / sqrt(1-2*mu+2*mu^2); % This normalization ensures that, if both mixedfilters and mixedsig have unit covariance, then so will sig_use
+0094 end
+0095 
+0096 % Perform ICA
+0097 [ica_A, numiter] = fpica_standardica(sig_use, nIC, ica_A_guess, termtol, maxrounds);
+0098 
+0099 % Sort ICs according to skewness of the temporal component
+0100 ica_W = ica_A';
+0101 
+0102 ica_sig = ica_W * mixedsig;
+0103 ica_filters = reshape((mixedfilters*diag(CovEvals.^(-1/2))*ica_A)', nIC, nx);  % This is the matrix of the generators of the ICs
+0104 ica_filters = ica_filters / npix^2;
+0105 
+0106 icskew = skewness(ica_sig');
+0107 [icskew, ICord] = sort(icskew, 'descend');
+0108 ica_A = ica_A(:,ICord);
+0109 ica_sig = ica_sig(ICord,:);
+0110 ica_filters = ica_filters(ICord,:);
+0111 ica_filters = reshape(ica_filters, nIC, pixw, pixh);
+0112 
+0113 % Note that with these definitions of ica_filters and ica_sig, we can decompose
+0114 % the sphered and original movie data matrices as:
+0115 %     mov_sphere ~ mixedfilters * mixedsig = ica_filters * ica_sig = (mixedfilters*ica_A') * (ica_A*mixedsig),
+0116 %     mov ~ mixedfilters * pca_D * mixedsig.
+0117 % This gives:
+0118 %     ica_filters = mixedfilters * ica_A' = mov * mixedsig' * inv(diag(pca_D.^(1/2)) * ica_A'
+0119 %     ica_sig = ica_A * mixedsig = ica_A * inv(diag(pca_D.^(1/2))) * mixedfilters' * mov
+0120 
+0121     function [B, iternum] = fpica_standardica(X, nIC, ica_A_guess, termtol, maxrounds)
+0122         
+0123         numSamples = size(X,2);
+0124         
+0125         B = ica_A_guess;
+0126         BOld = zeros(size(B));
+0127         
+0128         iternum = 0;
+0129         minAbsCos = 0;
+0130         
+0131         errvec = zeros(maxrounds,1);
+0132         while (iternum < maxrounds) && ((1 - minAbsCos)>termtol)
+0133             iternum = iternum + 1;
+0134             
+0135             % Symmetric orthogonalization.
+0136             B = (X * ((X' * B) .^ 2)) / numSamples;
+0137             B = B * real(inv(B' * B)^(1/2));
+0138             
+0139             % Test for termination condition.
+0140             minAbsCos = min(abs(diag(B' * BOld)));
+0141             
+0142             BOld = B;
+0143             errvec(iternum) = (1 - minAbsCos);
+0144         end
+0145         
+0146         if iternum<maxrounds
+0147             fprintf('Convergence in %d rounds.\n', iternum)
+0148         else
+0149             fprintf('Failed to converge; terminating after %d rounds, current change in estimate %3.3g.\n', ...
+0150                 iternum, 1-minAbsCos)
+0151         end
+0152     end
+0153 end
+0154
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortICAplot.html b/doc/CellsortICAplot.html new file mode 100644 index 0000000..f2a6471 --- /dev/null +++ b/doc/CellsortICAplot.html @@ -0,0 +1,489 @@ + + + + Description of CellsortICAplot + + + + + + + + + +
CellSort 1.0 > CellsortICAplot.m
+ + + +

CellsortICAplot +

+ +

PURPOSE ^

+
CellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
+ +

SYNOPSIS ^

+
function CellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
+ +

DESCRIPTION ^

+
 CellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
+
+ Display the results of ICA analysis in the form of paired spatial filters
+ and signal time courses
+
+ Inputs:
+     mode - 'series' shows each spatial filter separately; 'contour'
+     displays a single plot with contours for all spatial filters
+     superimposed on the mean fluorescence image
+     ica_filters - nIC x X x Y array of ICA spatial filters
+     ica_sig - nIC x T matrix of ICA temporal signals
+     f0 - mean fluorescence image
+     tlims - 2-element vector specifying the range of times to be displayed
+     dt - time step corresponding to individual movie time frames
+     ratebin - size of time bins for spike rate computation
+     plottype - type of spike plot to use:
+         plottype = 1: plot cellular signals
+         plottype = 2: plot cellular signals together with spikes
+         plottype = 3: plot spikes only
+         plottype = 4: plot spike rate over time
+     ICuse - vector of indices of cells to be plotted
+     spt - vector of spike times (needed for plottype > 1)
+     spc - vector of spike indices (which cell) (needed for plottype > 1)
+
+ Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+ Email: eran@post.harvard.edu, mschnitz@stanford.edu
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function CellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
+0002 % CellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
+0003 %
+0004 % Display the results of ICA analysis in the form of paired spatial filters
+0005 % and signal time courses
+0006 %
+0007 % Inputs:
+0008 %     mode - 'series' shows each spatial filter separately; 'contour'
+0009 %     displays a single plot with contours for all spatial filters
+0010 %     superimposed on the mean fluorescence image
+0011 %     ica_filters - nIC x X x Y array of ICA spatial filters
+0012 %     ica_sig - nIC x T matrix of ICA temporal signals
+0013 %     f0 - mean fluorescence image
+0014 %     tlims - 2-element vector specifying the range of times to be displayed
+0015 %     dt - time step corresponding to individual movie time frames
+0016 %     ratebin - size of time bins for spike rate computation
+0017 %     plottype - type of spike plot to use:
+0018 %         plottype = 1: plot cellular signals
+0019 %         plottype = 2: plot cellular signals together with spikes
+0020 %         plottype = 3: plot spikes only
+0021 %         plottype = 4: plot spike rate over time
+0022 %     ICuse - vector of indices of cells to be plotted
+0023 %     spt - vector of spike times (needed for plottype > 1)
+0024 %     spc - vector of spike indices (which cell) (needed for plottype > 1)
+0025 %
+0026 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0027 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0028 %
+0029 
+0030 colord=[         0         0    1.0000
+0031     0    0.4000         0
+0032     1.0000         0         0
+0033     0    0.7500    0.7500
+0034     0.7500         0    0.7500
+0035     0.8, 0.5, 0
+0036     0         0    0.5
+0037     0         0.85      0];
+0038 
+0039 % Check input arguments
+0040 nIC = size(ica_sig,1);
+0041 if nargin<5 || isempty(tlims)
+0042     tlims = [0, size(ica_sig,2)*dt]; % seconds
+0043 end
+0044 if nargin<8 || isempty(plottype)
+0045     plottype = 1;
+0046 end
+0047 if nargin<9 || isempty(ICuse)
+0048     ICuse = [1:nIC];
+0049 end
+0050 if size(ICuse,2)==1
+0051     ICuse = ICuse';
+0052 end
+0053 
+0054 % Reshape the filters
+0055 [pixw,pixh] = size(f0);
+0056 if size(ica_filters,1)==nIC
+0057     ica_filters = reshape(ica_filters,nIC,pixw,pixh);
+0058 elseif size(ica_filters,2)==nIC
+0059     ica_filters = reshape(ica_filters,nIC,pixw,pixh);
+0060 end
+0061 
+0062 switch mode
+0063     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+0064     case {'series'}
+0065         colmax = 20; % Maximum # of ICs in one column
+0066         ncols = ceil(length(ICuse)/colmax);
+0067         if plottype==4
+0068             ncols=1;
+0069         end
+0070         nrows = ceil(length(ICuse)/ncols);
+0071 
+0072         if size(ica_filters(:,:,1))==size(f0(:,:,1))
+0073             ica_filters = permute(ica_filters,[3,1,2]);
+0074         end
+0075 
+0076         subplot(1,3*ncols,[2:3])
+0077         tlims(2) = min(tlims(2),size(ica_sig,2)*dt);
+0078         tlims(1) = max(tlims(1),0);
+0079 
+0080         clf
+0081         f_pos = get(gcf,'Position');
+0082         f_pos(4) = max([f_pos(4),500,50*nrows]);
+0083         f_pos(3) = max(400*ncols,0.9*f_pos(4));
+0084 
+0085         colormap(hot)
+0086         colord=get(gca,'ColorOrder');
+0087         ll=0;
+0088         filtax = [];
+0089         if ~isempty(ica_filters)
+0090             for k=0:ncols-1
+0091                 jj=3*k;
+0092                 nrows_curr = min(nrows,length(ICuse)-k*nrows);
+0093                 for j=1:nrows_curr
+0094                     filtax= [filtax,subplot(nrows_curr + (plottype==4),3*ncols, jj+1)];
+0095                     jj=jj+3*ncols;
+0096                     ll=ll+1;
+0097                     imagesc(squeeze(ica_filters(ICuse(ll),:,:)))
+0098                     axis image tight off
+0099                 end
+0100             end
+0101         end
+0102 
+0103         ax = [];
+0104         for j=0:ncols-1
+0105             ax(j+1)=subplot(1,3*ncols,3*j+[2:3]);
+0106             ICuseuse = ICuse([1+j*nrows:min(length(ICuse),(j+1)*nrows)]);
+0107             if plottype<=2
+0108                 complot(ica_sig, ICuseuse, dt)
+0109             end
+0110             formataxes
+0111             if plottype>=2
+0112                 [spcICuse,spcord] = ismember(spc,ICuseuse);
+0113                 spc_curr = spcord(spcICuse&(spt<=tlims(2))&(spt>=tlims(1)));
+0114                 spt_curr = spt(spcICuse&(spt<=tlims(2))&(spt>=tlims(1)));
+0115                 alpha = diff(ylim)/length(ICuseuse);
+0116                 switch plottype
+0117                     case 2
+0118                         hold on
+0119                         scatter(spt_curr,-alpha*(spc_curr-1)+2,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled')
+0120                         hold off
+0121                     case 3
+0122                         cla
+0123                         scatter(spt_curr,-alpha*(spc_curr-1)+0.4*alpha,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled')
+0124                         yticks=sort(-alpha*([1:length(ICuseuse)]-1)+0.4*alpha);
+0125                         yl = sort(-alpha*[length(ICuseuse)-1,-1]);
+0126                         set(gca,'YTick',yticks)
+0127                         set(gca,'YTicklabel',num2str(fliplr(ICuseuse)'),'YLim',yl)
+0128                         set(gca,'YLim',sort(-alpha*[length(ICuseuse)-1,-0.6]))
+0129 
+0130                     case 4
+0131                         ax4=[];
+0132                         tvec = [ratebin/2:ratebin:size(ica_sig,2)*dt];
+0133                         binsize = ratebin*ones(size(tvec));
+0134                         binsize(end) = size(ica_sig,2)*dt-tvec(end-1)-ratebin/2;
+0135                         fprintf('IC mean rate\n')
+0136                         sprm = []; sprs = [];
+0137                         for jj=1:length(ICuseuse)
+0138                             sptj = spt_curr(spc_curr==jj);
+0139                             spn = hist(sptj,tvec);
+0140                             spr = spn./binsize;
+0141                             sperr = sqrt(spn)./binsize;
+0142 
+0143                             ax4(jj) = subplot(length(ICuseuse)+1,3,[3*jj-1,3*jj]);
+0144                             errorbar(tvec,spr,sperr,'Color',colord(mod(jj-1,size(colord,1))+1,:),'LineWidth',1)
+0145                             hold on
+0146                             bar(tvec,spr)
+0147                             hold off
+0148                             axis tight
+0149                             formataxes
+0150                             ylabel({['IC',int2str(ICuseuse(jj))],[num2str(mean(spr),2),'']})
+0151                             sprm(jj) = sum(spn)/sum(binsize);
+0152                             fprintf('%3.3f\n', sprm(jj))
+0153                             box off
+0154 
+0155                         end
+0156                         histax = subplot(length(ICuseuse)+1,3, 3*length(ICuseuse)+1);
+0157                         hist(sprm,[0:0.1:1])
+0158                         formataxes
+0159                         xlabel('Spike rate (Hz)')
+0160                         ylabel('# Cells')
+0161                         xlim([0,1])
+0162 
+0163                         spn = hist(spt_curr,tvec);
+0164                         spr = spn./(binsize*length(ICuseuse));
+0165                         sperr = sqrt(spn)./(binsize*length(ICuseuse));  %Standard error for measurement of mean from Poisson process
+0166 
+0167                         ax4(jj+1) = subplot(length(ICuseuse)+1,3,3*length(ICuseuse)+[2,3]);
+0168                         errorbar(tvec,spr,sperr,'k','LineWidth',3)
+0169                         formataxes
+0170                         set(ax4,'XLim',tlims)
+0171                         set(ax4(1:jj-1),'XTick',[])
+0172                         ylabel({'Pop. mean',[num2str(mean(spr),2),' Hz']},'FontAngle','i')
+0173                         fprintf('Pop. ave\tPop. SD\n')
+0174                         fprintf('%3.3f\t\t%3.3f\n', mean(sprm), std(sprm))
+0175                         xlabel('Time (s)','FontAngle','i')
+0176                         box off
+0177                         ax(j+1) = ax4(1);
+0178                     case 5
+0179                         % Scatter plot of spikes, showing synchrony
+0180                         cla
+0181                         [nsp, nsp_t] = hist(spt_curr, unique(spt_curr));
+0182                         nsp_t = nsp_t(nsp>1);
+0183                         nsp = nsp(nsp>1);
+0184                         for jj = 1:length(nsp_t)
+0185                             plot( nsp_t(jj)*ones(nsp(jj),1), spc_curr(spt_curr==nsp_t(jj)), 'o-', ...
+0186                                 'Color', colord(mod(nsp(jj)-2,size(colord,1))+1,:), 'LineWidth', 2, 'MarkerSize',10)
+0187                             hold on
+0188                         end
+0189                         scatter(spt_curr,-spc_curr,20,colord(mod(spc_curr-1,size(colord,1))+1,:),'o','filled')
+0190                         hold off
+0191                         yticks=sort([1:max(spc_curr)]);
+0192                         yl = sort([max(spc_curr)+1,min(spc_curr)-1]);
+0193                         set(gca,'YTick',yticks)
+0194                         set(gca,'YTicklabel',num2str(fliplr(ICuse)'),'YLim',yl)
+0195                 end
+0196             end
+0197             formataxes
+0198             xlabel('Time (s)')
+0199             xlim(tlims)
+0200             yl = ylim;
+0201             drawnow
+0202         end
+0203         set(gcf,'Color','w','PaperPositionMode','auto')
+0204 
+0205         %%%%
+0206         % Resize plots to appropriate size
+0207         if (plottype<4)&(length(ICuse)>=3)
+0208             bigpos = get(ax(1),'Position');
+0209             aheight = 0.9*bigpos(4)/nrows;
+0210             for k=1:length(filtax)
+0211                 axpos = get(filtax(k),'Position');
+0212                 axpos(3) = aheight;
+0213                 axpos(4) = aheight;
+0214                 set(filtax(k),'Position',axpos)
+0215             end
+0216 
+0217             set(gcf,'Units','normalized')
+0218             fpos = get(gcf,'Position');
+0219             for j=1:ncols
+0220                 axpos = get(ax(j),'OuterPosition');
+0221                 filtpos = get(filtax(1+(j-1)*nrows),'Position');
+0222                 axpos(1) = filtpos(1) + filtpos(3)*1.1;
+0223                 set(ax(j),'OuterPosition',axpos,'ActivePositionProperty','outerposition')
+0224                 axpos = get(ax(j),'Position');
+0225             end
+0226             set(gcf,'Resize','on','Units','characters')
+0227         end
+0228 
+0229         for j=1:ncols
+0230             ax = [ax,axes('Position',get(ax(j),'Position'),'XAxisLocation','top','Color','none')];
+0231             if plottype==4
+0232                 xt = get(ax4(end),'XTick');
+0233             else
+0234                 xt = get(ax(j),'XTick');
+0235             end
+0236             xlim(tlims)
+0237             formataxes
+0238             set(gca,'YTick',[],'XTick',xt,'XTickLabel',num2str(xt'/dt, '%15.0f'))
+0239             xlabel('Frame number')
+0240             axes(ax(j))
+0241             box on
+0242         end
+0243         linkaxes(ax,'xy')
+0244 
+0245 
+0246         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+0247     case {'contour'}
+0248         f_pos = get(gcf,'Position');
+0249         if f_pos(4)>f_pos(3)
+0250             f_pos(4) = 0.5*f_pos(3);
+0251             set(gcf,'Position',f_pos);
+0252         end
+0253         set(gcf,'Renderer','zbuffer','RendererMode','manual')
+0254 
+0255         subplot(1,2,2)
+0256 
+0257         clf
+0258         colormap(gray)
+0259 
+0260         set(gcf,'DefaultAxesColorOrder',colord)
+0261         subplot(1,2,1)
+0262 
+0263         crange = [min(min(min(ica_filters(ICuse,:,:)))),max(max(max(ica_filters(ICuse,:,:))))];
+0264         contourlevel = crange(2) - diff(crange)*[1,1]*0.8;
+0265 
+0266         cla
+0267         if ndims(f0)==2
+0268             imagesc(f0)
+0269         else
+0270             image(f0)
+0271         end
+0272         cax = caxis;
+0273         sigsm = 1;
+0274         shading flat
+0275         hold on
+0276         for j=1:length(ICuse)
+0277             ica_filtersuse = gaussblur(squeeze(ica_filters(ICuse(j),:,:)), sigsm);
+0278             contour(ica_filtersuse, [1,1]*(mean(ica_filtersuse(:))+4*std(ica_filtersuse(:))), ...
+0279                 'Color',colord(mod(j-1,size(colord,1))+1,:),'LineWidth',2)
+0280         end
+0281         for j=1:length(ICuse)
+0282             ica_filtersuse = gaussblur(squeeze(ica_filters(ICuse(j),:,:)), sigsm);
+0283 
+0284             % Write the number at the cell center
+0285             [ypeak, xpeak] = find(ica_filtersuse == max(max(ica_filtersuse)),1);
+0286             text(xpeak,ypeak,num2str(j), 'horizontalalignment','c','verticalalignment','m','color','y')
+0287         end
+0288         hold off
+0289         caxis(cax)
+0290         formataxes
+0291         axis image tight off
+0292         title('Avg of movie, with contours of ICs')
+0293 
+0294         ax = subplot(1,2,2);
+0295         if plottype<=2
+0296             complot(ica_sig, ICuse, dt)
+0297         end
+0298         set(gca,'ColorOrder',colord)
+0299         if plottype>=2
+0300             [spcICuse,spcord] = ismember(spc,ICuse);
+0301             spc = spcord(spcICuse&(spt<=tlims(2))&(spt>=tlims(1)));
+0302             spt = spt(spcICuse&(spt<=tlims(2))&(spt>=tlims(1)));
+0303             alph = diff(ylim)/length(ICuse);
+0304             switch plottype
+0305                 case 2
+0306                     % Traces with spikes
+0307                     hold on
+0308                     scatter(spt,-alph*(spc-1)+0.4*alph,20,colord(mod(spc-1,size(colord,1))+1,:),'o','filled')
+0309                     hold off
+0310                 case 3
+0311                     % Spike raster only
+0312                     cla
+0313                     scatter(spt,-alph*(spc-1)+0.4*alph,20,colord(mod(spc-1,size(colord,1))+1,:),'o','filled')
+0314                     yticks=sort(-alph*([1:length(ICuse)]-1)+0.4*alph);
+0315                     yl = sort(-alph*[length(ICuse)-1,-1]);
+0316                     set(gca,'YTick',yticks)
+0317                     set(gca,'YTicklabel',num2str(fliplr(ICuse)'),'YLim',yl)
+0318                     set(gca,'YLim',sort(-alph*[length(ICuse)-1,-0.6]))
+0319                 case 4
+0320                     % Binned spike rate
+0321                     tvec = [ratebin/2:ratebin:size(ica_sig,2)*dt];
+0322                     binsize = ratebin*ones(size(tvec));
+0323                     binsize(end) = size(ica_sig,2)*dt-tvec(end-1)-ratebin/2;
+0324                     for j=1:length(ICuse)
+0325                         sptj = spt(spc==j);
+0326                         spn = hist(sptj,tvec);
+0327                         spr(j,:) = spn./binsize;
+0328                         sperr(j,:) = sqrt(spn)./binsize;
+0329                     end
+0330                     spn = hist(spt,tvec);
+0331                     mspr = spn./(binsize*length(ICuse));
+0332                     msperr = sqrt(spn)./(binsize*length(ICuse));
+0333 
+0334                     hold on
+0335                     plot(tvec,mspr,'k','LineWidth',3);
+0336                     patch([tvec(end:-1:1),tvec],[mspr(end:-1:1)+msperr(end:-1:1),mspr-msperr],0.5*[1,1,1],'FaceAlpha',1)
+0337                     hold off
+0338                     axis tight
+0339                     formataxes
+0340             end
+0341         end
+0342         formataxes
+0343         xlim(tlims)
+0344         xlabel('Time (s)','FontAngle','i')
+0345         if plottype<=3
+0346             ylabel('IC #','FontAngle','i')
+0347         else
+0348             ylabel('Spike rate (Hz)','FontAngle','i')
+0349         end
+0350         set(gcf,'Color','w','PaperPositionMode','auto')
+0351         set(gca,'yticklabel',num2str(fliplr([1:length(ICuse)])'))
+0352 
+0353         axes('Position',get(ax,'Position'),'XAxisLocation','top','Color','none')
+0354         xt = get(ax,'XTick');
+0355         xlim(tlims)
+0356         formataxes
+0357         set(gca,'YTick',[], ...
+0358             'XTick',xt,'XTickLabel',num2str(xt'/dt, '%15.0f'))
+0359         xlabel('Frame number')
+0360         axes(ax)
+0361         box on
+0362 
+0363 end
+0364 
+0365 %%%%%%%%%%%%%%%%%%%%%
+0366 function complot(sig, ICuse, dt)
+0367 
+0368 for i = 1:length(ICuse)
+0369     zsig(i, :) = zscore(sig(ICuse(i),:));
+0370 end
+0371 
+0372 alpha = mean(max(zsig')-min(zsig'));
+0373 if islogical(zsig)
+0374     alpha = 1.5*alpha;
+0375 end
+0376 
+0377 zsig2 = zsig;
+0378 for i = 1:size(ICuse,2)
+0379     zsig2(i,:) = zsig(i,:) - alpha*(i-1)*ones(size(zsig(1,:)));
+0380 end
+0381 
+0382 tvec = [1:size(zsig,2)]*dt;
+0383 if islogical(zsig)
+0384     plot(tvec, zsig2','LineWidth',1)
+0385 else
+0386     plot(tvec, zsig2','LineWidth',1)
+0387 end
+0388 axis tight
+0389 
+0390 set(gca,'YTick',(-size(zsig,1)+1)*alpha:alpha:0);
+0391 set(gca,'YTicklabel',fliplr(ICuse));
+0392 
+0393 
+0394 function formataxes
+0395 
+0396 set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out')
+0397 set(gcf,'Color','w','PaperPositionMode','auto')
+0398 
+0399 function fout = gaussblur(fin, smpix)
+0400 %
+0401 % Blur an image with a Gaussian kernel of s.d. smpix
+0402 %
+0403 
+0404 if ndims(fin)==2
+0405     [x,y] = meshgrid([-ceil(3*smpix):ceil(3*smpix)]);
+0406     smfilt = exp(-(x.^2+y.^2)/(2*smpix^2));
+0407     smfilt = smfilt/sum(smfilt(:));
+0408 
+0409     fout = imfilter(fin, smfilt, 'replicate', 'same');
+0410 else
+0411     [x,y] = meshgrid([-ceil(smpix):ceil(smpix)]);
+0412     smfilt = exp(-(x.^2+y.^2)/(2*smpix^2));
+0413     smfilt = smfilt/sum(smfilt(:));
+0414 
+0415     fout = imfilter(fin, smfilt, 'replicate', 'same');
+0416 end
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortPCA.html b/doc/CellsortPCA.html new file mode 100644 index 0000000..78e226e --- /dev/null +++ b/doc/CellsortPCA.html @@ -0,0 +1,401 @@ + + + + Description of CellsortPCA + + + + + + + + + +
CellSort 1.0 > CellsortPCA.m
+ + + +

CellsortPCA +

+ +

PURPOSE ^

+
[mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
+ +

SYNOPSIS ^

+
function [mixedsig, mixedfilters, CovEvals, covtrace, movm,movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
+ +

DESCRIPTION ^

+
 [mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
+
+ CELLSORT
+ Read TIFF movie data and perform singular-value decomposition (SVD)
+ dimensional reduction.
+
+ Inputs:
+   fn - movie file name. Must be in TIFF format.
+   flims - 2-element vector specifying the endpoints of the range of
+   frames to be analyzed. If empty, default is to analyze all movie
+   frames.
+   nPCs - number of principal components to be returned
+   dsamp - optional downsampling factor. If scalar, specifies temporal
+   downsampling factor. If two-element vector, entries specify temporal
+   and spatial downsampling, respectively.
+   outputdir - directory in which to store output .mat files
+   badframes - optional list of indices of movie frames to be excluded
+   from analysis
+
+ Outputs:
+   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
+   points.
+   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
+   X x Y spatial points.
+   CovEvals - largest eigenvalues of the covariance matrix
+   covtrace - trace of covariance matrix, corresponding to the sum of all
+   eigenvalues (not just the largest few)
+   movm - average of all movie time frames at each pixel
+   movtm - average of all movie pixels at each time frame, after
+   normalizing each pixel deltaF/F
+
+ Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+ Email: eran@post.harvard.edu, mschnitz@stanford.edu
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function [mixedsig, mixedfilters, CovEvals, covtrace, movm, ...
+0002     movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
+0003 % [mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
+0004 %
+0005 % CELLSORT
+0006 % Read TIFF movie data and perform singular-value decomposition (SVD)
+0007 % dimensional reduction.
+0008 %
+0009 % Inputs:
+0010 %   fn - movie file name. Must be in TIFF format.
+0011 %   flims - 2-element vector specifying the endpoints of the range of
+0012 %   frames to be analyzed. If empty, default is to analyze all movie
+0013 %   frames.
+0014 %   nPCs - number of principal components to be returned
+0015 %   dsamp - optional downsampling factor. If scalar, specifies temporal
+0016 %   downsampling factor. If two-element vector, entries specify temporal
+0017 %   and spatial downsampling, respectively.
+0018 %   outputdir - directory in which to store output .mat files
+0019 %   badframes - optional list of indices of movie frames to be excluded
+0020 %   from analysis
+0021 %
+0022 % Outputs:
+0023 %   mixedsig - N x T matrix of N temporal signal mixtures sampled at T
+0024 %   points.
+0025 %   mixedfilters - N x X x Y array of N spatial signal mixtures sampled at
+0026 %   X x Y spatial points.
+0027 %   CovEvals - largest eigenvalues of the covariance matrix
+0028 %   covtrace - trace of covariance matrix, corresponding to the sum of all
+0029 %   eigenvalues (not just the largest few)
+0030 %   movm - average of all movie time frames at each pixel
+0031 %   movtm - average of all movie pixels at each time frame, after
+0032 %   normalizing each pixel deltaF/F
+0033 %
+0034 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0035 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0036 
+0037 tic
+0038 fprintf('-------------- CellsortPCA %s: %s -------------- \n', date, fn)
+0039 
+0040 %-----------------------
+0041 % Check inputs
+0042 if isempty(dir(fn))
+0043     error('Invalid input file name.')
+0044 end
+0045 if (nargin<2)||(isempty(flims))
+0046     nt_full = tiff_frames(fn);
+0047     flims = [1,nt_full];
+0048 end
+0049 
+0050 useframes = setdiff((flims(1):flims(2)), badframes);
+0051 nt = length(useframes);
+0052 
+0053 if nargin<3 || isempty(nPCs)
+0054     nPCs = min(150, nt);
+0055 end
+0056 if nargin<4 || isempty(dsamp)
+0057     dsamp = [1,1];
+0058 end
+0059 if nargin<5 || isempty(outputdir)
+0060     outputdir = [pwd,'/cellsort_preprocessed_data/'];
+0061 end
+0062 if nargin<6
+0063     badframes = [];
+0064 end
+0065 if isempty(dir(outputdir))
+0066     mkdir(pwd, '/cellsort_preprocessed_data/')
+0067 end
+0068 if outputdir(end)~='/';
+0069     outputdir = [outputdir, '/'];
+0070 end
+0071 
+0072 % Downsampling
+0073 if length(dsamp)==1
+0074     dsamp_time = dsamp(1);
+0075     dsamp_space = 1;
+0076 else
+0077     dsamp_time = dsamp(1);
+0078     dsamp_space = dsamp(2); % Spatial downsample
+0079 end
+0080 
+0081 [fpath, fname] = fileparts(fn);
+0082 if isempty(badframes)
+0083     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)), '_', date,'.mat'];
+0084 else
+0085     fnmat = [outputdir, fname, '_',num2str(flims(1)),',',num2str(flims(2)),'_selframes_', date,'.mat'];
+0086 end
+0087 if ~isempty(dir(fnmat))
+0088     fprintf('CELLSORT: Movie %s already processed;', ...
+0089         fn)
+0090     forceload = input(' Re-load data? [0-no/1-yes] ');
+0091     if isempty(forceload) || forceload==0
+0092         load(fnmat)
+0093         return
+0094     end
+0095 end
+0096 
+0097 fncovmat = [outputdir, fname, '_cov_', num2str(flims(1)), ',', num2str(flims(2)), '_', date,'.mat'];
+0098 
+0099 [pixw,pixh] = size(imread(fn,1));
+0100 npix = pixw*pixh;
+0101 
+0102 fprintf('   %d pixels x %d time frames;', npix, nt)
+0103 if nt<npix
+0104     fprintf(' using temporal covariance matrix.\n')
+0105 else
+0106     fprintf(' using spatial covariance matrix.\n')
+0107 end
+0108 
+0109 % Create covariance matrix
+0110 if nt < npix
+0111     [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp);
+0112 else
+0113     [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp);
+0114 end
+0115 covtrace = trace(covmat) / npix;
+0116 movm = reshape(movm, pixw, pixh);
+0117 
+0118 if nt < npix
+0119     % Perform SVD on temporal covariance
+0120     [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix);
+0121 
+0122     % Load the other set of principal components
+0123     [mixedfilters] = reload_moviedata(pixw*pixh, mov, mixedsig, CovEvals);
+0124 else
+0125     % Perform SVD on spatial components
+0126     [mixedfilters, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix);
+0127 
+0128     % Load the other set of principal components
+0129     [mixedsig] = reload_moviedata(nt, mov', mixedfilters, CovEvals);
+0130 end
+0131 mixedfilters = reshape(mixedfilters, pixw,pixh,nPCs);
+0132 
+0133 firstframe_full = imread(fn,1);
+0134 firstframe = firstframe_full;
+0135 if dsamp_space>1
+0136     firstframe = imresize(firstframe, size(mov(:,:,1)),'bilinear');
+0137 end
+0138 
+0139 %------------
+0140 % Save the output data
+0141 save(fnmat,'mixedfilters','CovEvals','mixedsig', ...
+0142     'movm','movtm','covtrace')
+0143 fprintf(' CellsortPCA: saving data and exiting; ')
+0144 toc
+0145 
+0146     function [covmat, mov, movm, movtm] = create_xcov(fn, pixw, pixh, useframes, nt, dsamp)
+0147         %-----------------------
+0148         % Load movie data to compute the spatial covariance matrix
+0149 
+0150         npix = pixw*pixh;
+0151 
+0152         % Downsampling
+0153         if length(dsamp)==1
+0154             dsamp_time = dsamp(1);
+0155             dsamp_space = 1;
+0156         else
+0157             dsamp_time = dsamp(1);
+0158             dsamp_space = dsamp(2); % Spatial downsample
+0159         end
+0160 
+0161         if (dsamp_space==1)
+0162             mov = zeros(pixw, pixh, nt);
+0163             for jjind=1:length(useframes)
+0164                 jj = useframes(jjind);
+0165                 mov(:,:,jjind) = imread(fn,jj);
+0166                 if mod(jjind,500)==1
+0167                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0168                     toc
+0169                 end
+0170             end
+0171         else
+0172             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
+0173             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
+0174             for jjind=1:length(useframes)
+0175                 jj = useframes(jjind);
+0176                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
+0177                 if mod(jjind,500)==1
+0178                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0179                     toc
+0180                 end
+0181             end
+0182         end
+0183 
+0184         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0185         toc
+0186         mov = reshape(mov, npix, nt);
+0187 
+0188         % DFoF normalization of each pixel
+0189         movm = mean(mov,2); % Average over time
+0190         movmzero = (movm==0);
+0191         movm(movmzero) = 1;
+0192         mov = mov ./ (movm * ones(1,nt)) - 1; % Compute Delta F/F
+0193         mov(movmzero, :) = 0;
+0194 
+0195         if dsamp_time>1
+0196             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
+0197             mov = downsample(mov', dsamp)';
+0198         end
+0199 
+0200         movtm = mean(mov,2); % Average over space
+0201         clear movmzeros
+0202 
+0203         c1 = (mov*mov')/size(mov,2);
+0204         toc
+0205         covmat = c1 - movtm*movtm';
+0206         clear c1
+0207     end
+0208 
+0209     function [covmat, mov, movm, movtm] = create_tcov(fn, pixw, pixh, useframes, nt, dsamp)
+0210         %-----------------------
+0211         % Load movie data to compute the temporal covariance matrix
+0212         npix = pixw*pixh;
+0213 
+0214         % Downsampling
+0215         if length(dsamp)==1
+0216             dsamp_time = dsamp(1);
+0217             dsamp_space = 1;
+0218         else
+0219             dsamp_time = dsamp(1);
+0220             dsamp_space = dsamp(2); % Spatial downsample
+0221         end
+0222 
+0223         if (dsamp_space==1)
+0224             mov = zeros(pixw, pixh, nt);
+0225             for jjind=1:length(useframes)
+0226                 jj = useframes(jjind);
+0227                 mov(:,:,jjind) = imread(fn,jj);
+0228                 if mod(jjind,500)==1
+0229                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0230                     toc
+0231                 end
+0232             end
+0233         else
+0234             [pixw_dsamp,pixh_dsamp] = size(imresize( imread(fn,1), 1/dsamp_space, 'bilinear' ));
+0235             mov = zeros(pixw_dsamp, pixh_dsamp, nt);
+0236             for jjind=1:length(useframes)
+0237                 jj = useframes(jjind);
+0238                 mov(:,:,jjind) = imresize( imread(fn,jj), 1/dsamp_space, 'bilinear' );
+0239                 if mod(jjind,500)==1
+0240                     fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0241                     toc
+0242                 end
+0243             end
+0244         end
+0245 
+0246         fprintf(' Read frame %4.0f out of %4.0f; ', jjind, nt)
+0247         toc
+0248         mov = reshape(mov, npix, nt);
+0249 
+0250         % DFoF normalization of each pixel
+0251         movm = mean(mov,2); % Average over time
+0252         movmzero = (movm==0); % Avoid dividing by zero
+0253         movm(movmzero) = 1;
+0254         mov = mov ./ (movm * ones(1,nt)) - 1;
+0255         mov(movmzero, :) = 0;
+0256 
+0257         if dsamp_time>1
+0258             mov = filter(ones(dsamp,1)/dsamp, 1, mov, [], 2);
+0259             mov = downsample(mov', dsamp)';
+0260         end
+0261 
+0262         c1 = (mov'*mov)/npix;
+0263         movtm = mean(mov,1); % Average over space
+0264         covmat = c1 - movtm'*movtm;
+0265         clear c1
+0266     end
+0267 
+0268     function [mixedsig, CovEvals, percentvar] = cellsort_svd(covmat, nPCs, nt, npix)
+0269         %-----------------------
+0270         % Perform SVD
+0271 
+0272         covtrace = trace(covmat) / npix;
+0273 
+0274         opts.disp = 0;
+0275         opts.issym = 'true';
+0276         if nPCs<size(covmat,1)
+0277             [mixedsig, CovEvals] = eigs(covmat, nPCs, 'LM', opts);  % pca_mixedsig are the temporal signals, mixedsig
+0278         else
+0279             [mixedsig, CovEvals] = eig(covmat);
+0280             CovEvals = diag( sort(diag(CovEvals), 'descend'));
+0281             nPCs = size(CovEvals,1);
+0282         end
+0283         CovEvals = diag(CovEvals);
+0284         if nnz(CovEvals<=0)
+0285             nPCs = nPCs - nnz(CovEvals<=0);
+0286             fprintf(['Throwing out ',num2str(nnz(CovEvals<0)),' negative eigenvalues; new # of PCs = ',num2str(nPCs),'. \n']);
+0287             mixedsig = mixedsig(:,CovEvals>0);
+0288             CovEvals = CovEvals(CovEvals>0);
+0289         end
+0290 
+0291         mixedsig = mixedsig' * nt;
+0292         CovEvals = CovEvals / npix;
+0293 
+0294         percentvar = 100*sum(CovEvals)/covtrace;
+0295         fprintf([' First ',num2str(nPCs),' PCs contain ',num2str(percentvar,3),'%% of the variance.\n'])
+0296     end
+0297 
+0298     function [mixedfilters] = reload_moviedata(npix, mov, mixedsig, CovEvals)
+0299         %-----------------------
+0300         % Re-load movie data
+0301         nPCs = size(mixedsig,1);
+0302 
+0303         Sinv = inv(diag(CovEvals.^(1/2)));
+0304 
+0305         movtm = mean(mov,1); % Average over space
+0306         movuse = mov - ones(npix,1) * movtm;
+0307         mixedfilters = reshape(movuse * mixedsig' * Sinv, npix, nPCs);
+0308     end
+0309
+0310 
+0311     function j = tiff_frames(fn)
+0312         %
+0313         % n = tiff_frames(filename)
+0314         %
+0315         % Returns the number of slices in a TIFF stack.
+0316         %
+0317         % Modified April 9, 2013 for compatibility with MATLAB 2012b
+0318 
+0319         j = length(imfinfo(fn));
+0320     end
+0321 end
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortPlotPCspectrum.html b/doc/CellsortPlotPCspectrum.html new file mode 100644 index 0000000..d9cb181 --- /dev/null +++ b/doc/CellsortPlotPCspectrum.html @@ -0,0 +1,142 @@ + + + + Description of CellsortPlotPCspectrum + + + + + + + + + +
CellSort 1.0 > CellsortPlotPCspectrum.m
+ + + +

CellsortPlotPCspectrum +

+ +

PURPOSE ^

+
function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
+ +

SYNOPSIS ^

+
function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
+ +

DESCRIPTION ^

+
CellsortPlotPCspectrum(fn, CovEvals, pcuse)
+
+ Plot the principal component (PC) spectrum and compare with the
+ corresponding random-matrix noise floor
+
+ Inputs:
+   fn - movie file name. Must be in TIFF format.
+   CovEvals - eigenvalues of the covariance matrix
+   pcuse - [optional] - indices of PCs included in dimensionally reduced
+   data set
+
+ Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+ Email: eran@post.harvard.edu, mschnitz@stanford.edu
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + +

SUBFUNCTIONS ^

+ +

SOURCE CODE ^

+
0001 function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
+0002 % function CellsortPlotPCspectrum(fn, CovEvals, pcuse)
+0003 %
+0004 % Plot the principal component (PC) spectrum and compare with the
+0005 % corresponding random-matrix noise floor
+0006 %
+0007 % Inputs:
+0008 %   fn - movie file name. Must be in TIFF format.
+0009 %   CovEvals - eigenvalues of the covariance matrix
+0010 %   pcuse - [optional] - indices of PCs included in dimensionally reduced
+0011 %   data set
+0012 %
+0013 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0014 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0015 
+0016 if nargin<3
+0017     pcuse = [];
+0018 end
+0019 
+0020 [pixw,pixh] = size(imread(fn,1));
+0021 npix = pixw*pixh;
+0022 nt = tiff_frames(fn);
+0023 
+0024 % Random matrix prediction (Sengupta & Mitra)
+0025 p1 = npix; % Number of pixels
+0026 q1 = nt; % Number of time frames
+0027 q = max(p1,q1);
+0028 p = min(p1,q1);
+0029 sigma = 1;
+0030 lmax = sigma*sqrt(p+q + 2*sqrt(p*q));
+0031 lmin = sigma*sqrt(p+q - 2*sqrt(p*q));
+0032 lambda = [lmin: (lmax-lmin)/100.0123423421: lmax];
+0033 rho = (1./(pi*lambda*(sigma^2))).*sqrt((lmax^2-lambda.^2).*(lambda.^2-lmin^2));
+0034 rho(isnan(rho)) = 0;
+0035 rhocdf = cumsum(rho)/sum(rho);
+0036 noiseigs = interp1(rhocdf, lambda, [p:-1:1]'/p, 'linear', 'extrap').^2 ;
+0037 
+0038 % Normalize the PC spectrum
+0039 normrank = min(nt-1,length(CovEvals));
+0040 pca_norm = CovEvals*noiseigs(normrank) / (CovEvals(normrank)*noiseigs(1));
+0041 
+0042 clf
+0043 plot(pca_norm, 'o-', 'Color', [1,1,1]*0.3, 'MarkerFaceColor', [1,1,1]*0.3, 'LineWidth',2)
+0044 hold on
+0045 plot(noiseigs / noiseigs(1), 'b-', 'LineWidth',2)
+0046 plot(2*noiseigs / noiseigs(1), 'b--', 'LineWidth',2)
+0047 if ~isempty(pcuse)
+0048     plot(pcuse, pca_norm(pcuse), 'rs', 'LineWidth',2)
+0049 end
+0050 hold off
+0051 formataxes
+0052 set(gca,'XScale','log','YScale','log', 'Color','none')
+0053 xlabel('PC rank')
+0054 ylabel('Normalized variance')
+0055 axis tight
+0056 if isempty(pcuse)
+0057     legend('Data variance','Noise floor','2 x Noise floor')
+0058 else
+0059     legend('Data variance','Noise floor','2 x Noise floor','Retained PCs')
+0060 end
+0061 
+0062 fntitle = fn;
+0063 fntitle(fn=='_') = ' ';
+0064 title(fntitle)
+0065 
+0066 function formataxes
+0067 
+0068 set(gca,'FontSize',12,'FontWeight','bold','FontName','Helvetica','LineWidth',2,'TickLength',[1,1]*.02,'tickdir','out')
+0069 set(gcf,'Color','w','PaperPositionMode','auto')
+0070 
+0071
+0072
+0073     function j = tiff_frames(fn)
+0074         %
+0075         % n = tiff_frames(filename)
+0076         %
+0077         % Returns the number of slices in a TIFF stack.
+0078         %
+0079         % Modified April 9, 2013 for compatibility with MATLAB 2012b
+0080 
+0081         j = length(imfinfo(fn));
+0082 end
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/CellsortSegmentation.html b/doc/CellsortSegmentation.html new file mode 100644 index 0000000..00955c7 --- /dev/null +++ b/doc/CellsortSegmentation.html @@ -0,0 +1,217 @@ + + + + Description of CellsortSegmentation + + + + + + + + + +
CellSort 1.0 > CellsortSegmentation.m
+ + + +

CellsortSegmentation +

+ +

PURPOSE ^

+
[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
+ +

SYNOPSIS ^

+
function [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
+ +

DESCRIPTION ^

+
 [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
+
+CellsortSegmentation
+ Segment spatial filters derived by ICA
+
+ Inputs:
+     ica_filters - X x Y x nIC matrix of ICA spatial filters
+     smwidth - standard deviation of Gaussian smoothing kernel (pixels)
+     thresh - threshold for spatial filters (standard deviations)
+     arealims - 2-element vector specifying the minimum and maximum area
+     (in pixels) of segments to be retained; if only one element is
+     specified, use this as the minimum area
+     plotting - [0,1] whether or not to show filters
+
+ Outputs:
+     ica_segments - segmented spatial filters
+     segmentabel - indices of the ICA filters from which each segment was derived
+     segcentroid - X,Y centroid, in pixels, of each segment
+
+ Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009 
+ Email: eran@post.harvard.edu, mschnitz@stanford.edu
+ + +

CROSS-REFERENCE INFORMATION ^

+This function calls: +
    +
+This function is called by: +
    +
+ + + +

SOURCE CODE ^

+
0001 function [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
+0002 % [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
+0003 %
+0004 %CellsortSegmentation
+0005 % Segment spatial filters derived by ICA
+0006 %
+0007 % Inputs:
+0008 %     ica_filters - X x Y x nIC matrix of ICA spatial filters
+0009 %     smwidth - standard deviation of Gaussian smoothing kernel (pixels)
+0010 %     thresh - threshold for spatial filters (standard deviations)
+0011 %     arealims - 2-element vector specifying the minimum and maximum area
+0012 %     (in pixels) of segments to be retained; if only one element is
+0013 %     specified, use this as the minimum area
+0014 %     plotting - [0,1] whether or not to show filters
+0015 %
+0016 % Outputs:
+0017 %     ica_segments - segmented spatial filters
+0018 %     segmentabel - indices of the ICA filters from which each segment was derived
+0019 %     segcentroid - X,Y centroid, in pixels, of each segment
+0020 %
+0021 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
+0022 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
+0023 
+0024 tic
+0025 fprintf('-------------- CellsortSegmentation %s -------------- \n', date)
+0026 
+0027 if (nargin<3)||isempty(thresh)
+0028     thresh = 2;
+0029 end
+0030 if (nargin<4)||isempty(arealims)
+0031     arealims = 200;
+0032 end
+0033 if (nargin<5)||isempty(plotting)
+0034     plotting = 0;
+0035 end
+0036 [nic,pixw,pixh] = size(ica_filters);
+0037 
+0038 ica_filtersorig = ica_filters / abs(std(ica_filters(:)));
+0039 ica_filters = (ica_filters - mean(ica_filters(:)))/abs(std(ica_filters(:)));
+0040 if smwidth>0
+0041     % Smooth mixing filter with a Gaussian of s.d. smwidth pixels
+0042     smrange = max(5,3*smwidth);
+0043     [x,y] = meshgrid([-smrange:smrange]);
+0044 
+0045     smy = 1; smx = 1;
+0046     ica_filtersfilt = exp(-((x/smx).^2 + (y/smy).^2)/(2*smwidth^2));
+0047     
+0048     ica_filtersfilt = ica_filtersfilt/sum(ica_filtersfilt(:));
+0049     ica_filtersbw = false(pixw,pixh,nic);
+0050     tic
+0051     for j = 1:size(ica_filters,1)
+0052         ica_filtersuse = ica_filters(j,:,:);
+0053         ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:)));
+0054         ica_filtersbw(:,:,j) = (imfilter(ica_filtersuse, ica_filtersfilt, 'replicate', 'same') > thresh);
+0055     end
+0056 else
+0057     ica_filtersbw = (permute(ica_filters,[2,3,1]) > thresh);
+0058     ica_filtersfilt = 1;
+0059 end
+0060 
+0061 tic
+0062 if plotting
+0063     clf
+0064     set(gcf,'Color','w')
+0065     colormap(gray)
+0066     
+0067     subplot(223)
+0068     imagesc(squeeze(sum(ica_filters,1)))
+0069     axis image off
+0070     hold on
+0071 end
+0072 ica_filterslabel = [];
+0073 ica_segments = [];
+0074 k=0;
+0075 L=[];
+0076 segmentlabel = [];
+0077 segcentroid = [];
+0078 [x,y] = meshgrid([1:pixh], [1:pixw]);
+0079 for j = 1:nic
+0080     % Label contiguous components
+0081     L = bwlabel(ica_filtersbw(:,:,j), 4);
+0082     Lu = 1:max(L(:));
+0083 
+0084     % Delete small components
+0085     Larea = struct2array(regionprops(L, 'area'));
+0086     Lcent = regionprops(L, 'Centroid');
+0087     
+0088     if length(arealims)==2
+0089         Lbig = Lu( (Larea >= arealims(1))&(Larea <= arealims(2)));
+0090         Lsmall = Lu((Larea < arealims(1))|(Larea > arealims(2)));
+0091     else
+0092         Lbig = Lu(Larea >= arealims(1));
+0093         Lsmall = Lu(Larea < arealims(1));
+0094     end
+0095     
+0096     L(ismember(L,Lsmall)) = 0;   
+0097         
+0098     for jj = 1:length(Lbig)
+0099         segcentroid(jj+k,:) = Lcent(Lbig(jj)).Centroid;
+0100     end
+0101     
+0102     ica_filtersuse = squeeze(ica_filtersorig(j,:,:));
+0103     for jj = 1:length(Lbig)
+0104         ica_segments(jj+k,:,:) = ica_filtersuse .* ( 0*(L==0) + (L==Lbig(jj)) );  % Exclude background
+0105     end
+0106     
+0107     if plotting && ~isempty(Lbig)
+0108         if smwidth>0
+0109             subplot(2,2,2)
+0110             ica_filtersuse = squeeze(ica_filters(j,:,:));
+0111             ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:)));
+0112             imagesc(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'),[-1,4])
+0113             hold on
+0114             contour(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'), [1,1]*thresh, 'k')
+0115             hold off
+0116             hc = colorbar('Position',[0.9189    0.6331    0.0331    0.2253]);
+0117             ylabel(hc,'Std. dev.')
+0118             title(['IC ',num2str(j),' smoothed'])
+0119             axis image off
+0120 
+0121             subplot(2,2,1)
+0122         else
+0123             subplot(211)
+0124         end
+0125         imagesc(squeeze(ica_filters(j,:,:)))
+0126         title(['IC ',num2str(j),' original'])
+0127         axis image off
+0128 
+0129         colord = lines(k+length(Lbig));
+0130         for jj = 1:length(Lbig)
+0131             subplot(223)
+0132             contour(ica_filtersbw(:,:,j), [1,1]*0.5, 'color',colord(jj+k,:),'linewidth',2)
+0133             hold on
+0134             text(segcentroid(jj+k,1), segcentroid(jj+k,2), num2str(jj+k), 'horizontalalignment','c', 'verticalalignment','m')
+0135             set(gca, 'ydir','reverse','tickdir','out')
+0136             axis image
+0137             xlim([0,pixw]); ylim([0,pixh])
+0138             
+0139             subplot(224)
+0140             imagesc(squeeze(ica_segments(jj+k,:,:)))
+0141             hold on
+0142             plot(segcentroid(jj+k,1), segcentroid(jj+k,2), 'bo')
+0143             hold off
+0144             axis image off
+0145             title(['Segment ',num2str(jj+k)])
+0146             drawnow
+0147         end
+0148     end
+0149     k = size(ica_segments,1);
+0150 end
+0151 toc
+
Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/images/Linking to Harvard Library E-Resources - Harvard College Library.webarchive b/doc/images/Linking to Harvard Library E-Resources - Harvard College Library.webarchive new file mode 100644 index 0000000..2d53eeb Binary files /dev/null and b/doc/images/Linking to Harvard Library E-Resources - Harvard College Library.webarchive differ diff --git a/doc/images/Mukamel1a.jpg b/doc/images/Mukamel1a.jpg new file mode 100644 index 0000000..1178fc0 Binary files /dev/null and b/doc/images/Mukamel1a.jpg differ diff --git a/doc/images/alpha.png b/doc/images/alpha.png new file mode 100644 index 0000000..c73de7b Binary files /dev/null and b/doc/images/alpha.png differ diff --git a/doc/images/c++.png b/doc/images/c++.png new file mode 100644 index 0000000..24f56e6 Binary files /dev/null and b/doc/images/c++.png differ diff --git a/doc/images/c.png b/doc/images/c.png new file mode 100644 index 0000000..c39fbf0 Binary files /dev/null and b/doc/images/c.png differ diff --git a/doc/images/demoicon.gif b/doc/images/demoicon.gif new file mode 100644 index 0000000..c89e7ac Binary files /dev/null and b/doc/images/demoicon.gif differ diff --git a/doc/images/down.png b/doc/images/down.png new file mode 100644 index 0000000..d41104a Binary files /dev/null and b/doc/images/down.png differ diff --git a/doc/images/fortran.png b/doc/images/fortran.png new file mode 100644 index 0000000..350c572 Binary files /dev/null and b/doc/images/fortran.png differ diff --git a/doc/images/hp.png b/doc/images/hp.png new file mode 100644 index 0000000..d09f988 Binary files /dev/null and b/doc/images/hp.png differ diff --git a/doc/images/left.png b/doc/images/left.png new file mode 100644 index 0000000..404df04 Binary files /dev/null and b/doc/images/left.png differ diff --git a/doc/images/linux.png b/doc/images/linux.png new file mode 100644 index 0000000..42c0c32 Binary files /dev/null and b/doc/images/linux.png differ diff --git a/doc/images/matlabicon.gif b/doc/images/matlabicon.gif new file mode 100644 index 0000000..7c2b5e6 Binary files /dev/null and b/doc/images/matlabicon.gif differ diff --git a/doc/images/mex.png b/doc/images/mex.png new file mode 100644 index 0000000..396f1bc Binary files /dev/null and b/doc/images/mex.png differ diff --git a/doc/images/right.png b/doc/images/right.png new file mode 100644 index 0000000..067c5ba Binary files /dev/null and b/doc/images/right.png differ diff --git a/doc/images/sgi.png b/doc/images/sgi.png new file mode 100644 index 0000000..20052bc Binary files /dev/null and b/doc/images/sgi.png differ diff --git a/doc/images/simulinkicon.gif b/doc/images/simulinkicon.gif new file mode 100644 index 0000000..1386ddd Binary files /dev/null and b/doc/images/simulinkicon.gif differ diff --git a/doc/images/solaris.png b/doc/images/solaris.png new file mode 100644 index 0000000..e31e8a2 Binary files /dev/null and b/doc/images/solaris.png differ diff --git a/doc/images/up.png b/doc/images/up.png new file mode 100644 index 0000000..b348b9a Binary files /dev/null and b/doc/images/up.png differ diff --git a/doc/images/windows.png b/doc/images/windows.png new file mode 100644 index 0000000..6cec95b Binary files /dev/null and b/doc/images/windows.png differ diff --git a/doc/index.html b/doc/index.html new file mode 100644 index 0000000..06c5e41 --- /dev/null +++ b/doc/index.html @@ -0,0 +1,52 @@ + + + + Index for Directory CellSort 1.0 + + + + + + + + + + +
Index for CellSort 1.0 >
+ + +

Index for CellSort 1.0

+ + + + + + + +

Matlab files in this directory:

+
CellSort is a MATLAB toolbox containing code that accompanies the manuscript, Automated analysis of cellular signals from large-scale calcium imaging data by Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, NEURON (2009). Please address comments and questions to Eran Mukamel (eran at post dot harvard dot edu) or Mark Schnitzer (mschnitz at stanford dot edu). +
+ + + + + + + + + + + + + + + + +
  Step 1: CellsortPCA[mixedsig, mixedfilters, CovEvals, covtrace, movm, movtm] = CellsortPCA(fn, flims, nPCs, dsamp, outputdir, badframes)
  Step 2a: CellsortChoosePCs[PCuse] = CellsortChoosePCs(fn, mixedfilters)
  Step 2b: CellsortPlotPCspectrumCellsortPlotPCspectrum(fn, CovEvals, pcuse)
  Step 3a: CellsortICA[ica_sig, ica_filters, ica_A, numiter] = CellsortICA(mixedsig, mixedfilters, PCuse, mu, nIC, ica_A_guess, termtol, maxrounds)
  Step 3b: CellsortICAplotCellsortICAplot(mode, ica_filters, ica_sig, f0, tlims, dt, ratebin, plottype, ICuse, spt, spc)
  Step 4a: CellsortSegmentation[ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
  Step 4b:CellsortApplyFiltercell_sig = CellsortApplyFilter(fn, ica_segments, flims, movm, subtractmean)
  Step 5:CellsortFindspikes[spmat, spt, spc, zsig] = CellsortFindspikes(ica_sig, thresh, dt, deconvtau, normalization)
+ +

+ +

Generated on Wed 29-Jul-2009 12:46:53 by m2html © 2003
+ + \ No newline at end of file diff --git a/doc/m2html.css b/doc/m2html.css new file mode 100644 index 0000000..9d88fb8 --- /dev/null +++ b/doc/m2html.css @@ -0,0 +1,75 @@ +body { + background: white; + color: black; + font-family: arial,sans-serif; + margin: 0; + padding: 1ex; +} + +div.fragment { + width: 98%; + border: 1px solid #CCCCCC; + background-color: #f5f5f5; + padding-left: 4px; + margin: 4px; +} + +div.box { + width: 98%; + background-color: #f5f5f5; + border: 1px solid #CCCCCC; + color: black; + padding: 4px; +} + +.comment { + color: #228B22; +} +.string { + color: #B20000; +} +.keyword { + color: #0000FF; +} + +.keywordtype { color: #604020; } +.keywordflow { color: #e08000; } +.preprocessor { color: #806020; } +.stringliteral { color: #002080; } +.charliteral { color: #008080; } + +a { + text-decoration: none; +} + +a:hover { + background-color: #006699; + color:#FFFFFF; +} + +a.code { + font-weight: normal; + color: #A020F0; +} + +a.code:hover { + background-color: #FF0000; + color: #FFFFFF; +} + +h1 { + background: transparent; + color: #006699; + font-size: x-large; + text-align: center; +} + +h2 { + background: transparent; + color: #006699; + font-size: large; +} + +address { + font-size:small; +}