Python MNE is an open source magnetoencephalogram (MEG) and electroencephalogram (EEG) analysis and visualization tool. It is highly versatile because it can be applied to the data formats of many devices. In this article, follow Most Basic Tutorials , MEG and EEG MNE analysis procedure is explained.

Execution environment

Mac OS 10.15.3
Python 3.6.5
mne 0.20.4

Installation

You can install it with either Anaconda or pip. On the Official page, Anaconda is recommended.

For Anaconda


conda install -c conda-forge mne

For pip


pip install -U mne

Also, GPU settings seems to be possible, so it takes time for filtering and sampling. If it takes, it may be better to set it.


$ conda install cupy
$ MNE_USE_CUDA=true python -c "import mne; mne.cuda.init_cuda(verbose=True)"
Enabling CUDA with 1.55 GB available memory

>>> mne.utils.set_config('MNE_USE_CUDA', 'true')  
# mne.io.Raw.filter()Etc. n_jobs='cuda'give

MEG / EEG analysis tutorial

Typical of EEG analysis according to this tutorial I will explain the basic processing and the basic data structure of MNE. I will explain it in a concise manner, so if you have knowledge of English and EEG analysis, you should read the original page.

Analysis procedure

[Load data](#Load data)
[Pre-processing](# Pre-processing)
[Event detection](#Event detection)
[Cut into epoch](#Cut into epoch)
[Time frequency analysis](#Time frequency analysis)
[Estimation of evoked response](# Estimate of evoked response)

Module import

import os
import numpy as np
import mne

Data loading

Load the data prepared as a sample in mne (List of all datasets). By the way, MEG and EEG have several file formats for each device, but there is a function to load those data (Supported Data Formats. stable / overview / implementation.html # data-formats))

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_filt-0-40_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file)

This time, we are using the data that has been low-pass filtered at 40Hz called sample_audvis_filt-0-40_raw.fif. It took about 1 minute to download the data (1.54GB).

The function mne.io.read_raw_fif () reads the data of the MNE data structure called Raw object from the file in fif format.

By the way, in the experiment of this data, the subject is randomly presented with visual stimulus and auditory stimulus every 750 ms. In addition, a laughing face is sometimes displayed on the screen, and the subject presses a button when he sees it. The presentation of these stimuli is called an event.

Display the loaded object.

print(raw)
print(raw.info)

# output
<Raw | sample_audvis_filt-0-40_raw.fif, 376 x 41700 (277.7 s), ~3.6 MB, data not loaded>
<Info | 15 non-empty values
 bads: 2 items (MEG 2443, EEG 053)
 ch_names: MEG 0113, MEG 0112, MEG 0111, MEG 0122, MEG 0123, MEG 0121, MEG ...
 chs: 204 GRAD, 102 MAG, 9 STIM, 60 EEG, 1 EOG
 custom_ref_applied: False
 dev_head_t: MEG device -> head transform
 dig: 146 items (3 Cardinal, 4 HPI, 61 EEG, 78 Extra)
 file_id: 4 items (dict)
 highpass: 0.1 Hz
 hpi_meas: 1 item (list)
 hpi_results: 1 item (list)
 lowpass: 40.0 Hz
 meas_date: 2002-12-03 19:01:10 UTC
 meas_id: 4 items (dict)
 nchan: 376
 projs: PCA-v1: off, PCA-v2: off, PCA-v3: off, Average EEG reference: off
 sfreq: 150.2 Hz
>

For MNE, Raw, Epochs /mne.Epochs.html#mne.Epochs), Evoked Data is managed by objects called objects. These are distinguished by the time division of the data, and Raw contains the data of all the time from the start to the stop of EEG recording, and it is good to think that it is the state before being separated by the event. I think. raw.info holds information such as the number of channels, names, and sampling rate of raw.

The power spectral density (PSD) of this Raw object and the waveform for each sensor are displayed. plot () can be scrolled and scaled in an interactive shell.

raw.plot_psd(fmax=50)
raw.plot(duration=5, n_channels=30)

Preprocessing

In the noise processing of EEG / MEG analysis, the standard independent component analysis (ICA) is performed. Simply put, by decomposing a multi-channel signal into independent components, processing that removes things that are independent of brain activity, such as power supply noise and blinking, is performed.

Here, for the signal of 364 channels, after decomposing into 20 independent components, [1, 2] and the index of the component to be removed are specified (actual removal is the following ʻica.apply). It is done by (raw) ). Display the information of independent components with ʻica_properties () .

# set up and fit the ICA
ica = mne.preprocessing.ICA(n_components=20, random_state=97, max_iter=800)
ica.fit(raw)
ica.exclude = [1, 2]  # details on how we picked these are omitted here
ica.plot_properties(raw, picks=ica.exclude)

# output
Fitting ICA to data using 364 channels (please be patient, this may take a while)
Inferring max_pca_components from picks
Selecting by number: 20 components
Fitting ICA took 2.4s.
    Using multitaper spectrum estimation with 7 DPSS windows
138 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
138 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped

In order to get the signal with the specified independent component removed, apply the ICA object to the Raw object.

orig_raw = raw.copy()
raw.load_data()
ica.apply(raw)

When the front channel is displayed, the noise is removed as shown in the following two figures.

# show some frontal channels to clearly illustrate the artifact removal
chs = ['MEG 0111', 'MEG 0121', 'MEG 0131', 'MEG 0211', 'MEG 0221', 'MEG 0231',
       'MEG 0311', 'MEG 0321', 'MEG 0331', 'MEG 1511', 'MEG 1521', 'MEG 1531',
       'EEG 001', 'EEG 002', 'EEG 003', 'EEG 004', 'EEG 005', 'EEG 006',
       'EEG 007', 'EEG 008']
chan_idxs = [raw.ch_names.index(ch) for ch in chs]
orig_raw.plot(order=chan_idxs, start=12, duration=4)
raw.plot(order=chan_idxs, start=12, duration=4)

# output
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
Transforming to ICA space (20 components)
Zeroing out 2 ICA components

Event detection

This sample data includes what is called the STIM channel, which contains data on the event type and its time. This time, the channel named STI 014 is applicable, so specify it and retrieve the event of the experiment.

events = mne.find_events(raw, stim_channel='STI 014')
print(events[:5])  # show the first 5

# output
319 events found
Event IDs: [ 1  2  3  4  5 32]

[[6994    0    2]
 [7086    0    3]
 [7192    0    1]
 [7304    0    4]
 [7413    0    2]]

ʻEventsis a numpy array, the first column represents the sampling time and the last column represents the event type (ID) (0 in the middle can be ignored). Here, it seems that there is an event with ID[1 2 3 4 5 32]`.

Event ID	Condition
1	auditory stimulus (tone) to the left ear
2	auditory stimulus (tone) to the right ear
3	visual stimulus (checkerboard) to the left visual field
4	visual stimulus (checkerboard) to the right visual field
5	smiley face (catch trial)
32	subject button press

The Event ID and the description of the event are stored in a dictionary type. By the way, the key here can be any character string.

event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4, 'smiley': 5, 'buttonpress': 32}

Let's plot the presentation timing of each event on the time axis. The ID defined earlier by ʻevent_id = even_dict` corresponds to the event name on the plot.

fig = mne.viz.plot_events(events, event_id=event_dict, sfreq=raw.info['sfreq'],
                          first_samp=raw.first_samp)

Cut out into an epoch

In MEG / EEG analysis, when investigating the response to a certain stimulus, the signal is divided into units called epochs and treated as one data. By cutting out a certain section before and after the stimulus, it is made into one data unit. This cutting process is called epoching. Raw objects are connected for the entire time during recording, but only the part of the response to the stimulus is cut out and stored in the data type called Epochs object.

In epoching, a process called threshold rejection is often performed. This means that if there is a signal that is too loud during the time, that part will have been heavily noisy and will be excluded from the data. Determine the criteria for the exclusion.

reject_criteria = dict(mag=4000e-15,     # 4000 fT
                       grad=4000e-13,    # 4000 fT/cm
                       eeg=150e-6,       # 150 µV
                       eog=250e-6)       # 250 µV

Next, create an Epochs object from the Raw object. The arguments of mne.Epochs here are as follows.

raw: Raw object containing signal data
ʻEvents`: numpy array with time and event ID information
ʻEvent_id`: Dictionary type that associates event IDs with their names
tmin: How many seconds of stimulus presentation to cut out (add a minus before stimulus presentation)
tmax: How many seconds after the stimulus is presented to cut out the data
reject: Data exclusion criteria

epochs = mne.Epochs(raw, events, event_id=event_dict, tmin=-0.2, tmax=0.5,
                    reject=reject_criteria, preload=True)

# output
319 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 4)
4 projection items activated
Loading data for 319 events and 106 original time points ...
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on MAG : ['MEG 1711']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on MAG : ['MEG 1711']
    Rejecting  epoch based on EEG : ['EEG 008']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
    Rejecting  epoch based on EOG : ['EOG 061']
10 bad epochs dropped

It seems that 10 epochs have been excluded this time.

Since the timing of pressing the button of the subject is not dealt with this time, only the audiovisual stimulation event is targeted. In addition, in order to prevent the number of data from being biased among the conditions, the number of each condition is made uniform by deleting the events with a large number of conditions with ʻepochs.equalize_event_counts () `.

conds_we_care_about = ['auditory/left', 'auditory/right',
                       'visual/left', 'visual/right']
epochs.equalize_event_counts(conds_we_care_about)  # this operates in-place

In ʻevent_dict, it was specified as'auditory / left' and ʻauditory / right, but by passing it to the Epochs object, if ʻauditory` is specified, it will point to both of them. Here, we divide it into data for auditory stimuli and data for visual stimuli.

aud_epochs = epochs['auditory']
vis_epochs = epochs['visual']
del raw, epochs  # free up memory

The channel is specified by the auditory data, and the averaging waveform with the same stimulus presentation time is displayed.

aud_epochs.plot_image(picks=['MEG 1332', 'EEG 021'])

# output
136 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
136 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped

As a supplement, since time series waveform data can be acquired in the form of numpy array by the get_data () method of the Epochs object, it is possible to input it to the machine learning model by taking an array for each condition. ..

aud_epochs.get_data().shape

(136, 376, 106) # (Epoch, channel, time)

The same applies to the following frequency processing that returns a numpy array.

Time frequency analysis

Time-frequency analysis is performed to investigate the relationship between time and frequency before and after stimulus presentation in MEG / EEG waveforms. Here, the response to auditory stimuli is calculated by wavelet analysis and plotted. As shown on the color bar, the darker the red, the greater the amplitude of that frequency at that time. (More detailed example)

frequencies = np.arange(7, 30, 3)
power = mne.time_frequency.tfr_morlet(aud_epochs, n_cycles=2, return_itc=False,
                                      freqs=frequencies, decim=3)
power.plot(['MEG 1332'])

# output
No baseline correction applied

Estimate of evoked response

The Evoked object is the data type of the MNE obtained by averaging the Epochs objects. MEG / EEG data has a lot of noise other than brain activity, and the difference between conditions is often investigated by performing averaging.

Here, we plot the difference between the added waveforms of the auditory data and the visual data.

aud_evoked = aud_epochs.average()
vis_evoked = vis_epochs.average()

mne.viz.plot_compare_evokeds(dict(auditory=aud_evoked, visual=vis_evoked),
                             legend='upper left', show_sensors='upper right')

Evoked objects have methods for further analysis. If you display the waveform for each channel, you can see the runout (P200) that becomes negative in 100ms (N100) in the front channel and has a positive peak around 200ms.

aud_evoked.plot_joint(picks='eeg')

You can also display the distribution on the scalp called topography. Here, it can be seen that brain waves with positive potentials for red and negative potentials for blue are observed on the scalp as follows.

aud_evoked.plot_topomap(times=[0., 0.08, 0.1, 0.12, 0.2], ch_type='eeg')

Also, to see the difference between the conditions, mne.combine_evoked () creates an Evoked object with the data of the two condition differences (at this time, one object is given a minus sign). If you view the waveform of this object, you can see the waveform of the difference between the two conditions.

evoked_diff = mne.combine_evoked([aud_evoked, -vis_evoked], weights='equal')
evoked_diff.pick_types('mag').plot_topo(color='r', legend=False)

in conclusion

Only the basic functions are introduced here, but there are many Tutorial and Example On the MNE official page. tools / stable / auto_examples / index.html) is available, so please refer to it as the first step for further detailed analysis.

Also, since Python has many libraries for statistics and machine learning, I think it will be easier to analyze by MEG / EEG machine learning.

EEG analysis in Python: Python MNE tutorial