Adding dentate spike and CSD analysis

Author: emilyasterjones

Filters/ejdsfilter4.mat

@@ -0,0 +1,76 @@
+function out = DSMUA(index, excludeperiods, mua, dspikes, chinfo, varargin)
+%parse the options
+win = [-0.4 0.4];
+binsize = 0.1;
+
+Fs = 1000;
+spikeFs = 30000;
+
+for option = 1:2:length(varargin)-1
+    if ischar(varargin{option})
+        switch(varargin{option})
+            case 'event_window'
+                win = varargin{option+1};
+            case 'binsize'
+                binsize = varargin{option+1};
+            case 'trig'
+                trigcrit = varargin{option+1};
+            case 'probe'
+                probecrit = varargin{option+1};
+            otherwise
+                error(['Option ',varargin{option},' unknown.']);
+        end
+    else
+        error('Options must be strings, followed by the variable');
+    end
+end
+
+% search through chinfo struct to find channels meeting trig and probe criteria
+trigindex = evaluatefilter(chinfo,trigcrit);
+probeindex = evaluatefilter(chinfo,probecrit) ;
+out.trigindex = trigindex;
+out.probeindex = probeindex;
+
+fulltimes = mua{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)}.starttime:1/Fs:mua{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)}.endtime;
+validtimes = ~isExcluded(fulltimes, excludeperiods);
+
+%iterate over trigger channels
+for t = 1:size(trigindex,1)
+    
+    %get dspike trigger times
+    d = dspikes{trigindex(t,1)}{trigindex(t,2)}{trigindex(t,3)};
+    valid = find(~isExcluded(d.starttime,excludeperiods) & ~isExcluded(d.endtime,excludeperiods));
+    
+    d.starttime = d.starttime(valid);
+    d.endtime = d.endtime(valid);
+    d.peaktime = mua{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)}.starttime + d.peakind / Fs;
+    triggers = d.peaktime(valid);
+    out.dstimes{t} = triggers;
+    
+    %iterate through probe channels
+    for p = 1:size(probeindex,1)
+        pmua = mua{probeindex(p,1)}{probeindex(p,2)}{probeindex(p,3)};
+        spiketimes = pmua.spiketimes/spikeFs;
+        bins = fulltimes(1):binsize:fulltimes(end);
+        validinds = ~isExcluded(spiketimes,excludeperiods);
+        
+        bininds = ~isExcluded(bins,excludeperiods);
+        out.baserate{t}{p} = sum(validinds)/sum(validtimes)*Fs; %FR over all immob
+        out.sdbaserate{t}{p} = std(histcounts(spiketimes(validinds),bins(bininds))/binsize);
+        
+        bins = win(1):binsize:win(2);
+        out.psth{t}{p} = NaN(length(triggers),(length(bins)-1)); %no 10x
+        out.peakrate{t}{p} = NaN(length(triggers),1);
+        out.DSrate{t}{p} = NaN(length(triggers),1);
+        out.normrate{t}{p} = NaN(length(triggers),1);
+        
+        %iterate over dspikes
+        for r = 1:length(triggers)
+            bins = win(1)+triggers(r):binsize:triggers(r)+win(2); %binsize
+            out.psth{t}{p}(r,:) = histcounts(spiketimes,bins)/binsize; %binsize
+            out.peakrate{t}{p}(r) = length(spiketimes(spiketimes>triggers(r) & spiketimes<triggers(r)+binsize))/binsize;
+            out.DSrate{t}{p}(r) = length(spiketimes(spiketimes>d.starttime(r) & spiketimes<d.endtime(r)))/(d.endtime(r) - d.starttime(r));
+            out.normrate{t}{p}(r) = (out.DSrate{t}{p}(r)-out.baserate{t}{p})/out.sdbaserate{t}{p}';
+        end
+    end
+end

Rest Analysis/DSMUA.m

@@ -0,0 +1,34 @@
+function out = DSdescript(index, excludeperiods, dspikes, eeg, chinfo, varargin)
+% out = ripdescript(index, excludeperiods, ripples, options)
+%  Computes basic characteristics of detected ripple events
+
+for option = 1:2:length(varargin)-1
+    if ischar(varargin{option})
+        switch(varargin{option})
+            case 'trig'
+                trigcrit = varargin{option+1};
+            otherwise
+                error(['Option ',varargin{option},' unknown.']);
+        end
+    else
+        error('Options must be strings, followed by the variable');
+    end
+end
+
+%open dspikes
+trigindex = evaluatefilter(chinfo,trigcrit);
+out.trigindex = trigindex;
+d = dspikes{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)};
+e = eeg{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)};
+valid = find(~isExcluded(d.starttime,excludeperiods) & ~isExcluded(d.endtime,excludeperiods));
+
+%calculate DS rate in 1min bins
+Fs = e.samprate;
+fulltimes = e.starttime:1/Fs:e.endtime;
+validtimeinds = find(~isExcluded(fulltimes, excludeperiods));
+timeseq = [validtimeinds(1:60000:end); validtimeinds(end)];
+timecounts = histcounts(d.startind(valid),timeseq)/60; %rips per 1min bin
+
+out.dspeakamps = d.peak(valid);
+out.dslength = d.endind(valid) - d.startind(valid);   %in ms
+out.timecounts = timecounts;

Rest Analysis/DSdescript.m

@@ -0,0 +1,95 @@
+function out = RTCSD(index, excludeperiods, eeg, ripples, chinfo, varargin)
+%parse the options
+win = [-0.2 0.2];
+minstd = 5;
+interripwin = 0;  %does not eliminate ripples that occur in doublets or triplets
+Fs = 1000;
+
+for option = 1:2:length(varargin)-1
+    if ischar(varargin{option})
+        switch(varargin{option})
+            case 'win'
+                win = varargin{option+1};
+            case 'minstd'
+                minstd = varargin{option+1};
+            case 'trig'
+                trigcrit = varargin{option+1};
+            case 'freqband'
+                freqs = varargin{option+1};
+            otherwise
+                error(['Option ',varargin{option},' unknown.']);
+        end
+    else
+        error('Options must be strings, followed by the variable');
+    end
+end
+
+out.win = win;
+
+% search through chinfo struct to find channels meeting trig criteria
+trigindex = evaluatefilter(chinfo,trigcrit);
+out.trigindex = trigindex;
+
+%get EEG & rest times
+epoch_eeg = eeg{trigindex(1,1)}{trigindex(1,2)};
+Fs_temp = length(epoch_eeg{1}.data)/(epoch_eeg{1}.endtime-epoch_eeg{1}.starttime);
+fulltimes = epoch_eeg{1}.starttime:1/Fs_temp:epoch_eeg{1}.endtime;
+validtimes = ~isExcluded(fulltimes, excludeperiods);
+
+%if freq band specified, filter trace
+if exist('freqs','var')
+    for i=1:length(epoch_eeg)
+        coeff = designeegfilt(Fs,freqs(1),freqs(2));
+        data = double(epoch_eeg{i}.data);
+        epoch_eeg{i}.data = int16(filtfilt(coeff,1,data));
+    end
+end
+
+%smoothing via triangular kernel
+for i = 2:length(epoch_eeg)-1
+    %smoothing via triangular kernel
+    b = epoch_eeg{i}.data;
+    a = epoch_eeg{i-1}.data;
+    c = epoch_eeg{i+1}.data;
+    smooth_eeg{i}.data = (a+2*b+c)/4;
+end
+%second spatial derivative
+csd = NaN(length(epoch_eeg),length(epoch_eeg{1}.data)); %channels x time
+for i = 3:length(epoch_eeg)-2
+    %smoothing via triangular kernel
+    b = smooth_eeg{i}.data;
+    a = smooth_eeg{i-1}.data;
+    c = smooth_eeg{i+1}.data;
+    csd(i,:) = (a-2*b+c)/4;
+end
+out.csd = csd;
+
+%mean & SD over rest trace
+csd_base = csd(:,validtimes(1:length(csd)));
+out.baseline = mean(csd_base,2);
+out.sd = std(csd_base,0,2);
+
+%get ripple trigger times
+rip = ripples{trigindex(1,1)}{trigindex(1,2)}{trigindex(1,3)};
+[valid, ~] = getvalidrips(ripples,index,trigindex, excludeperiods, win, minstd, interripwin);
+rip.starttime = rip.starttime(valid);
+rip.endtime = rip.endtime(valid);
+triggers = rip.threshind(valid);
+out.ripinds = triggers;
+
+%bin by ripple triggers
+ripcsd = NaN(length(triggers),length(epoch_eeg),length(win(1):1/Fs:win(2))); %rips x channels x time
+for r = 1:length(triggers)
+    ripcsd(r,:,:) = csd(:,win(1)*Fs+triggers(r):triggers(r)+win(2)*Fs);
+end
+out.ripcsd = ripcsd;
+
+%calculate peaks
+out.peakripsource = max(ripcsd,[],3); %rips x channels
+out.peakripsink = min(ripcsd,[],3); %rips x channels
+
+%Z-score
+out.zripsource = (out.peakripsource - out.baseline') ./ out.sd';
+out.zripsink = (out.peakripsink - out.baseline') ./ out.sd';
+
+end