From 8182109466470ce670b89fbc3650952bf3ca1f35 Mon Sep 17 00:00:00 2001 From: Demetris Date: Thu, 24 Oct 2019 20:45:40 -0700 Subject: [PATCH] initial attempt at coactive z swr burst v noburst --- DFAnalysis/dfa_perripspikingcorr.m | 3 +- DFAnalysis/dfa_suCoactiveXcorr.m | 178 ++++++++++++++++++++++ DFFunctions/dfa_lickBoutSpikeCorr.m | 8 +- DFFunctions/getLickBout.m | 5 +- DFFunctions/load_filter_params.m | 31 +++- notebooks/coactiveBurstNoBurst_20191024.m | 65 ++++++++ 6 files changed, 274 insertions(+), 16 deletions(-) create mode 100644 DFAnalysis/dfa_suCoactiveXcorr.m create mode 100644 notebooks/coactiveBurstNoBurst_20191024.m diff --git a/DFAnalysis/dfa_perripspikingcorr.m b/DFAnalysis/dfa_perripspikingcorr.m index 576a7e37..e9778905 100644 --- a/DFAnalysis/dfa_perripspikingcorr.m +++ b/DFAnalysis/dfa_perripspikingcorr.m @@ -4,7 +4,7 @@ function out = dfa_perripspikingcorr(index, excludeperiods, spikes, eventscons, varargin) % Calculate per rip corr between the spikes from ntA, ntB % exclude periods and other filters apply to the rips - +% COMBINES SPIKES FROM ALL CLUSTERS ON EACH NTRODE %%% Inputs % index [day epoch tetrode cell tetrode cell] @@ -64,6 +64,7 @@ end % ones = included, zeros = excluded + includetimes = ~isExcluded(times, excludeperiods); includetimes = includetimes(:); diff --git a/DFAnalysis/dfa_suCoactiveXcorr.m b/DFAnalysis/dfa_suCoactiveXcorr.m new file mode 100644 index 00000000..1d82381d --- /dev/null +++ b/DFAnalysis/dfa_suCoactiveXcorr.m @@ -0,0 +1,178 @@ +%{ +meant to be an updated version on calcxcorrmeasures +that is compatible with my changes to singlecellanal (which is actually +feeding in multiple su).. +.. i.e. the way data is fed in as a varargin + +also see work i had done for mu xcorr in dfa_perripspikingcorr + +also looks like i did something similar with: dfa_lickBoutSpikeCorr.. +except that i took out the coactive z stuff. also it looks like i was using +that to run MU + +%} + +function out = dfa_suCoactiveXcorr(idx, excludetimes, varargin) +%function out = calcxcorrmeasures(index, excludetimes, spikes, varargin) +% Calculates the excess correlation and RMS time lag for the specified cell +% pairs using only spikes not excluded by the excludetimes +% +% +% Options: +% 'bin', n binsize in sec. Default 0.002 (2 ms) +% 'tmax', n maximum time for cross correlation. Default 1 sec. +% 'edgespikes', 0 or 1 +% 0 will enforce the exclude periods and will thus +% exclude all spikes outside of the exclude windows. (default) +% 1 will include all spikes where at least one of the spikes is +% not excluded +% 'sw1', n The smaller smoothing width for the excess correlation measure +% Default 0.005 (5 ms) +% 'sw2', n The larger smoothing width for the excess correlation measure +% Default 0.250 (250 ms) +% 'rmstmax', n The maximum absolute time to be used for the rms calculation +% Default 0.1 (100 ms) +% 'rmsmincounts', n The minimum number of counts in the histogram for the rms +% measurement to be carried out. +% Default 50. +% 'calclinfields', 1 or 0 -- set to 1 to calculate linearized place fields +% 'calctrajxcorr', 1 or 0 -- set to 1 to calculate one cross correlation per +% trajectory (default 0); +% + +calclinfields = 0; +calctrajxcorr = 0; +appendindex = 0; +bin = 0.002; +tmax = 1; +edgespikes = 0; +sw1 = 0.005; +sw2 = 0.250; +rmstmax = 0.1; +rmsmincounts = 50; + +if ~isempty(varargin) + assign(varargin{:}) +end + +% set the defaults +out.time = []; +out.c1vsc2 = []; +out.ac1 = []; +out.ac2 = []; +out.index = idx; +out.ec = NaN; +out.rms = NaN; +out.probcoa = NaN; +out.coactivez_numsp = NaN; +out.coactivez = NaN; +out.lf1 = []; +out.lf2 = []; + +% for each cell we calculate the cross correlation +try + t1 = spikes{idx(1)}{idx(2)}{idx(3)}{idx(4)}.data; + t2 = spikes{idx(1)}{idx(2)}{idx(5)}{idx(6)}.data; +catch + % if either of those produced an error, we return NaNs + return; +end + +%apply the exclude rule +t1inc = []; +t2inc = []; +if (length(t1)) + t1inc = t1(find(~isExcluded(t1(:,1), excludetimes)),1); +else + return +end +if (length(t2)) + t2inc = t2(find(~isExcluded(t2(:,1), excludetimes)),1); +else + return +end + +ac1 = spikexcorr(t1inc, t1inc, bin, tmax); +ac2 = spikexcorr(t2inc, t2inc, bin, tmax); +out.ac1 = ac1.c1vsc2; +out.ac2 = ac2.c1vsc2; +out.time = ac1.time; + +out.index = idx; +if (~calctrajxcorr) + xc = spikexcorr(t1inc, t2inc, bin, tmax); + % if we want to include edge spikes, we need to add in the correlation of + % the excluded t1 spikes with the included t2spikes + if ((edgespikes) & (~isempty(xc.time))) + t1ex = t1(find(isExcluded(t1(:,1), excludetimes))); + if (~isempty(t1ex)) + tmpxc = spikexcorr(t1ex, t2inc, bin, tmax); + % add these values to the original histogram + xc.c1vsc2 = xc.c1vsc2 + tmpxc.c1vsc2; + end + end + out.c1vsc2 = xc.c1vsc2; + % compute the excess correlation at 0 lag + out.ec = excesscorr(xc.time, xc.c1vsc2, xc.nspikes1, xc.nspikes2, sw1, sw2); + out.rms = xcorrrms(xc.time, xc.c1vsc2, rmstmax, rmsmincounts); + + % get the probability that both were coactive in the included intervals + n1 = nInInterval(t1inc, excludetimes, 0); + n2 = nInInterval(t2inc, excludetimes, 0); + out.probcoa = mean(logical(n1) & logical(n2)); + + % get the zscore for this probabliity of coactivity + out.coactivez = coactivezscore(n1, n2); + out.coactivez_numsp = coactivezscore(n1, n2, 'anyspikes', 0); +else + % calculated individual trajectory cross correlations + statematrix = linpos{idx(1)}{idx(2)}.statematrix; + statevector = statematrix.traj; + statevector(find(isExcluded(statematrix.time, excludetimes))) = -1; + % look up the state for each spike + posindexfield = 7; + t1 = t1(:,[1 posindexfield]); + t2 = t2(:,[1 posindexfield]); + t1(:,3) = statevector(t1(:,2)); + t2(:,3) = statevector(t2(:,2)); + % get rid of -1 trajectories and excluded times + t1 = t1(find(t1(:,3) ~= -1),:); + t2 = t2(find(t2(:,3) ~= -1),:); + t2 = t2(find(~isExcluded(t2(:,1), excludetimes)), :); + for i = 1:max(statevector) + t1tmp = t1(find(t1(:,3) == i), 1); + t1inc = t1tmp(find(~isExcluded(t1tmp, excludetimes))); + t2tmp = t2(find(t2(:,3) == i), 1); + t2inc = t2tmp(find(~isExcluded(t2tmp, excludetimes))); + xc = spikexcorr(t1inc, t2inc, bin, tmax); + if ((edgespikes) & (~isempty(xc.time))) + t1ex = t1tmp(find(isExcluded(t1tmp, excludetimes))); + if (~isempty(t1ex)) + tmpxc = spikexcorr(t1ex, t2inc, bin, tmax); + % add these values to the original histogram + xc.c1vsc2 = xc.c1vsc2 + tmpxc.c1vsc2; + end + end + out.c1vsc2{i} = xc.c1vsc2; + % compute the excess correlation at 0 lag + out.ec(i) = excesscorr(xc.time, xc.c1vsc2, xc.nspikes1, xc.nspikes2, sw1, sw2); + out.rms(i) = xcorrrms(xc.time, xc.c1vsc2, rmstmax, rmsmincounts); + + % get the probability that both were coactive in the included intervals + n1 = nInInterval(t1inc, excludetimes, 0); + n2 = nInInterval(t2inc, excludetimes, 0); + out.probcoa(i) = mean(logical(n1) & logical(n2)); + + % get the zscore for this probabliity of coactivity + out.coactivez(i) = coactivezscore(n1, n2); + out.coactivez_numsp(i) = coactivezscore(n1, n2, 'anyspikes', 0); + end +end + +if (calclinfields) + % calculate the linear fields + out.lf1 = filtercalclinfields(idx(1:4), excludetimes, spikes, linpos, 'phasedist', 1); + out.lf2 = filtercalclinfields([idx(1:2) idx(5:6)], excludetimes, spikes, linpos, 'phasedist', 1); +end + + diff --git a/DFFunctions/dfa_lickBoutSpikeCorr.m b/DFFunctions/dfa_lickBoutSpikeCorr.m index 5541701c..57d7de4c 100644 --- a/DFFunctions/dfa_lickBoutSpikeCorr.m +++ b/DFFunctions/dfa_lickBoutSpikeCorr.m @@ -11,16 +11,12 @@ end pconf = paramconfig; -bin = 0.01; % 10 ms bins -sw1 = bin*3; -% sw1 = 0.005; +bin = 0.01; % seconds +sw1 = bin*3; % seconds sw2 = 0.250; rmstmax = 0.1; rmsmincounts = 50; - tmax = .5 ; -numshuff = 1000; -shuff = 100; % ms if ~isempty(varargin) assign(varargin{:}); end diff --git a/DFFunctions/getLickBout.m b/DFFunctions/getLickBout.m index aa6f57f2..0c1fba5a 100644 --- a/DFFunctions/getLickBout.m +++ b/DFFunctions/getLickBout.m @@ -7,9 +7,12 @@ % dio times along with.. /home/droumis/Src/Matlab/filterframework_dr/Functions/get_licks.m +'datadir' arg1 currently needs to be there bc of the way setfiltertime +calls the time funcs... but i don't use it so outside of setfiltertime i +should just pass an empty arg there.. %} -function out = getLickBout(animal, epochs, varargin) +function out = getLickBout(datadir, animal, epochs, varargin) maxIntraBurstILI = .25; % max burst ili threshold in seconds minBoutLicks = 3; %filter out bouts with less than boutNum licks diff --git a/DFFunctions/load_filter_params.m b/DFFunctions/load_filter_params.m index 0b4b9ff5..7101430c 100644 --- a/DFFunctions/load_filter_params.m +++ b/DFFunctions/load_filter_params.m @@ -162,19 +162,21 @@ options = {'savefigas', 'png', 'eventName',eventName}; case 'lickbouts' - lickGap = 0.5; - boutNum = 10; + maxIntraBurstILI = 0.25; % max burst ili threshold in seconds + minBoutLicks = 3; %filter out bouts with less than boutNum licks % timefilt1: lick bouts, - timefilter{end+1} = {'getLickBout', '($lickBout == 1)', 'lickGap', lickGap, ... - 'boutNum', boutNum}; + timefilter{end+1} = {'getLickBout', '($lickBout == 1)', ... + 'maxIntraBurstILI', maxIntraBurstILI, ... + 'minBoutLicks', minBoutLicks}; case 'nolickbouts' % timefilt: immoble not lick bout - lickGap = 0.5; - boutNum = 10; + maxIntraBurstILI = 0.25; % max burst ili threshold in seconds + minBoutLicks = 3; % filter out bouts with less than boutNum licks + minTimeFromLick = .5; % time duration from closest lick timefilter{end+1} = {'getLickBout', ... - '(($lickBout == 0) & ($velocity < 1) & ($timeFromLick > .5))', ... - 'lickGap', lickGap, 'boutNum', boutNum}; + sprintf('($lickBout == 0) & ($timeFromLick > %d)',minTimeFromLick), ... + 'maxIntraBurstILI', maxIntraBurstILI, 'minBoutLicks', minBoutLicks}; case 'dfa_lickphaseSUclustering' eventName = 'lick'; @@ -366,6 +368,19 @@ iterator = 'singlecellanal'; filtfunction = 'calcxcorrmeasures'; datatypes = {'spikes', 'ca1rippleskons', 'pos', 'task'}; + + case 'dfa_suCoactiveXcorr' + iterator = 'singlecellanal'; + filtfunction = 'dfa_suCoactiveXcorr'; + + cellpairfilter = {'allcomb', ... + '($numspikes > 100) && (isequal($area, ''ca1'')) && (all(cellfun(''isempty'',(arrayfun(@(x) strfind(x,''mua''), $tags, ''un'', 0)))))', ... + '($numspikes > 100) && (isequal($area, ''ca1'')) && (all(cellfun(''isempty'',(arrayfun(@(x) strfind(x,''mua''), $tags, ''un'', 0)))))'}; + + eventName = 'swr'; + iterator = 'singlecellanal'; + datatypes = {'spikes'}; + options = {'savefigas', 'png', 'eventName',eventName}; case 'mua_calcxcorrmeasures' % requires 'tetrodepairs', Fp.tetpairfilter, diff --git a/notebooks/coactiveBurstNoBurst_20191024.m b/notebooks/coactiveBurstNoBurst_20191024.m new file mode 100644 index 00000000..a83a2dfb --- /dev/null +++ b/notebooks/coactiveBurstNoBurst_20191024.m @@ -0,0 +1,65 @@ + + +create_filter = 0; +run_ff = 0; +save_ffdata = 0; +load_ffdata = 1; + +%% get the lickphase for each swr +% data filter params +Fp.animals = {'D10'}; %, 'D12', 'D13', 'JZ1', 'JZ2', 'JZ3', 'JZ4'}; %, 'JZ2', 'JZ4'}; +Fp.filtfunction = 'dfa_suCoactiveXcorr'; +% condition = 'lickburst'; +condition = 'noburst'; +% Fp.params = {'savefigs', 'wtrackdays', 'ripples', 'lickbouts', Fp.filtfunction}; +Fp.params = {'savefigs', 'wtrackdays', 'ripples', 'nolickbouts', Fp.filtfunction}; +% FF +Fp = load_filter_params(Fp); +if create_filter + F = createfilter('animal', Fp.animals, 'epochs', Fp.epochfilter, ... + 'excludetime', Fp.timefilter,'cellpairs', Fp.cellpairfilter, 'iterator',Fp.iterator); + F = setfilterfunction(F, Fp.filtfunction, Fp.datatypes, Fp.options{:}); +end +if run_ff; F = runfilter(F); + for d = 1:length(F); F(d).datafilter_params = Fp; end +end +if save_ffdata + save_data(F, Fp.paths.filtOutputDirectory, Fp.paths.filenamesave, ... + 'filetail', sprintf('_%s_%s', Fp.epochEnvironment, condition)); +end +if load_ffdata + Flb = load_data(Fp.paths.filtOutputDirectory, Fp.paths.filenamesave, Fp.animals, ... + 'filetail', sprintf('_%s_%s', Fp.epochEnvironment, 'lickburst')); + Fno = load_data(Fp.paths.filtOutputDirectory, Fp.paths.filenamesave, Fp.animals, ... + 'filetail', sprintf('_%s_%s', Fp.epochEnvironment, 'noburst')); +end +%% +% i need to combine epoch results for each cell pair (take mean?) +% plot the cumulative coactive z for each condition. run kstest? +ani = 1; +Flb_struct = cell2mat(Flb(ani).output{1}'); +Flb_coZ = [Flb_struct.coactivez]'; + +Fno_struct = cell2mat(Fno(ani).output{1}'); +Fno_coZ = [Fno_struct.coactivez]'; +%% +[lickh, lickb] = histcounts(Flb_coZ,1000,'Normalization', 'cdf'); +plot(lickb(1:end-1)+diff(lickb) / 2, lickh, 'k') +hold on; +[noh, nob] = histcounts(Fno_coZ,1000,'Normalization', 'cdf'); +plot(nob(1:end-1)+diff(nob) / 2, noh, 'b') +hold off +[H, pValue, KSstatistic] = kstest2(Fno_coZ, Flb_coZ); +%% + +histogram(Flb_coZ, 40, 'Normalization', 'probability') +hold on; +histogram(Fno_coZ, 40, 'Normalization', 'probability', 'DisplayStyle','stairs', 'linewidth', 2) +hold off +legend({'lickburst swr', 'nonburst swr'}) +title('coactiveZ ca1') +save_figure([pconf.andef{4} '/coactiveZburstNoburst/'], 'coactiveZburstNoburst_ca1'); + +% then for each pair get diff of coactiveZ for burst-noburst... and then +% generate a condition-shuffled distribution of diff coactiveZ.. to +% compare that to?