From fe582fa3d3956cc306cd803e989c1127c824ec09 Mon Sep 17 00:00:00 2001 From: droumis Date: Tue, 13 Oct 2020 21:46:51 +0200 Subject: [PATCH] dfa for XPprepostSWR XPtrigAvgRip pctILB rewTrigSWRXP ripPos template --- DFAnalysis/dfa_XPprepostSWR.m | 127 ++++++++++++++++++++++++++++++++++ DFAnalysis/dfa_XPtrigAvgRip.m | 93 +++++++++++++++++++++++++ DFAnalysis/dfa_pctILB.m | 97 ++++++++++++++++++++++++++ DFAnalysis/dfa_rewTrigSWRXP.m | 101 +++++++++++++++++++++++++++ DFAnalysis/dfa_ripPos.m | 107 ++++++++++++++++++++++++++++ DFAnalysis/dfa_template.m | 33 +++++++++ 6 files changed, 558 insertions(+) create mode 100644 DFAnalysis/dfa_XPprepostSWR.m create mode 100644 DFAnalysis/dfa_XPtrigAvgRip.m create mode 100644 DFAnalysis/dfa_pctILB.m create mode 100644 DFAnalysis/dfa_rewTrigSWRXP.m create mode 100644 DFAnalysis/dfa_ripPos.m create mode 100644 DFAnalysis/dfa_template.m diff --git a/DFAnalysis/dfa_XPprepostSWR.m b/DFAnalysis/dfa_XPprepostSWR.m new file mode 100644 index 00000000..950ce2d0 --- /dev/null +++ b/DFAnalysis/dfa_XPprepostSWR.m @@ -0,0 +1,127 @@ + + +function out = dfa_XPprepostSWR(idx, timeFilter, varargin) +% get XP offset, pre and post SWR +% use with singledayanal +% DR_2020-09-30 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'lick'}; +check_required(reqData, varargin) + +run_shuffle = 1; +eventType = 'ca1rippleskons'; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); % init output +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); +%% Get SWRs in timeFilter + +evid = find(contains(varargin(1:2:end), eventType)); +o = [1:2:length(varargin)]+1; +swr = varargin{o(evid)}; +swrTimes = []; +maxthresh = []; +swrEnd = []; +for e = 1:length(eps) + ep = eps(e); + try + swrTimes = [swrTimes; swr{day}{ep}{1}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}{1}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}{1}.endtime]; + catch + try + swrTimes = [swrTimes; swr{day}{ep}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}.endtime]; + catch + fprintf('no swr events detected for day%d ep%d\n', day, ep) + return + end + end +end +% apply timefilter to swrs +incSWRs = ~isIncluded(swrTimes(:,1), timeFilter); +fprintf('%d of %d swr events discarded bc excluded periods in timefilter: d%d\n',... + sum(incSWRs), length(incSWRs), day) +swrTimes = swrTimes(incSWRs,:); +% out.maxthresh = maxthresh; +% out.swrEnd = swrEnd; + +if isempty(swrTimes) + return +end + +%% Get licks in lick-burst intervals +[intraBoutXP, boutTimes] = getLickBoutLicks(animal, [repmat(day,length(eps),1) eps'], ... + varargin{:}); +boutTimes = cell2mat({boutTimes{day}{eps}}'); +intraBoutXP = cell2mat({intraBoutXP{day}{eps}}'); +intraBoutXP = intraBoutXP(:,1); +fprintf('%d XP within %d bursts \n', numel(intraBoutXP), size(boutTimes,1)) + + +%% get xp pre post swr +% get swr events within lick-burst intervals +burstSWRTimes = swrTimes(isIncluded(swrTimes, boutTimes)); +burstSWRTimes = sort(burstSWRTimes(~isnan(burstSWRTimes))); +fprintf(' %d swrs within %d lickbouts. (%.02f pct swrs) \n', size(boutTimes,1), ... + numel(burstSWRTimes), length(burstSWRTimes)/length(swrTimes)*100) +if isempty(burstSWRTimes) + fprintf('no swrs in lick bouts for %d %d.. skipping\n', day, ep) + return +end +% Get the swr-containing licks +% histc given licks as edges: get bin idx. lickedges(binidx) is the prior edge of bin, +% aka the swr-preceding lick +% [~,~,swrLickBinidx] = histcounts(burstSWRTimes, intraBoutXP); +% +% % time since last lick +% swrTimeSinceLick = burstSWRTimes - intraBoutXP(swrLickBinidx); +% if any(swrTimeSinceLick > .250) +% error('wut') +% end +% swr-containing ILI +for b = 1:length(burstSWRTimes) + XPpreSWR(b,1) = max(intraBoutXP(intraBoutXP=burstSWRTimes(b))); +end +XPpreSWR_offset = burstSWRTimes - XPpreSWR; +XPpostSWR_offset = burstSWRTimes - XPpostSWR; + +if run_shuffle + r = rand(length(XPpreSWR),1); + ILI = XPpostSWR - XPpreSWR; % ILI + SWRtimes_shuf = XPpreSWR + (ILI .* r);% rand fraction of each ILI as shuf swr + XPpreSWR_offset_shuf = SWRtimes_shuf - XPpreSWR; + XPpostSWR_offset_shuf = SWRtimes_shuf - XPpostSWR; +end + +out.XPpreSWR = XPpreSWR; +out.XPpostSWR = XPpostSWR; +out.XPpreSWR_offset = XPpreSWR_offset; +out.XPpostSWR_offset = XPpostSWR_offset; +out.XPpreSWR_offset_shuf = XPpreSWR_offset_shuf; +out.XPpostSWR_offset_shuf = XPpostSWR_offset_shuf; +out.SWRTimes_iLB = burstSWRTimes; +out.XPTimes_iLB = intraBoutXP; +end + +function out = init_out() +out.index = []; +out.animal = []; + +out.XPpreSWR = []; +out.XPpostSWR = []; +out.XPpreSWR_offset = []; +out.XPpostSWR_offset = []; +out.XPpreSWR_offset_shuf = []; +out.XPpostSWR_offset_shuf = []; +out.SWRTimes_iLB = []; +out.XPTimes_iLB = []; +end \ No newline at end of file diff --git a/DFAnalysis/dfa_XPtrigAvgRip.m b/DFAnalysis/dfa_XPtrigAvgRip.m new file mode 100644 index 00000000..6406e2b6 --- /dev/null +++ b/DFAnalysis/dfa_XPtrigAvgRip.m @@ -0,0 +1,93 @@ + +function out = dfa_XPtrigAvgRip(idx, timeFilter, varargin) +% get XP trig avg of rip power +% use with singledayanal +% DR_2020-09-30 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'lick'}; +check_required(reqData, varargin) + +win = [-.5 .5]; +eventType = 'ca1rippleskons'; +srate = 1500; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); % init output +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); + +%% get rip pwr trace +evid = find(contains(varargin(1:2:end), eventType)); +o = [1:2:length(varargin)]+1; +swr = varargin{o(evid)}; + +rippwr = []; +time = []; + +for e = 1:eps + rippwr = [rippwr; swr{day}{eps(e)}{1}.powertrace']; + time = [time; swr{day}{eps(e)}{1}.eegtimesvec_ref']; +end + +% zscore +rippwr = rippwr'; +zrippwr = ((rippwr - nanmean(rippwr))./nanstd(rippwr)); + +%% Get licks in lick-burst intervals +[intraBoutXP, boutTimes] = getLickBoutLicks(animal, [repmat(day,length(eps),1) eps'], ... + varargin{:}); +boutTimes = cell2mat({boutTimes{day}{eps}}'); +intraBoutXP = cell2mat({intraBoutXP{day}{eps}}'); +intraBoutXP = intraBoutXP(:,1); +fprintf('%d XP within %d bursts \n', numel(intraBoutXP), size(boutTimes,1)) + +% for each iLB XP, get a zrippwr trace within window +% return this stack +% also return the nanmean and nanstd of the stack +zrippwr_XPtrig = []; +rippwr_XPtrig = []; +for iXP = 1:length(intraBoutXP) + xptimeIdx = min(find(time>=intraBoutXP(iXP))); + if diff([intraBoutXP(iXP) time(xptimeIdx)]) > .001 + sprintf('%.03f diff, skipping',diff([intraBoutXP(iXP) time(xptimeIdx)])) + continue + end + startIdx = xptimeIdx-abs(win(1)*srate)-1; + endIdx = xptimeIdx+abs(win(2)*srate); + try + zrippwr_XPtrig = [zrippwr_XPtrig; zrippwr(startIdx:endIdx)]; + rippwr_XPtrig = [rippwr_XPtrig; rippwr(startIdx:endIdx)]; + catch + disp('window exceeds bounds, skipping') + continue + end +end + +out.zrippwr_XPtrig = zrippwr_XPtrig; +out.mean_zrippwr_XPtrig = nanmean(zrippwr_XPtrig,1); +out.sem_zrippwr_XPtrig = nanstd(zrippwr_XPtrig,1)/sqrt(size(zrippwr_XPtrig,1)); +out.rippwr_XPtrig = rippwr_XPtrig; +out.mean_rippwr_XPtrig = nanmean(rippwr_XPtrig,1); +out.sem_rippwr_XPtrig = nanstd(rippwr_XPtrig,1)/sqrt(size(rippwr_XPtrig,1)); +out.win = win; +out.srate = srate; +out.time = linspace(win(1),win(2),length(out.mean_zrippwr_XPtrig)); +end + +function out = init_out() +out.index = []; +out.animal = []; +out.time = []; + +out.zrippwr_XPtrig = []; +out.mean_zrippwr_XPtrig = []; +out.sem_zrippwr_XPtrig = []; +out.rippwr_XPtrig = []; +out.mean_rippwr_XPtrig = []; +out.sem_rippwr_XPtrig = []; +end \ No newline at end of file diff --git a/DFAnalysis/dfa_pctILB.m b/DFAnalysis/dfa_pctILB.m new file mode 100644 index 00000000..ba7be875 --- /dev/null +++ b/DFAnalysis/dfa_pctILB.m @@ -0,0 +1,97 @@ +function out = dfa_pctILB(idx, timeFilter, varargin) +% get pct swr iLB vs eLB +% use with singledayanal +% notebook: pctILB_2020_10_02.m +% DR_2020-10-02 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'lick'}; +check_required(reqData, varargin) + +stdthresh = [2 4 6 8]; +eventType = 'ca1rippleskons'; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); % init output +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); + +%% Get SWRs in timeFilter + +evid = find(contains(varargin(1:2:end), eventType)); +o = [1:2:length(varargin)]+1; +swr = varargin{o(evid)}; +swrTimes = []; +maxthresh = []; +swrEnd = []; +for e = 1:length(eps) + ep = eps(e); + try + swrTimes = [swrTimes; swr{day}{ep}{1}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}{1}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}{1}.endtime]; + catch + try + swrTimes = [swrTimes; swr{day}{ep}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}.endtime]; + catch + fprintf('no swr events detected for day%d ep%d\n', day, ep) + return + end + end +end +% apply timefilter to swrs +incSWRs = ~isIncluded(swrTimes(:,1), timeFilter); +fprintf('%d of %d swr events discarded bc excluded periods in timefilter: d%d\n',... + sum(incSWRs), length(incSWRs), day) +swrTimes = swrTimes(incSWRs,:); + +maxthresh = maxthresh(incSWRs); +if isempty(swrTimes) + return +end + +%% Get licks in lick-burst intervals +[intraBoutXP, boutTimes] = getLickBoutLicks(animal, [repmat(day,length(eps),1) eps'], ... + varargin{:}); +boutTimes = cell2mat({boutTimes{day}{eps}}'); +% fprintf('%d lick bursts \n', size(boutTimes,1)) + +%% get iLB swrs + +for istd = 1:length(stdthresh) + stdidx = maxthresh > stdthresh(istd); + swrTimes_std = swrTimes(stdidx); + burstSWRTimes = swrTimes_std(isIncluded(swrTimes_std, boutTimes)); +% burstSWRTimes = sort(burstSWRTimes(~isnan(burstSWRTimes))); + fprintf('::>%d Thresh:: %d of %d swrs within %d lickbouts. (%.02f pct swrs) \n', ... + stdthresh(istd), numel(burstSWRTimes), numel(swrTimes_std), size(boutTimes,1),... + length(burstSWRTimes)/length(swrTimes_std)*100) +% if isempty(burstSWRTimes) +% fprintf('no swrs in lick bouts for %d %d.. \n', day, ep) +% keyboard +% end + + out.numSWR(istd) = length(swrTimes_std); + out.numSWRiLB(istd) = length(burstSWRTimes); + out.burstSWRTimes{istd} = burstSWRTimes; + out.pctSWRiLB(istd) = length(burstSWRTimes)/length(swrTimes_std)*100; +% out.maxthresh = maxthresh; +end +end + + +function out = init_out() +out.index = []; +out.animal = []; +out.numSWR = []; +out.numSWRiLB = []; +out.burstSWRTimes = []; +out.pctSWRiLB = []; +% out.maxthresh = []; +end \ No newline at end of file diff --git a/DFAnalysis/dfa_rewTrigSWRXP.m b/DFAnalysis/dfa_rewTrigSWRXP.m new file mode 100644 index 00000000..7a4c06e6 --- /dev/null +++ b/DFAnalysis/dfa_rewTrigSWRXP.m @@ -0,0 +1,101 @@ +function out = dfa_rewTrigSWRXP(idx, timeFilter, varargin) + +% - reward triggered swr and licks +% - for use with singledayanal +% DR_2020-10-13 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'DIO','lick', 'task'}; +check_required(reqData, varargin) +win = [-2 10]; +bin = .001; +eventType = 'ca1rippleskons'; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); + +%% get DIO reward output + +[dioOut, dioOutfields] = get_dio_output(DIO, task, day, eps); + +%% Get SWRs in timeFilter +evid = find(contains(varargin(1:2:end), eventType)); +o = [1:2:length(varargin)]+1; +swr = varargin{o(evid)}; +swrTimes = []; +maxthresh = []; +swrEnd = []; +for e = 1:length(eps) + ep = eps(e); + try + swrTimes = [swrTimes; swr{day}{ep}{1}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}{1}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}{1}.endtime]; + catch + try + swrTimes = [swrTimes; swr{day}{ep}.starttime]; + maxthresh = [maxthresh; swr{day}{ep}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}.endtime]; + catch + fprintf('no swr events detected for day%d ep%d\n', day, ep) + return + end + end +end + +% apply timefilter to swrs +incSWRs = ~isIncluded(swrTimes(:,1), timeFilter); +fprintf('%d of %d swr events discarded bc excluded periods in timefilter: d%d\n',... + sum(incSWRs), length(incSWRs), day) +swrTimes = swrTimes(incSWRs,:); +% out.maxthresh = maxthresh; +% out.swrEnd = swrEnd; + +if isempty(swrTimes) + return +end + +%% Get licks in lick-burst intervals +% [intraBoutXP, boutTimes] = getLickBoutLicks(animal, [repmat(day,length(eps),1) eps'], ... +% varargin{:}); +% boutTimes = cell2mat({boutTimes{day}{eps}}'); +% intraBoutXP = cell2mat({intraBoutXP{day}{eps}}'); +% intraBoutXP = intraBoutXP(:,1); +% fprintf('%d XP within %d bursts \n', numel(intraBoutXP), size(boutTimes,1)) +xp = []; +for e = 1:length(eps) + ep = eps(e); + try + xp = [xp; lick{day}{ep}.eventtime]; + catch + xp = [xp; lick{day}{ep}.starttime]; % legacy name + end +end + +time = -abs(win(1))-(bin/2) : bin : win(2)+(bin/2); + +swr_prth = cell2mat(arrayfun(@(r) histc(swrTimes, r + time), dioOut(:,5), 'un', 0)')'; +xp_prth = cell2mat(arrayfun(@(r) histc(xp, r + time), dioOut(:,5), 'un', 0)')'; + +%% output +out.rewTrigSWR = swr_prth; +out.rewTrigXP = xp_prth; +out.time = time; + +end + +function out = init_out() +out.index = []; +out.animal = []; + +out.rewTrigSWR = []; +out.rewTrigXP = []; +out.time = []; +end + diff --git a/DFAnalysis/dfa_ripPos.m b/DFAnalysis/dfa_ripPos.m new file mode 100644 index 00000000..7d32d009 --- /dev/null +++ b/DFAnalysis/dfa_ripPos.m @@ -0,0 +1,107 @@ +function out = dfa_ripPos(idx, timeFilter, varargin) + +% get rip position +% use with singledayanal +% DR_2020-10-12 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'pos','lick'}; +check_required(reqData, varargin) + +eventType = 'ca1rippleskons'; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); % init output +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); + +%% Get SWRs in timeFilter + +evid = find(contains(varargin(1:2:end), eventType)); +o = [1:2:length(varargin)]+1; +swr = varargin{o(evid)}; +SWRTimes = []; +maxthresh = []; +swrEnd = []; +swr_ep = []; +for e = 1:length(eps) + ep = eps(e); + try + SWRTimes = [SWRTimes; swr{day}{ep}{1}.starttime]; + swr_ep = [swr_ep; repmat(ep,length(SWRTimes),1)]; + maxthresh = [maxthresh; swr{day}{ep}{1}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}{1}.endtime]; + catch + try + SWRTimes = [SWRTimes; swr{day}{ep}.starttime]; + swr_ep = [swr_ep; repmat(ep,length(SWRTimes),1)]; + maxthresh = [maxthresh; swr{day}{ep}.maxthresh]; + swrEnd = [swrEnd; swr{day}{ep}.endtime]; + catch + fprintf('no swr events detected for day%d ep%d\n', day, ep) + return + end + end +end +% apply timefilter to swrs +incSWRs = ~isIncluded(SWRTimes(:,1), timeFilter); +fprintf('%d of %d swr events discarded bc excluded periods in timefilter: d%d\n',... + sum(incSWRs), length(incSWRs), day) +SWRTimes = SWRTimes(incSWRs,:); +SWRTimes_maxThresh = maxthresh(incSWRs,:); +swr_ep = swr_ep(incSWRs); + +if isempty(SWRTimes) + return +end + +%% Get pos +posdata = []; +for e = 1:length(eps) + ep = eps(e); + posdata = [posdata; pos{day}{eps(e)}.data]; + posfields = pos{day}{eps(e)}.fields; +end +xstring = 'x-loess'; +xcol = find(cell2mat(cellfun(@(x) strcmp(x,xstring), strsplit(posfields, ' '), ... + 'UniformOutput', false))); +ystring = 'y-loess'; +ycol = find(cell2mat(cellfun(@(x) strcmp(x,ystring), strsplit(posfields, ' '), ... + 'UniformOutput', false))); + +swr_posidx = knnsearch(posdata(:,1), SWRTimes(:,1)); +SWRTimes_xpos = posdata(swr_posidx,xcol); +SWRTimes_ypos = posdata(swr_posidx,ycol); + +%% Get lick-burst intervals +[intraBoutXP, boutTimes] = getLickBoutLicks(animal, [repmat(day,length(eps),1) eps'], ... + varargin{:}); +boutTimes = cell2mat({boutTimes{day}{eps}}'); + +%% get swr events within lick-burst intervals +SWRTimes_iLB = isIncluded(SWRTimes, boutTimes); + +%% output +out.txy = posdata(:,[1 xcol ycol]); +out.SWR_time = SWRTimes; +out.SWR_iLB = SWRTimes_iLB; +out.SWR_xpos = SWRTimes_xpos; +out.SWR_ypos = SWRTimes_ypos; +out.SWR_maxThresh = SWRTimes_maxThresh; +end + +function out = init_out() +out.index = []; +out.animal = []; + +out.txy = []; +out.SWR_time = []; +out.SWR_iLB = []; +out.SWR_xpos = []; +out.SWR_ypos = []; +out.SWR_maxThresh = []; +end \ No newline at end of file diff --git a/DFAnalysis/dfa_template.m b/DFAnalysis/dfa_template.m new file mode 100644 index 00000000..0d8b50a2 --- /dev/null +++ b/DFAnalysis/dfa_template.m @@ -0,0 +1,33 @@ +function out = dfa_template(idx, timeFilter, varargin) + +% reward triggered swr and licks +% use with singledayanal +% DR_2020-10-13 + +fprintf('day %d \n',idx(1)) +reqData = {'ca1rippleskons', 'pos','lick'}; +check_required(reqData, varargin) + +eventType = 'ca1rippleskons'; +if ~isempty(varargin) + assign(varargin{:}); +end + +out = init_out(); +out.index = idx; +out.animal = animal; +day = idx(1); +eps = idx(2:end); + + + + +end + +function out = init_out() +out.index = []; +out.animal = []; + + +end +