Initial commit

Dentate Spikes/ag_extractdspikes.m

@@ -0,0 +1,137 @@
+function ag_extractdspikes(directoryname, fileprefix, day, min_suprathresh_duration, nstd, varargin)
+%based on extractripples
+%extract dentate spikes
+%with no filtering, extracts events that are large - mindur 10ms, 5SD above
+%baseline
+%exclude any events that coincide with ripples, since the deflection is
+%likely to just be the sharp wave
+%save one file per session, just like ripples
+%
+% automatically runs through all eps per day
+
+% Outputs:
+%dspikes 	- structue with various fields, including the following which
+%			describe each dentate spike.
+%	starttime - time of beginning of dspike
+%	endtime	  - time of end of dspike
+%	midtime   - time of midpoint of energy of event
+%	peak	  - peak height of waveform)
+%	maxthresh - the largest threshold in stdev units at which this dspike
+%			would still be detected.
+%	energy	  - total sum squared energy of waveform
+%	startind  - index of start time in eeg structure
+%	endind    - index of end time in eeg structure
+%	midind    - index of middle time in eeg structure
+%
+%
+%defaults
+stdev = 0;
+baseline = 0;
+maxpeakval = 1000;
+for option = 1:2:length(varargin)-1
+    if isstr(varargin{option})
+        switch(varargin{option})
+            case 'stdev'
+                stdev = varargin{option+1};
+            case 'baseline'
+                baseline = varargin{option+1};
+            case 'maxpeakval'
+                maxpeakval = varargin{option+1};
+            otherwise
+                error(['Option ',varargin{option},' unknown.']);
+        end
+    else
+        error('Options must be strings, followed by the variable');
+    end
+end
+
+% define the standard deviation for the Gaussian smoother which we
+% apply before thresholding (this reduces sensitivity to spurious
+% flucutations in the ripple envelope)
+smoothing_width = 0.004; % 4 ms
+
+d = day;
+
+% move to the directory
+cd(directoryname);
+
+tmpflist = dir(sprintf('*eeg%02d-*.mat', day));
+%get rip list from pyr1 chan
+load(sprintf('%s/%schinfo.mat',directoryname,fileprefix))
+location = '(isequal($area,''ca1'')&& contains($layer,''pyr 1'')  )';%
+pyr1chan = evaluatefilter(chinfo(day),location);
+pyr1chan = unique(pyr1chan(:,3));
+
+for i = 1:length(tmpflist) %iterate through all eps, all chans
+    % load the ripple file
+    load(tmpflist(i).name);
+    % get the epoch number
+    dash = find(tmpflist(i).name == '-');
+    e = str2num(tmpflist(i).name((dash(1)+1):(dash(2)-1)));
+    t = str2num(tmpflist(i).name((dash(2)+1:dash(2)+2)));
+    % raw trace envelope
+    renv = abs(hilbert(double(eeg{d}{e}{t}.data)));
+
+    % smooth the envelope:
+    samprate = eeg{d}{e}{t}.samprate;
+    kernel = gaussian(smoothing_width*samprate, ceil(8*smoothing_width*samprate));
+    renv = smoothvect(renv, kernel);
+    % calculate the threshold in uV units
+    baseline = mean(renv);
+    stdev = std(renv);
+    thresh = baseline + nstd * stdev;
+    % find the events
+    % calculate the duration in terms of samples
+    mindur = round(min_suprathresh_duration * samprate);
+
+    % extract the events if this is a valid trace
+    if (thresh > 0) & any(find(renv<baseline))
+        tmpevents = extracteventsnew(renv, thresh, baseline, 0, mindur, 0)';
+
+        %eliminate any events within 100ms of ripples; they are
+        %sharpwaves
+        load(sprintf('%s/%sripples%02d.mat',directoryname, fileprefix,day'))
+        ripinds = [ripples{d}{e}{pyr1chan}.startind - 100, ripples{d}{e}{pyr1chan}.endind + 100];
+        valids = ~isExcluded(tmpevents(:,8),ripinds); %midind anywhere in window of rips
+        % Assign the fields
+        % start and end indeces
+        ds.startind = tmpevents(valids,1);
+        ds.endind = tmpevents(valids,2);
+        % middle of energy index
+        ds.midind = tmpevents(valids,8);
+
+        %convert the samples to times for the first three fields
+        ds.starttime = eeg{d}{e}{t}.starttime + ds.startind / samprate;
+        ds.endtime = eeg{d}{e}{t}.starttime + ds.endind / samprate;
+        ds.midtime = eeg{d}{e}{t}.starttime + ds.midind / samprate;
+        ds.peak = tmpevents(valids,3);
+        ds.energy = tmpevents(valids,7);
+        ds.maxthresh = (tmpevents(valids,9) - baseline) / stdev;
+        ds.threshind = tmpevents(valids,10);
+        ds.threshtime = eeg{d}{e}{t}.starttime + tmpevents(valids,10) / samprate;
+    else
+        ds.startind = [];
+        ds.endind = [];
+        ds.midind = [];
+        ds.starttime = [];
+        ds.endtime = [];
+        ds.midtime = [];
+        ds.peak = [];
+        ds.energy = [];
+    end
+    if any(find(renv<baseline))==0
+        warning(['No below baseline values in data.  Fields left blank, ',tmpflist(i).name])
+    end
+
+    ds.timerange = [0 length(renv)/samprate] + eeg{d}{e}{t}.starttime;
+    ds.samprate = eeg{d}{e}{t}.samprate;
+    ds.threshold = thresh;
+    ds.baseline = baseline;
+    ds.std = stdev;
+    ds.minimum_duration = min_suprathresh_duration;
+
+    dspikes{d}{e}{t} = ds;
+    clear eeg ds ripples;
+end
+save(sprintf('%s/%sdspikes%02d.mat', directoryname, fileprefix, d), 'dspikes');
+end

Dentate Spikes/dayprocess_dspikes.m

@@ -0,0 +1,22 @@
+E3 = {'r1_01_', 'r1_12_', 'r1_18_', 'r1_21_', 'r1_29_', 'r1_33_', 'r1_34_',  'r1_39_',...
+    'r2_06_', 'r2_17_', 'r2_28_',...
+    'r3_02_', 'r3_11_'};
+E4 = {'r1_05_', 'r1_24_', 'r1_31_', 'r1_35_', ...
+    'r2_08_', 'r2_13_', 'r2_14_', 'r2_26_', 'r2_27_', 'r2_36_', 'r2_38_',...
+    'r3_04_', 'r3_09_', 'r3_20_', 'r3_32_', 'r3_40_'};
+animals = [E3 E4];
+
+animalbasedir = 'E:\\LFP Data\CRCNS\Cohort 1\Preprocessed Data';
+days = 1;
+
+for a = animals
+    animaldir = sprintf('%s/%s/',animalbasedir,a{1});
+    prefix = a{1};
+
+    for d = days
+        ag_extractdspikes(animaldir, prefix, d, .015, 3)
+        plotdspikes
+        pause
+        close all
+    end
+end

Dentate Spikes/plotdspikes.m

@@ -0,0 +1,54 @@
+
+Fs=1000;  %sampling freq (/sec)
+s = d;
+e = 1;
+
+%identify ca1 pyr channel
+location = '(isequal($area,''dg'')&& (contains($layer,''**gc**'')||contains($layer,''**hil**'')))';  %search terms to search for 'pyr 1' (best pyr chan) in 'ca1' (to be searched inside 'chinfo')
+load(sprintf('%s/%schinfo.mat',animaldir,prefix)) %loads chinfo
+temp = evaluatefilter(chinfo,location); %Opens chinfo returns matches from 'location' in format [session, epoch, channel] (some repeated values because of chinfo)
+if(isempty(temp))
+    error('Chinfo must contain at least one gc channel');
+end
+ctarget = unique(temp(:,3)); %(':,3) looks specifically in col3 for the unique entries
+%ctarget = 10
+
+window = 500; %time before and after in ms (1ms=1sample); 1sec total  %500
+x_shift = 700; %spacing btwn traces for plot x-direction   %700
+y_shift=1500; %spacing btwn traces for plot y-direction
+
+load(sprintf('%s/%sdspikes%02d.mat', animaldir, prefix, s)) %loads ripples structure for session 's' (from extract_ripples.m contains information from all events classified as ripples, all chan)
+numrips = 10;  %how many ripples to print
+starts = dspikes{s}{e}{ctarget(1)}.startind(dspikes{s}{e}{ctarget(1)}.startind<900000);
+totalripnum = length(starts);
+rips = randi(totalripnum,1,numrips); %gives 1xnumrips random integers indicating selected ripple events
+
+for r = rips
+    figure %generate one
+    %%get ripple event times
+
+    startindex=dspikes{s}{e}{ctarget(1)}.midind(r)-window;
+    endindex = dspikes{s}{e}{ctarget(1)}.midind(r)+window;
+    timevec = startindex:(startindex+2*window); %build time vector, of the window before and after start index
+
+    for c=1:length(dspikes{s}{e}) %To figure out how many channels you have (some are 32, 31, etc)
+
+        file = sprintf('%s/%seeg%02d-%d-%02d.mat', animaldir, ...
+            prefix, s, e, c); %specifies raw data for the session, epoch, channel
+        load(file);
+
+        if c<=8 %first shank of plots
+            hold on
+            axis off
+            plot(timevec,eeg{s}{e}{c}.data(startindex:endindex)-c*y_shift)
+        elseif c>8 && c<=16 %second shank of plots
+            plot((timevec+window+x_shift),eeg{s}{e}{c}.data(startindex:endindex)-(c-8)*y_shift)
+        elseif c>16 && c<=24 %third shank of plots
+            plot((timevec+2*window+2*x_shift),eeg{s}{e}{c}.data(startindex:endindex)-(c-16)*y_shift)
+        elseif c>24 %fourth shank of plots
+            plot((timevec+3*window+3*x_shift),eeg{s}{e}{c}.data(startindex:endindex)-(c-24)*y_shift)
+        end
+    end
+
+    set(gcf, 'Position', [1 38 1920 487])
+end

MUA/RTMUA.m

@@ -0,0 +1,86 @@
+function out = RTMUA(index, excludeperiods, mua, ripples, chinfo, varargin)
+%parse the options
+win = [-0.4 0.4];
+binsize = 0.1;
+minstd = 5;
+interripwin = 0;  %does not eliminate ripples that occur in doublets or triplets
+
+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 'minstd'
+                minstd = 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;
+
+if ~isempty(trigindex)
+    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);
+end
+
+%iterate through probe channels
+for t = 1:size(trigindex,1)
+
+    %get ripple trigger times
+    rip = ripples{trigindex(t,1)}{trigindex(t,2)}{trigindex(t,3)};
+    [valid, ~] = getvalidrips(ripples,index,trigindex, excludeperiods, win, minstd, interripwin);
+    rip.starttime = rip.starttime(valid);
+    rip.endtime = rip.endtime(valid);
+    triggers = rip.threshtime(valid);
+    out.riptimes{t} = triggers;
+    out.ripsizes{t} = rip.maxthresh(valid);
+    out.ripstarts{t} = rip.starttime(valid);
+    out.ripends{t} = rip.endtime(valid);
+
+    %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.wave{t}{p} = NaN(length(triggers),40);
+        out.peakrate{t}{p} = NaN(length(triggers),1);
+        out.SWRrate{t}{p} = NaN(length(triggers),1);
+        out.normrate{t}{p} = NaN(length(triggers),1);
+        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
+            if sum(out.psth{t}{p}(r,:))>0 %if there are any spikes in this SWR
+                out.wave{t}{p}(r,:) = nanmean(pmua.waveform(spiketimes>triggers(r) & spiketimes<triggers(r)+binsize,:),1);
+            end
+            out.peakrate{t}{p}(r) = length(spiketimes(spiketimes>triggers(r) & spiketimes<triggers(r)+binsize))/binsize;
+            out.SWRrate{t}{p}(r) = length(spiketimes(spiketimes>rip.starttime(r) & spiketimes<rip.endtime(r)))/(rip.endtime(r) - rip.starttime(r));
+            out.normrate{t}{p}(r) = (out.SWRrate{t}{p}(r)-out.baserate{t}{p})/out.sdbaserate{t}{p}';
+        end
+    end
+end

MUA/dayprocess_MUA.m

@@ -0,0 +1,26 @@
+E3 = {'r1_01_', 'r1_12_', 'r1_18_', 'r1_21_', 'r1_29_', 'r1_33_', 'r1_34_',  'r1_39_',...
+    'r2_06_', 'r2_17_', 'r2_28_',...
+    'r3_02_', 'r3_11_'};
+E4 = {'r1_05_', 'r1_24_', 'r1_31_', 'r1_35_', ...
+    'r2_08_', 'r2_13_', 'r2_14_', 'r2_26_', 'r2_27_', 'r2_36_', 'r2_38_',...
+    'r3_04_', 'r3_09_', 'r3_20_', 'r3_32_', 'r3_40_'};
+animals = [E3 E4];
+
+datadir = 'E:\\LFP Data\CRCNS\Cohort 1\Trodes Raw Data';
+animalbasedir = 'E:\\LFP Data\CRCNS\Cohort 1\Preprocessed Data';
+
+for a = 1:length(animals)
+    prefix = animals{a};
+    animaldir = [animalbasedir '/' prefix];
+
+    taskname = sprintf('%s/%stask.mat',animaldir,prefix);
+    load(taskname);
+    for d = 1:length(task)
+        for e = 1:length(task{d})
+            if(~contains(lower(task{d}{e}.env), '*fail*') && ...
+                    ~contains(lower(task{d}{e}.descript), '*fail*'))
+                extract_spikes(datadir, animaldir, prefix, d, e);
+            end
+        end
+    end
+end