123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384 |
- function Analyze_Spatial_Contrast_Adaptation_Stimulus(varargin)
- %
- %
- %%% Analyze_Spatial_Contrast_Adaptation_Stimulus %%%
- %
- % This function analyzes the response of the ganglion cells to the spatial
- % contrast adaptation stimulus. This is a supplement code for the Khani & Gollisch
- % (2017) study about the diversity of spatial contrast adaptation in the mouse retina.
- % This function accompanys dataset of Khani and Gollisch (2017). The dataset is available at:
- % https://gin.g-node.org/gollischlab/Khani_and_Gollisch_2017_RGC_spiketrains_for_spatial_contrast_adaptation
- % This function first loads the data in a format that is provided in by the
- % supplementary database of the study. It then generate PSTHs, spike-triggered
- % average (STA) and nonlinearity functions for the selected cell.
- %
- % ===============================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====================================
- %
- % PSTH, STAs and nonlinearities for the selected cell.
- %
- %--------------------------------------------------------------------------------------------------%
- %---------- main-function ----------%
- %--------------------------------------------------------------------------------------------------%
- % first load the data for an example cell
- [dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin{:});
- % analyze the data
- ca_data = sca_anylsis(dp, stimulusnumber, cellnumber, clusternumber);
- % plotting
- plot_sca(ca_data, cellnumber, clusternumber);
- end
- %--------------------------------------------------------------------------------------------------%
- %---------- sub-functions ----------%
- %--------------------------------------------------------------------------------------------------%
- function expdata = sca_anylsis(dp, stimulusnum, cellnum, clusternum,varargin)
- %
- expdata = load_data_for_selected_cell(dp, stimulusnum, cellnum, clusternum);
- % set all the required parameters
- para = expdata.stimpara;
- para.refreshrate = 60;
- para.upsampleratio = ceil(1e3/para.refreshrate);
- para.numNLbins = 40;
- para.meanIntensity = 0.5;
- para.statime = fliplr(linspace(0,(1e3/60*para.Nblinks)*(para.numNLbins/para.upsampleratio),para.numNLbins));
- Nruns = floor((length(expdata.frametimes)-1)/(para.Nblinks*para.switchFrames));
- para.Nframes = Nruns*para.Nblinks*para.switchFrames;
- para.Nrepeats = floor(Nruns/2); % number of runs per trial;
- para.binlength = 250 / 1e3; % 250 ms each bin
- % for the PSTH
- para.low_contduration = para.switchFrames ./ para.refreshrate;
- para.high_contduration = para.switchFrames ./ para.refreshrate;
- para.psthtime_low = linspace(0, para.low_contduration ,(para.low_contduration / para.binlength)+1);
- para.psthtime_high = linspace(0, para.high_contduration ,(para.high_contduration / para.binlength)+1);
- para.trialduration = para.low_contduration + para.high_contduration;
- para.psthtimelh = linspace(0, para.trialduration,(para.trialduration / para.binlength));
- % loading clusters and frametimings
- ft = expdata.frametimes; % work in seconds
- ftimes = ft(1:para.Nframes);
- % get frame times per trial
- ftimes = reshape(ftimes,para.switchFrames,Nruns)';
- ftimes_high = ftimes(2:2:end,:); ftimes_high = ftimes_high(1:para.Nrepeats,:);
- ftimes_low = ftimes(1:2:end,:); ftimes_low = ftimes_low(1:para.Nrepeats,:);
- % getting the stimulus right!
- stimulus = Spatial_Contrast_Adaptation_Stimulus(para.seed,para.Nframes,para.switchFrames,Nruns, para.meanIntensity);
- % binning spikes
- spkbins = (histc(expdata.spiketimes,ft(1:para.Nframes))); %#ok
- spkbins = reshape(spkbins,para.switchFrames,Nruns)';
- spkbins_high = spkbins(2:2:end,:); spkbins_high = spkbins_high(1:para.Nrepeats,:);
- spkbins_low = spkbins(1:2:end,:); spkbins_low = spkbins_low(1:para.Nrepeats,:);
- % setting output
- expdata.ftimes.high = ftimes_high;
- expdata.ftimes.low = ftimes_low;
- expdata.stimulus = stimulus;
- expdata.stimpara = para;
- % calculating rasters and psth
- expdata = calc_rasters_psth(expdata);
- % calculating STAs and nonlinearities for each contrast and location
- [locX_high.sta,locX_high.nlx,locX_high.nly,locX_high.time] = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locX_high, para);
- [locX_low.sta,locX_low.nlx,locX_low.nly,locX_low.time] = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locX_low, para);
- [locY_high.sta,locY_high.nlx,locY_high.nly,locY_high.time] = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locY_high, para);
- [locY_low.sta,locY_low.nlx,locY_low.nly,locY_low.time] = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locY_low, para);
- % setting output
- expdata.locXhigh = locX_high;
- expdata.locX_low = locX_low;
- expdata.locY_high = locY_high;
- expdata.locY_low = locY_low;
- end
- %--------------------------------------------------------------------------------------------------%
- function expdata = calc_rasters_psth(expdata)
- fthigh = expdata.ftimes.high;
- ftlow = expdata.ftimes.low;
- spks = expdata.spiketimes;
- para = expdata.stimpara;
- [ras_high, ras_low] = deal(cell(para.Nrepeats,1));
- delay = 50/1e3; % ignore first spikes after contrast switch
- for ii = 1:para.Nrepeats
-
- ras_high{ii} = spks(and(spks >= fthigh(ii,1),spks <= fthigh(ii,end))) - fthigh(ii,1);
- ras_low{ii} = spks(and(spks >= (ftlow(ii,1)+delay),spks <= ftlow(ii,end))) - (ftlow(ii,1)+delay);
-
- end
- rasters_high = CelltoMatUE(ras_high);
- rasters_low = CelltoMatUE(ras_low);
- psth_high = mean(histc(rasters_high',para.psthtime_high),2)*(1/para.binlength); %#ok
- psth_high = psth_high(1:end-1); psth_high(end) = psth_high(end-1);
- para.psthtime_high = para.psthtime_high(1:end-1);
- psth_low = mean(histc(rasters_low',para.psthtime_low),2)*(1/para.binlength); %#ok
- psth_low = psth_low(1:end-1); psth_low(end) = psth_low(end-1);
- para.psthtime_low = para.psthtime_low(1:end-1);
- expdata.rasters_high = rasters_high;
- expdata.rasters_low = rasters_low;
- expdata.psth_high = psth_high';
- expdata.psth_low = psth_low';
- expdata.psthlh = [psth_low; psth_high]';
- end
- function [staNorm,nlx,nly,timeVec] = calculate_sta_nonlinearity(runningbin, stimulus, para)
- toHzratio = para.Nblinks/para.refreshrate; % normalize the response to Hz
- rnstimEn = zeros(para.Nrepeats,para.numNLbins, para.switchFrames-para.numNLbins+1);
- for itrial = 1 : para.Nrepeats
- stimtrial = stimulus(itrial,:);
- hankelMat = hankel(1:para.numNLbins,para.numNLbins:para.switchFrames);
- rnstimEn(itrial,:,:) = stimtrial(hankelMat);
- end
- runningbin = reshape(permute(runningbin(:,para.numNLbins:end,:),[1 3 2]),1,[]);
- rnstimEn = reshape(permute(rnstimEn, [2 1 3]), para.numNLbins,[]);
- sta = bsxfun(@rdivide,runningbin*rnstimEn',sum(runningbin,2));
- staNorm = bsxfun(@times,sta,1./sqrt(sum(sta.^2,2)));
- timeVec = (-(para.numNLbins-1):1:0)*toHzratio; %in seconds
- runningGenerators = staNorm * rnstimEn;
- [nly,nlx] = calculate_nonlinearity(runningGenerators', runningbin', para.numNLbins, toHzratio);
- end
- %--------------------------------------------------------------------------------------------------%
- function [ values, centers,sevalues] = calculate_nonlinearity(linearOutput, spikes, nbins, dt)
- qtStep = 1/nbins;
- qtEdges = 0:qtStep:1;
- genEdges = quantile(linearOutput, qtEdges);
- if all(isnan(genEdges)) % this is for cases without spikes where everything is NaNs
- values = zeros(size(genEdges));
- sevalues = zeros(size(genEdges));
- else
- [a, spikebins] = histc(linearOutput,genEdges); %#ok
-
- values = accumarray(spikebins, spikes', [nbins+1 1], @sum);
- sdvalues = accumarray(spikebins, spikes', [nbins+1 1], @std);
-
- values = values./a/dt;
- sevalues = sdvalues./sqrt(a)/dt;
-
- end
- values(end) = [];
- sevalues(end) = [];
- centers = genEdges(1:nbins)+diff(genEdges)/2;
- values = transpose(values);
- sevalues = transpose(sevalues);
- 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']);
- canames = stimnames(~(cellfun('isempty',(strfind(stimnames,'Contrast_Adaptation')))));
-
- if numel(canames) > 1
- idx = listdlg('PromptString',{'Select an example dataset:'},...
- 'SelectionMode','single','ListString',canames,'listsize',[450 100]);
- canames = canames{idx};
- end
- stimulusnumber = str2double(extractBefore(canames,'_'));
- 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 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_sca(ca_data, cellnum, clusternum)
- % putting data together to make plotting easier
- para = ca_data.stimpara;
- sta = [ca_data.locXhigh.sta;ca_data.locX_low.sta;ca_data.locY_high.sta;ca_data.locY_low.sta]';
- statime = - [ca_data.locXhigh.time;ca_data.locX_low.time;ca_data.locY_high.time;ca_data.locY_low.time]'*1e3;
- nlx = [ca_data.locXhigh.nlx;ca_data.locX_low.nlx;ca_data.locY_high.nlx;ca_data.locY_low.nlx]';
- nly = [ca_data.locXhigh.nly;ca_data.locX_low.nly;ca_data.locY_high.nly;ca_data.locY_low.nly]';
- figure('pos',[600 50 900 950],'color','w');
- % get plots title and legend labels
- titr = {'Location X', 'Location Y'};
- if para.twofixed
- l = {'high-low','low-low'};
- s = 'Spatial contrast adaptation';
- else
- l = {'high-low','low-high'};
- s = 'Spatial contrast adaptation (alternating stimulus)';
- end
- % first plotting psths
- subplot(3,2,[1,2])
- bar(para.low_contduration+para.psthtime_high(1:end-1),ca_data.psth_high);
- hold on;
- bar(para.psthtime_low(1:end-1),ca_data.psth_low);
- pbaspect([3 1 1]); box off;
- xlabel('Time (sec)'); ylabel('Rate (Hz)');
- xlim([0 para.low_contduration+para.high_contduration]);
- xticks([0 para.low_contduration para.low_contduration+para.high_contduration]);
- yAx = get(gca,'YLim'); yAx = ceil(yAx(2)/2)*2; yticks([0 yAx/2 yAx]);
- legend(l,'Location','northwest'); legend boxoff;
- % plotting STAs and nonlinearities
- for ii = 1:2
-
- subplot(3,2,ii+2)
- plot(statime(:,2*(ii-1)+1),sta(:,2*(ii-1)+1),statime(:,2*ii),sta(:,2*ii),'LineWidth',2);
- pbaspect([1.2 1 1]); box off;
- axis([0 600 -0.5 0.5]);
- title(['STA ',titr{ii}]);
- xlabel('Delay (sec)'); ylabel('Filter (a.u.)');
- if ii==2, legend(l,'Location','southeast'); legend boxoff; end
-
- subplot(3,2,ii+4)
- plot(nlx(:,2*(ii-1)+1),nly(:,2*(ii-1)+1),nlx(:,2*ii),nly(:,2*ii),'LineWidth',2);
- pbaspect([1.2 1 1]); box off; ax = gca;
- yAx = ceil(ax.YLim(2)/2)*2;
- axis([-3.1 3.1 0 yAx]);
- title(['Nonlinearity ',titr{ii}]);
- xlabel('Input'); ylabel('\Delta Rate (Hz)');
- %ax.XAxisLocation = 'origin'; ax.YAxisLocation = 'origin';
- yticks(0:yAx/2:yAx);
- xticks(-3:1:3);
-
- end
- % mega title on top
- annotation('textbox', [0.1, 0.95, 0.8, 0], 'string', [s, ' properties of cell ',...
- num2str(cellnum),', cluster ',num2str(clusternum)],'FontSize',16,'EdgeColor','none',...
- 'HorizontalAlignment','center','VerticalAlignment','middle');
- end
- %--------------------------------------------------------------------------------------------------%
|