Content of the download_raw_ephys_snippets.py script
# %%# This script demonstrates raw Neuropixels AP-band data processing:# 1. Raw data — straight off the probe, unfiltered# 2. Butterworth — high-pass filter to remove LFP/DC drift# 3. Destriped — full IBL pre-processing pipeline (butter + ADC correction + k-filter)## viewephys opens an interactive viewer for each processing stage so you can compare them side by side.import scipy.signalimport ibldsp.voltagefrom one.api import ONEfrom brainbox.io.one import SpikeSortingLoaderfrom viewephys.gui import viewephys# ONE = Open Neurophysiology Environment, the IBL data access layer.# The public endpoint lets anyone download data with those public credentials.one_kwargs =dict( base_url='https://openalyx.internationalbrainlab.org', username='intbrainlab', password='international', silent=True,)one = ONE(**one_kwargs)# Probe insertion IDs (pids) from the IBL Benchmarks dataset — a curated set of high-quality recordings.pids = ['1a276285-8b0e-4cc9-9f0a-a3a002978724', # 00 - Benchmark PIDS start'1e104bf4-7a24-4624-a5b2-c2c8289c0de7','6638cfb3-3831-4fc2-9327-194b76cf22e1','749cb2b7-e57e-4453-a794-f6230e4d0226','d7ec0892-0a6c-4f4f-9d8f-72083692af5c','da8dfec1-d265-44e8-84ce-6ae9c109b8bd','dab512bd-a02d-4c1f-8dbc-9155a163efc0','dc7e9403-19f7-409f-9240-05ee57cb7aea','e8f9fba4-d151-4b00-bee7-447f0f3e752c','eebcaf65-7fa4-4118-869d-a084e84530e2','fe380793-8035-414e-b000-09bfe5ece92a', # Benchmark PIDS stop]T0 =500# start time in seconds — pick a quiet moment well into the recordingfor pid in pids:# SpikeSortingLoader is a convenience object that bundles the probe geometry,# spike-sorting results, and raw data access for a single insertion. ss = SpikeSortingLoader(pid=pid, one=one)# get the histology information channels = ss.load_channels()# stream=True fetches only the requested chunk on-the-fly — no need to download the full recording.# band='ap' selects the action-potential band (300 Hz high-pass, 30 kHz sample rate). sr = ss.raw_electrophysiology(band='ap', stream=True)# sr[samples, channels] — slice to get 0.5 s of data.# We transpose (.T) so the result is (n_channels, n_samples): the convention expected by ibldsp.# The last channel is the sync/TTL channel; dropping it keeps only the 384 neural channels. raw = sr[int(T0 * sr.fs):int((T0 +0.5) * sr.fs), :-1].T
5.1 Display snippets of raw data at different stages of pre-processing
Run the following code in a interactive Python session to visualize raw data at different stages of pre-processing. If using Ipython, use the magic command %gui qt
# %%# This script demonstrates raw Neuropixels AP-band data processing:# 1. Raw data — straight off the probe, unfiltered# 2. Butterworth — high-pass filter to remove LFP/DC drift# 3. Destriped — full IBL pre-processing pipeline (butter + ADC correction + k-filter)## viewephys opens an interactive viewer for each processing stage so you can compare them side by side.import scipy.signalimport ibldsp.voltagefrom one.api import ONEfrom brainbox.io.one import SpikeSortingLoaderfrom viewephys.gui import viewephysfrom iblatlas.atlas import BrainRegions# ONE = Open Neurophysiology Environment, the IBL data access layer.# The public endpoint lets anyone download data with those public credentials.one_kwargs =dict( base_url='https://openalyx.internationalbrainlab.org', username='intbrainlab', password='international', silent=True,)one = ONE(**one_kwargs)# Probe insertion IDs (pids) from the IBL Benchmarks dataset — a curated set of high-quality recordings.pids = ['1a276285-8b0e-4cc9-9f0a-a3a002978724', # 00 - Benchmark PIDS start'1e104bf4-7a24-4624-a5b2-c2c8289c0de7','6638cfb3-3831-4fc2-9327-194b76cf22e1','749cb2b7-e57e-4453-a794-f6230e4d0226','d7ec0892-0a6c-4f4f-9d8f-72083692af5c','da8dfec1-d265-44e8-84ce-6ae9c109b8bd','dab512bd-a02d-4c1f-8dbc-9155a163efc0','dc7e9403-19f7-409f-9240-05ee57cb7aea','e8f9fba4-d151-4b00-bee7-447f0f3e752c','eebcaf65-7fa4-4118-869d-a084e84530e2','fe380793-8035-414e-b000-09bfe5ece92a', # Benchmark PIDS stop]pid = pids[-1]# SpikeSortingLoader is a convenience object that bundles the probe geometry,# spike-sorting results, and raw data access for a single insertion.ss = SpikeSortingLoader(pid=pid, one=one)# get the histology informationchannels = ss.load_channels()# stream=True fetches only the requested chunk on-the-fly — no need to download the full recording.# band='ap' selects the action-potential band (300 Hz high-pass, 30 kHz sample rate).sr = ss.raw_electrophysiology(band='ap', stream=True)sr_lf = ss.raw_electrophysiology(band='lf', stream=True)# %%# ---------------------------------------------------------------------------# Load a short snippet of raw voltage data# ---------------------------------------------------------------------------T0 =500# start time in seconds — pick a quiet moment well into the recording# sr[samples, channels] — slice to get 0.5 s of data.# We transpose (.T) so the result is (n_channels, n_samples): the convention expected by ibldsp.# The last channel is the sync/TTL channel; dropping it keeps only the 384 neural channels.raw = sr[int(T0 * sr.fs):int((T0 +0.5) * sr.fs), :-1].T# %%# ---------------------------------------------------------------------------# Step 1 — Butterworth high-pass filter# ---------------------------------------------------------------------------# A 3rd-order Butterworth filter at 300 Hz removes slow LFP oscillations and# DC drift while preserving the sharp transients caused by action potentials.# sosfiltfilt applies the filter twice (forward + backward) for zero phase distortion.butter_kwargs = {'N': 3, 'Wn': 300, 'btype': 'highpass', 'fs': sr.fs}sos = scipy.signal.butter(**butter_kwargs, output='sos')butt = scipy.signal.sosfiltfilt(sos, raw)# ---------------------------------------------------------------------------# Step 2 — Full IBL destriping pipeline# ---------------------------------------------------------------------------# ibldsp.voltage.destripe runs the complete IBL pre-processing chain:# • High-pass Butterworth (same as above)# • ADC sample-shift correction (Neuropixels channels are not sampled simultaneously)# • Spatial k-filter (removes stripe artefacts correlated across channels)## channel_labels=True lets the function detect bad channels automatically from this# snippet — convenient for exploration, but in production you would pre-compute# labels from a longer stretch of data for robustness.destriped = ibldsp.voltage.destripe(raw, fs=sr.fs, h=sr.geometry, neuropixel_version=sr.major_version, channel_labels=channels['labels'])# ---------------------------------------------------------------------------# Display all three stages in viewephys for side-by-side comparison# ---------------------------------------------------------------------------# viewephys is an IBL interactive viewer: scroll with the mouse wheel,# adjust gain with +/-, and pan by clicking and dragging.regions = BrainRegions()eqcs = {}eqcs['raw'] = viewephys(raw, sr.fs, channels=channels, title='Raw (unfiltered)', t0=T0, br=regions)eqcs['butt'] = viewephys(butt, sr.fs, channels=channels, title='Butterworth high-pass 300 Hz', t0=T0, br=regions)eqcs['destriped'] = viewephys(destriped, sr.fs, channels=channels, title='Destriped (butter + ADC + k-filter)', t0=T0, br=regions)
The viewephys viewer showing the destriped raw data snippet
5.1.1 Navigating the viewer
Pan: left-click and drag
Zoom: scroll the mouse wheel (zooms on the cursor position), or right-click and drag
Pan with the keyboard: arrow keys (hold Ctrl for a bigger step)
Gain: Page Up / Ctrl+A decreases the gain (amplifies the signal), Page Down / Ctrl+Z increases it (you can also type a value directly in the gain box)
Propagate: Ctrl+P applies the current gain, colormap, sort order and axis links to all other open viewer windows; Shift+P does the same and additionally tiles the windows side by side
Copy the current view: Ctrl+S copies the plot to the clipboard as an image
Right-click on the plot for a context menu (view a single trace, spectrum or spectrogram)
This is handy for comparing the raw, Butterworth-filtered and destriped windows side by side: open the three viewers, arrange/zoom one of them as needed, then hit Ctrl+P (or Shift+P) in that window to propagate the view to the other two.
5.1.2 Cycling through the different pre-processing steps
Once you did Propagate: Ctrl+P, all of the 3 windows should now be synchronized. - Cycle through using Alt+Tab. How different are they ? - On the raw display, right-click on the data, and look at individual traces. What can be said about the DC offset and grounding ? The units are uV, and those probes saturate at 1200 uV, can you spot any saturation events on the raw data ?
5.1.3 Looking at different recordings in different brain regions
You can select different insertions (probe ids), by setting the pid = pids[i] index to visualize different brain regions.
In a shell, with the environment active type in viewephys. Select raw, and LFP, then open the File > Open menu and point to the .cbin file in the folder above.
How close is the raw voltage from saturation ? What is the shape of the spectrum if you right click and compute the PSD ?
6 Preprocessing (coding)
In this chapter we will explore in more detail the preprocessing steps used in extracellular electrophysiology.
The purpose of preprocessing is to isolate as best possible our signal of interest (action potentials) from background noise. In electrophysiology the noise sources are typically electrical (e.g. broken recording electrodes) or mechanical (e.g. physical movement of the probe in the brain).
This chapter will focus on the practical implementation of preprocessing steps using SpikeInterface. See the previous chapter for details on the underlying theory.
The preprocessing steps we will explore in this chapter are:
saturation detection: Detect periods in the recording where the electronics were saturated (i.e. the voltage exceeds the range of the recording electronics). We detect saturation periods on the raw recording, and set any detected periods to zero after filtering.
bandpass filtering: Filters out low-frequency oscillations (e.g. 50-60 Hz noise) and high-frequency oscillations (e.g. noise spikes) from the data. Applied separately per channel.
bad channel detection : Automatically detect channels that are dead, noisy or out-of-brain. These bad channels can be removed or interpolated (i.e ‘filled in’ using the data from surrounding channels).
common median referencing: Removes noise that takes the form of stripes through the data (i.e. occurring on all channels for a short period of time). This is applied across all channels, separately per time point.
motion correction : A step to account for physical drift of the probe in the brain. Attempts to estimate the drift and interpolate the channels (having the effect of shifting the signal up or down the probe) to account for it.
whitening : Removes correlated noise across channels and helps isolate spikes on a subset of channels. Has significant effects on the data scaling and waveform shape.
Warning
We are applying all of these preprocessing steps to explore their effects — In general you would not apply all of these steps to your recordings. The preprocessing steps you choose will depend on your data, e.g. your recording device, whether your implant is chronic or acute and the sorter you use. We will discuss this throughout the chapter.
6.1 Set up
We will preprocess a new dataset recorded on a 64-channel Cambridge Neurotech probe. Note that we will not apply the phase_shift step that was used in introduction preprocessing section, as that is only relevant for Neuropixel probes.
Before starting, we will import all required packages for the section:
from pathlib import Pathimport numpy as npimport matplotlib.pyplot as pltimport spikeinterface.core as siimport spikeinterface.preprocessing as spreimport spikeinterface.widgets as swfrom course_ephys_osss.data.utils import download_and_slice_nwb_from_dandi
6.1.1 Downloading the data
This is the first chapter in which we will use a core dataset used in this course, a postsubiculum recording on a 64-channel Cambridge Neurotech probe used in an experiment on the neural tuning of head direction cells (see here for details).
We will download the first 90 seconds of a recording from this dataset from the open access DANDI archive repository (an excellent resource)
# (optionally) change this to where you would like to download the data to.RAW_DATA_FOLDER = Path("data/raw-dandiset-000939")if RAW_DATA_FOLDER.exists():# if we have already downloaded the data, reload it. recording = si.load(RAW_DATA_FOLDER)else:# Otherwise, download the first 90 seconds from DANDI download_and_slice_nwb_from_dandi( s3_url="https://dandiarchive.s3.amazonaws.com/blobs/a8f/800/a8f8003e-4483-4b50-8a45-91ac5971f5d5", duration_s=90, local_folder=RAW_DATA_FOLDER ) recording = si.load(RAW_DATA_FOLDER)
Saved https://dandiarchive.s3.amazonaws.com/blobs/a8f/800/a8f8003e-4483-4b50-8a45-91ac5971f5d5 to data/raw-dandiset-000939!
Alternatively, the data can be downloaded directly from the course data Zenodo, with folder name dandiset-000939-90s.zip (will need unzipping before use).
6.1.2 Creating a plotting function
Next, we will define a function to plot and compare multiple preprocessing steps. We will reuse this function through the rest of the tutorial.
def plot_prepro(recordings, time_range, return_in_uV=True, clim=None):""" Plot multiple recordings in subplots. Parameters ---------- recordings : A dictionary of recordings. The recordings will be plot and the corresponding keys used as plot titles. time_range: A tuple with the (start, stop) times of the data to plot. clim: A tuple with the (min, max) colour limits for the heatmap. """ fig, axes = plt.subplots(1, len(recordings), figsize=(6*len(recordings), 5), squeeze=False, layout="constrained", )# Grab the first (and only) row of axes axes = axes[0, :]for idx, (name, recording) inenumerate(recordings.items()): ax = axes[idx] sw.plot_traces( recording, time_range=time_range, mode="map", # plot as a heatmap ax=ax, with_colorbar=True, return_in_uV=True, clim=clim, ) ax.set_title(name) ax.tick_params(axis="x", rotation=45)
6.2 Detecting saturation periods
The electronics that record the extracellular voltage have dynamic ranges in which they operate. If the voltage exceeds this range, the equipment is said to ‘saturate’ and the recorded values are uninformative.
The working ranges are associated with the specific recording equipment, most commonly the amplifier and analog-digital converters. For Neuropixels, these electronics are contained in the probe and the probe manual specifies the input voltage range (IBL recommend a threshold of 1200 μV for NP1; NP2 probes have a wider range and are harder to saturate).
The data we are working with was recorded on a Cambridge Neurotech probe, using an Intan RHD2000-series amplifier/headstage with range -5 mV to 5 mV. Therefore, we can use SpikeInterface to detect any periods of the recording that fall outside ±5 mV (5000 μV).
We need to detect saturation events before any preprocessing has occurred, as we need the raw voltage values to determine saturation. However, we must remove these events by zeroing them out after filtering, as sharp transitions to zero can cause filter artifacts. Therefore, we will detect the saturation periods now but actually zero out the data later.
We can see that the detected saturation event indeed falls outside the ±5 mV range. The data recorded here is thus uninformative and we will zero it out later.
First, let’s have a look at the effect of bandpass filtering the data, using the bandpass_filter() function. We will plot the first second of the data:
recording_bp = spre.bandpass_filter( recording, # commonly used values to filter noise # while not affecting spikes freq_min=300, freq_max=6000) plot_prepro( {"raw": recording, "filtered": recording_bp}, time_range=(0, 1))
The filtering operation removes a lot of contaminating low-frequency noise, realizing a significant improvement to the data. Electrical noise at 50 or 60 Hz (i.e. oscillating 50 or 60 times per second) is common as mains power electricity oscillates at those frequencies in many countries.
Bandpass filtering helps suppress this slow periodic component (although in practice good grounding and shielding are also important for reducing mains noise at the source).
Plotting a shorter time window demonstrates how this improves resolution of individual spikes:
Individual recording electrodes on a probe sometimes break and become ‘dead’ (no signal), ‘noise’ (excessive electrical noise) or may be ‘out of brain’ (i.e. a contiguous set of channels at the end of the probe with no signal, indicating they are not in brain tissue). In our case, there is a small number of channels at the top of the previous plot whose traces look different than all the rest of the signal. These are likely out of the brain.
First, let’s use the detect_bad_channels() function on our data, printing a list of the channel label (“good”, “out”, “dead”, “noise”), removing the bad channels and visualizing the traces:
It didn’t work! We can still see the dead channels on our recording. Why is this?
The IBL algorithm was developed for Neuropixel probes (with 384 recording channels) but we are using Cambridge Neurotech 64-channel probes. As such, the default settings are not appropriate for our data, and we need to adjust them.
The new values were found through trial and error. n_neighbors default is 11 so we reduced it for the 64-channel Cambridge Neurotech probe. We also set outside_channels_location=“bottom”. By convention the algorithm expects out-of-brain channels at the top of the depth-ordered channels (as on Neuropixels), but by this probe’s channel-ordering convention they sit at the bottom (index 0). Finally, we reduced the thresholds through trial and error to ensure the dead and out-of-brain channels were flagged correctly (once set, this configuration should work for all recordings in this dataset, but of course we should always check).
# Let's try again with different settingsbad_channel_ids, channel_labels = spre.detect_bad_channels( recording_bp, dead_channel_threshold=-0.15, outside_channels_location="bottom", outside_channel_threshold=-0.25, n_neighbors=5)print("second attempt labels:\n", channel_labels)rec_clean = recording_bp.remove_channels(bad_channel_ids)plot_prepro( {"filtered": recording_bp, "bad channels removed": rec_clean}, TIME_RANGE)
Now the dead channels, identifiable by visual inspection, are properly detected by the algorithm.
Warning
This example demonstrates why it is critically important to always visualize your data to ensure the default arguments of the algorithms we are using are suitable for the recording.
6.5 Common average referencing
Electrophysiological recordings are often contaminated by noise artifacts that appear on all channels at once for a short period of time. Striping artifacts are typically caused by electrical noise that is shared across many recording channels, for example fluctuations in the reference electrode or electromagnetic interference. These appear as vertical ‘stripes’ on the heatmap view of the data, as can be seen below.
Common average referencing (CAR) takes the average signal across all of the channels, and subtracts it. Sometimes you will see this called ‘common median referencing’ (CMR) when the median is used instead of the mean.
We can use SpikeInterface’s common_reference() function to perform common median referencing:
We can see that contaminating noise spikes that run vertically across the heatmap are removed through this processing method.
The implicit assumption of CAR or CMR is that the variability in the noise stripe is fairly small across channels (and so can be approximated well by removing only the central tendency measure). However, in practice this assumption is often broken, with stripes appearing only on a subset of channels (e.g. the top half). As such, the common_reference function contains a number of arguments allowing the operation to be performed on subsets of channels, and these can be explored in the ‘things to try’ section below.
6.6 Zeroing the saturation artifacts
We do not want to zero the saturation periods before filtering because the sharp change will create artifacts during filtering. Therefore, we apply the zeroing after the filtering steps:
There’s no need to get through all of these — pick whatever looks interesting.
The time window we plot is important. Each spike is around 1-5 ms long, so if we plot long time windows we don’t have enough pixels on the screen to resolve spikes. Nonetheless, longer windows can reveal other patterns in the data (e.g. 50 Hz noise). Play around with plotting different time windows (e.g. 10 ms, 200 ms, 1 s) to get a feel for this.
Plot the data using the mode="line" argument. Same data — a completely different perspective!
What happens if you set the bandpass filter minimum cutoff to zero? What frequencies are included now? Or set the maximum cutoff to 15000?
Try common median referencing (CMR) particular subsets of channels, using the reference, groups arguments as well as use different channels as reference using the ref_channel_ids argument. Recordings have a get_channel_ids() function that returns a list of ids.
Compute the Fourier transform of a data before and after bandpass filtering (using the .get_traces() function and scipy’s fourier transform functions) and plot the frequency spectrum.
6.7 Motion correction
The recording probe can move during an experiment. The major axis of movement is typically the direction along the probe’s insertion path, as this offers less physical resistance than moving into the surrounding brain. Movement along this axis has the effect of shifting the recorded waveforms up or down the probe.
The majority of motion correction algorithms are based around the movement of detected spike waveforms. The pattern of detected spikes is used as a ‘barcode’ across time, allowing motion to be tracked through the recording (see this infographic for details - the grey dots in the plots are individual spikes which have been detected). Most motion correction algorithms have three main steps:
Spike detection and localization: This step detects spikes in the data and estimates their physical location in space (x, y, z coordinates) from the waveform.
Motion estimation: This estimates the probe movement based on how the detected spikes move up or down along the probe during the recording. Estimation may be rigid (capturing movement of the entire probe) or non-rigid (capturing movement that differs along the probe length, such as that caused by movement of brain tissue).
Interpolation: Given the motion trace and recorded data, this interpolates the data (i.e. estimates the true data at a point given the collected data at surrounding points) back to what it would have been if there was no motion.
The key drift correction algorithms out there include Kilosort’s drift correction (kilosort-like in SpikeInterface), DREDge and MEDiCINe. SpikeInterface exposes these through the correct_motion() function that uses convenient presets which define how spike detection, motion estimation and interpolation is run.
Today we will use dredge_fast and set output_motion_info=True so we can inspect the results of motion correction.
Typically, non-rigid estimation algorithms break up the probe into sections and perform rigid alignment within each. Because this Cambridge Neurotech probe has only 64 channels, non-rigid estimation does not perform as well, and so we will turn it off (a full list of these arguments can be found here):
An InterpolateMotionRecording is returned. This recording object applies motion correction to the requested chunk of data based on the estimated motion.
Once motion correction is run, we can plot the outputs to inspect the results:
In the top two plots, we see a spike raster plot before (left) and after (right) motion correction. Each dot indicates a spike, with its depth on the y-axis and recording time on the x-axis.
The bottom left plot shows the estimated physical motion in μm over time.
Because this is only a 90-second recording, we see very little drift (note: the y-spacing between the electrodes is around 10 microns, while the drift is less than 1 micron.). On longer recordings, the drift and its correction will be more apparent, for example on this example recording with injected drift (from this useful explainer):
Figure 6.1: Example DREDGE drift-correction on a longer recording with injected drift.
In our case, if we compute the drift across the entire recording we get the below plot:
Figure 6.2
We see that there is little drift in the recording, and that motion correction has in fact introduced an artefact around the 4000 second mark visible in the top-right plot.
This highlights an important point: if the drift plots for your recordings indicate little drift, you may want to consider turning off motion correction. This is because the interpolation step will be performed based on small, but non-zero estimated motion and can add noise to the recording.
Similarly, erroneous motion correction can result in a vertical ‘smearing’ of spikes in the raster plot. This is especially true for low-channel count probes or probes with low firing rate neurons (see here for a discussion on the phenomenon) e.g:
Figure 6.3: Erroneous motion estimation from the kilosort_like preset on a low-channel-count probe. A large error, around 4700 seconds, can be seen in the motion vector which ‘smears’ the estimated peak locations in the corrected raster plot (top right).
The takeaway from this section is to always plot your estimation motion (even if you don’t eventually apply motion correction), using it to make the decision whether to apply this to your sorted data or not. You can check drift maps in SpikeInterface, through Kilosort outputs or using the interactive driftplots package.
6.8 Whitening
Finally, we will explore the whitening step. The whitening step reduces noise correlated across channels, and can help localize spikes onto a smaller number of channels.
Whitening is quite an aggressive preprocessing step, significantly rescaling the data and changing spike waveforms, and is sensitive to how the covariance between channels is estimated.
In SpikeInterface, whitening is done by the whiten() function:
This is all we need to do to get a whitened version of our recording, however due to the scaling that whitening performs it is difficult to see what has happened.
To get a better understanding of how whitening changes the data, let’s rescale the whitened data back to its original scaling (roughly) for visualization purposes. Note we would not ever do this for recordings used in downstream analysis, this is just for visualization purposes.
NoteHow the rescaling works.
The correlation matrix is the key structure of the whitening step. It is computed by correlating all channels with each other:
Figure 6.4
Here \(c_i\) represents a channel’s data indexed by some integer \(i\). \(c_ic_j\) represents the sum of squares of the data from channels \(i\) and \(j\) (i.e. multiply each value in the vector with the corresponding value and add these all together). For example \(c_1c_1\) is the variance for channel \(c_1\), and \(c_1c_2\) represents the covariance between channels 1 and 2.
For those interested, you can see how this comes from the covariance definition when the mean of the channel data is zero and we ignore the scaling factor \(N\).
The goal of the whitening step is to remove all cross-channel correlation from the data. As we can see from the above matrix, this equates to setting all the off-diagonal entries to zero. In fact it goes a little further, setting the variance of each channel to 1 as well (making the covariance matrix the ‘identity’ matrix):
Figure 6.5
Rescaling the variances is why this step has such a large effect on the data scaling.
We can estimate the original channel scaling by multiplying the data by the average standard deviation (i.e. noise) of the data. When we compute the covariance matrix of the scaled data, we should see this approximation of the original channel variances along the diagonal, instead of ones (with zeros still on the off-diagonals, at least in theory).
We can zoom in on some spikes, to really see how the waveforms are changed (below, we extract the raw data and slice it to visualize some representative spikes that display the phenomenon clearly):
# Select the times we want to display and convert them into samplesstart_frame =int(TIME_RANGE[0] * rec_corrected.get_sampling_frequency())end_frame =int(TIME_RANGE[1] * rec_corrected.get_sampling_frequency())# Get the data from this time range, and further slice it around a few spikesorig_data = rec_corrected.get_traces( start_frame=start_frame, end_frame=end_frame, return_in_uV=True)[200:400, :4]whitened_data = scaled_white_recording.get_traces( start_frame=start_frame, end_frame=end_frame, return_in_uV=True,)[200:400, :4]# Plot the spikes, adding an offset between traces so they # are displayed vertically stacked.plt.figure(figsize=(8, 5))plt.plot(orig_data + np.arange(orig_data.shape[1]) *750, color="k")plt.plot(whitened_data + np.arange(whitened_data.shape[1]) *750, color="r")plt.xlabel("sample number")plt.show()
Here, we can see that spikes are ‘squashed’ onto a smaller number of channels (which can aid spike sorting). The second spike is now only found on 4 channels instead of 3.
However, the 1st and 3rd spikes, which was likely a single waveform distributed over 4 channels, has now been split into two separate spikes. This could lead to a single unit being split into two during sorting.
6.9 Saving the data
SpikeInterface makes it very easy to save preprocessed data, which is useful as applying preprocessing to an entire long recording can take some time. Here we will save the data preprocessed with bandpass filter, bad channel removal and common referencing for use later in the course.
In general, we can perform these lighter preprocessing steps in SpikeInterface and let the sorter use its internal preprocessing for other steps.
Finally, when we sort our preprocessed data we must remember to turn off the relevant steps in Kilosort so we do not apply them twice! (for filtering, it is not possible to fully turn it off so we set the cutoff to a very low value so that the step does little to the data).
For Kilosort this is an example to make this point clear. In practice, most people who run Kilosort use Kilosort’s own whitening and motion correction.
Warning
Take care when turning off preprocessing steps that a sorter applies. For example, Kilosort runs whitening prior to motion correction, while in SpikeInterface it is recommended to perform motion correction first. There may be subtle differences in how these steps are applied, and sorters may make assumptions internally that the data will be preprocessed in a certain way.
Typically it is safe to perform filtering and common referencing in SpikeInterface. When overwriting more aggressive preprocessing steps (e.g. motion correction, whitening), it is worth checking the sorting outputs are improved compared to the sorter’s own defaults.
import spikeinterface.sorters as sssorting = ss.run_sorter("kilosort4", recording_whiten, folder="_tmp/ks4_out", remove_existing_folder=True, verbose=True, skip_kilosort_preprocessing=True, # turn off whitening and other steps do_CAR=False, # turn off CAR and filtering to be sure highpass_cutoff=0.01, do_correction=False, # turn off motion correction)print(sorting)
