Analyse des ondes cérébrales avec Python: tutoriel Python MNE

Python MNE est un outil open source pour analyser et visualiser des diagrammes magnétiques cérébraux (MEG) et des ondes cérébrales (EEG). On peut dire qu'elle est très polyvalente car elle peut être appliquée aux formats de données de nombreux appareils. Dans cet article, suivez Most Basic Tutorials , La procédure d'analyse MEG et EEG MNE est expliquée.

Environnement d'exécution

Installation

Vous pouvez l'installer avec Anaconda ou pip. Sur la page officielle, Anaconda est recommandé.


conda install -c conda-forge mne


pip install -U mne

De plus, les paramètres GPU semblent être possibles, donc cela prend du temps pour le filtrage et l'échantillonnage. Si cela prend, il vaut peut-être mieux le régler.


$ 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'donner

Tutoriel d'analyse MEG / EEG

Typique de l'analyse des ondes cérébrales selon ce tutoriel J'expliquerai le traitement de base et la structure des données de base du MNE. Je vais l'expliquer de manière concise, donc si vous avez des connaissances en anglais et en analyse des ondes cérébrales, vous devriez lire la page originale.

Procédure d'analyse

Importation de module

import os
import numpy as np
import mne

Chargement des données

Chargez les exemples de données fournis dans mne (Liste de tous les ensembles de données). À propos, MEG et EEG ont plusieurs formats de fichiers pour chaque appareil, mais il existe une fonction pour charger ces données (Formats de données pris en charge. 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)

Cette fois, nous utilisons les données avec un filtre passe-bas à 40Hz appelé sample_audvis_filt-0-40_raw.fif. Le téléchargement des données a pris environ 1 minute (1,54 Go). download.png

La fonction mne.io.read_raw_fif () lit les données de la structure de données MNE appelée objet Raw à partir du fichier au format fif.

À propos, dans l'expérience de ces données, le sujet est présenté avec un stimulus visuel et un stimulus auditif au hasard à gauche et à droite toutes les 750 ms. Aussi, parfois, un visage riant est affiché à l'écran et le sujet appuie sur un bouton lorsqu'il le voit. La présentation de ces stimuli s'appelle un événement.

Affichez l'objet chargé.

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
>

Pour MNE, Raw, Epochs /mne.Epochs.html#mne.Epochs), Evoked Les données sont gérées par des objets appelés objets. Ceux-ci se distinguent par la division temporelle des données, et Raw inclut les données de tout le temps depuis le début jusqu'à l'arrêt de l'enregistrement des ondes cérébrales, et il est bon de penser que c'est l'état avant d'être séparé par l'événement. Je pense. raw.info contient des informations telles que le nombre de canaux, les noms et les taux d'échantillonnage de raw.

La densité spectrale de puissance (PSD) de cet objet brut et la forme d'onde de chaque capteur sont affichées. plot () peut être déroulé et mis à l'échelle dans un shell interactif.

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

Prétraitement

Dans le traitement du bruit de l'analyse EEG / MEG, l'analyse standard des composants indépendants (ICA) est effectuée. En termes simples, en décomposant un signal multicanal en composants indépendants, un traitement qui supprime les éléments indépendants de l'activité cérébrale, tels que le bruit de l'alimentation et le clignotement, est effectué.

Ici, pour le signal de 364 canaux, après décomposition en 20 composantes indépendantes, «[1, 2]» et l'index du composant à supprimer sont spécifiés (le retrait réel est le suivant ʻica.apply). Cela se fait par (raw) ). Affichez les informations des composants indépendants avec ʻ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
ica_1.png ica_2.png

Appliquez l'objet ICA à l'objet Raw pour obtenir le signal avec le composant indépendant spécifié supprimé.

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

Lorsque le canal avant est affiché, le bruit est supprimé comme indiqué dans les deux figures suivantes.

# 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
orig_raw.png prune_raw.png

Détection d'événement

Ces exemples de données incluent ce que l'on appelle le canal STIM, qui contient des données sur le type d'événement et son heure. Cette fois, le canal nommé «STI 014» est applicable, alors spécifiez-le et récupérez l'événement de l'expérience.

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]]

ʻEvents` est un tableau numpy, la première colonne représente le temps d'échantillonnage et la dernière colonne représente le type d'événement (ID) (0 au milieu peut être ignoré). Ici, il semble y avoir un événement avec l'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

L'ID d'événement et la description de l'événement sont stockés dans un type de dictionnaire. À propos, la clé ici peut être n'importe quelle chaîne de caractères.

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

Tracons le timing de présentation de chaque événement sur l'axe des temps. L'ID défini précédemment par ʻevent_id = even_dict` correspond au nom de l'événement sur le tracé.

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

Découpé dans une époque

Dans l'analyse MEG / EEG, en particulier lors de l'étude de la réponse à un certain stimulus, le signal est divisé en unités appelées époques et traité comme une seule donnée. En découpant une certaine section avant et après le stimulus, il est transformé en une unité de données. Ce processus de découpe est appelé époque. Les objets bruts sont connectés pendant toute la durée de l'enregistrement, mais seule la partie de la réponse au stimulus est découpée et stockée dans le type de données appelé objet Epochs.

À l'époque, un processus appelé rejet de seuil est souvent effectué. Cela signifie que s'il y a un signal trop fort pendant le temps concerné, cette partie aura été très bruyante et sera exclue des données. Déterminez les critères d'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

Ensuite, créez un objet Epochs à partir de l'objet Raw. Les arguments de mne.Epochs ici sont les suivants.

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

Il semble que 10 époques aient été exclues cette fois.

Le moment d'appuyer sur le bouton du sujet n'étant pas traité cette fois, seul l'événement de stimulation audiovisuelle est ciblé. De plus, afin d'éviter que le nombre de données ne soit biaisé entre les conditions, le nombre de chaque condition est uniformisé en supprimant les événements d'un grand nombre de conditions avec ʻ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

Dans ʻevent_dict, il a été spécifié comme'auditory / left' et ʻauditory / right, mais en le passant à l'objet Epochs, si ʻauditory` est spécifié, il pointera vers les deux. Ici, nous le divisons en données pour les stimuli auditifs et en données pour les stimuli visuels.

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

Un canal est spécifié par des données auditives, et une forme d'onde de moyenne avec le même temps de présentation du stimulus est affichée.

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
meg_aud.png meg_viz.png

En complément, étant donné que les données de forme d'onde de séries temporelles peuvent être acquises sous forme de tableau numpy par la méthode get_data () de l'objet Epochs, il est possible de les saisir dans le modèle d'apprentissage automatique en prenant un tableau pour chaque condition. ..

aud_epochs.get_data().shape

(136, 376, 106) # (Époque, canal, heure)

La même chose s'applique au traitement de fréquence suivant qui renvoie un tableau numpy.

Analyse temps-fréquence

Une analyse temps-fréquence est effectuée pour étudier la relation entre le temps et la fréquence avant et après la présentation du stimulus dans la forme d'onde MEG / EEG. Ici, la réponse aux stimuli auditifs est calculée par analyse d'ondelettes et tracée. Comme indiqué dans la barre de couleur, plus le rouge est foncé, plus l'amplitude de cette fréquence est élevée à ce moment-là. (Exemple plus détaillé)

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
tf.png

Estimation de la réponse évoquée

En tant que type de données MNE, l'objet évoqué est obtenu en calculant la moyenne de l'objet Epochs. Les données MEG / EEG ont beaucoup de bruit autre que l'activité cérébrale, et la différence entre les conditions est souvent étudiée en effectuant une moyenne d'addition.

Ici, la différence entre les formes d'onde ajoutées des données auditives et les données visuelles est tracée.

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')
meg_evoked.png eeg_evoked.png grad_evoked.png

Les objets évoqués ont des méthodes pour une analyse plus approfondie. Si vous affichez la forme d'onde pour chaque canal, vous pouvez voir le voile (P200) qui devient négatif en 100 ms (N100) dans le canal avant et a un pic positif autour de 200 ms.

aud_evoked.plot_joint(picks='eeg')
aud_eeg_evoked.png

Vous pouvez également afficher la distribution sur le cuir chevelu appelée topographie. Ici, on peut voir que des ondes cérébrales avec des potentiels positifs pour le rouge et des potentiels négatifs pour le bleu sont observées sur le cuir chevelu comme suit.

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

Aussi, pour voir la différence entre les conditions, mne.combine_evoked () crée un objet Evoked avec les données des deux différences de condition (à ce moment, un objet reçoit un signe moins). Si vous visualisez la forme d'onde de cet objet, vous pouvez voir la forme d'onde de la différence entre les deux conditions.

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

en conclusion

Seules les fonctions de base ont été introduites ici, mais il existe de nombreux Tutoriel et [Exemple](https: // mne.) Sur la page officielle du MNE. tools / stable / auto_examples / index.html) est disponible, veuillez donc vous y référer comme première étape pour une analyse plus détaillée.

De plus, étant donné que Python dispose de nombreuses bibliothèques pour les statistiques et l'apprentissage automatique, je pense qu'il sera plus facile à analyser par l'apprentissage automatique MEG / EEG.

Recommended Posts

Analyse des ondes cérébrales avec Python: tutoriel Python MNE
Analyse d'association en Python
Analyse de régression avec Python
Analyse des contraintes symétriques axiales avec Python
Analyse de régression simple avec Python
Première analyse de régression simple en Python
Tutoriel Python
Analyse du squelette planaire dans Python (2) Hotfix
[Didacticiel d'analyse Python en base de données avec SQL Server 2017]
Tutoriel de recommandation utilisant l'analyse d'association (implémentation python)
Analyse résiduelle en Python (Supplément: règles Cochrane)
Quadtree en Python --2
Python en optimisation
CURL en Python
Métaprogrammation avec Python
Python 3.3 avec Anaconda
Géocodage en python
SendKeys en Python
Méta-analyse en Python
Tutoriel Python Django (5)
Tutoriel Python Django (2)
Unittest en Python
Résumé du didacticiel Python
Analyse de données python
Époque en Python
Discord en Python
Allemand en Python
DCI en Python
tri rapide en python
nCr en python
N-Gram en Python
Programmation avec Python
Tutoriel Python Django (8)
Tutoriel Python Django (6)
Plink en Python
Constante en Python
FizzBuzz en Python
Sqlite en Python
Étape AIC en Python
LINE-Bot [0] en Python
CSV en Python
Assemblage inversé avec Python
Réflexion en Python
Constante en Python
nCr en Python.
format en python
Scons en Python 3
Tutoriel Python Django (7)
Tutoriel Python Django (1)
Puyopuyo en python
python dans virtualenv
PPAP en Python
Tutoriel du didacticiel Python Django
Quad-tree en Python
Réflexion en Python
Chimie avec Python
Hashable en Python
DirectLiNGAM en Python
Tutoriel Python Django (3)