latency analysis stats

end of June 2024

Author: dshimaoka

fitPSTH_pop.m

@@ -38,7 +38,7 @@
 successTr = find(dd.successTrials==1);
 % wrongSuccTr = find(dd.successTrials==0&~isnan(dd.tOnset)&succDist>=distTh);
 % quiescentTr = find(dd.successTrials==0&~isnan(dd.tOnset)&succDist<distTh);
-wrongSuccTr = find(dd.successTrials==0&~isnan(dd.tOnset)&~isnan(dd.srt));
+wrongSuccTr = find(dd.successTrials==0&~isnan(dd.tOnset)&~isnan(dd.srt)); %ALLWAYS EMPTY
 quiescentTr = find(dd.successTrials==0&~isnan(dd.tOnset)&isnan(dd.srt));
 
 choiceOutcome = nan(nTrials,1);

fitPSTH_pop_limpredictor.m

@@ -0,0 +1,72 @@
+%return exemplery cells 2024 March
+
+[saveServer, rootFolder] = getReady();
+load(fullfile(saveServer,'param20240625.mat'),'param');
+
+animal = 'hugo';%'ollie';%'andy';% 'andy' '
+tgtModality = 'all';%'eyeSpeed';
+
+dataDir = '/mnt/syncitium/Daisuke/cuesaccade_data OBS/figPSTH_pop20231026hugo/';
+
+
+%% load gain info <> avg tgt resp [HACK]
+%from fitPSTH_pop_limpredictor.m (?)
+load(fullfile(dataDir,'fitPSTH_pop20231026hugo.mat'),...
+    'id_pop','kernel_pop','tlags');
+
+%preferred direction & amplitude of each kernel modality
+tgtRange = [0.05 0.15; 0.03 0.25; -0.1 0.1];
+prefDirOption = 0;
+[kernelPrefDir, kernelAmp] = getKernelPrefDirAmp(kernel_pop, tlags, tgtRange, param.cardinalDir, prefDirOption);
+
+
+%% load kernel info from tgt units
+% not used any more??
+% gainKernelName = fullfile(dataDir,'gainkernel_splt_all.mat');
+% if exist(gainKernelName, 'file')
+%     load(gainKernelName);
+% end
+
+%% exclude NG units
+%from pickUnitsByClass
+load(fullfile(dataDir,'pickUnitsByClass.mat'),"funcClass",'nUnits');
+
+%% tgt units
+highAmp = true;
+if highAmp
+    suffix = '_highAmp';
+    param.ampTh = .5;
+else
+    suffix = '';
+end
+for itgtModality = 1:3
+    switch itgtModality
+        case 1
+            tgtModality = 'vision';
+        case 2
+            tgtModality = 'eyeSpeed';
+        case 3
+            tgtModality = 'eyePosition';
+    end
+    for iprefDir =1%:8
+
+        prefDir_q = quantizeDir(kernelPrefDir(:,itgtModality), param.cardinalDir);
+        if ~highAmp
+            tgtID = id_pop(prefDir_q== param.cardinalDir(iprefDir));
+        else
+            tgtID = id_pop(prefDir_q== param.cardinalDir(iprefDir) & ...
+                kernelAmp(:,itgtModality) > param.ampTh);
+        end
+
+        tgtID = intersect(funcClass.id_all, tgtID);%exclude NG units
+
+        [~,tgtIDidx] = intersect(id_all, tgtID);
+
+
+        id_tgtModality{itgtModality, iprefDir} = id_pop(tgtIDidx);
+    
+
+    end
+end
+
+save(fullfile(dataDir,'exampleCells.mat'), 'id_tgtModality');

getChoiceOutcome.m

@@ -5,7 +5,7 @@
 % created from showTonsetResp.m
 
 [avgOnsetResp, winSamps, singleOnsetResp, ...
-    sortedOnsetLabels, uniqueOnsetLabels] ...
+    tgtDir_valid, uniqueOnsetLabels] ...
     = eventLockedAvg(y_r', t_r, onsetTimes, tgtDir, param.figTWin);
 
 tonsetRespAmp = characteriseResp(singleOnsetResp, ...
@@ -15,47 +15,9 @@
 
 %prefDir = 180/pi*circ_mean(tgtDir'/180*pi,tonsetRespAmp_p(:,psthIdx)); %NG
 %psthIdx = find(strcmp(psthNames, 'psth'));
-[fitPars, fitErr, r] = fitResponse(tgtDir, tonsetRespAmp, param.cardinalDir);
+[fitPars, fitErr, r] = fitResponse(tgtDir_valid, tonsetRespAmp, param.cardinalDir);
 prefDir = fitPars(1); %[deg]
 [~, prefDirIdx] = min(abs(circ_dist(pi/180*param.cardinalDir, pi/180*prefDir)));
 prefDir_quantized = param.cardinalDir(prefDirIdx);
-prefDirTrials =  find(tgtDir == prefDir_quantized);%trials with cell's preferred direction
+prefDirTrials =  find(tgtDir_valid == prefDir_quantized);%trials with cell's preferred direction
 
-
-
-% % created from dirTuinig_pupil.m
-% % UNDER CONSTRUCTION
-% 
-% %% direction tuning
-% dirs = unique(dd.targetloc);
-% mresp = zeros(length(dirs),3);
-% mparea = [];
-% for idir = 1:length(dirs)
-% 
-% 
-%     theseTr = intersect(find(dd.successTrials), find(dd.targetloc==dirs(idir)));
-% 
-%     resp_c = [];
-%     parea_c = [];
-%     for itr = 1:length(theseTr)
-%         thisTr = theseTr(itr);
-%         theseTimes  = intersect(find(eyeData(thisTr).t - dd.tOnset(thisTr) > 1e-3*respWin(1)), ...
-%             find(eyeData(thisTr).t - dd.tOnset(thisTr) < 1e-3*respWin(2)));
-%         resp_c(itr) = sum(psth_tr{thisTr}(theseTimes))/diff(respWin)/1e-3; %[spikes/s]
-% 
-%         theseTimes_pre  = intersect(find(eyeData(thisTr).t - dd.tOnset(thisTr) > 1e-3*preWin(1)), ...
-%             find(eyeData(thisTr).t - dd.tOnset(thisTr) < 1e-3*preWin(2)));
-% 
-%        parea_c(itr) = mean(eyeData_rmblk_tr(thisTr).parea(theseTimes_pre));
-%     end
-%     mresp(idir,1) = mean(resp_c);
-%     mresp(idir,2) = nanmean(resp_c(parea_c<parea_avg));
-%     mresp(idir,3) = nanmean(resp_c(parea_c>parea_avg));
-%     mparea(idir) = mean(parea_c);
-%     length(find(parea_c>parea_avg))
-% end
-% 
-% plot(dirs, mresp);
-% legend('all trials','parea<mean','parea>mean');
-% xlabel('tgt direction [deg]');
-% %polarplot(dirs/180*pi, mresp);

getExampleCells.m

@@ -0,0 +1,10 @@
+function [prevDir_quantized, prevDirTrials] = getPrevDir(tgtDir, param)
+%[prevDir_quantized, prevDirTrials] = getPrevDir(tgtDir, param)
+% return the most prevalent stimulus direction
+% Note this direction is irrelevant of neural response
+% see also: getPrefDir
+
+prevDir = mode(tgtDir);
+[~, prevDirIdx] = min(abs(circ_dist(pi/180*param.cardinalDir, pi/180*prevDir)));
+prevDir_quantized = param.cardinalDir(prevDirIdx);
+prevDirTrials =  find(tgtDir == prevDir_quantized);%trials with cell's preferred direction

getPrefDir.m

@@ -14,6 +14,7 @@
         %addpath(genpath('/home/daisuke/Documents/git'));
         rmpath(genpath('/home/daisuke/Documents/git/VariabilityEarlyVisualCortex/matlab/'))
         rmpath(genpath('/home/daisuke/Documents/git/dsbox/chronux_2_11/'));
+        rmpath(genpath('/home/daisuke/Documents/git/chronux_2_11'));
         %saveServer = '~/Documents/cuesaccade_data';
         %saveServer = '/mnt/syncitium/Daisuke/cuesaccade_data';
         saveServer = '/media/daisuke/cuesaccade_data';

getPrevDir.m

@@ -1,5 +1,5 @@
 function [trialEndTimes, successTimes, trialOutcome] = getRewardTimes(dd)
-%[rewardTimes, unishTimes, successTimes, trialOutcome] = getRewardTimes(dd)
+%[trialEndTimes, successTimes, trialOutcome] = getRewardTimes(dd)
 %returns times of reward delivery detected as item2delivered
 %
 %cf. concatenate_eye.m

getReady.m

@@ -1,24 +1,29 @@
-function [latency_neuro, thresh_neuro] = getSingleTrialLatency(singleResp_t, winSamps_t, tWin_t, threshOption, Thresh)
+function [latency_neuro, thresh_neuro] = getSingleTrialLatency(singleResp_t, winSamps_t, ...
+    tWin_t, threshParam)
 %[latency_neuro, thresh_neuro] = getSingleTrialLatency(singleResp_t, winSamps_t, tWin_t, threshOption)
 % from     latencyV1MT/secondary/findSingleTrialLatency.m
 % singleResp_t: trial x time
 
 baseline = 'tonset'; %used to have an option of fonset but turned out to be noisier
 
-if nargin < 5
-    Thresh = 4; %s.d.
-end
 if nargin < 4
-    threshOption = 'individual'; %india's choice
+    threshParam.dur = 0;
+    threshParam.thresh = 4; %s.d.
+    threshParam.option = 'individual'; %india's choice
 end
 
+assert(strcmp(threshParam.option, 'individual') + strcmp(threshParam.option, 'uniform') )
+
 latencyLimit = [0 tWin_t(2)];
 x = 1e-2*round(1e2*winSamps_t);
 sdf_t = singleResp_t; %trial x time
 
+dt = median(diff(x));
+tidx_dur = round(threshParam.dur/dt)-1;
+
 %% define threshold
 thresh_neuro = nan(size(sdf_t,1),1);
-if strcmp(threshOption, 'individual')
+if strcmp(threshParam.option, 'individual')
     for itr = 1:size(sdf_t,1)
 
         if strcmp(baseline, 'tonset')
@@ -27,17 +32,17 @@
 
         sd = std(sdf_pre(:));
         mu = mean(sdf_pre(:));
-        thresh = mu + Thresh*sd; %compute each trail or across all trials?
+        thresh = mu + threshParam.thresh*sd; %compute each trail or across all trials?
         thresh_neuro(itr) = thresh;
     end
-elseif strcmp(threshOption, 'uniform');
+elseif strcmp(threshParam.option, 'uniform');
     if strcmp(baseline, 'tonset')
         sdf_pre = sdf_t(:,x<0);
     end
 
     sd = std(sdf_pre(:));
     mu = mean(sdf_pre(:));
-    thresh = mu + Thresh*sd; %compute each trail or across all trials?
+    thresh = mu + threshParam.thresh*sd; %compute each trail or across all trials?
     thresh_neuro(:) = thresh;
 end
 
@@ -46,8 +51,23 @@
 latency_neuro = nan(size(sdf_t,1),1);
 for itr = 1:size(sdf_t,1)
     j = find(x >= 0,1);
-    if j ~= size(sdf_t,2)
-        latency = x(find(sdf_t(itr,j:end) >= thresh_neuro(itr),1) + j-1); %Calculates latency as first time point after stimulus where sdf is >= threshold
+    tidx_a = find(sdf_t(itr,j:end) >= thresh_neuro(itr)) + j-1;
+
+    if j ~= size(sdf_t,2) && ~isempty(tidx_a)
+        %Calculates latency as first time point after stimulus where sdf is >= threshold
+        
+        for xxx = 1:numel(tidx_a)
+            tidx = tidx_a(xxx);
+
+            %minimal duration surpassing the threshold
+            if isempty(find(tidx+tidx_dur<=numel(x), 1))
+                latency = NaN; continue;
+            elseif sum(sdf_t(itr, tidx:tidx+tidx_dur) >= thresh_neuro(itr)) == tidx_dur+1
+                latency = x(tidx); break;
+            else 
+                latency = NaN;
+            end
+        end
     else
         latency = NaN;
     end
@@ -60,6 +80,7 @@
     if isempty(latency)
         latency = NaN;
     end
+    
 
     latency_neuro(itr) = latency;
 end

getRewardTimes.m

@@ -1,16 +1,14 @@
 function [latency_neuro, latency_bhv, latcorr, f, fneuro] = getTgtLatencyCorr(PSTH_f, t_r, onsets_cat, ...
-    catEvTimes, tWin_t, Thresh, param, dd )
+    catEvTimes, tWin_t, param, dd, validEvents )
 % [latency_neuro, latency_bhv, r, f, fneuro] = getTgtLatencyCorr(PSTH_f, t_r, onsets_cat, ...
 % catEvTimes, tWin_t, Thresh, param, dd)
 % computes neural and behavioural latencies with or without cue
 
-%validEvents = intersect(find(~isnan(onset)), find(dd.cueOn==icue-1));
-%< this condition only includes all trials irrespective of the trial outcome
-validEvents = find(~isnan(catEvTimes.tOnset) .* ~isnan(catEvTimes.fOnset)); %with or without cue
+corrOption = 'Spearman'; %'Pearson'
 
 %% neural latency
 [latency_neuro, ~,~,fneuro] = getTgtNeuroLatency(PSTH_f,  t_r,  onsets_cat,  ...
-    catEvTimes, tWin_t, Thresh, param, dd, validEvents);
+    catEvTimes, tWin_t, param.threshParam, param, dd, validEvents);
 
 %% behavioural latency
 option = 1;
@@ -26,37 +24,80 @@
 fail = choiceOutcome(validEvents) >= 2;
 
 
-%% select trials with preferred direction (showTonsetResp)
 theseEvents = validEvents;%(success);%better to use all conditions??
 onsetTimes = catEvTimes.tOnset; 
 onsetTimes = onsetTimes(theseEvents);
 tgtDir = getTgtDir(dd.targetloc, param.cardinalDir);
 tgtDir = tgtDir(theseEvents);
+
+%% select trials with preferred direction (showTonsetResp)
 [prefDir, prefDirTrials] = getPrefDir(PSTH_f, t_r, onsetTimes, tgtDir, param);
-tgt_in = zeros(numel(validEvents),1);
-tgt_in(prefDirTrials) = 1;
+tgt_pref = zeros(numel(validEvents),1);
+tgt_pref(prefDirTrials) = 1;
 
+%% select trials with most frequent direction
+[prevDir, prevDirTrials] = getPrevDir(tgtDir, param);
+tgt_prev = zeros(numel(validEvents),1);
+tgt_prev(prevDirTrials) = 1;
 
 validLatency = ~isnan(latency_neuro);
 
 %correlation beteween behaviour and neuro using only successful trials
-latcorr.r = corr(latency_bhv(logical(validLatency.*~isnan(latency_bhv))), latency_neuro(logical(validLatency.*~isnan(latency_bhv))));
-latcorr.r_success = corr(latency_bhv(logical(success.*validLatency)), latency_neuro(logical(success.*validLatency)));
-latcorr.r_pref = corr(latency_bhv(logical(validLatency.*~isnan(latency_bhv).*tgt_in)), latency_neuro(logical(validLatency.*~isnan(latency_bhv).*tgt_in)));
-latcorr.r_success_pref = corr(latency_bhv(logical(success.*validLatency.*tgt_in)), latency_neuro(logical(success.*validLatency.*tgt_in)));
+try
+    [latcorr.r, latcorr.p] = corr(latency_bhv(logical(validLatency.*~isnan(latency_bhv))), ...
+        latency_neuro(logical(validLatency.*~isnan(latency_bhv))), 'type', corrOption);
+catch err
+    latcorr.r = nan;
+    latcorr.p = nan;
+end
+try
+    [latcorr.r_success, latcorr.p_success] = corr(latency_bhv(logical(success.*validLatency)), ...
+        latency_neuro(logical(success.*validLatency)), 'type', corrOption);
+catch err
+    latcorr.r_success = nan;
+    latcorr.p_success = nan;
+end
+try
+    [latcorr.r_pref, latcorr.p_pref] = corr(latency_bhv(logical(validLatency.*~isnan(latency_bhv).*tgt_pref)), ...
+        latency_neuro(logical(validLatency.*~isnan(latency_bhv).*tgt_pref)), 'type', corrOption);
+catch err
+    latcorr.r_pref = nan;
+    latcorr.p_pref = nan;
+end
+try
+    [latcorr.r_success_pref, latcorr.p_success_pref] = corr(latency_bhv(logical(success.*validLatency.*tgt_pref)), ...
+        latency_neuro(logical(success.*validLatency.*tgt_pref)), 'type', corrOption);
+catch err
+    latcorr.r_success_pref  = nan;
+    latcorr.p_success_pref = nan;
+end
+try
+    [latcorr.r_success_prev, latcorr.p_success_prev] = corr(latency_bhv(logical(success.*validLatency.*tgt_prev)), ...
+        latency_neuro(logical(success.*validLatency.*tgt_prev)), 'type', corrOption);
+catch err
+    latcorr.r_success_prev  = nan;
+    latcorr.p_success_prev = nan;
+end
 
 if nargout > 3
     f = figure;
     plot(latency_bhv(logical(wCue.*success)), latency_neuro(logical(wCue.*success)),'ro','DisplayName','success w cue'); hold on;
     plot(latency_bhv(logical(woCue.*success)), latency_neuro(logical(woCue.*success)),'bo','DisplayName','succes wo cue'); hold on;
-    plot(latency_bhv(logical(wCue.*success.*tgt_in)), latency_neuro(logical(wCue.*success.*tgt_in)),'o','MarkerFaceColor','r','DisplayName','success w cue prefDir'); hold on;
-    plot(latency_bhv(logical(woCue.*success.*tgt_in)), latency_neuro(logical(woCue.*success.*tgt_in)),'o','MarkerFaceColor','b','DisplayName','succes wo cue prefDir'); hold on;
+    plot(latency_bhv(logical(wCue.*success.*tgt_pref)), latency_neuro(logical(wCue.*success.*tgt_pref)),'o','MarkerFaceColor','r','DisplayName','success w cue prefDir'); hold on;
+    plot(latency_bhv(logical(woCue.*success.*tgt_pref)), latency_neuro(logical(woCue.*success.*tgt_pref)),'o','MarkerFaceColor','b','DisplayName','succes wo cue prefDir'); hold on;
     plot(latency_bhv(logical(wCue.*fail)), latency_neuro(logical(wCue.*fail)),'rx','DisplayName','failure w cue'); hold on;
     plot(latency_bhv(logical(woCue.*fail)), latency_neuro(logical(woCue.*fail)),'bx','DisplayName','failure wo cue'); hold on;
     xlabel('behavioural latency');
     ylabel('neural latency');
     legend('location','eastoutside');
-    tname = sprintf('all dir: %.2f, success all dir: %.2f \n success pref dir: %.2f', latcorr.r, latcorr.r_success,  latcorr.r_success_pref);
+    tname = sprintf(['all dir: %.2f (p=%.2f), \n' ...
+        'success all dir: %.2f (p=%.2f) \n' ...
+        'success pref dir: %.2f (p=%.2f)\n' ...
+         'success prev dir: %.2f (p=%.2f)'], ...
+        latcorr.r, latcorr.p, ...
+        latcorr.r_success, latcorr.p_success, ...
+        latcorr.r_success_pref, latcorr.p_success_pref,...
+        latcorr.r_success_prev, latcorr.p_success_prev);
     title(tname); %title(['r=' num2str(r) ', r(success)=' num2str(r_success)]);
     squareplot;
     axis padded; set(gca,'TickDir','out');