Skip to content

Commit

Permalink
updating calc su mod
Browse files Browse the repository at this point in the history
  • Loading branch information
@droumis
droumis committed Apr 23, 2020
1 parent 487944a commit bd58135
Show file tree
Hide file tree
Showing 4 changed files with 161 additions and 155 deletions.
125 changes: 63 additions & 62 deletions Functions/calcPhaseMod.m
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

function out = calcPhaseMod(F, dmat, varagin)
% [out] = calcSUmod(F, varargin)
% TODO: rename func to calcPhaseModSU
% calculate event based phase modulation of spiking
%
% args:
Expand All @@ -10,11 +11,6 @@
% - dayeps: [day ep;...] per ev
% - dm: binary. (event x set)
% - expvars: dmat labels as cell array of strings
%{
Saw
- barn:rat:beer:saw
@DKR
%}

pconf = paramconfig;
comuteShuf = 1;
Expand All @@ -23,20 +19,19 @@
minILIthresh = .06; % seconds
maxILIthresh = .250; % seconds
numPhaseBins = 36;
try
assign(varagin{:})
catch
end
minILIspikes = 100; % TODO: implement min ILI spikes cumulative per unit

try assign(varagin{:}); catch; end

for a = 1:length(F)
animal = F(a).animal{3};
out(a).animal = F(a).animal;
fprintf('=========== %s ===========\n', animal);
out(a).output = {};
out(a).outputIdx = dmat(a).expvars;
for iv = 1:size(dmat(a).expvars,2)
for iv = 1:size(dmat(a).expvars,2) % for each experimental condition
out(a).output{iv} = [];
for c = 1:length(F(a).output{1})
for c = 1:length(F(a).output{1}) % could parfor this for each unit
cF = init_out();
cF.animal = animal;
% collect event design mat for this cluster (per day)
Expand All @@ -58,25 +53,25 @@
end
end
cF.index = idx;

% for each condition, get events' psth,
% for each condition, get events' psth
eT = F(a).output{1}(c).eventTimes(dayDM(:,iv));
time = F(a).output{1}(c).time;
psth = F(a).output{1}(c).psth(dayDM(:,iv),:);
% then collect spikes into inter event bins and compute their pct of interval
ILI = diff(eT);
cIdx = knnsearch(time', 0);
ILIidx = knnsearch(time', ILI);
% collect and transform spikes into inter event bins and
% compute spike bin pct of interval
ILI = diff(eT); % diff gives duration from each lick to the next
cIdx = knnsearch(time', 0); % find center (event time index)
ILIidx = knnsearch(time', ILI); % get index into time for each ILI duration
pSinceLick = [];
spkOffsetAll = [];
for e = 1:length(ILI) % per event interval
if ILI(e) < maxILIthresh && ILI(e) > minILIthresh
% get spikes between eventA and eventB
if ILI(e) < maxILIthresh && ILI(e) > minILIthresh % if valid ILI
% get spikes within current ILI
spIliIdx = find(psth(e,cIdx:ILIidx(e)));
if ~isempty(spIliIdx)
% get event interval percent of spikes
spkOffset = spIliIdx*bin; % idx distance from center (event), scaled to time
pSinceLick{e,1} = [spkOffset ./ ILI(e)]';
pSinceLick{e,1} = [spkOffset ./ ILI(e)]'; % spiketime pct of ILI
spkOffsetAll{e,1} = spkOffset;
end
else
Expand All @@ -85,50 +80,57 @@
end
spikePctSinceLick = cell2mat(pSinceLick);
if ~isempty(spikePctSinceLick)
% convert percent into phase
% convert percent into phase. since pct is 0-1, just 2*pi*pct
spikeIEIphase = 2*pi*spikePctSinceLick;
% get mean resultant vec mag and angle
% get mean resultant vec mag and angle *****
meanvec = mean(exp(1i*spikeIEIphase));
meanMRVmag = abs(meanvec);
vecang = angle(meanvec);
meanMRVmag = abs(meanvec); % magnitude
vecang = angle(meanvec); % angle of mean vec
% get log normalized Raleigh Z as 'phasemod'
[pval, z] = circ_rtest(spikeIEIphase);
phasemod = log(z);

cF.pval = pval;
cF.ILI = ILI;
cF.spkOffsetAll = spkOffsetAll;
cF.spikePctSinceLick = spikePctSinceLick;
cF.spikeIEIPhase = spikeIEIphase;
cF.meanMRVmag = meanMRVmag;
cF.vecang = vecang;
cF.phasemod = phasemod;
cF.area = idata.area;
cF.subarea = idata.subarea;

% get phaseHist
phasemod = log(z); % log normalize
% make phase histogram edges
binphaseEdges = linspace(0, 2*pi, numPhaseBins);
cf.spikeIEIPhaseHistProb = histcounts(spikeIEIphase, binphaseEdges, ...
% compute phase histogram *****
spikeIEIPhaseHistProb = histcounts(spikeIEIphase, binphaseEdges, ...
'Normalization', 'probability');
%
%% Shuffle
if comuteShuf
% shuffle percent of interval for all spikes
% rand shuffle percent of interval for all spikes
% TODO: shuffle spikeIEIphase instead of spikePctSinceLick
r = rand(length(spikePctSinceLick), numShuf);
cF.phasemodSh = [];
cF.vecangSh = [];
cF.meanMRVmagSh = [];
meanMRVmagSh = nan(numShuf);
vecangSh = nan(numShuf);
phasemodSh = nan(numShuf);
for s = 1:numShuf
spikeIEIphase = 2*pi*r(:,s);
meanvec = mean(exp(1i*spikeIEIphase));
cF.meanMRVmagSh = [cF.meanMRVmagSh; abs(meanvec)];
cF.vecangSh = [cF.vecangSh; angle(meanvec)];
meanMRVmagSh(s) = abs(meanvec);
vecangSh(s) = angle(meanvec);
[~, z] = circ_rtest(spikeIEIphase);
cF.phasemodSh = [cF.phasemodSh; log(z)];
phasemodSh(s) = log(z);
end
% real-mod shuff-pct
cF.modPctRank = 100*(1-(sum(cF.phasemodSh > cF.phasemod)./numShuf));
fprintf('iv%d mod > %.02f pctShufs.\n', iv, cF.modPctRank)
% real-phasemod vs shuff-phasemod distribution
modPctRank = 100*(1-(sum(phasemodSh > phasemod)./numShuf));
fprintf('iv%d mod > %.02f pctShufs.\n', iv, modPctRank)
% Shuff output
cF.meanMRVmagSh = meanMRVmagSh;
cF.vecangSh = vecangSh;
cF.phasemodSh = phasemodSh;
cF.modPctRank = modPctRank;
end
%% unit output
cF.area = idata.area;
cF.subarea = idata.subarea;
cF.ILI = ILI;
cF.spkOffsetAll = spkOffsetAll;
cF.spikePctSinceLick = spikePctSinceLick;
cF.spikeIEIPhase = spikeIEIphase;
cF.meanMRVmag = meanMRVmag;
cF.vecang = vecang;
cF.pval = pval; % circ_r
cF.phasemod = phasemod;
cF.spikeIEIPhaseHistProb = spikeIEIPhaseHistProb;
end
out(a).output{iv} = [out(a).output{iv}; cF];
end
Expand All @@ -137,21 +139,20 @@
end

function cf = init_out()
cf.pval = nan;
cf.ILI = nan;
cf.spkOffsetAll = nan;
cf.spikePctSinceLick = nan;
cf.spikeIEIPhase = nan;
% are nans the right init for all these? timemod uses mostly []
cf.area = '';
cf.subarea = '';
cf.ILI = [];
cf.spkOffsetAll = [];
cf.spikePctSinceLick = [];
cf.spikeIEIPhase = [];
cf.meanMRVmag = nan;
cf.vecang = nan;
cf.phasemod = nan;
cf.modPctRank = nan;

cf.area = '';
cf.subarea = '';

cf.spikeIEIPhaseHistProb = [];
cf.meanMRVmagSh = nan;
cf.vecangSh = nan;
cf.phasemodSh = nan;
cF.phasemodSh = [];
cF.vecangSh = [];
cF.meanMRVmagSh = [];
cf.pval = nan;
end
99 changes: 43 additions & 56 deletions Functions/calcSUmod.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

function out = calcSUmod(F, dmat, varargin)
% [out] = calcSUmod(F, varargin)
% TODO: rename func to calcTimeModSU
% calculate event-triggered modulation of SU spiking
%
% args:
Expand All @@ -12,39 +13,28 @@
% - dm: binary. (event x set)
% - expvars: dmat labels as cell array of strings
% varargs:
%
%
%{
Wheelbarrow
- barn:rat:beer:wheelbarrow
Notes;
need to add this func to dfa_eventTrigSpiking..
or add that to this and make this a dfa..
FFPhy
@DKR
%}

pconf = paramconfig;
% filtF filter to select animal, epochs, ntrode, clusers
rwin = [0 .2]; % response period in s rel to swr on
bwin = [-.7 -.2]; % baseline period in s rel to swr on
% minNumEvents = 10;
% minSpikesResp = 10;
% minSpikesBase = 10;
% minNumSwrSpikes = 10;
minSpikes = 50;
% runShuffle = 1;
minSpikes = 50; % TODO.. need to implement
computeShuf = 0;
nshuffs = 1000;
shuffbyMax = .7; % s to rand shuff by
saveResult = 1;
filetail = '';
% dmat = [];
% dmatIdx = {'all'};
use_psfr = 1; % use firing rate. otherwise uses instantaneous firing rate (see rat)

if ~isempty(varargin)
assign(varargin{:})
end
try assign(varargin{:}); catch; end

out = struct;
for a = 1:length(F)
Expand All @@ -55,7 +45,7 @@
out(a).dmatIdx = dmat(a).expvars;
for iv = 1:size(dmat(a).expvars,2)
out(a).output{iv} = [];
for c = 1:length(F(a).output{1}) % for each cluster
for c = 1:length(F(a).output{1}) % could parfor this for each cluster
cF = init_out();
cF.animal = animal;
% collect event design mat for this cluster (per day)
Expand All @@ -77,7 +67,6 @@
end
end
cF.index = idx;
use_psfr = 1;
if use_psfr
time = idata.wtime;
psfr = idata.psfr;
Expand All @@ -87,9 +76,8 @@
end
cF.time = time;
%% meanmod score for this cluster, per condition

ivIdx = find(dayDM(:,iv));
evIFR = psfr(ivIdx,:);
evIFR = psfr(ivIdx,:);
cF.evMean = nanmean(evIFR,1); % full mean per DM.
cF.evMeanZ = zscore(cF.evMean);

Expand All @@ -114,10 +102,10 @@
numBSpikes = sum(sum(idata.psth(dayDM(:,iv), psthbIdx(1):psthbIdx(2))));

fprintf('spikes in baseline period: %d \n', numBSpikes);
% if (numBSpikes + numRSpikes) < minSpikes
% fprintf('skipping\n')
% continue
% end
if (numBSpikes + numRSpikes) < minSpikes
fprintf('skipping\n')
continue
end

% take the diff of base and resp period for each event
% mean pct change from baseline
Expand All @@ -128,37 +116,37 @@

evIFRz = zscore(evIFR, [], 2);
mZResp = nanmean(nanmean(evIFRz(:,rIdx(1):rIdx(2)),2));
% OP.numRSpikes{iv} = numRSpikes;
% OP.numBSpikes{iv} = numBSpikes;

%% Shuffle
numEv = size(evIFR,1);
% binsize = diff(time(1:2));
shuffbyMaxBins = knnsearch(time', shuffbyMax)-knnsearch(time', 0);
r = randi([-shuffbyMaxBins shuffbyMaxBins], numEv, nshuffs);
mPctChangeSh = nan(nshuffs,1);
mZChangeSh = nan(nshuffs,1);
parfor i = 1:nshuffs % can use parfor
% for each event, select a shuf rand period as r and b
% time-mean fr in response period
evRespMSh = nanmean(cell2mat(arrayfun(@(x,y) evIFR(x,rIdx(1)+y:rIdx(2)+y), ...
1:size(r,1), r(:,i)', 'uni',0)'),2);
evBaseMSh = nanmean(cell2mat(arrayfun(@(x,y) ...
evIFR(x,bIdx(1)+y:bIdx(2)+y), 1:size(r,1), r(:,i)', 'uni',0)'),2);
evBaseSTDSh = nanstd(cell2mat(arrayfun(@(x,y) ...
evIFR(x,bIdx(1)+y:bIdx(2)+y), 1:size(r,1), r(:,i)', 'uni',0)'),[],2);
% mean pct change from baseline
useB = ~(evBaseMSh == 0);
mPctChangeSh(i) = nanmean((evRespMSh(useB)-evBaseMSh(useB))./evBaseMSh(useB))*100;
useBs = ~(evBaseSTDSh == 0);
mZChangeSh(i) = nanmean((evRespMSh(useBs)-evBaseMSh(useBs))./evBaseSTDSh(useBs));
if computeShuf
numEv = size(evIFR,1);
shuffbyMaxBins = knnsearch(time', shuffbyMax)-knnsearch(time', 0);
r = randi([-shuffbyMaxBins shuffbyMaxBins], numEv, nshuffs);
mPctChangeSh = nan(nshuffs,1);
mZChangeSh = nan(nshuffs,1);
% use parfor at unit level above to reduce worker data copies
for i = 1:nshuffs % can use parfor.
% for each event, select a shuf rand period as r and b
% time-mean fr in response period
evRespMSh = nanmean(cell2mat(arrayfun(@(x,y) evIFR(x,rIdx(1)+y:rIdx(2)+y), ...
1:size(r,1), r(:,i)', 'uni',0)'),2);
evBaseMSh = nanmean(cell2mat(arrayfun(@(x,y) ...
evIFR(x,bIdx(1)+y:bIdx(2)+y), 1:size(r,1), r(:,i)', 'uni',0)'),2);
evBaseSTDSh = nanstd(cell2mat(arrayfun(@(x,y) ...
evIFR(x,bIdx(1)+y:bIdx(2)+y), 1:size(r,1), r(:,i)', 'uni',0)'),[],2);
% mean pct change from baseline
useB = ~(evBaseMSh == 0);
mPctChangeSh(i) = nanmean((evRespMSh(useB)-evBaseMSh(useB))./evBaseMSh(useB))*100;
useBs = ~(evBaseSTDSh == 0);
mZChangeSh(i) = nanmean((evRespMSh(useBs)-evBaseMSh(useBs))./evBaseSTDSh(useBs));
end
% real-mod shuff-pct
modPctRank = 100*(1-(sum(mPctChangeSh > mPctChange)./nshuffs));
modZRank = 100*(1-(sum(mZChangeSh > mZChange)./nshuffs));
fprintf('iv%d Pctmod > %.02f pctShuffs.\n', iv, modPctRank)
fprintf('iv%d Zmod > %.02f pctShuffs.\n', iv, modZRank)
end
%%
% real-mod shuff-pct
modPctRank = 100*(1-(sum(mPctChangeSh > mPctChange)./nshuffs));
modZRank = 100*(1-(sum(mZChangeSh > mZChange)./nshuffs));
fprintf('iv%d Pctmod > %.02f pctShuffs.\n', iv, modPctRank)
fprintf('iv%d Zmod > %.02f pctShuffs.\n', iv, modZRank)

%% unit output
cF.mPctChange = mPctChange;
cF.mPctChangeSh = mPctChangeSh;
cF.modPctRank = modPctRank;
Expand All @@ -169,15 +157,14 @@
cF.area = idata.area;
cF.subarea = idata.subarea;
cF.cellInfo = idata.cellInfo;

out(a).output{iv} = [out(a).output{iv}; cF];
end

out(a).output{iv} = [out(a).output{iv}; cF];
end
end
end

function cF = init_out()
% are empty lists the right init for all these? phasmod uses mostly nan
cF.mZChange = [];
cF.animal = '';
cF.dmat = [];
Expand Down
Loading

0 comments on commit bd58135

Please sign in to comment.