{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Working with CTF data: the Brainstorm auditory dataset\n\nHere we compute the evoked from raw for the auditory Brainstorm\ntutorial dataset. For comparison, see :footcite:`TadelEtAl2011` and the\nassociated [brainstorm site](https://neuroimage.usc.edu/brainstorm/Tutorials/Auditory).\n\nExperiment:\n\n - One subject, 2 acquisition runs 6 minutes each.\n - Each run contains 200 regular beeps and 40 easy deviant beeps.\n - Random ISI: between 0.7s and 1.7s seconds, uniformly distributed.\n - Button pressed when detecting a deviant with the right index finger.\n\nThe specifications of this dataset were discussed initially on the\n[FieldTrip bug tracker](http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2300)_.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Mainak Jas \n# Eric Larson \n# Jaakko Leppakangas \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport pandas as pd\n\nimport mne\nfrom mne import combine_evoked\nfrom mne.datasets.brainstorm import bst_auditory\nfrom mne.io import read_raw_ctf\nfrom mne.minimum_norm import apply_inverse" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To reduce memory consumption and running time, some of the steps are\nprecomputed. To run everything from scratch change ``use_precomputed`` to\n``False``. With ``use_precomputed = False`` running time of this script can\nbe several minutes even on a fast computer.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "use_precomputed = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data was collected with a CTF 275 system at 2400 Hz and low-pass\nfiltered at 600 Hz. Here the data and empty room data files are read to\nconstruct instances of :class:`mne.io.Raw`.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = bst_auditory.data_path()\n\nsubject = \"bst_auditory\"\nsubjects_dir = data_path / \"subjects\"\n\nraw_fname1 = data_path / \"MEG\" / subject / \"S01_AEF_20131218_01.ds\"\nraw_fname2 = data_path / \"MEG\" / subject / \"S01_AEF_20131218_02.ds\"\nerm_fname = data_path / \"MEG\" / subject / \"S01_Noise_20131218_01.ds\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the memory saving mode we use ``preload=False`` and use the memory\nefficient IO which loads the data on demand. However, filtering and some\nother functions require the data to be preloaded into memory.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw = read_raw_ctf(raw_fname1)\nn_times_run1 = raw.n_times\n\n# Here we ignore that these have different device<->head transforms\nmne.io.concatenate_raws([raw, read_raw_ctf(raw_fname2)], on_mismatch=\"ignore\")\nraw_erm = read_raw_ctf(erm_fname)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data array consists of 274 MEG axial gradiometers, 26 MEG reference\nsensors and 2 EEG electrodes (Cz and Pz). In addition:\n\n - 1 stim channel for marking presentation times for the stimuli\n - 1 audio channel for the sent signal\n - 1 response channel for recording the button presses\n - 1 ECG bipolar\n - 2 EOG bipolar (vertical and horizontal)\n - 12 head tracking channels\n - 20 unused channels\n\nNotice also that the digitized electrode positions (stored in a .pos file)\nwere automatically loaded and added to the `~mne.io.Raw` object.\n\nThe head tracking channels and the unused channels are marked as misc\nchannels. Here we define the EOG and ECG channels.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.set_channel_types({\"HEOG\": \"eog\", \"VEOG\": \"eog\", \"ECG\": \"ecg\"})\nif not use_precomputed:\n # Leave out the two EEG channels for easier computation of forward.\n raw.pick([\"meg\", \"stim\", \"misc\", \"eog\", \"ecg\"]).load_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For noise reduction, a set of bad segments have been identified and stored\nin csv files. The bad segments are later used to reject epochs that overlap\nwith them.\nThe file for the second run also contains some saccades. The saccades are\nremoved by using SSP. We use pandas to read the data from the csv files. You\ncan also view the files with your favorite text editor.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "annotations_df = pd.DataFrame()\noffset = n_times_run1\nfor idx in [1, 2]:\n csv_fname = data_path / \"MEG\" / \"bst_auditory\" / f\"events_bad_0{idx}.csv\"\n df = pd.read_csv(csv_fname, header=None, names=[\"onset\", \"duration\", \"id\", \"label\"])\n print(f\"Events from run {idx}:\")\n print(df)\n\n df[\"onset\"] += offset * (idx - 1)\n annotations_df = pd.concat([annotations_df, df], axis=0)\n\nsaccades_events = df[df[\"label\"] == \"saccade\"].values[:, :3].astype(int)\n\n# Conversion from samples to times:\nonsets = annotations_df[\"onset\"].values / raw.info[\"sfreq\"]\ndurations = annotations_df[\"duration\"].values / raw.info[\"sfreq\"]\ndescriptions = annotations_df[\"label\"].values\n\nannotations = mne.Annotations(onsets, durations, descriptions)\nraw.set_annotations(annotations)\ndel onsets, durations, descriptions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we compute the saccade and EOG projectors for magnetometers and add\nthem to the raw data. The projectors are added to both runs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "saccade_epochs = mne.Epochs(\n raw,\n saccades_events,\n 1,\n 0.0,\n 0.5,\n preload=True,\n baseline=(None, None),\n reject_by_annotation=False,\n)\n\nprojs_saccade = mne.compute_proj_epochs(\n saccade_epochs, n_mag=1, n_eeg=0, desc_prefix=\"saccade\"\n)\nif use_precomputed:\n proj_fname = data_path / \"MEG\" / \"bst_auditory\" / \"bst_auditory-eog-proj.fif\"\n projs_eog = mne.read_proj(proj_fname)[0]\nelse:\n projs_eog, _ = mne.preprocessing.compute_proj_eog(raw.load_data(), n_mag=1, n_eeg=0)\nraw.add_proj(projs_saccade)\nraw.add_proj(projs_eog)\ndel saccade_epochs, saccades_events, projs_eog, projs_saccade # To save memory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visually inspect the effects of projections. Click on 'proj' button at the\nbottom right corner to toggle the projectors on/off. EOG events can be\nplotted by adding the event list as a keyword argument. As the bad segments\nand saccades were added as annotations to the raw data, they are plotted as\nwell.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Typical preprocessing step is the removal of power line artifact (50 Hz or\n60 Hz). Here we notch filter the data at 60, 120 and 180 to remove the\noriginal 60 Hz artifact and the harmonics. The power spectra are plotted\nbefore and after the filtering to show the effect. The drop after 600 Hz\nappears because the data was filtered during the acquisition. In memory\nsaving mode we do the filtering at evoked stage, which is not something you\nusually would do.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if not use_precomputed:\n raw.compute_psd(tmax=np.inf, picks=\"meg\").plot(\n picks=\"data\", exclude=\"bads\", amplitude=False\n )\n notches = np.arange(60, 181, 60)\n raw.notch_filter(notches, phase=\"zero-double\", fir_design=\"firwin2\")\n raw.compute_psd(tmax=np.inf, picks=\"meg\").plot(\n picks=\"data\", exclude=\"bads\", amplitude=False\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also lowpass filter the data at 100 Hz to remove the hf components.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if not use_precomputed:\n raw.filter(\n None,\n 100.0,\n h_trans_bandwidth=0.5,\n filter_length=\"10s\",\n phase=\"zero-double\",\n fir_design=\"firwin2\",\n )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Epoching and averaging.\nFirst some parameters are defined and events extracted from the stimulus\nchannel (UPPT001). The rejection thresholds are defined as peak-to-peak\nvalues and are in T / m for gradiometers, T for magnetometers and\nV for EOG and EEG channels.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmin, tmax = -0.1, 0.5\nevent_id = dict(standard=1, deviant=2)\nreject = dict(mag=4e-12, eog=250e-6)\n# find events\nevents = mne.find_events(raw, stim_channel=\"UPPT001\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The event timing is adjusted by comparing the trigger times on detected\nsound onsets on channel UADC001-4408.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sound_data = raw[raw.ch_names.index(\"UADC001-4408\")][0][0]\nonsets = np.where(np.abs(sound_data) > 2.0 * np.std(sound_data))[0]\nmin_diff = int(0.5 * raw.info[\"sfreq\"])\ndiffs = np.concatenate([[min_diff + 1], np.diff(onsets)])\nonsets = onsets[diffs > min_diff]\nassert len(onsets) == len(events)\ndiffs = 1000.0 * (events[:, 0] - onsets) / raw.info[\"sfreq\"]\nprint(f\"Trigger delay removed (\u03bc \u00b1 \u03c3): {np.mean(diffs):0.1f} \u00b1 {np.std(diffs):0.1f} ms\")\nevents[:, 0] = onsets\ndel sound_data, diffs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We mark a set of bad channels that seem noisier than others. This can also\nbe done interactively with ``raw.plot`` by clicking the channel name\n(or the line). The marked channels are added as bad when the browser window\nis closed.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.info[\"bads\"] = [\"MLO52-4408\", \"MRT51-4408\", \"MLO42-4408\", \"MLO43-4408\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The epochs (trials) are created for MEG channels. First we find the picks\nfor MEG and EOG channels. Then the epochs are constructed using these picks.\nThe epochs overlapping with annotated bad segments are also rejected by\ndefault. To turn off rejection by bad segments (as was done earlier with\nsaccades) you can use keyword ``reject_by_annotation=False``.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epochs = mne.Epochs(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n picks=[\"meg\", \"eog\"],\n baseline=(None, 0),\n reject=reject,\n preload=False,\n proj=True,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We only use first 40 good epochs from each run. Since we first drop the bad\nepochs, the indices of the epochs are no longer same as in the original\nepochs collection. Investigation of the event timings reveals that first\nepoch from the second run corresponds to index 182.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epochs.drop_bad()\n\n# avoid warning about concatenating with annotations\nepochs.set_annotations(None)\n\nepochs_standard = mne.concatenate_epochs(\n [epochs[\"standard\"][range(40)], epochs[\"standard\"][182:222]]\n)\nepochs_standard.load_data() # Resampling to save memory.\nepochs_standard.resample(600, npad=\"auto\")\nepochs_deviant = epochs[\"deviant\"].load_data()\nepochs_deviant.resample(600, npad=\"auto\")\ndel epochs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The averages for each conditions are computed.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked_std = epochs_standard.average()\nevoked_dev = epochs_deviant.average()\ndel epochs_standard, epochs_deviant" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Typical preprocessing step is the removal of power line artifact (50 Hz or\n60 Hz). Here we lowpass filter the data at 40 Hz, which will remove all\nline artifacts (and high frequency information). Normally this would be done\nto raw data (with :func:`mne.io.Raw.filter`), but to reduce memory\nconsumption of this tutorial, we do it at evoked stage. (At the raw stage,\nyou could alternatively notch filter with :func:`mne.io.Raw.notch_filter`.)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for evoked in (evoked_std, evoked_dev):\n evoked.filter(l_freq=None, h_freq=40.0, fir_design=\"firwin\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we plot the ERF of standard and deviant conditions. In both conditions\nwe can see the P50 and N100 responses. The mismatch negativity is visible\nonly in the deviant condition around 100-200 ms. P200 is also visible around\n170 ms in both conditions but much stronger in the standard condition. P300\nis visible in deviant condition only (decision making in preparation of the\nbutton press). You can view the topographies from a certain time span by\npainting an area with clicking and holding the left mouse button.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked_std.plot(window_title=\"Standard\", gfp=True, time_unit=\"s\")\nevoked_dev.plot(window_title=\"Deviant\", gfp=True, time_unit=\"s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show activations as topography figures.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "times = np.arange(0.05, 0.301, 0.025)\nfig = evoked_std.plot_topomap(times=times)\nfig.suptitle(\"Standard\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig = evoked_dev.plot_topomap(times=times)\nfig.suptitle(\"Deviant\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the MMN effect more clearly by looking at the difference between\nthe two conditions. P50 and N100 are no longer visible, but MMN/P200 and\nP300 are emphasised.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "evoked_difference = combine_evoked([evoked_dev, evoked_std], weights=[1, -1])\nevoked_difference.plot(window_title=\"Difference\", gfp=True, time_unit=\"s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Source estimation.\nWe compute the noise covariance matrix from the empty room measurement\nand use it for the other runs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "reject = dict(mag=4e-12)\ncov = mne.compute_raw_covariance(raw_erm, reject=reject)\ncov.plot(raw_erm.info)\ndel raw_erm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The transformation is read from a file:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "trans_fname = data_path / \"MEG\" / \"bst_auditory\" / \"bst_auditory-trans.fif\"\ntrans = mne.read_trans(trans_fname)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To save time and memory, the forward solution is read from a file. Set\n``use_precomputed=False`` in the beginning of this script to build the\nforward solution from scratch. The head surfaces for constructing a BEM\nsolution are read from a file. Since the data only contains MEG channels, we\nonly need the inner skull surface for making the forward solution. For more\ninformation: `CHDBBCEJ`, :func:`mne.setup_source_space`,\n`bem-model`, :func:`mne.bem.make_watershed_bem`.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if use_precomputed:\n fwd_fname = data_path / \"MEG\" / \"bst_auditory\" / \"bst_auditory-meg-oct-6-fwd.fif\"\n fwd = mne.read_forward_solution(fwd_fname)\nelse:\n src = mne.setup_source_space(\n subject, spacing=\"ico4\", subjects_dir=subjects_dir, overwrite=True\n )\n model = mne.make_bem_model(\n subject=subject, ico=4, conductivity=[0.3], subjects_dir=subjects_dir\n )\n bem = mne.make_bem_solution(model)\n fwd = mne.make_forward_solution(evoked_std.info, trans=trans, src=src, bem=bem)\n\ninv = mne.minimum_norm.make_inverse_operator(evoked_std.info, fwd, cov)\nsnr = 3.0\nlambda2 = 1.0 / snr**2\ndel fwd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sources are computed using dSPM method and plotted on an inflated brain\nsurface. For interactive controls over the image, use keyword\n``time_viewer=True``.\nStandard condition.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "stc_standard = mne.minimum_norm.apply_inverse(evoked_std, inv, lambda2, \"dSPM\")\nbrain = stc_standard.plot(\n subjects_dir=subjects_dir,\n subject=subject,\n surface=\"inflated\",\n time_viewer=False,\n hemi=\"lh\",\n initial_time=0.1,\n time_unit=\"s\",\n)\ndel stc_standard, brain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Deviant condition.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "stc_deviant = mne.minimum_norm.apply_inverse(evoked_dev, inv, lambda2, \"dSPM\")\nbrain = stc_deviant.plot(\n subjects_dir=subjects_dir,\n subject=subject,\n surface=\"inflated\",\n time_viewer=False,\n hemi=\"lh\",\n initial_time=0.1,\n time_unit=\"s\",\n)\ndel stc_deviant, brain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Difference.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "stc_difference = apply_inverse(evoked_difference, inv, lambda2, \"dSPM\")\nbrain = stc_difference.plot(\n subjects_dir=subjects_dir,\n subject=subject,\n surface=\"inflated\",\n time_viewer=False,\n hemi=\"lh\",\n initial_time=0.15,\n time_unit=\"s\",\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n.. footbibliography::\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }