Complete pipeline

Introductory notes:

This notebook presents minimal functionality needed to go through the cleaning, ICA, spectral, and event analysis stages.

  • For the cleaning part, the functionality consists of resampling, filtering, bad channels and bad data spans annotation, and bad channels interpolation.

  • For the ICA part, it is fitting and selecting components you want to exclude.

  • For the spectral analyses part, it is spectrogram+hypnogram, PSD per sleep stage, and topomap per sleep stage per frequency band.

  • For the detection of the events, it is spindles, slow waves, and rapid eye movements detection and analysis.

For the extended functionality check out the corresponding notebooks.

Import module

import os
import sys
import sleepeegpy

Run the following code

###### This code contains widgets that will be previewed later in the notebook. ######
####### Run this code and continue. ######
import ipywidgets as widgets
from IPython.display import display

use_example_set = widgets.RadioButtons(
    options=[True, False], description="Use example set?", value=True
)
bad_channels_selection = widgets.RadioButtons(
    options=["automatically", "manually", "from-file"],
    description="Clean Bad Channels:",
    value="from-file",
)
annotations_selection = widgets.RadioButtons(
    options=["automatically", "manually", "from-file"],
    description="Clean bad epochs Annotations:",
    value="from-file",
)

hypno_selection = widgets.RadioButtons(
    options=["automatically", "from-file"], description="Hypnogram:", value="from-file"
)
import pooch
from sleepeegpy.pipeline import (
    CleaningPipe,
    ICAPipe,
    SpectralPipe,
    SpindlesPipe,
    SlowWavesPipe,
    RapidEyeMovementsPipe,
)
from sleepeegpy.dashboard import create_dashboard

Setup Input Files

By default, all the input files are assumed to be saved in input_files, which will be created (if not already exist) in the notebook path. Change the following strings to use another path. Changing the output directory is also optional.

from os import makedirs

output_dir = "output_folder"  # Output path and name can be changed here
input_dir = "input_files"  # input files dir can be changed here
makedirs(input_dir, exist_ok=True)
makedirs(output_dir, exist_ok=True)

Pipeline preference:

Change the values according to how you want to use the pipeline.

display(use_example_set)
display(bad_channels_selection)
display(annotations_selection)
display(hypno_selection)

Add requested files

Run the code below, add all requested files (according to your preference) in the input file, and write their names.

Instructions:

  • Put your files in the input directory.

  • EEG file must include a montage that works with mne.read_raw.

  • For more information about the mne.read_raw supported formats, see mne documentation

  • If your file has no montage, but there is a channel mapping, you can add montage in the cleaning part (there is an example below of one way to add it).

#### If you use your own files (not the example set), write the relevant filenames below. ###

subject_code = "DK8"  # The subject code that will appear in the dashboard
eeg_file = "resampled_raw.fif"
bad_channels = "bad_channels.txt"
annotations = "annotations.txt"

If you are using hypnogram file:

  • Modify your hypnogram file name (Point-per-row type of hypnogram) below.

  • If needed, change the Hypnogram’s sampling frequency

  • If you don’t have a hypnogram file, follow the next notebook instructions.

hypnogram_filename = "staging.txt"
hypno_freq = 1  # If required, change the Hypnogram's sampling frequency (Visbrain's hypnograms default to 1)

Modify hypnogram prediction arguments

If you don’t have hypnogram file, choose prediction arguments. These values will use YASA’s algorithm to create a new hypnogram.

Make sure the selected channel names exist in your montage

if hypno_selection.value == "automatically":
    hypnogram_path = "predict"
    #### If you selected automatic hypnogram, select prediction arguments. ####
    hypno_predict_arguments = {
        "eeg_name": "E183",  # (C4)
        "eog_name": "E252",  # (EOG1)
        "emg_name": "E247",
        "ref_name": "E26",
        "save": False,
    }
else:
    hypno_predict_arguments = None

Adjust variables

  • If required, change n_components - should equal or less than the number of channels. see more information

  • Modify the picked channel. picked_channel represents the EEG channel selected for plotting and computations of hypnospectrogram, PSDs, and TFRs.

  • Modify the variable in the next notebook cells, to have separate channels for each plot.

  • Modify loc_chname and roc_chname (EOG). For more information see: Left and Right Ocular Canthi.

  • Make sure all channel names are part of the montage.

n_components = 30
picked_channel = "E101"  # (Pz)
loc_chname = "E252"  # (EOG1)
roc_chname = "E226"  # (EOG2)

You can now run the notebook (you can change values according to your needs)

if use_example_set.value:
    cache_dir = pooch.os_cache("sleepeegpy_dataset")
    doi = "10.5281/zenodo.10362189"
    odie = pooch.create(
        path=cache_dir,
        base_url=f"doi:{doi}",
    )
    odie.load_registry_from_doi()
if use_example_set.value:
    ## Nap dataset files: ##
    bad_channels = odie.fetch("nap_bad_channels.txt")
    annotations = odie.fetch("nap_annotations.txt")
    path_to_eeg = odie.fetch("nap_resampled_raw.fif", progressbar=True)
    if hypno_selection.value == "from-file":
        hypnogram_path = odie.fetch("nap_staging.txt", progressbar=True)
else:
    path_to_eeg = os.path.join(input_dir, eeg_file)
    if hypno_selection.value == "from-file":
        hypnogram_path = os.path.join(input_dir, hypnogram_filename)
    bad_channels = (
        None
        if bad_channels_selection == "automatically"
        else os.path.join(input_dir, bad_channels)
    )
    annotations = (
        None
        if annotations_selection == "automatically"
        else os.path.join(input_dir, annotations)
    )

Cleaning

Initialize CleaningPipe object by providing it with path to eeg file and output directory in which you want the data to be saved.

pipe = CleaningPipe(
    path_to_eeg=path_to_eeg,
    output_dir=output_dir,
)
##### Adding montage example ####
# import mne
# channels_map = {"E59": "C3", "E183": "C4", "E36": "F3", "E224": "F4", "E47": "F7", "E2": "F8","E37": "Fp1", "E18": "Fp2", "E21": "Fz", "E116": "O1", "E150": "O2", "E87": "P3", "E153": "P4", "E101": "Pz","E69": "T3", "E202": "T4", "E96": "T5", "E170": "T6", "E94": "A1", "E190": "A2"}
# current_channel_names = set(pipe.mne_raw.ch_names)
# channels_to_drop = list(current_channel_names - set(channels_map.keys()))
# pipe.mne_raw.drop_channels(channels_to_drop)
# mne.rename_channels(pipe.mne_raw.info, channels_map)
# montage = mne.channels.make_standard_montage('standard_1020')
# pipe.mne_raw.set_montage(montage)

Resampling

This can take more than an hour depending on eeg signal size and specs of the computer you’re running the analysis on.

pipe.resample(sfreq=250)
pipe.filter(l_freq=0.75, h_freq=40)
pipe.notch(freqs="50s")

Select bad channels and epochs

If manually bad channels were selected, select bad channels in the pop-up window. Note that automatically bad channels selection takes time.

if bad_channels_selection.value == "manually":
    pipe.plot(save_bad_channels=True)
elif bad_channels_selection.value == "automatically":
    bad_channels = pipe.auto_detect_bad_channels()
pipe.read_bad_channels(
    path=None if bad_channels_selection.value == "Manually" else bad_channels
)
pipe.interpolate_bads(reset_bads=True)
Interpolated channels: ['E9', 'E45', 'E54', 'E122', 'E200', 'E256', 'VREF']

Select bad epochs

Click “a” -> “Add description” -> Enter BAD_EPOCH -> Annotate bad data spans

if annotations_selection.value == "manually":
    pipe.plot(butterfly=True, save_annotations=True, overwrite=True)
    pipe.read_annotations()
elif annotations_selection.value == "from-file":
    pipe.read_annotations(path=annotations)
elif annotations_selection.value == "automatically":
    pipe.auto_set_annotations()
fig = create_dashboard(
    subject_code=subject_code,
    prec_pipe=pipe,
    hypno_psd_pick=picked_channel,
    hypnogram=hypnogram_path,
    predict_hypno_args=hypno_predict_arguments,
    hypno_freq=hypno_freq,
    reference="average",
)
average reference has been applied
Hypnogram is SHORTER than data by 0.0 seconds. Padding hypnogram with last value to match data.size.
../_images/946106439de8e172d27e738d1066f2fbd24fe092dac5926dccdb3ec3f229ee1a.png

ICA

Pass the preceding (cleaning) pipe to the ICAPipe.

ica_pipe = ICAPipe(prec_pipe=pipe, n_components=n_components)

Fit the ICA on the 1 Hz high-pass filtered data.

ica_pipe.fit()

Visually inspect ICA components.

ica_pipe.plot_sources()
2024-10-15 11:30:57,649 - OpenGL.acceleratesupport - INFO - No OpenGL_accelerate module loaded: No module named 'OpenGL_accelerate'
<mne_qt_browser._pg_figure.MNEQtBrowser at 0x3a73945e0>

Pass to the exclude argument indices of components you want to remove from the raw signal.

ica_pipe.apply()
Excluded ICA components: []

Spectral

Pass the preceding (cleaning or ICA) pipe to the SpectralPipe. Also provide pass to the hypnogram and don’t forget to pass its frequency to the corresponding parameter.

spectral_pipe = SpectralPipe(
    prec_pipe=ica_pipe,
    path_to_hypno=hypnogram_path,
    hypno_freq=hypno_freq,
)
Hypnogram is SHORTER than data by 0.0 seconds. Padding hypnogram with last value to match data.size.

If you don’t have a hypnogram, predict_hypno will use YASA’s algorithm.

if hypnogram_path is None:
    spectral_pipe.predict_hypno(
        eeg_name=hypno_predict_arguments["eeg_name"],
        eog_name=hypno_predict_arguments["eog_name"],
        emg_name=hypno_predict_arguments["emg_name"],
        ref_name=hypno_predict_arguments["ref_name"],
        save=False,
    )

By default, picked_channel will be used to calculate spectrogram. You can pass another electrode name (make sure it exists)

spectral_pipe.plot_hypnospectrogram(picks=[picked_channel])
../_images/3cb0a583e686957aac6d050520d37d4fc1d276910664368974f2d749803bc64f.png
spectral_pipe.compute_psd(
    sleep_stages={"Wake": 0, "N1": 1, "N2/3": (2, 3), "REM": 4},
    reference="average",
    # Additional arguments passed to the Welch method:
    n_fft=1024,
    n_per_seg=1024,
    n_overlap=512,
    window="hamming",
    verbose=False,
)
spectral_pipe.plot_psds(picks=[picked_channel], psd_range=(-30, 30))
../_images/d23c74c3774e1c970cfcfc526cfd18857b3e97ed86a163524891565a4d40904f.png

Create a collage with rows for sleep stages and columns for bands.

spectral_pipe.plot_topomap_collage()
../_images/316762981c18eee9341e03e385f22d06582c32284b87c43d57c1f2eef3097a30.png

Events

Pass the preceding (cleaning or ICA or spectral) pipe to one of the SpindlesPipe, SlowWavesPipe, or RapidEyeMovementsPipe. If the preceding is cleaning or ICA - provide a path to the hypnogram and don’t forget to pass its frequency to the corresponding parameter.

spindles_pipe = SpindlesPipe(prec_pipe=spectral_pipe)

spindles_pipe.detect()
spindles_pipe.plot_average(
    center="Peak",
    hue="Stage",
    time_before=1,
    time_after=1,
)
../_images/fc9d014fcc397ada6c161e5318dc2da5850afae854eae060a87048153ab0bd4f.png
spindles_pipe.results.summary(grp_chan=False, grp_stage=True)
Count Density Duration Amplitude RMS AbsPower RelPower Frequency Oscillations Symmetry
Stage
1 6938 365.157895 0.869961 35.703243 7.572310 1.669308 0.370737 13.764517 11.291150 0.485610
2 20299 738.145455 0.866909 40.815780 8.733241 1.786240 0.363001 13.773122 11.369723 0.495506
3 907 106.705882 0.724123 41.968714 9.486414 1.902562 0.334536 13.590072 9.345094 0.459087
spindles_pipe.compute_tfr(freqs=(10, 20), n_freqs=100, time_before=1, time_after=1)
spindles_pipe.tfrs["N2"].plot([picked_channel])
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 257/257 [00:15<00:00, 16.48it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 257/257 [00:28<00:00,  9.13it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 213/213 [00:03<00:00, 56.33it/s]
../_images/bd5ff733ba5eb1a5750e07d1d5216cc7140edf122ba6eacfb16e0a3458f1dd42.png
[<Figure size 640x480 with 2 Axes>]
slow_waves_pipe = SlowWavesPipe(prec_pipe=spindles_pipe)
slow_waves_pipe.detect()

slow_waves_pipe.plot_average(
    center="NegPeak",
    hue="Stage",
    time_before=0.4,
    time_after=0.8,
)
Hypnogram is SHORTER than data by 0.0 seconds. Padding hypnogram with last value to match data.size.
Setting up band-pass filter from 0.3 - 1.5 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.30
- Lower transition bandwidth: 0.20 Hz (-6 dB cutoff frequency: 0.20 Hz)
- Upper passband edge: 1.50 Hz
- Upper transition bandwidth: 0.20 Hz (-6 dB cutoff frequency: 1.60 Hz)
- Filter length: 4125 samples (16.500 s)
../_images/6f64acd8ca17470aa210631dea5f42b06ab0c65c477bf2ea28f8f79d8107290b.png
slow_waves_pipe.compute_tfr(
    freqs=(0.5, 5), n_freqs=100, time_before=4, time_after=4, n_cycles=2
)
ploted_channel = "VREF"  # Cz
slow_waves_pipe.tfrs["N3"].plot([ploted_channel])
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 109/109 [00:05<00:00, 19.18it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 203/203 [00:10<00:00, 18.78it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 168/168 [00:08<00:00, 18.87it/s]
../_images/c0f61f5c9ddf0edc04c402237bc24dd7552f744b9f977e182526ee5d4bb9594e.png
[<Figure size 640x480 with 2 Axes>]
rems_pipe = RapidEyeMovementsPipe(prec_pipe=slow_waves_pipe)

rems_pipe.detect(
    loc_chname=loc_chname,
    roc_chname=roc_chname,
)

rems_pipe.plot_average(
    center="Peak",
    time_before=0.5,
    time_after=0.5,
    filt=(None, None),
    mask=None,
)
Hypnogram is SHORTER than data by 0.0 seconds. Padding hypnogram with last value to match data.size.
Setting up band-pass filter from 0.5 - 5 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 0.50
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 0.25 Hz)
- Upper passband edge: 5.00 Hz
- Upper transition bandwidth: 2.00 Hz (-6 dB cutoff frequency: 6.00 Hz)
- Filter length: 1651 samples (6.604 s)
../_images/bfb238c26733030cc3bd85150561cc4c6c3b776d45819e3ae1747f6f81a560e0.png