From 49184ec1de30e138c1c4075d35348dce0865949f Mon Sep 17 00:00:00 2001 From: Ryan Harvey Date: Sun, 5 Jan 2025 15:03:33 -0500 Subject: [PATCH] Create detect_spindles.m --- lfp/detect_spindles.m | 324 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 324 insertions(+) create mode 100644 lfp/detect_spindles.m diff --git a/lfp/detect_spindles.m b/lfp/detect_spindles.m new file mode 100644 index 0000000..8413501 --- /dev/null +++ b/lfp/detect_spindles.m @@ -0,0 +1,324 @@ +function detect_spindles(varargin) +% detect_spindles: Detect spindles in Local Field Potential (LFP) data. +% +% This function detects spindles in LFP data by preprocessing the signal, +% selecting appropriate channels based on spindle power, and running the +% detection algorithm. It uses multiple parameters to customize the process, +% such as the brain region to analyze, EMG thresholds, and whether to restrict +% detection to specific intervals (e.g., NREM sleep or immobility). +% +% Usage: +% detect_spindles('ParameterName', ParameterValue, ...) +% +% Parameters: +% 'basepath' (string, default: current directory) +% Path to the session directory containing the data. +% +% 'channel' (numeric array, default: []) +% Channel(s) to use for spindle detection (1-indexed). If empty, a +% channel is automatically selected based on spindle power. +% +% 'peak_threshold' (numeric, default: 5) +% Minimum z-score of the peak amplitude for spindle detection. +% +% 'tail_threshold' (numeric, default: 2.5) +% Minimum z-score for the start and end of a spindle. +% +% 'brainRegion' (string, char, or cell array of strings, default: {'PFC', 'MEC', 'LEC', 'EC', 'ILA', 'PL', 'Cortex', 'CTX'}) +% Specifies brain regions for channel selection. If no channel is +% provided, the function selects channels from these regions. +% +% 'manual_pick_channel' (logical, default: false) +% If true, opens Neuroscope2 for manual channel selection if no channel +% is provided. The user must tag the selected channel as 'spindle'. +% +% 'EMG_threshold' (numeric, default: 0.6) +% Threshold for EMG (electromyography) data. Periods with EMG below +% this value are considered immobility periods for spindle detection. +% Set to Inf to disable the immobility restriction. +% +% 'NREM_restrict' (logical, default: true) +% If true, restrict spindle detection to previously detected NREM +% sleep intervals. +% +% 'passband' (numeric array, default: [9, 17]) +% Frequency range (in Hz) used for filtering the signal when detecting +% spindles. +% +% 'ignore_intervals' (numeric array, default: []) +% Intervals to ignore during spindle detection (e.g., optogenetic +% stimulation periods). +% +% 'forceDetect' (logical, default: false) +% If true, forces spindle detection even if a detection file already +% exists. Overwrites outdated detection files. +% +% 'save_filename' (string, default: 'spindles') +% Filename for saving spindle detection results. +% +% Outputs: +% Results are saved to a file in the session directory in CellExplorer +% event format. The output file includes spindle start/stop times, peak +% amplitudes, durations, and metadata about the detection process. +% +% Example: +% detect_spindles('basepath', '/path/to/data', 'brainRegion', 'PFC'); +% +% Notes: +% - If no channel is provided, the function automatically selects a channel +% based on spindle power within the specified brain regions. +% - If 'manual_pick_channel' is true, Neuroscope2 is launched to allow the +% user to select a channel. +% - EMG and NREM restrictions help ensure that detected spindles occur during +% immobility or specific sleep states. +% +% Dependencies: +% - Neuroscope2 (optional for manual channel selection). +% - CellExplorer event file format for saving results. + + +% parse inputs +p = inputParser; +p.PartialMatching = true; % allows the use of abbreviations (e.g. "EMG_thres" instead of "EMG_threshold" +addParameter(p, 'basepath', pwd, @(x) any([isfolder(x), iscell(x)])); +addParameter(p, 'channel', [], @isnumeric); +addParameter(p, 'peak_threshold', 5, @isnumeric); +addParameter(p, 'tail_threshold', 2.5, @isnumeric); +addParameter(p, 'brainRegion', {'PFC', 'MEC', 'LEC', 'EC', 'ILA', 'PL', 'Cortex', 'CTX'}, ... + @(x) ischar(x) || isstring(x) || (iscell(x) && all(cellfun(@(y) ischar(y) || isstring(y), x)))); +addParameter(p, 'manual_pick_channel', false, @islogical); +addParameter(p, 'EMG_threshold', 0.6, @isnumeric); +addParameter(p, 'NREM_restrict', true, @islogical); +addParameter(p, 'passband', [9, 17], @isnumeric); +addParameter(p, 'ignore_intervals', [], @isnumeric); +addParameter(p, 'forceDetect', false, @islogical); +addParameter(p, 'save_filename', 'spindles', @ischar); + +parse(p, varargin{:}) +basepath = p.Results.basepath; + +% if single basepath, package in cell array +if ~iscell(basepath) + basepath = {basepath}; +end +% iterate over basepaths +for i = 1:length(basepath) + run(basepath{i}, p) +end +end + +function run(basepath, p) +% local function to carry out detection + +brainRegion = p.Results.brainRegion; +EMG_threshold = p.Results.EMG_threshold; +peak_threshold = p.Results.peak_threshold; +tail_threshold = p.Results.tail_threshold; +passband = p.Results.passband; +channel = p.Results.channel; +NREM_restrict = p.Results.NREM_restrict; +manual_pick_channel = p.Results.manual_pick_channel; +ignore_intervals = p.Results.ignore_intervals; +forceDetect = p.Results.forceDetect; +save_filename = p.Results.save_filename; + +disp(basepath) + +basename = basenameFromBasepath(basepath); + +% check if file aready exists +event_file = fullfile(basepath, [basename, '.', save_filename, '.events.mat']); +if exist(event_file, "file") && ~forceDetect + % verify file contains minimum fields (timestamps, peaks) + load(event_file, 'spindle') + if isfield(spindle, 'timestamps') && isfield(spindle, 'peaks') + return + else + warning('existing file is out of date, and will be overwritten') + end +end + +% load session meta data +load(fullfile(basepath, [basename, '.session.mat']), 'session') + +if NREM_restrict + load(fullfile(basepath, [basename, '.SleepState.states.mat']), 'SleepState') + lfp_restrict = SleepState.ints.NREMstate; + if isempty(lfp_restrict) + warning("no NREM sleep, skipping session") + return + end +else + lfp_restrict = [-inf, inf]; +end + +% checks for manually added spindle and if doesn't exist will open neuroscope +if manual_pick_channel + if isfield(session.channelTags, 'spindle') + channel = session.channelTags.spindle.channels; + warning('using spindle channel that was already tagged') + else + disp("Pick spindle channel from neuroscope") + disp("add channel tag as 'spindle'") + NeuroScope2('basepath', basepath) + % load session meta data + load(fullfile(basepath, [basename, '.session.mat']), 'session') + if isfield(session.channelTags, 'spindle') + channel = session.channelTags.spindle.channels; + else + warning("you did not pick a channel, or did not name it 'spindle'") + end + end +end + + +if isempty(channel) || length(channel) > 1 + + if isempty(channel) + % find nearest region name + regions = fields(session.brainRegions); + region_idx = contains(regions, brainRegion); + if isempty(region_idx) + warning("no channels for brainRegion") + return + end + region_name = regions(region_idx); + + % find channels for this brain region + channels = []; + for region_i = 1:length(region_name) + channels = [channels, session.brainRegions.(region_name{region_i}).channels]; + end + + % find and remove bad channels + channels = channels(~ismember(channels, session.channelTags.Bad.channels)); + + if isempty(channels) + warning("no channels for brainRegion") + return + end + else + channels = channel; + end + + + [lfp, timestamps] = load_lfp(basepath, session, channels, lfp_restrict); + + % get spindle power to choose channel + try + pBand = bandpower(double(lfp), ... + session.extracellular.srLfp, passband); + + pTot = bandpower(double(lfp), ... + session.extracellular.srLfp, ... + [1, (session.extracellular.srLfp / 2) - 1]); + catch + for c = 1:size(lfp, 2) + pBand(c) = bandpower(double(lfp(:, c)), ... + session.extracellular.srLfp, passband); + + pTot(c) = bandpower(double(lfp(:, c)), ... + session.extracellular.srLfp, ... + [1, (session.extracellular.srLfp / 2) - 1]); + end + end + % find max theta power, normalized by wide band + [~, c_idx] = max(pBand./pTot); + + % only leave channel and load lfp + channel = channels(:, c_idx); + lfp = lfp(:, c_idx); +else + [lfp, timestamps] = load_lfp(basepath, session, channel, lfp_restrict); +end + +% remove artifacts from lfp +try + [clean_lfp, ~, ~] = CleanLFP([timestamps, double(lfp)], ... + 'thresholds', [6, 10], 'manual', false); +catch + % if CleanLFP fails, use original... CleanLFP needs testing + clean_lfp = [timestamps, double(lfp)]; +end +% remove lfp struct because we now have clean lfp +clear lfp + +if EMG_threshold < Inf + % keep waves during low emg + EMG = getStruct(basepath, 'EMG'); + immobility = EMG.timestamps(FindInterval(EMG.data < EMG_threshold)); + immobility(diff(immobility, [], 2) < 1, :) = []; + % reshape interval if needed + if size(immobility, 2) == 1 + immobility = immobility'; + end +else + immobility = [-inf, inf]; +end + +% filter in spindle band +filtered = FilterLFP(clean_lfp, 'passband', passband); + +% detect spindles, Nx4 matrix [start peak_t end peak_z] +spindles_ = FindSpindles(Restrict(filtered, immobility), ... + 'threshold', tail_threshold, 'peak', peak_threshold); + +% put into table (like pandas but much worse) +spindles = array2table(spindles_, 'VariableNames', {'start', 'peak_t', 'stop', 'peak_z'}); + +% remove spindles that intersect with ignore intervals +if ~isempty(ignore_intervals) + keep_intervals = ~IntervalsIntersect(spindles{:, {'start', 'stop'}}, ignore_intervals); + spindles = spindles(keep_intervals, :); +end + +% store to cell explorer event file +spindles_ = struct(); +spindles_.timestamps = spindles{:, {'start', 'stop'}}; +spindles_.peaks = spindles.peak_t; +spindles_.amplitude = spindles.peak_z; +spindles_.amplitudeUnits = 'zscore'; +spindles_.eventID = []; +spindles_.eventIDlabels = []; +spindles_.eventIDbinary = false; +spindles_.center = median(spindles_.timestamps, 2); +spindles_.duration = spindles_.timestamps(:, 2) - spindles_.timestamps(:, 1); +spindles_.detectorinfo.detectorname = 'detect_spindles.m'; +spindles_.detectorinfo.detectiondate = datetime("today"); +spindles_.detectorinfo.detectionintervals = [clean_lfp(1, 1), clean_lfp(end, 1)]; +spindles_.detectorinfo.detectionparms = p.Results; +try + spindles_.detectorinfo.detectionchannel1 = channel; + spindles_.detectorinfo.detectionchannel = channel - 1; +catch + +end +spindles_.detectorinfo.threshold = [tail_threshold, peak_threshold]; + +spindles = spindles_; +save(event_file, 'spindles'); +end + +function [lfp, restricted_timestamps] = load_lfp(basepath, session, channels, restrict_intervals) + +fileName = fullfile(basepath, [basenameFromBasepath(basepath), '.lfp']); +if ~exist(fileName, 'file') + fileName = fullfile(basepath, [basenameFromBasepath(basepath), '.eeg']); +end +nChannels = session.extracellular.nChannels; + +filenamestruct = dir(fileName); +% determine number of bytes per sample +dataTypeNBytes = numel(typecast(cast(0, 'int16'), 'uint8')); +% Number of samples per channel +nSamp = filenamestruct.bytes / (nChannels * dataTypeNBytes); +% memory map file +mmf = memmapfile(fileName, 'Format', {'int16', [nChannels, nSamp], 'data'}, ... + 'Writable', false); + +timestamps = (1:nSamp)' / session.extracellular.srLfp; + +[restricted_timestamps, idx] = Restrict(timestamps, restrict_intervals); + +lfp = mmf.Data.data(channels, idx)'; +end