123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649 |
- function Analyze_Chromatic_Integration_Stimulus(varargin)
- %
- %
- %%% Analyze_Chromatic_Integration_Stimulus %%%
- %
- %
- % This function analyzes the response of the ganglion cells to the chromatic
- % integration stimulus. This is a supplement code for the Khani & Gollisch
- % (2021) study of chromatic integration in the mouse retina.
- % This function first loads the data in a format that is provided in by the
- % supplementary database of the study. It then generate rasters and PSTHs for
- % all the 22 contrast combinations that was shown to each retina. The order
- % of the contrasts were randomized and it is recreated using Fisher-Yates
- % random permutation algorithem. Check fisher_Yates_shuffle_order and
- % reorder_contrast_spikes functions below for more information.
- % The contrast values of both colors range from -20 to 20 in steps of 2% change.
- % They are categorized in two groups of Green-On-UV-Off and Green-Off-UV-On.
- %-------------------------------Green-On-UV-Off-----------------------------
- % Green => 20 18 16 14 12 10 8 6 4 2 0
- % UV => 0 -2 -4 -6 -8 -10 -12 -14 -16 -18 -20
- %-------------------------------Green-Off-UV-On-----------------------------
- % Green => -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0
- % UV => 0 2 4 6 8 10 12 14 16 18 20
- %
- %
- % ===============================Inputs====================================
- %
- % datapath : path to the folder in the database.
- % stimulusnumber : number to select correct stimulus.
- % cellnumber : an example cell number from the list of good cells.
- % clusternumber : a cluster number for the selected cell.
- %
- %================================Output====================================
- %
- % chromatic integration curves and PSTHs for the selected cell.
- %
- % written by Mohammad, 12.02.2021
- %--------------------------------------------------------------------------------------------------%
- %---------- main-function ----------%
- %--------------------------------------------------------------------------------------------------%
- % first load the data for an example cell
- [dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin{:});
- % get the parameters for the selected example cell
- expdata = ci_parameters(dp, stimulusnumber, cellnumber, clusternumber);
- % calculate the rasters and PSTH for the stimulus and the gray frame before the stimulus
- ci_data = ci_analysis(expdata);
- % plot chromatic integration curves and all the psths
- plot_chromatic_integration(ci_data, cellnumber, clusternumber, expdata.stimpara);
- end
- %--------------------------------------------------------------------------------------------------%
- %---------- sub-functions ----------%
- %--------------------------------------------------------------------------------------------------%
- function expdata = ci_parameters(dp, stimulusnum, cellnum, clusternum,varargin)
- %
- expdata = load_data_for_selected_cell(dp, stimulusnum, cellnum, clusternum);
- para = expdata.stimpara;
- % loading clusters and frametimings
- ft = expdata.frametimes; % work in seconds
- % getting contrasts of the experiment (-20:2:20)
- contgr = round([para.mincontrast:para.contrastdiff:0,0:para.contrastdiff:para.maxcontrast]*100);
- contuv = round([0:para.contrastdiff:para.maxcontrast,para.mincontrast:para.contrastdiff:0]*100);
- % re-ordering stimulus for plotting and presentation purpose, here we start
- % with on green and reverse off green instead of off to on green as in stimulus program.
- para.contrastIdx = [length(contgr) : -1 : 1 + length(contgr)/2, 1 : length(contuv)/2];
- para.contrastgreen = contgr(para.contrastIdx);
- para.contrastuv = contuv(para.contrastIdx);
- % generating contrast values
- para.numtrials = floor( length( ft(2 : 2 : end)) / length(contgr) );
- seed = para.seed;
- stimorder = nan(para.numtrials,length(contgr));
- % getting the order and location of each contrast from the psudo random sequence
- for ii = 1 : para.numtrials
- [stimorder(ii,:),seed,~] = fisher_Yates_shuffle_order(seed,contgr);
- end
- % experiment parameters
- refreshrate = 60;
- para.binlength = 10 / 1e3; % 10 ms each bin
- para.delay = 0; % 0 ms delay
- para.activeregion = [0.025, 0.25]; % region used to measure response (50 to 250 ms)
- para.refreshrate = refreshrate;
- % make a bin vector for stimulus, used for making PSTH
- para.stimduration = (para.stimduration / refreshrate); % in sec
- stimBin = linspace(0, para.stimduration, (para.stimduration / para.binlength)+1); % +1 here to avoid issues with histc last bin.
- % make a bin vector for preframe, used for making PSTH
- para.pfrduration = (para.preframes / refreshrate); % in sec
- pfrBin = linspace(0, para.pfrduration,(para.pfrduration / para.binlength) + 1);
- % seting output
- expdata.stimorder = stimorder;
- expdata.stimulusbinvector = stimBin;
- expdata.preframebinvector = pfrBin;
- expdata.stimpara = para;
- end
- %--------------------------------------------------------------------------------------------------%
- function ci_data = ci_analysis(expdata)
- %
- ft = expdata.frametimes;
- spks = expdata.spiketimes;
- para = expdata.stimpara;
- % for whole experiment generally has 2000 ms preframe and 500 ms stimulus
- [~, allpreframespikes] = spikes_per_frametimes(ft, spks, para.stimduration, para.pfrduration,para);
- % this is to get rasters and psth for the graybackground before the stimulus
- [ci_data.preframe.allpfrRasters, ci_data.preframe.allpfrPSTH] = reorder_contrast_spikes(allpreframespikes,...
- expdata.stimorder,expdata.preframebinvector,para);
- %%% for 500 ms before and after stimulus
- [stimulusspikes, preframespikes] = spikes_per_frametimes(ft,spks,para.stimduration,para.stimduration,para);
- % rasters and pstrh for the 500 ms stimulus
- [ci_data.stimulusRasters, ci_data.stimulusPSTH] = reorder_contrast_spikes(stimulusspikes,...
- expdata.stimorder,expdata.stimulusbinvector,para);
- % rasters and pstrh for the 500 ms gray background before the onset of the stimulus
- [ci_data.preframe.preframeRasters, ci_data.preframe.preframePSTH] = reorder_contrast_spikes(preframespikes,...
- expdata.stimorder,expdata.stimulusbinvector,para);
- ci_data.preframePSTH = ci_data.preframe.preframePSTH;
- end
- %--------------------------------------------------------------------------------------------------%
- function varargout = fisher_Yates_shuffle_order(seed,inputVec,varargin)
- %
- %%% fisher_Yates_shuffle_order %%%
- %
- %
- % This function generate psudorandom permution similar to randperm in MATLAB
- % but works with psudorandom number generator ran1. It also gives back the
- % most recent seed value to continue the permuation in case of repeated trials.
- % note that the direction of permutation is along x-axix or for columns of
- % MATALB not for the rows.
- %
- %
- % ===============================Inputs====================================
- %
- % seed : seed value for random number generation.
- % inputVec : input vector used for permutation.
- %
- %================================Output====================================
- %
- % testOrder : vector of permuted indices for the inputVec.
- % newSeed : the recent seed that used in the ran1 function.
- % outputVec : the permuted input vector along x-axis
- %
- % Note that the permution algorithem is based on Fisher-Yates shuffle
- % algorithem identical to what is used in the stimulus program.
- % for more info check :
- % https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle
- %
- % written by Mohammad, 01.02.2016
- newSeed = seed;
- testOrder = zeros(1,size(inputVec,2));
- testOrder(1) = 1;
- for i = 2:length(inputVec)-1
-
- [randVal,newSeed] = ran1(newSeed);
-
- j = ceil(i*randVal); % based on Fischer-Yates algorithem
- testOrder(i) = testOrder(j);
- testOrder(j) = i;
-
- end
- testOrder = [testOrder(end),testOrder(1:end-1)]+1; % to match MATLAB indexing
- varargout{1} = testOrder;
- varargout{2} = newSeed;
- for j = 1:size(inputVec,1)
- varargout{3}(j,:) = inputVec(j,testOrder);
- end
- end
- %--------------------------------------------------------------------------------------------------%
- function [stimSpk, pfrSpk] = spikes_per_frametimes(ft,spk,stimdur,pfrdur,para,varargin)
- %
- delay = para.delay / 1e3; % change delay time to second to match ft and spk
- [stimSpk,pfrSpk] = deal( cell( ceil( numel( ft )/ 2), 1));
- iter = 1;
- for ii = 2:2:numel(ft)
-
- stimSpk{iter} = spk(and(spk > ft(ii)+delay,spk <= ft(ii)+(stimdur+delay))) - (ft(ii)+delay);
- % before stimulus onset to get preframe
- pfrSpk{iter} = spk(and(spk > ft(ii)-(pfrdur+delay),spk <= ft(ii))) - (ft(ii)-(pfrdur+delay));
-
- if isempty(stimSpk{iter}), stimSpk{iter} = NaN; end
- if isempty(pfrSpk{iter}), pfrSpk{iter} = NaN; end
- iter = iter+1;
- end
- end
- %--------------------------------------------------------------------------------------------------%
- function varargout = reorder_contrast_spikes(spk,order,binVec,para,varargin)
- %
- neworder = reshape(order',1,[]); % reshape stimulus order
- spk = spk(1:length(neworder))'; % cutting the last unfinished stimulus series
- % pre-allocating memory
- spkorder = cell(size(order));
- contspk = cell(1,max(neworder));
- psth = nan(length(binVec),max(neworder));
- for ii = 1 : max(neworder)
-
- thisorder = neworder == ii; % indexing the orders from 1 to end
- spkorder(:,ii) = spk( thisorder ); % sorting spikes based on previous order
- contspk{ii} = CelltoMatUE( spkorder(:,ii));
- if size(contspk{ii},2) < 2
- psth(:,ii) = (histc(contspk{ii}',binVec)/size(contspk{ii},1))*(1/para.binlength); %#ok to fix situation with 1 spike
- else
- psth(:,ii) = mean(histc(contspk{ii}',binVec),2)*(1/para.binlength); %#ok must transpose contspk for histc to work properly.
- end
- end
- psth = psth(1:length(binVec)-1,:); % delete last zero bin.
- % Here we reorder the psth and spikes in a way that the opposing stimuli
- % face each other. the indices 22:-1:11 are set to 1:11, these are green on
- % to uv off. The indices 1:1:11 are set to 12 to 22, these are green off uv
- % on. for the rest of the function no index correction is made.
- % the idices for -20:2:20 are as follows
- % for first 11 or Green-ON, UV-OFF
- % Green => 20 18 16 14 12 10 8 6 4 2 0
- % UV => 0 -2 -4 -6 -8 -10 -12 -14 -16 -18 -20
- % for second 11 or UV-ON, Green-OFF
- % Green => -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0
- % UV => 0 2 4 6 8 10 12 14 16 18 20
- contidx = [length(para.contrastgreen) : -1:1 + length(para.contrastgreen)/2, 1: length(para.contrastuv)/2];
- varargout{1} = contspk(contidx);
- varargout{2} = psth(:,contidx);
- varargout{3} = spkorder(:,contidx);
- end
- %--------------------------------------------------------------------------------------------------%
- function [ci, ciang] = chromatic_integration_curves(stimpsth, pfrpsth, para, varargin)
- %
- contnum = size(stimpsth, 2) / 2;
- timevec = linspace(0, para.stimduration, para.stimduration / para.binlength);
- activtime = timevec >= para.activeregion(1) & timevec <= para.activeregion(2);
- % only between 50-250 ms is counted!
- spkdiff = stimpsth(activtime,:) - pfrpsth(1+end-sum(activtime):end,:);
- respdiff = mean(spkdiff,1); % difference in response
- respsem = nan_sem(spkdiff); % standard error of mean
- grONuvOFF = respdiff(1 : contnum);
- grOFFuvON = respdiff(1 + contnum : contnum*2);
- grONuvOFFsem = respsem(1 : contnum);
- grOFFuvONsem = respsem(1 + contnum : contnum*2);
- gNuFspk = spkdiff(:, 1:contnum);
- gFuNspk = spkdiff(:, 1+contnum : contnum*2);
- ci.grONuvOFF = grONuvOFF;
- ci.grONuvOFFsem = grONuvOFFsem;
- ci.grONuvOFFspkdiff = gNuFspk;
- ci.grOFFuvON = grOFFuvON;
- ci.grOFFuvONsem = grOFFuvONsem;
- ci.grOFFuvONspkdiff = gFuNspk;
- if nargout > 1
- gNuFangles = [para.contrastgreen(1:contnum);para.contrastuv(1:contnum)];
- gFuNangles = [para.contrastgreen(1 + contnum:contnum*2); para.contrastuv(1 + contnum:contnum*2)];
- ciang.gNuFangles = gNuFangles;
- ciang.gFuNangles = gFuNangles;
- end
- end
- %--------------------------------------------------------------------------------------------------%
- function [dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin)
- % first parse the user inputs
- p = inputParser(); % check the user options.
- p.addParameter('datapath', [], @(x) ischar(x) || isstring(x));
- p.addParameter('stimulusnumber', [], @isnumeric);
- p.addParameter('cellnumber', [], @isnumeric);
- p.addParameter('clusternumber', [], @isnumeric);
- p.addParameter('help', false, @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));
- p.parse(varargin{:});
- % defualt parameters
- defpara = p.Results;
- if isempty(defpara.datapath)
- dp = uigetdir('Select an experiment folder:');
- else
- dp = defpara.datapath;
- end
- if isempty(defpara.cellnumber) || isempty(defpara.clusternumber)
- listofgoodcells = load([dp,filesep,'list_of_good_cells.txt'],'-ascii');
- cellabels = sprintfc('cell %g, cluster %g',listofgoodcells);
- if numel(cellabels) > 1
- idx = listdlg('PromptString',{'Select an example ganglion cell:'},...
- 'SelectionMode','single','ListString',cellabels,'listsize',[250 250]);
- else
- idx = 1;
- end
- cellnumber = listofgoodcells(idx,1);
- clusternumber = listofgoodcells(idx,2);
- else
- cellnumber = defpara.cellnumber;
- clusternumber = defpara.clusternumber;
- end
- if isempty(defpara.stimulusnumber)
- stimnames = importdata([dp,filesep,'stimuli_names.txt']);
- cinames = stimnames(~(cellfun('isempty',(strfind(stimnames,'_Chromatic_Integration_')))));
- cinames = cinames((cellfun('isempty',(strfind(cinames,'Local')))));
- cinames = cinames((cellfun('isempty',(strfind(cinames,'Grating')))));
-
- if numel(cinames) > 1
- idx = listdlg('PromptString',{'Select an example dataset:'},...
- 'SelectionMode','single','ListString',cinames,'listsize',[450 100]);
- cinames = cinames{idx};
- end
- stimulusnumber = str2double(extractBefore(cinames,'_Chromatic'));
- else
- stimulusnumber = defpara.stimulusnumber;
- end
- end
- %--------------------------------------------------------------------------------------------------%
- function expdata = load_data_for_selected_cell(dp, stimnum, cellnum, clusternum)
- %
- ftfolder = [dp,filesep,'frametimes'];
- if ~exist(ftfolder,'dir'), ftfolder = dp; end
- rasfolder = [dp,filesep,'spiketimes'];
- if ~exist(rasfolder,'dir'), rasfolder = dp; end
- paramfolder = [dp,filesep,'stimulusparameters'];
- if ~exist(paramfolder,'dir'), paramfolder = dp; end
- ftnames = dir([ftfolder, filesep, num2str(stimnum,'%02g'),'*_frametimings.txt']);
- rastername = [rasfolder, filesep, sprintf('%d_SP_C%d%02d.txt', stimnum, cellnum, clusternum)];
- warning('off','all'); % this is supress the warning about long names, hopefully nothing happen in between these lines!
- expdata.spiketimes = load(rastername, '-ascii');
- expdata.frametimes = load([ftfolder,filesep,ftnames.name], '-ascii');
- stimpath = dir([paramfolder,filesep,num2str(stimnum,'%02g'),'*_parameters.txt']);
- % now reading the txt file with parameters
- fid = fopen([paramfolder,filesep,stimpath.name]);
- tline = fgetl(fid);
- while ischar(tline)
- tline = fgetl(fid);
- if tline == -1, break; end
- fn = extractBefore(tline,' = ');
- if isempty(fn), continue; end
- val = extractAfter(tline,' = ');
- % to convert to double
- if ~isnan(str2double(val))
- val = str2double(val);
- end
- % to convert to logical
- if strcmp(val,'true'), val = true; end
- if strcmp(val,'false'), val = false; end
- % store the values in a structure
- stimpara.(fn) = val;
- end
- fclose(fid);
- expdata.stimpara = stimpara;
- expdata.name = extractBefore(ftnames.name,'_frametimings.txt');
- expdata.paraname = extractBefore(stimpath.name,'_parameters.txt');
- if strcmp(expdata.name , expdata.paraname)
- expdata = rmfield(expdata,'paraname');
- end
- warning('on','all');
- end
- %--------------------------------------------------------------------------------------------------%
- function y = nan_sem(x,dim)
- %
- if isempty(x), y = NaN; return; end
- if nargin < 2, dim = find(size(x)~=1, 1 ); if isempty(dim), dim = 1; end; end
- % Find NaNs in x and nanmean(x)
- nans = isnan(x);
- count = size(x,dim) - sum(nans,dim);
- % Protect against a all NaNs in one dimension
- i = find(count==0);
- count(i) = 1;
- y = nanstd(x,[],dim)./sqrt(count);
- y(i) = i + NaN;
- end
- %--------------------------------------------------------------------------------------------------%
- function OutputMat = CelltoMatUE (input_cell)
- %
- inptsize = size(input_cell);
- maxLength = max(cellfun('length',input_cell));
- if inptsize(2) > inptsize(1)
- OutputMat = cell2mat(cellfun(@(x)cat(1,x,nan(maxLength-length(x),1)),input_cell,'un',0));
- else
- OutputMat = cell2mat(cellfun(@(x)cat(2,x',nan(1,maxLength-length(x))),input_cell,'un',0));
- end
- end
- %--------------------------------------------------------------------------------------------------%
- %---------- plotting modules ----------%
- %--------------------------------------------------------------------------------------------------%
- function plot_chromatic_integration(ci_data, cellnum, clusternum, para)
- %
- figure('pos',[400 150 1100 750],'color','w');
- % first plot chromatic integration curve
- subplot(1,2,1)
- plot_chromatic_integration_curves(ci_data, para);
- % then plot the psths
- sp = subplot(1,2,2);
- plot_chromatic_integration_psths(ci_data, para);
- % massive title for both subplots
- annotation('textbox', [0.1, 0.95, 0.8, 0], 'string', ['Chromatic integration responses for cell ',...
- num2str(cellnum),', cluster ',num2str(clusternum)],'FontSize',16,'EdgeColor','none',...
- 'HorizontalAlignment','center','VerticalAlignment','middle');
- % stretch the second plot a bit.
- sp.Position = [sp.Position(1)-0.05, sp.Position(2)-0.07 sp.Position(3:4)+0.05];
- %annotation('textbox', [0.5, 0.06, 0, 0.8], 'string',sprintf('%2g\t\n\n',1:11),'VerticalAlignment','cap');
- end
- %--------------------------------------------------------------------------------------------------%
- function [ci, datplt]= plot_chromatic_integration_curves(ci_data, para)
- %
- stimpsth = ci_data.stimulusPSTH;
- pfrpsth = ci_data.preframePSTH;
- ptch = @(x,y,ysem,col)(patch([x,fliplr(x)],[y-ysem,fliplr(y+ysem)],1,...
- 'facecolor',col,'edgecolor','none','facealpha',0.3));
- pltline = @(x,y,col)plot(x, y,'-o','color',col,'markerfacecolor',col,...
- 'markeredgecolor',col,'markersize',7,'linewidth',2);
- gNuFcol = [1 0.0781 0.5742]; % green on, uv off color
- gFuNcol = [0 0.7461 1]; % green off, uv on color
- noaxis = false;
- contnum = size(stimpsth,2)/2;
- ci = chromatic_integration_curves(stimpsth, pfrpsth, para);
- % plotting...plotting...plotting
- hold on;
- ptch(1 : contnum, ci.grONuvOFF, ci.grONuvOFFsem, gNuFcol);
- ptch(1 : contnum, ci.grOFFuvON, ci.grOFFuvONsem, gFuNcol);
- respdiff = [ci.grONuvOFF+ci.grONuvOFFsem, ci.grONuvOFF-ci.grONuvOFFsem, ...
- ci.grOFFuvON+ci.grOFFuvONsem, ci.grOFFuvON-ci.grOFFuvONsem];
- datplt.plot1 = pltline(1:contnum, ci.grONuvOFF,gNuFcol);
- datplt.plot2 = pltline(1:contnum, ci.grOFFuvON,gFuNcol);
- datplt.zeroline = plot(linspace(0.5,contnum + 0.5,contnum),zeros(1,contnum),'--k','linewidth',0.5);
- yAx = round(linspace(floor(min(respdiff)/2)*2,ceil(max(respdiff)/2)*2, 4));
- if all(yAx == 0), yAx = [0,10]; end
- if all(yAx > 0), yAx = [0,yAx]; end
- if any(yAx <= 0), yAx = unique(sort([0,yAx])); end
- axis([0, 1+contnum, yAx(1)-1, yAx(end)+1]); axis square; box off;
- yAxneg = yAx(yAx<0); yAxpos = yAx(yAx>0);
- if not(isempty(yAxneg))
- yAxneg = [ceil(yAxneg(1)/5)*5, floor(yAxneg(2:end)/5)*5];
- if all(yAxneg == 0), yAxneg = yAx(1); end
- end
- yAxpos = floor(yAxpos/5)*5; yAxpos = yAxpos(yAxpos~=0);
- ylab = unique([yAxneg,0,yAxpos]);
- xticks(1 : 1 : contnum); yticks(ylab);
- ax = gca; ax.TickDir = 'out';
- if noaxis
- yrn = (yAx(end)+1)- (yAx(1)-1); %#ok
- yrn = ceil((yrn/10)/5)*5;
- if not(ishold), hold on; end
- plot([1,1],[ylab(end)/4, yrn + ylab(end)/4],'k','linewidth',1);
- text(1.2, ylab(end)/4 + yrn/2, [num2str(yrn),' Hz'], 'fontsize', 8);
- box off; ax.YTick = []; ax.YColor = 'w';
- end
- xlabel('Stimulus index');
- ylabel('Rate (Hz)');
- title('Chromatic Integration Curve');
- legend([datplt.plot1, datplt.plot2],'Green-On-UV-Off','Green-Off-UV-On','location','southoutside');
- legend boxoff;
- end
- %--------------------------------------------------------------------------------------------------%
- function plot_chromatic_integration_psths(ci_data, para)
- %
- stimpsth = ci_data.stimulusPSTH;
- pfrpsth = ci_data.preframePSTH;
- gNuFcol = [1 0.0781 0.5742]; % green on, uv off color
- gFuNcol = [0 0.7461 1]; % green off, uv on color
- contnum = size(stimpsth,2)/2;
- ygap = 1.15;
- xgap = 0.2;
- pfrcol = 0.75*[1 1 1]; % silver gray color for preframes
- % re-arrange input data
- cipsth = ([pfrpsth(:,1:contnum);stimpsth(:,1:contnum);...
- pfrpsth(:,1+contnum:2*contnum);stimpsth(:,1+contnum:2*contnum)]);
-
- para.pfrduration = para.stimduration;
- pfrnumbin = para.pfrduration / para.binlength;
- stimnumbin = para.stimduration / para.binlength;
- timevec = [linspace(0,para.pfrduration,pfrnumbin);...
- linspace(para.pfrduration,para.pfrduration+para.stimduration,stimnumbin);...
- linspace(para.pfrduration+para.stimduration,2*para.pfrduration+para.stimduration,pfrnumbin);...
- linspace(2*para.pfrduration+para.stimduration,2*(para.pfrduration+para.stimduration),stimnumbin)];
- tIdx = reshape(1:length(timevec(:)),[],4)';
- ygap = ceil((max(cipsth(:))* ygap)/10)*10;
- if (ygap == 0), ygap = 10; end
- axis([-xgap, (xgap+para.stimduration+para.pfrduration)*2 ,-ygap/4-15, ygap*contnum]);
- hold on; box off;
- % add shades
- shadefun(xgap);
- % draw bars first
- for ii = 1 : contnum % start from below
-
- % plotting bars
- barfun(timevec(1,:), cipsth(tIdx(1,:),ii), ((ii-1)* ygap), pfrcol,0.5/2);
- barfun(timevec(2,:), cipsth(tIdx(2,:),ii), ((ii-1)* ygap), gNuFcol,0.5);
- barfun(xgap + timevec(3,:), cipsth(tIdx(3,:),ii), ((ii-1)* ygap),pfrcol,0.5/2);
- barfun(xgap + timevec(4,:), cipsth(tIdx(4,:),ii), ((ii-1)* ygap),gFuNcol,0.5);
- scalebarfun(xgap, ygap, para, ii);
- end
- stairsfun([timevec(1,:),timevec(2,:)], cipsth(1:size(cipsth,1)/2,:), ygap, contnum, 0, 0.5);
- stairsfun([xgap + timevec(3,:),xgap + timevec(4,:)],cipsth(1+size(cipsth,1)/2:size(cipsth,1),:),...
- ygap,contnum, 0, 0.5);
- datascalefun(xgap+0.05, ygap, para);
- % getting rid of labels
- set(gca,'xtick', [],'ytick', [],'xcolor', 'w','ycolor', 'w');
- pbaspect([1 1.6 1]);
- end
- %--------------------------------------------------------------------------------------------------%
- function stairsfun(x,y,pltYgap,contnum,drawbaseline,lw)
- %
- gaps = 0 : pltYgap : pltYgap * (contnum - 1);
- yAx = y + repmat(gaps, size(y,1),1);
- yAx = [gaps ; yAx ; yAx(end,:); gaps]; % this is to close the ending of stairs function
- xAx = [x(1), x , x(end), x(end)]; % this is to make x axis match the ending
- stairs( xAx, yAx, 'k', 'linewidth', lw);
- if drawbaseline
- if ~ishold, hold on; end
- line([xAx(1), xAx(end)],[yAx(1,:) ; yAx(1,:)], 'color', 'k', 'linewidth', lw);
- end
- end
- %--------------------------------------------------------------------------------------------------%
- function b = barfun(x,y,g,col,transP,varargin)
- %
- if iscolumn(y), y = y'; end
- x = [x;x];
- y = [y;y];
- xplt = x([2:end end]);
- yplt = y(1:end);
- b = patch([xplt, fliplr(xplt)],[yplt + g, fliplr( repmat(g, size(yplt)))], 0 ,...
- 'FaceColor', col, 'EdgeColor', 'none', 'FaceAlpha', transP);
- end
- %--------------------------------------------------------------------------------------------------%
- function shadefun(xgap, shaderegion, col)
- %
- if nargin <2, shaderegion = [0.55,0.75,1.55,1.75]; end
- if nargin <3, col = [0.9792 0.9792 0.8203]; end
- patch([shaderegion(1), shaderegion(1), shaderegion(2), shaderegion(2)],...
- [get(gca,'YLim'), fliplr(get(gca,'YLim'))], col, 'edgecolor','none');
- patch(xgap + [shaderegion(3), shaderegion(3), shaderegion(4),shaderegion(4)],...
- [get(gca,'YLim'), fliplr(get(gca,'YLim'))], col, 'edgecolor','none');
- end
- %--------------------------------------------------------------------------------------------------%
- function scalebarfun(xgap, ygap, para, iter )
- %
- gr = para.contrastgreen;
- uv = para.contrastuv;
- grcol = [0 0.8008 0];
- uvcol = [0.5391 0.1680 0.8828];
- tranp = 0.5;
- sc = (ygap*0.3)/max(gr); % scaling ratio for scale bars
- tsc = (sc*-8)/2; % text distance from border of each scale bar
- barfun([-0.03,0.015], [sc*gr(iter), sc*gr(iter)], ((iter-1)* ygap) + ygap/2, grcol, tranp);
- barfun([-0.03,0.015], [sc*uv(iter), sc*uv(iter)], (((iter-1)* ygap) + ygap/2), uvcol, tranp);
- text(-0.06,tsc+((((iter-1)* ygap)+ygap/2) + sc*gr(iter)),num2str(gr(iter)),'horiz','right','fontsize',6);
- text(-0.06,((((iter-1)* ygap)+ygap/2) + sc*uv(iter))-tsc,num2str(uv(iter)),'horiz','right','fontsize',6);
- barfun([xgap+(1.1-0.03),xgap+(1.115)],[-sc*gr(iter),-sc*gr(iter)],((iter-1)* ygap)+ygap/2,grcol,tranp);
- barfun([xgap+(1.1-0.03),xgap+(1.115)],[-sc*uv(iter),-sc*uv(iter)],(((iter-1)* ygap)+ygap/2),uvcol,tranp);
- text(xgap+(1.1-0.06),((((iter-1)* ygap)+ygap/2) - sc*gr(iter))-tsc,num2str(-gr(iter)),'horiz','right','fontsize',6);
- text(xgap+(1.1-0.06),tsc+((((iter-1)* ygap)+ygap/2) - sc*uv(iter)),num2str(-uv(iter)),'horiz','right','fontsize',6);
- end
- %--------------------------------------------------------------------------------------------------%
- function datascalefun(xgap, ygap, para)
- % adding scale bar for firing rate
- scalloc = -ygap/3.5;
- tvec = (para.stimduration + para.pfrduration) * 2;
- scalresp = ceil(ygap/50)*10;
- plot(xgap + [(0.1/tvec)+tvec tvec tvec],[scalloc scalloc (scalloc + scalresp)],'k-','linewidth',0.5);
- text(xgap + tvec + (0.05/tvec),scalloc-5,'100 ms','horiz','center','vert','top','fontsize',7);
- text(xgap + tvec - 0.02,(scalloc + scalloc/2)/3,[num2str(scalresp),' spikes'],...
- 'horiz','right','vert','middle','fontsize',7);
- end
- %--------------------------------------------------------------------------------------------------%
|