[1]:
import os
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
import openneuro
from ndspflow.workflows.workflow import WorkFlow
from neurodsp.filt import filter_signal
from neurodsp.spectral import compute_spectrum
from fooof import FOOOF, FOOOFGroup
from bycycle import Bycycle
BIDS WorkFlow¶
Below, we fetch a BIDS dataset using the openneuro cli. Then define and execute the following workflow:
Load the BIDS directory using MNE BIDS
Apply custom functions to select channels and epoch the signal
Low-pass filter with neurodsp
Fit FOOOF
Fit Bycycle
[ ]:
%%capture
dataset = 'ds003844'
include = ['sub-RESP0059/ses-SITUATION1A', 'sub-RESP0280/ses-SITUATION1A']
# Create BIDs folder
bids_path = f'{os.getcwd()}/bids'
if not os.path.isdir(bids_path):
os.makedirs(bids_path)
openneuro.download(dataset=dataset, target_dir=bids_path, include=include)
[3]:
def select_channels(y_array, channels):
"""Select sub-set of channels."""
inds = [ind for ind, ch in enumerate(channels)
if ch in channels]
return y_array[inds]
def epoch(y_array, n_epochs, epoch_len):
"""Epoch function."""
y_array = y_array[:, :int(n_epochs) * epoch_len]
y_array = y_array.reshape(-1, epoch_len)
return y_array
[4]:
wf = WorkFlow(bids_path=bids_path, session='SITUATION1A', task='acute')
wf.read_bids(subject='sub-RESP0059', allow_ragged=True)
# Custom transforms
wf.transform(select_channels, ['Gr16', 'Gr17'])
wf.transform(epoch, 2, 5000)
wf.transform(filter_signal, 'self.fs', 'lowpass', 200, remove_edges=False)
wf.fork(0)
wf.transform(compute_spectrum, 'self.fs')
wf.fit(FOOOF(verbose=False, max_n_peaks=5), (1, 100), axis=-1)
wf.fork(0)
wf.fit(Bycycle(), 'self.fs', (50, 100), axis=-1)
wf.run(n_jobs=1, progress=tqdm)
[5]:
wf.results.tolist()
[5]:
[[[<fooof.objs.fit.FOOOF at 0x7f9618bf3af0>,
<fooof.objs.fit.FOOOF at 0x7f961bc99a60>,
<fooof.objs.fit.FOOOF at 0x7f9619c689d0>,
<fooof.objs.fit.FOOOF at 0x7f9619c681c0>],
[<bycycle.objs.fit.Bycycle at 0x7f961bc95970>,
<bycycle.objs.fit.Bycycle at 0x7f961bc95d30>,
<bycycle.objs.fit.Bycycle at 0x7f961bc95310>,
<bycycle.objs.fit.Bycycle at 0x7f961bdb3790>]],
[[<fooof.objs.fit.FOOOF at 0x7f9618b9f7f0>,
<fooof.objs.fit.FOOOF at 0x7f961bc997c0>,
<fooof.objs.fit.FOOOF at 0x7f9619c688e0>,
<fooof.objs.fit.FOOOF at 0x7f9619c68310>],
[<bycycle.objs.fit.Bycycle at 0x7f961bc8f610>,
<bycycle.objs.fit.Bycycle at 0x7f961bc8f2e0>,
<bycycle.objs.fit.Bycycle at 0x7f961bc8f070>,
<bycycle.objs.fit.Bycycle at 0x7f961bc8f9d0>]]]
Inspecting Workflows¶
Preprocessing workflows may be inspected by calling .run()
with no .fit()
defined in the workflow. The resulting array (.y_array
) is collected below, from a series of .transform()
calls.
[6]:
wf = WorkFlow(bids_path=bids_path, session='SITUATION1A', task='acute')
wf.read_bids(subject='sub-RESP0059', allow_ragged=True)
wf.transform(select_channels, ['Gr16', 'Gr17'])
wf.transform(epoch, 2, 5000)
wf.transform(filter_signal, 'self.fs', 'lowpass', 200, remove_edges=False)
wf.run()
[7]:
# Plot example channels
plt.figure(figsize=(12, 2))
for subj in wf.y_array:
for ch in subj[2:]:
plt.plot(ch)
plt.ylabel('Voltage')
plt.xlabel('Time (samples)');
