Skip to content

Commit 55d85ae

Browse files
dshimaokadshimaoka
committed
population analysis w/wo success trials for fitting
only 150 units in 2021 March was analyzed
1 parent ee3784e commit 55d85ae

11 files changed

+373
-189
lines changed

‎fitPSTH_pop.m ‎fitPSTH_pop_limpredictor.m

+103-55
Original file line numberDiff line numberDiff line change
@@ -3,16 +3,16 @@
33
%align number of rows of population data
44

55
%% get ready
6+
limSuffix = '';%'_woSuccess';%_woSuccess';
7+
68
[saveServer, rootFolder] = getReady();
7-
saveSuffix_p = 'fitPSTH_pop20231102';
9+
saveSuffix_p = ['fitPSTH_pop' limSuffix];
810

911
%% recorded data
1012
animal = 'hugo';
1113
dataType = 0;%0: each channel, 1: all channels per day
1214
tWin = [0 0.5];%[s]
1315

14-
load(fullfile(saveServer,'param20231102.mat'),'param');
15-
1616
%alldata = [];
1717
mFiringRate_pop = [];
1818
kernel_pop = [];
@@ -30,8 +30,10 @@
3030
id_pop = [];
3131
Rsqadj_pop = [];
3232
gainInfo_pop = [];
33+
expval_trig_pop = [];
34+
corr_trig_pop = [];
3335
errorIDs= cell(1);
34-
for yy = 1:3
36+
for yy = 1%:3
3537
switch yy
3638
case 1
3739
year = '2021';
@@ -54,12 +56,12 @@
5456
%dataByYear = [];%struct(length(channels),1);
5557
%datech_pop = nan(length(channels),1);
5658
%corrcoef_pred_spk_pop = nan(length(channels),1);
57-
for idata = 1:length(channels)
59+
for idata = 1:150%length(channels)
5860
datech = [months{idata} '/' dates{idata} '/' num2str(channels{idata})];
5961
thisid = [animal '/' year '/' datech];
6062
disp(thisid);
6163

62-
saveSuffix = [animal replace(datech,'/','_') '_linear_rReg'];
64+
saveSuffix = [animal replace(datech,'/','_') '_linear_rReg' limSuffix];
6365

6466
thisDate = [months{idata} '_' dates{idata}];
6567

@@ -73,8 +75,7 @@
7375
S = load(saveName, 'PSTH_f','predicted_all', 'predicted', ...
7476
'kernelInfo','t_r','cellclassInfo','param','mFiringRate','t_cat',...
7577
'dds');
76-
77-
78+
7879
if isfield(S,'kernelInfo')
7980
mFiringRate_pop = cat(2,mFiringRate_pop,S.mFiringRate);
8081
kernel_pop = cat(3,kernel_pop,S.kernelInfo.kernel);
@@ -88,39 +89,39 @@
8889

8990
eyeName = fullfile(saveFolder,[
9091
'eyeCat_' animal thisDate '.mat']);
91-
load(fullfile(saveFolder,['predictorInfo_' animal thisDate '.mat']), ...
92-
'predictorInfo');
93-
load(eyeName,'eyeData_rmotl_cat','catEvTimes');
94-
95-
if ~isfield(S,'dds')
96-
load(loadNames{idata},'dd'); %slow to read
97-
dds = getCueSaccadeSmall(dd);
98-
save(saveName, 'dds','-append');
99-
S.dds = dd;
100-
end
101-
dd = S.dds;
92+
% load(fullfile(saveFolder,['predictorInfo_' animal thisDate '.mat']), ...
93+
% 'predictorInfo');
94+
load(eyeName,'eyeData_rmotl_cat','catEvTimes','startSaccNoTask');
10295

96+
% if ~isfield(S,'dds')
97+
load(loadNames{idata},'dd'); %slow to read
98+
% dds = getCueSaccadeSmall(dd);
99+
% save(saveName, 'dds','-append');
100+
% S.dds = dd;
101+
% end
102+
% dd = S.dds;
103103

104-
%% Rsq adj of subjset of variables
105-
nsub=4;
106-
Rsqadj = zeros(nsub,1);
107-
for jj = 1:nsub
108-
switch jj
109-
case 1 %full model
110-
tgtGroups = 1:5;
111-
case 2 %omit vision %added 26/10/2023
112-
tgtGroups = setxor(1:5, 1);
113-
case 3 %omit eye speed
114-
tgtGroups = setxor(1:5, 2);
115-
case 4 %omit eye position
116-
tgtGroups = setxor(1:5, 3);
117-
end
118-
[Rsqadjusted,rr,r0] = fitSubset(S.PSTH_f, predictorInfo, ...
119-
tgtGroups, param);%, idxTgtOnsets);
120104

121-
Rsqadj(jj) = Rsqadjusted;
122-
end
123-
Rsqadj_pop = [Rsqadj_pop Rsqadj];
105+
%% Rsq adj of subjset of variables
106+
% nsub=4;
107+
% Rsqadj = zeros(nsub,1);
108+
% for jj = 1:nsub
109+
% switch jj
110+
% case 1 %full model
111+
% tgtGroups = 1:5;
112+
% case 2 %omit vision %added 26/10/2023
113+
% tgtGroups = setxor(1:5, 1);
114+
% case 3 %omit eye speed
115+
% tgtGroups = setxor(1:5, 2);
116+
% case 4 %omit eye position
117+
% tgtGroups = setxor(1:5, 3);
118+
% end
119+
% [Rsqadjusted,rr,r0] = fitSubset(S.PSTH_f, predictorInfo, ...
120+
% tgtGroups, param);%, idxTgtOnsets);
121+
%
122+
% Rsqadj(jj) = Rsqadjusted;
123+
% end
124+
% Rsqadj_pop = [Rsqadj_pop Rsqadj];
124125

125126
% %Rsq_adjusted computed from pre-target
126127
% periods - NG. can be < 0
@@ -135,10 +136,10 @@
135136
% Rsqadj_tgt_pop = [Rsqadj_tgt_pop Rsqadjusted_tgt];
136137

137138
%% response to target
138-
PtonsetResp_pop = [PtonsetResp_pop S.cellclassInfo.PtonsetResp];
139-
140-
%% response to saccade
141-
PsaccResp_pop = [PsaccResp_pop S.cellclassInfo.PsaccResp];
139+
%PtonsetResp_pop = [PtonsetResp_pop S.cellclassInfo.PtonsetResp];
140+
141+
% %% response to saccade
142+
% PsaccResp_pop = [PsaccResp_pop S.cellclassInfo.PsaccResp];
142143

143144
%% explained variance per kernel
144145
expval = zeros(size(S.predicted,2)+1,1);
@@ -152,7 +153,7 @@
152153

153154
%% explained variance for target response
154155
expval_tgt(1,1) = getExpVal_tgt(S.PSTH_f, S.predicted_all, catEvTimes, S.t_r, tWin);
155-
expval_tgt(2:6,1) = getExpVal_tgt(S.PSTH_f, S.predicted, catEvTimes, S.t_r, tWin);
156+
expval_tgt(2:7,1) = getExpVal_tgt(S.PSTH_f, S.predicted, catEvTimes, S.t_r, tWin);
156157

157158
expval_tgt_pop = [expval_tgt_pop expval_tgt];
158159

@@ -164,6 +165,31 @@
164165
% expval_avgtgt_pop = [expval_avgtgt_pop expval_avgtgt];
165166
% corr_avgtgt_pop = [corr_avgtgt_pop corr_avgtgt];
166167

168+
%% compute time-resolved explained variance
169+
[choiceOutcome] = getChoiceOutcome(dd);
170+
171+
expval_trig = []; corr_trig = [];
172+
for ievtype = 1:4
173+
% 1: success trial
174+
% 2: failed quiescent trial
175+
% 3: failed wrong saccade direction
176+
% 4: saccade outside trials
177+
178+
if ievtype <= 3
179+
theseEvents = find(choiceOutcome==ievtype);
180+
eventTimes = catEvTimes.tOnset(theseEvents);
181+
elseif ievtype == 4
182+
eventTimes = startSaccNoTask;
183+
end
184+
[expval_trig(:,:,ievtype), corr_trig(:,:,ievtype), winSamps] = ...
185+
getExpVal_trig(S.PSTH_f, [S.predicted_all S.predicted], S.t_r, eventTimes, [-0.5 0.5]);
186+
187+
%ax_expval_trig(ievtype) = subplot(4,1,ievtype);
188+
%plot(winSamps, squeeze(expval_trig(:,:,ievtype))');
189+
end
190+
expval_trig_pop = cat(4, expval_trig_pop, expval_trig);
191+
corr_trig_pop = cat(4, corr_trig_pop, expval_trig);
192+
167193
%% number of target trials
168194
ntargetTrials_pop = [ntargetTrials_pop sum(~isnan(catEvTimes.tOnset))];
169195

@@ -175,17 +201,20 @@
175201
onlySuccess = 0;
176202
respWin = [0.05 0.35]; %[s]
177203
y_r = cat(2,S.PSTH_f,S.predicted_all);
178-
gainInfo = getGainInfo(S.t_r, y_r, param.cardinalDir, catEvTimes, ...
204+
gainInfo = getGainInfo(S.t_r, y_r, S.param.cardinalDir, catEvTimes, ...
179205
dd, figTWin, onlySuccess, respWin);
180206
gainInfo_pop = [gainInfo_pop gainInfo];
181207

182208
%retrieve just once
183209
tlags = S.kernelInfo.tlags;
210+
param = S.param;
211+
184212
clear S
185213
end
186214
catch err
187215
errorIDs{numel(errorIDs)+1} = thisid;
188-
continue;
216+
disp('ERROR');
217+
%continue;
189218
end
190219
end
191220
end
@@ -195,24 +224,31 @@
195224

196225
% save('fitPSTH_pop20220202','avgPupilResp_pop', '-append');
197226
kernel_pop = squeeze(kernel_pop);
198-
save(fullfile(saveServer,saveSuffix_p,[saveSuffix_p animal]),'mFiringRate_pop','kernel_pop','expval_pop','corrcoef_pop',...
227+
save(fullfile(saveServer,[saveSuffix_p animal '.mat']),'mFiringRate_pop','kernel_pop','expval_pop','corrcoef_pop',...
199228
'corrcoef_pred_spk_pop','id_pop','ntotTrials_pop','ntargetTrials_pop','param',...
200229
'PsaccResp_pop','PtonsetResp_pop','expval_ind_pop','expval_tgt_pop','tlags',...
201-
'corr_avgtgt_pop','expval_avgtgt_pop','Rsqadj_pop','gainInfo_pop');
230+
'corr_avgtgt_pop','expval_avgtgt_pop','gainInfo_pop','expval_trig_pop','corr_trig_pop','winSamps'); %'Rsqadj_pop',
231+
232+
% from syncitium/Daisuke/cuesaccade_data OBS/param20231102.mat
233+
param.mfiringRateTh = 5;
234+
param.expvalTh = 3;
235+
param.ntargetTrTh = 200;
236+
param.ptonsetRespTh = 0.05;
202237

203238
%% apply inclusion critetia
204239
% [okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
205240
% = inclusionCriteria(mFiringRate_pop, expval_pop, ntargetTrials_pop, PtonsetResp_pop, param);
206241
% [okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
207242
% = inclusionCriteria(mFiringRate_pop, expval_tgt_pop(1,:), ntargetTrials_pop, PtonsetResp_pop, param);
208-
[okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
209-
= inclusionCriteria(mFiringRate_pop, expval_ind_pop(1,:), ntargetTrials_pop, PtonsetResp_pop, param);
243+
% [okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
244+
% = inclusionCriteria(mFiringRate_pop, expval_ind_pop(1,:), ntargetTrials_pop, PtonsetResp_pop, param);
245+
246+
%TEMP as PtonsetResp_pop is unavailable:
247+
okunits = find((mFiringRate_pop>param.mfiringRateTh) ...
248+
+ (expval_ind_pop(1,:)>param.expvalTh) ...
249+
+ (ntargetTrials_pop> param.ntargetTrTh) == 3);
250+
save(fullfile(saveServer,[saveSuffix_p animal '.mat']), 'okunits','-append');
210251

211-
%% hack exclude redundant data
212-
[~, mfiringRate_u] = unique(mFiringRate_pop);
213-
[~, expval_u] = unique(expval_ind_pop(1,:));
214-
okunits_u = intersect(mfiringRate_u, expval_u);
215-
okunits = intersect(okunits, okunits_u);
216252

217253
kernel_pop = kernel_pop(:,okunits);
218254
expval_ind_pop = expval_ind_pop(:,okunits);
@@ -221,7 +257,19 @@
221257
%corr_avgtgt_pop = corr_avgtgt_pop(:,okunits);
222258
id_pop = id_pop(okunits);
223259
mFiringRate_pop = mFiringRate_pop(okunits);
224-
Rsqadj_pop = Rsqadj_pop(:,okunits);
260+
%Rsqadj_pop = Rsqadj_pop(:,okunits);
261+
expval_trig_pop = expval_trig_pop(:,:,:,okunits);
262+
263+
mexpval_trig_pop = squeeze(mean(expval_trig_pop,4));
264+
265+
for ievtype = 1:4
266+
ax_expval_trig(ievtype) = subplot(4,1,ievtype);
267+
plot(winSamps, squeeze(mexpval_trig_pop(:,:,ievtype))');
268+
end
269+
linkaxes(ax_expval_trig);
270+
ylim([-10 40])
271+
legend(['all',param.predictorNames],'location','west');
272+
screen2png(fullfile(saveServer,[saveSuffix_p animal '_expval_trig' limSuffix '.png']));
225273

226274

227275
%% selected units

‎getExpVal.m

+2
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
function [expval, mse, corr] = getExpVal(observed, predicted)
22
%[expval, mse, corr] = getExpVal(observed, predicted)
33

4+
predicted = predicted - mean(predicted) + mean(observed); %13/6/24
5+
46
mse = mean((observed - predicted).^2);
57
expval = 100*(1 - mse / mean((observed - mean(observed)).^2));
68
R = corrcoef(observed, predicted);

‎getExpVal_trig.m

+21
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
function [expval_trig, corr_trig, winSamps] = getExpVal_trig(PSTH_f, predicted, t_r, eventTimes, tWin)
2+
%[expval_trig, corr_trig] = getExpVal_trig(PSTH_f, predicted, t_r, eventTimes, tWin)
3+
% returns explained variance of triggered events, each moment in time from the event
4+
%
5+
% output:
6+
% expval_trig: [kernel type x time]
7+
% corr_trig: [kernel type x time]
8+
% 13/6/24 created from getExpVal_tgt
9+
10+
%% event-triggered traces
11+
[~, ~, resp] = eventLockedAvg(PSTH_f', t_r, eventTimes, [], tWin);
12+
[~, winSamps, resp_predicted] = eventLockedAvg(predicted', t_r, eventTimes, [], tWin);
13+
14+
%% explained variance for each moment in time
15+
expval_trig = zeros(size(predicted,2),numel(winSamps));
16+
corr_trig = zeros(size(predicted,2),numel(winSamps));
17+
for ivar = 1:size(predicted,2)
18+
for it = 1:numel(winSamps)
19+
[expval_trig(ivar, it), ~, corr_trig(ivar, it)] = getExpVal(resp(:,1,it), resp_predicted(:,ivar,it));
20+
end
21+
end

‎inclusionCriteria.m

+9
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,14 @@
11
function [okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
22
= inclusionCriteria(mFiringRate_pop, expval_pop, ntargetTrials_pop, PtonsetResp_pop, param)
3+
%[okunits, mfiringRateOK, expvalOK, ntargetTrOK, ptonsetRespOK] ...
4+
% = inclusionCriteria(mFiringRate_pop, expval_pop, ntargetTrials_pop, PtonsetResp_pop, param)
5+
%
6+
% param requires following fields:
7+
% mfiringRateTh
8+
% expvalTh
9+
% ntargetTrTh
10+
% ptonsetRespTh
11+
312
mfiringRateOK = (mFiringRate_pop > param.mfiringRateTh);
413
expvalOK = (expval_pop>param.expvalTh);
514
ntargetTrOK = (ntargetTrials_pop>param.ntargetTrTh );

‎limitPredictor.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
eventTimes = [];
1111

1212
for itr = 1:dd.numTrials
13-
%only register trials with reward
13+
%only register trials with reward from fixation behaviour till reward or punishment
1414
fOnset = onsets_cat.fOnset(itr);
1515
cOnset = onsets_cat.cOnset(itr);
1616
if isnan(cOnset) || isnan(fOnset) || trialOutcome(itr) == -1

0 commit comments

Comments
 (0)