Computing ERFs from HCPΒΆ

In this tutorial we compare different ways of arriving at event related fields (ERF) starting from different HCP outputs. We will first reprocess the HCP dat from scratch, then read the preprocessed epochs, finally read the ERF files. Subsequently we will compare these outputs.

# Author: Denis A. Enegemann
# License: BSD 3 clause

import os.path as op

import numpy as np
import matplotlib.pyplot as plt
import mne
import hcp
import hcp.preprocessing as preproc

mne.set_log_level('WARNING')

# we assume our data is inside its designated folder under $HOME
storage_dir = op.expanduser('~')
hcp_params = dict(
    hcp_path=op.join(storage_dir, 'mne-hcp-data', 'HCP'),
    subject='105923',
    data_type='task_working_memory')

We first reprocess the data from scratch

That is, almost from scratch. We’re relying on the ICA solutions and data annotations.

In order to arrive at the final ERF we need to pool over two runs. for each run we need to read the raw data, all annotations, apply the reference sensor compensation, the ICA, bandpass filter, baseline correction and decimation (downsampling)

# these values are looked up from the HCP manual
tmin, tmax = -1.5, 2.5
decim = 4
event_id = dict(face=1)
baseline = (-0.5, 0)

# we first collect events
trial_infos = list()
for run_index in [0, 1]:
    hcp_params['run_index'] = run_index
    trial_info = hcp.read_trial_info(**hcp_params)
    trial_infos.append(trial_info)


# trial_info is a dict
# it contains a 'comments' vector that maps on the columns of 'codes'
# 'codes is a matrix with its length corresponding to the number of trials
print(trial_info['stim']['comments'][:10])  # which column?
print(set(trial_info['stim']['codes'][:, 3]))  # check values

# so according to this we need to use the column 7 (index 6)
# for the time sample and column 4 (index 3) to get the image types
# with this information we can construct our event vectors

all_events = list()
for trial_info in trial_infos:
    events = np.c_[
        trial_info['stim']['codes'][:, 6] - 1,  # time sample
        np.zeros(len(trial_info['stim']['codes'])),
        trial_info['stim']['codes'][:, 3]  # event codes
    ].astype(int)

    # for some reason in the HCP data the time events may not always be unique
    unique_subset = np.nonzero(np.r_[1, np.diff(events[:, 0])])[0]
    events = events[unique_subset]  # use diff to find first unique events

    all_events.append(events)

# now we can go ahead
evokeds = list()
for run_index, events in zip([0, 1], all_events):

    hcp_params['run_index'] = run_index

    raw = hcp.read_raw(**hcp_params)
    raw.load_data()
    # apply ref channel correction and drop ref channels
    preproc.apply_ref_correction(raw)

    annots = hcp.read_annot(**hcp_params)
    # construct MNE annotations
    bad_seg = (annots['segments']['all']) / raw.info['sfreq']
    annotations = mne.Annotations(
        bad_seg[:, 0], (bad_seg[:, 1] - bad_seg[:, 0]),
        description='bad')

    raw.annotations = annotations
    raw.info['bads'].extend(annots['channels']['all'])
    raw.pick_types(meg=True, ref_meg=False)

    # Note: MNE complains on Python 2.7
    raw.filter(0.50, None, method='iir',
               iir_params=dict(order=4, ftype='butter'), n_jobs=1)
    raw.filter(None, 60, method='iir',
               iir_params=dict(order=4, ftype='butter'), n_jobs=1)

    # read ICA and remove EOG ECG
    # note that the HCP ICA assumes that bad channels have already been removed
    ica_mat = hcp.read_ica(**hcp_params)

    # We will select the brain ICs only
    exclude = [ii for ii in range(annots['ica']['total_ic_number'][0])
               if ii not in annots['ica']['brain_ic_vs']]
    preproc.apply_ica_hcp(raw, ica_mat=ica_mat, exclude=exclude)

    # now we can epoch
    events = np.sort(events, 0)
    epochs = mne.Epochs(raw, events=events[events[:, 2] == 1],
                        event_id=event_id, tmin=tmin, tmax=tmax,
                        reject=None, baseline=baseline, decim=decim,
                        preload=True)

    evoked = epochs.average()
    # now we need to add back out channels for comparison across runs.
    evoked = preproc.interpolate_missing(evoked, **hcp_params)
    evokeds.append(evoked)
    del epochs, raw

Out:

[u'1. Run Number' u'2. Block Number within run'
 u'3. Nan.  ( This column has been reserved to contain the image ID number which is not encoded in the trigger values. This is not yet implemented.)'
 u'4. imgType : 1- Face, 2- Tools  0- Fixation'
 u'5. memoryType :  1: 0-Back   2: 2-Back'
 u'6. targetType : 1- target, 2- nontarget,  3: lure '
 u'7. Trial trigger onset Sample ' u'8. Trial trigger offset Sample'
 u'9. Sequence of image in the block'
 u'10. isPressed :  0- user did not press any response button, 1- user pressed a response button']
set([0.0, 1.0, 2.0])
The default output type is "ba" in 0.13 but will change to "sos" in 0.14
The default output type is "ba" in 0.13 but will change to "sos" in 0.14
The default output type is "ba" in 0.13 but will change to "sos" in 0.14
The default output type is "ba" in 0.13 but will change to "sos" in 0.14

Now we can compute the same ERF based on the preprocessed epochs

These are obtained from the ‘tmegpreproc’ pipeline. Things are pythonized and simplified however, so

evokeds_from_epochs_hcp = list()

for run_index, events in zip([0, 1], all_events):
    hcp_params['run_index'] = run_index

    epochs_hcp = hcp.read_epochs(**hcp_params)
    # for some reason in the HCP data the time events may not always be unique
    unique_subset = np.nonzero(np.r_[1, np.diff(events[:, 0])])[0]
    evoked = epochs_hcp[unique_subset][events[:, 2] == 1].average()

    del epochs_hcp
    # These epochs have different channels.
    # We use a designated function to re-apply the channels and interpolate
    # them.

    evoked.baseline = baseline
    evoked.apply_baseline()
    evoked = preproc.interpolate_missing(evoked, **hcp_params)

    evokeds_from_epochs_hcp.append(evoked)

Finally we can read the actual official ERF file

These are obtained from the ‘eravg’ pipelines. We read the matlab file, MNE-HCP is doing some conversions, and then we search our condition of interest. Here we’re looking at the image as onset. and we want the average, not the standard deviation.

evoked_hcp = None
del hcp_params['run_index']
hcp_evokeds = hcp.read_evokeds(onset='stim', **hcp_params)

for ev in hcp_evokeds:
    if not ev.comment == 'Wrkmem_LM-TIM-face_BT-diff_MODE-mag':
        continue

# Once more we add and interpolate missing channels
evoked_hcp = preproc.interpolate_missing(ev, **hcp_params)

Time to compare the outputs

evoked = mne.combine_evoked(evokeds, weights='equal')
evoked_from_epochs_hcp = mne.combine_evoked(
    evokeds_from_epochs_hcp, weights='equal')

fig1, axes = plt.subplots(3, 1, figsize=(12, 8))

evoked.plot(axes=axes[0], show=False)
axes[0].set_title('MNE-HCP')

evoked_from_epochs_hcp.plot(axes=axes[1], show=False)
axes[1].set_title('HCP epochs')

evoked_hcp.plot(axes=axes[2], show=False)
axes[2].set_title('HCP evoked')
fig1.canvas.draw()

plt.show()

# now some correlations

plt.figure()
r1 = np.corrcoef(evoked_from_epochs_hcp.data.ravel(),
                 evoked_hcp.data.ravel())[0][1]
plt.plot(evoked_from_epochs_hcp.data.ravel()[::10] * 1e15,
         evoked_hcp.data.ravel()[::10] * 1e15,
         linestyle='None', marker='o', alpha=0.1,
         mec='orange', color='orange')
plt.annotate("r=%0.3f" % r1, xy=(-300, 250))
plt.ylabel('evoked from HCP epochs')
plt.xlabel('evoked from HCP evoked')
plt.show()

plt.figure()
r1 = np.corrcoef(evoked.data.ravel(), evoked_hcp.data.ravel())[0][1]
plt.plot(evoked.data.ravel()[::10] * 1e15,
         evoked_hcp.data.ravel()[::10] * 1e15,
         linestyle='None', marker='o', alpha=0.1,
         mec='orange', color='orange')
plt.annotate("r=%0.3f" % r1, xy=(-300, 250))
plt.ylabel('evoked from scratch with MNE-HCP')
plt.xlabel('evoked from HCP evoked file')
plt.show()
  • ../_images/sphx_glr_plot_reproduce_erf_001.png
  • ../_images/sphx_glr_plot_reproduce_erf_002.png
  • ../_images/sphx_glr_plot_reproduce_erf_003.png

Total running time of the script: ( 2 minutes 38.477 seconds)

Generated by Sphinx-Gallery