Description of CellsortICAplot

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
CellsortICAplot

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
