MNE-phantom-KIT-data BIDS-ification script

bidsify_kit_phantom.py

"""BIDSify and trim the KIT phantom data."""
from pathlib import Path

import numpy as np
from scipy.signal import find_peaks
import matplotlib.pyplot as plt

import mne
import mne_bids
from mne.io.constants import FIFF

data_path = Path(__file__).parent  # should be run from old dir
raw = mne.io.read_raw_kit(
    data_path / "002_phantom_11Hz_100uA.con",
    stim=[166, *range(176, 190)],
    slope="+",
    stim_code="channel",
    stimthresh=2,
).load_data()

# To find the events, we can look at the MISC channel that recorded the activations.
# Here we use a very simple thresholding approach to find the events.
# The MISC 017 channel holds the dipole activations, which are 2-cycle 11 Hz sinusoidal
# bursts with the initial sinusoidal deflection downward, so we do a little bit of
# signal manipulation to help :func:`~scipy.signal.find_peaks`.

# cut from ~800 to ~300s for speed, and also at convenient dip stim boundaries
# chosen by examining MISC 017 by eye.
raw.crop(0, 302.9).load_data()  # also could crop at 829.7

# Figure out events
dip_act, dip_t = raw["MISC 017"]
dip_act = dip_act[0]  # 2D to 1D array
dip_act -= dip_act.mean()  # remove DC offset
dip_act *= -1  # invert so first deflection is positive
thresh = np.percentile(dip_act, 90)
dip_freq = 11.0
min_dist = raw.info["sfreq"] / dip_freq * 0.9  # 90% of period, to be safe
peaks = find_peaks(dip_act, height=thresh, distance=min_dist)[0]
assert len(peaks) % 2 == 0  # 2-cycle modulations
peaks = peaks[::2]  # take only first peaks of each 2-cycle burst
# From looking at ax.plot(dip_t, dip_act) and zooming in, we need to cut out the
# first stimulation interval and the last because they are incomplete
t_window = (11.7, 829.7)  # from entire file
peaks = peaks[(dip_t[peaks] > t_window[0]) & (dip_t[peaks] < t_window[1])]

fig, ax = plt.subplots(layout="constrained", figsize=(12, 4))
stop = int(30 * raw.info["sfreq"])
ax.plot(dip_t[:stop], dip_act[:stop], color="k", lw=1)
ax.axhline(thresh, color="r", ls="--", lw=1)
peak_idx = peaks[peaks < stop]
ax.plot(dip_t[peak_idx], dip_act[peak_idx], "ro", zorder=5, ms=5)
ax.set(xlabel="Time (s)", ylabel="Dipole activation (AU)\n(MISC 017 adjusted)")
ax.set(xlim=dip_t[[0, stop - 1]])

# We know that there are 32 dipoles, so mark the first ones as well
n_dip = 49
# shift to start
onsets = peaks - np.round(raw.info["sfreq"] / dip_freq / 4.0).astype(int)
onset_idx = onsets[onsets < stop]
n_rep = len(peaks) // n_dip
events = np.zeros((len(peaks), 3), int)
events[:, 0] = onsets + raw.first_samp
events[:, 2] = np.tile(np.arange(1, n_dip + 1), n_rep)
raw.plot(events=events)
event_id = {f"dip{ii:02d}": ii for ii in range(1, n_dip + 1)}

# Sanity check and determine epoching params
deltas = np.diff(events[:, 0], axis=0)
group_deltas = deltas[n_dip - 1 :: n_dip] / raw.info["sfreq"]  # gap between 49 and 1
assert (group_deltas > 0.8).all()
assert (group_deltas < 0.9).all()
others = np.delete(deltas, np.arange(n_dip - 1, len(deltas), n_dip))  # remove 49->1
others = others / raw.info["sfreq"]
assert (others > 0.25).all()
assert (others < 0.3).all()
tmax = 1 / dip_freq * 2.0  # 2 cycles
tmin = tmax - others.min()
assert tmin < 0

# We need to correct the ``dev_head_t`` because it's incorrect for these data!
# relative to the helmet: hleft, forward, up
translation = mne.transforms.translation(x=0.01, y=-0.015, z=-0.088)
# pitch down (rot about x/R), roll left (rot about y/A), yaw left (rot about z/S)
rotation = mne.transforms.rotation(
    x=np.deg2rad(5),
    y=np.deg2rad(-1),
    z=np.deg2rad(-3),
)
raw.info["dev_head_t"]["trans"][:] = translation @ rotation

# Drop some channels to reduce file size
raw.drop_channels([
    ch_name for ch_name in raw.ch_names
    if not (
        ch_name.startswith("MEG") or
        ch_name in ("MISC 001", "MISC 002", "MISC 003", "MISC 017")
    )
])
# Avoid MNE-BIDS read warning
for ch in raw.info['chs'][-4:]:
    assert ch['ch_name'].startswith('MISC')
    ch["unit"] = FIFF.FIFF_UNIT_NONE

# Finally let's write the results
bids_root = Path("..") / "MNE-phantom-KIT-data"
bids_root.mkdir(exist_ok=True)
mne_bids.make_dataset_description(
    path=bids_root,
    name="KIT phantom data",
    authors=["Paul Sowman", "Judy Zhu", "Eric Larson"],
    how_to_acknowledge='If you use this data, please cite MNE-Python and MNE-BIDS.',
    data_license='CC-BY-SA',
    references_and_links=[],
    overwrite=True,
)
(bids_root / "version.txt").write_text("0.2")  # for mne.datasets use
(bids_root / 'README').write_text("""
KIT Phantom Data
================

This dataset contains MEG data recorded on a KIT system using a 49-dipole phantom.
It is a cropped version of the original dataset in https://osf.io/svnt3.

The data were converted with MNE-BIDS and MNE-Python:

- Appelhoff, S., Sanderson, M., Brooks, T., Vliet, M., Quentin, R., Holdgraf, C., Chaumon, M., Mikulan, E., Tavabi, K., Höchenberger, R., Welke, D., Brunner, C., Rockhill, A., Larson, E., Gramfort, A. and Jas, M. (2019). MNE-BIDS: Organizing electrophysiological data into the BIDS format and facilitating their analysis. Journal of Open Source Software 4: (1896). https://doi.org/10.21105/joss.01896
- Niso, G., Gorgolewski, K. J., Bock, E., Brooks, T. L., Flandin, G., Gramfort, A., Henson, R. N., Jas, M., Litvak, V., Moreau, J., Oostenveld, R., Schoffelen, J., Tadel, F., Wexler, J., Baillet, S. (2018). MEG-BIDS, the brain imaging data structure extended to magnetoencephalography. Scientific Data, 5, 180110. https://doi.org/10.1038/sdata.2018.110
- Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., Goj, R., Jas, M., Brooks, T., Parkkonen, L., & Hämäläinen, M. (2013). MEG and EEG data analysis with MNE-Python. Frontiers in Neuroscience, 7. https://doi.org/10.3389/fnins.2013.00267

For BIDS-ification code, see:

https://gist.github.com/larsoner/e1f856d8b98a45a8c8896170609c59f9
""".strip())
bids_path = mne_bids.BIDSPath(
    subject="01",
    task="phantom",
    run="01",
    suffix="meg",
    datatype="meg",
    root=bids_root,
)
mne_bids.write_raw_bids(
    raw,
    bids_path,
    events=events,
    event_id=event_id,
    overwrite=True,
    allow_preload=True,
    format="FIF",
)