Skip to content

Commit

Permalink
initial attempt at coactive z swr burst v noburst
Browse files Browse the repository at this point in the history
  • Loading branch information
droumis committed Oct 25, 2019
1 parent d8ca332 commit 8182109
Show file tree
Hide file tree
Showing 6 changed files with 274 additions and 16 deletions.
3 changes: 2 additions & 1 deletion DFAnalysis/dfa_perripspikingcorr.m
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -64,6 +64,7 @@
end

% ones = included, zeros = excluded

includetimes = ~isExcluded(times, excludeperiods);
includetimes = includetimes(:);

Expand Down
178 changes: 178 additions & 0 deletions DFAnalysis/dfa_suCoactiveXcorr.m
Original file line number Diff line number Diff line change
@@ -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


8 changes: 2 additions & 6 deletions DFFunctions/dfa_lickBoutSpikeCorr.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion DFFunctions/getLickBout.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
31 changes: 23 additions & 8 deletions DFFunctions/load_filter_params.m
Original file line number Diff line number Diff line change
Expand Up @@ -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';
Expand Down Expand Up @@ -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,
Expand Down
65 changes: 65 additions & 0 deletions notebooks/coactiveBurstNoBurst_20191024.m
Original file line number Diff line number Diff line change
@@ -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?

0 comments on commit 8182109

Please sign in to comment.