.. _sphx_glr_auto_tutorials_plot_reproduce_erf.py: .. _tut_reproduce_erf: ======================= 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. .. code-block:: python # 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) .. code-block:: python # 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 .. rst-class:: sphx-glr-script-out 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 .. code-block:: python 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. .. code-block:: python 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 .. code-block:: python 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() .. rst-class:: sphx-glr-horizontal * .. image:: /auto_tutorials/images/sphx_glr_plot_reproduce_erf_001.png :scale: 47 * .. image:: /auto_tutorials/images/sphx_glr_plot_reproduce_erf_002.png :scale: 47 * .. image:: /auto_tutorials/images/sphx_glr_plot_reproduce_erf_003.png :scale: 47 **Total running time of the script:** ( 2 minutes 38.477 seconds) .. container:: sphx-glr-download :download:`Download Python source code: plot_reproduce_erf.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: plot_reproduce_erf.ipynb ` .. rst-class:: sphx-glr-signature `Generated by Sphinx-Gallery `_