Closed AngCamp closed 10 months ago
Would something like this work?
import numpy as np
from scipy.signal import hilbert, butter, filtfilt
def detect_ripples(lfp_data, sampling_rate):
# Define the frequency range for ripple detection (e.g., 80-250 Hz)
ripple_band = [80, 250]
# Apply Hilbert transform to get the analytic signal
analytic_signal = hilbert(lfp_data)
# Compute the instantaneous amplitude
amplitude = np.abs(analytic_signal)
# Design a Butterworth filter for ripple bandpass filtering
nyquist = 0.5 * sampling_rate
low_cutoff = ripple_band[0] / nyquist
high_cutoff = ripple_band[1] / nyquist
b, a = butter(4, [low_cutoff, high_cutoff], btype='band')
# Apply the Butterworth filter to the amplitude
filtered_amplitude = filtfilt(b, a, amplitude)
return filtered_amplitude
# Load your LFP data (channels by samples)
lfp_data = data # put CA1 neuropixels data here as a channels by sample array
# Sampling rate (replace with your actual sampling rate)
sampling_rate = 2500 #the standard LFP frequency for neuropixels probes
# Detect ripples
filtered_amplitude = detect_ripples(lfp_data, sampling_rate)
Hi, To adapt the code to other cases is not trivial, are you using a Allen Institute neuropixel dataset? The heart of the detection is here https://github.com/RobertoDF/De-Filippo-et-al-2022/blob/7585ff1b1e3f9013e01d4d8774ecde28b197c98d/Calculations/Calculate_ripples.py#L51-L95
I am using Allen Brain data set and the IBL. What you called the UCL dataset.
I have also been using this library from the Eden Kramer Lab. I'm not sure there is much difference is the way they process the ripples. It may be identical to yours. I'm mainly interested in how you selected channels on the probes. You seemed to select the channel based on the highest ripple strength as well as skewness of the lfp?
"Ripples were detected on CA1 LFP traces, the best channel (higher ripple strength) was selected by looking at the SD of the envelope of the filtered trace, if multiple SD peaks were present across space (possibly caused by sharp waves in stratum radiatum and ripple activity in stratum pyramidale) we subsequently looked at the channel with higher skewness, in this way we could reliably identify the best ripple channel." So were you selecting one channel for each recording or was each channel selected for each ripple?
I selected one channel per area (one for CA1, one for CA3.....) in each session. I picked the channel with higher STD of the amplitude envelope https://github.com/RobertoDF/De-Filippo-et-al-2022/blob/7585ff1b1e3f9013e01d4d8774ecde28b197c98d/Utils/Utils.py#L766-L773 If there were multiple channels showing higher STD I looked at the skewness. Skewness could help picking channels with ripples and not the ones with LFP deflections not ripple-related Here is the code: https://github.com/RobertoDF/De-Filippo-et-al-2022/blob/7585ff1b1e3f9013e01d4d8774ecde28b197c98d/Utils/Process_lfp.py#L117-L141
Thanks that's very helpful, did you then manually determine by visual inspection which events were SWRs or did you trust your detector was working reasonably well and just accept there would be false positives? I know its standard to sort through events and rely on expert inspection but it seems like a ton of data to sort through visually.
In the process of developing the detection algorithm, I conducted extensive visual validations to ensure its accuracy. I wrote a little app to help with that https://github.com/RobertoDF/vizs/blob/main/lfp_explorer.gif. The thresholds for detection are quite high compared to what I saw in other papers so I´m quite confident the ripples detected are genuine. Given the size of the dataset I preferred to minimize false positive and accept some false negatives.
Hi I am trying to use your method to detect SWR in other data but these scripts seem tailor made to looping through the ABI visual coding dataset. I'm really just interested in the butterworth filter you've implemented and cleaning the lfp data. It seems straight forward but there are so many dependencies in the script it's hard to pull the function out for my use case. Could you explain what it's doing and how I can run that in numpy myself if I have the lfp as a 2d numpy array (channel by sample)? I have already multiplid the LFP data I'm using by it's scale factor and subtracted the median, I have identified the probes that are in CA1 but I want to try running your filter on it and it's hard because it calls so many other functions.
I apologize I am a statistician and not familiar with the intricacies of the electrical physics and signal processing going on and be able to pull out the relevant numpy functions you are using.