Seminar 1. EEG analysis

# For Colab only
!pip install mne

import mne
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns

%matplotlib notebook

!wget "https://www.physionet.org/files/eegmmidb/1.0.0/S003/S003R03.edf"

!wget "https://www.physionet.org/files/eegmmidb/1.0.0/S003/S003R03.edf.event"

!ls

sample = mne.io.read_raw_edf('S003R03.edf', verbose=False, preload=True)

sample.info

# Sampling frequency
sample.info['sfreq']

# Length in seconds
len(sample) / sample.info['sfreq']

# Number of channels
len(sample.ch_names)

sample.ch_names[:3]

# fix trailing dots in channel names
# use sample.rename_channels(map)

# YOUR CODE HERE

sample.ch_names[:3]

# 19 channels from 10-20. no A1 and A2 here
# Be careful. Pure 10-20 labeling differs from high-resolution montages
# In MNE, 10-20 montage is actually an extended high-resulution version of 10-20
# FYI, mapping from pure 10-20 to high-resolution versions
# T3 = T7
# T4 = T8
# T5 = P7
# T6 = P8

channels_to_use = [
    # prefrontal
    'Fp1',
    'Fp2',
    # frontal
    'F7',
    'F3',
    'F4',
    'Fz',
    'F8',
    # central and temporal
    'T7',
    'C3',
    'Cz',
    'C4',
    'T8',
    # parietal
    'P7',
    'P3',
    'Pz',
    'P4',
    'P8',
    # occipital
    'O1',
    'O2',
]

sample_1020 = sample.copy().pick_channels(channels_to_use)

# check that everything is OK
assert len(channels_to_use) == len(sample_1020.ch_names)

ch_map = {ch.lower(): ch for ch in sample_1020.ch_names}

ten_twenty_montage = mne.channels.make_standard_montage('standard_1020')
len(ten_twenty_montage.ch_names)

ten_twenty_montage.ch_names = [ch_map[ch.lower()] if ch.lower() in ch_map else ch 
                               for ch in ten_twenty_montage.ch_names]

sample_1020.set_montage(ten_twenty_montage)

sample_1020.plot_sensors(show_names=True);

sample_1020.plot_psd();

sample_1020.filter(l_freq=1, h_freq=50, method='iir')

# Plot psd after filtering

# YOUR CODE HERE

sample_1020.plot(n_channels=8, duration=20);

# Plot in better scale. Use 'scaling' argument

# YOUR CODE HERE

events, events_dict = mne.events_from_annotations(sample_1020)

events_dict

events[:5]

epochs = mne.Epochs(sample_1020, events,  tmin=-0.5, tmax=0.8, preload=True)

pd.DataFrame(epochs.events, columns=['_', '__', 'event_id'])['event_id'].value_counts()

for epoch in epochs:
    break
epoch.shape

epoch.shape[1] / sample_1020.info['sfreq']

sample_1020.to_data_frame().shape

df = epochs.to_data_frame()
df.head(3).iloc[:, :10]

df[sample_1020.ch_names + ['epoch']].groupby('epoch').agg(lambda arr: arr.max() - arr.min()).hist(figsize=[10, 10]);
plt.tight_layout()

epochs = mne.Epochs(sample_1020, events,  tmin=-0.5, tmax=0.8, reject={'eeg': 600e-6}, preload=True, baseline=(-.1, 0))

epochs.plot_psd();

epochs.plot(n_channels=8, scalings={'eeg':1e-4});

epochs.event_id

# check number of events of each type
# use epochs.events

# Your code here

pd.DataFrame(epochs.events, columns=['_', '__', 'event_id'])['event_id'].value_counts()

evoked_T0 = epochs['1'].average()
evoked_T1 = epochs['2'].average()
evoked_T2 = epochs['3'].average()

evoked_T0.plot(spatial_colors=True);

evoked_T1.plot(spatial_colors=True);

evoked_T2.plot(spatial_colors=True);

ica = mne.preprocessing.ICA(n_components=10, random_state=42)

ica.fit(sample_1020)

ica.plot_sources(sample_1020);

ica.plot_components();

ica.plot_properties(sample_1020, picks=[4]);

ica.plot_overlay(sample_1020, exclude=[0, 1, 4, 5, 8, 9], picks=['F3']);

ica.exclude = [0, 1]

sample_1020_clr = sample_1020.copy()

ica.apply(sample_1020_clr)

# plot channels

# YOUR CODE HERE

epochs = mne.Epochs(sample_1020_clr, events,  tmin=-0.5, tmax=0.8, reject={'eeg': 600e-6}, preload=True, baseline=(-.1, 0))

evoked_T0 = epochs['1'].average()
evoked_T1 = epochs['2'].average()
evoked_T2 = epochs['3'].average()

evoked_T0.plot(spatial_colors=True);

evoked_T1.plot(spatial_colors=True);

evoked_T2.plot(spatial_colors=True);

evoked_T0.plot_topomap(times=[0, .2, .4, .6, .8], vmin=-50, vmax=50);

evoked_T1.plot_topomap(times=[0, .2, .4, .6, .8], vmin=-50, vmax=50);

evoked_T2.plot_topomap(times=[0, .2, .4, .6, .8], vmin=-50, vmax=50);

evoked_T0_alpha = evoked_T0.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T1_alpha = evoked_T1.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T2_alpha = evoked_T2.copy().filter(l_freq=7, h_freq=13, method='iir', verbose=False).apply_hilbert(envelope=True)

evoked_T1_alpha.plot(spatial_colors=True);

evoked_T0_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vmin=0, vmax=30);

evoked_T1_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vmin=0, vmax=30);

evoked_T2_alpha.plot_topomap(times=[0, .1, .2, .3, .4, .6], vmin=0, vmax=30);

evoked_T0_beta_low = evoked_T0.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T1_beta_low = evoked_T1.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)
evoked_T2_beta_low = evoked_T2.copy().filter(l_freq=13, h_freq=20, method='iir', verbose=False).apply_hilbert(envelope=True)

evoked_T0_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vmin=0, vmax=30);

evoked_T1_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vmin=0, vmax=30);

evoked_T2_beta_low.plot_topomap(times=[0, .2, .4, .6, .8], vmin=0, vmax=30);

conn_T1, freqs, times, _, _ = mne.connectivity.spectral_connectivity(epochs['2'], method='coh')

def plot_topomap_connectivity_2d(info, con, picks=None, pairs=None, vmin=None, vmax=None, cm=None, show_values=False, show_names=True):
    """
    Plots connectivity-like data in 2d
    
    Drawing every pair of channels will likely make a mess
    There are two options to avoid it:
    - provide picks
    - provide specific pairs of channels to draw
    """
    
    # get positions
    _, pos, _, _, _, _, _ = mne.viz.topomap._prepare_topomap_plot(info, 'eeg');
    
#     if picks is None and pairs is None:
#         picks = info.ch_names
    
    ch_names_lower = [ch.lower() for ch in info.ch_names]
    if picks:
        picks_lower = [ch.lower() for ch in picks]
    if pairs:
        pairs_lower = [tuple(sorted([ch1.lower(), ch2.lower()])) for ch1, ch2 in pairs]
    
    rows = []
    for idx1, ch1 in enumerate(ch_names_lower):
        for idx2, ch2 in enumerate(ch_names_lower):
            if ch1 >= ch2:
                continue
            if picks and (ch1 not in picks_lower or ch2 not in picks_lower):
                    continue
            if pairs and (ch1, ch2) not in pairs_lower:
                    continue
            rows.append((
                pos[idx1],
                pos[idx2],
                con[idx1, idx2]
            ))
    
    if not len(rows):
        raise ValueError('No pairs to plot')
    
    con_to_plot = np.array([row[2] for row in rows])
    if vmin is None:
        vmin = np.percentile(con_to_plot, 2)
    if vmax is None:
        vmax = np.percentile(con_to_plot, 98)
    norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
    
    if cm is None:
        cm = sns.diverging_palette(240, 10, as_cmap=True)
    
    fig, ax = plt.subplots(figsize=[5, 5])
    mne.viz.utils.plot_sensors(info, show_names=show_names, show=False, axes=ax);
    for row in rows:
        rgba_color = cm(norm(row[2]))
        plt.plot([row[0][0], row[1][0]], [row[0][1], row[1][1]], color=rgba_color)
        if show_values:
            plt.text((row[0][0] + row[1][0]) / 2, 
                     (row[0][1] + row[1][1]) / 2, 
                     '{:.2f}'.format(row[2]))

conn_T0, freqs, times, _, _ = mne.connectivity.spectral_connectivity(epochs['1'], method='coh', verbose=False);
conn_T1, freqs, times, _, _ = mne.connectivity.spectral_connectivity(epochs['2'], method='coh', verbose=False);
conn_T2, freqs, times, _, _ = mne.connectivity.spectral_connectivity(epochs['3'], method='coh', verbose=False);

freqs

conn_T0_beta = conn_T0[:, :, 12:27].mean(axis=2)
conn_T0_beta = conn_T0_beta + conn_T0_beta.T

conn_T1_beta = conn_T1[:, :, 12:27].mean(axis=2)
conn_T1_beta = conn_T1_beta + conn_T1_beta.T

conn_T2_beta = conn_T2[:, :, 12:27].mean(axis=2)
conn_T2_beta = conn_T2_beta + conn_T2_beta.T

plot_topomap_connectivity_2d(epochs.info, conn_T1_beta, picks=epochs.ch_names);

plot_topomap_connectivity_2d(epochs.info, conn_T0_beta, 
                             pairs=[('F7', 'F4'), ('O2', 'T7'), ('C3', 'C4'), ('P7', 'P8'), ('F8', 'T8'), ('O1', 'O2'), ('O1', 'P4')],
                             show_values=True,
                             show_names=False
                            
                            );

# calculate coherence in alpha band
# plot 5-10 pairs that you are interested in

# YOUR CODE HERE