Preproc - Basic source space pipeline

This example script shows how to perform basic preprocessing manually. The main input to this pipeline is an MEEG object that has already had the coregistration and forward model run - refer to the relevant practicals for examples (in particular, the CTF data import tutorial). This pipeline goes from a raw sensorspace coregistered MEEG to a source reconstruction.

Contents

Overview

Having read the raw files into an SPM MEEG object, we are now ready to implement a basic source reconstruction pipeline. The main steps are

although you may not need or want to include all of these steps, depending on the analysis you are running.

First, we will load in the MEEG file. This example is a continuous resting state CTF recording

data_dir = fullfile(osldir,'example_data','preproc_simple_example');
spm_file = fullfile(data_dir,'3006.mat');
D = spm_eeg_load(spm_file);

Note that at any point, we can find an MEEG object on disk by checking the fullfile property, so you don't really need to keep track of the file name separately.

D.fullfile

ans =

    '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006.mat'

Update D.inv paths

One issue you might encounter is that the paths to structural files in D.inv may be different if coregistration was performed on someone else's machine, or even on the same machine but in a different directory.

D.inv{1}.mesh

ans = 

  struct with fields:

       template: 0
           sMRI: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/3006_CRG.nii'
            def: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/y_3006_CRG.nii'
         Affine: [4×4 double]
          Msize: 2
       tess_mni: [1×1 struct]
       tess_ctx: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/3006_CRGcortex_8196.surf.gii'
     tess_scalp: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/3006_CRGscalp_2562.surf.gii'
    tess_oskull: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/3006_CRGoskull_2562.surf.gii'
    tess_iskull: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_Data/ctf_preprocessing/3006/3006_CRGiskull_2562.surf.gii'
            fid: [1×1 struct]
     tess_rhino: [1×1 struct]

Notice how the structural files for this example are stored in the preproc_simple_example folder, but this is not the location contained in the MEEG object. This can be a problem if you want to display the coregistration e.g. with rhino_display. You can update these paths using osl_update_inv_dir(). The inputs to this file are the MEEG object, and the new path of the folder containing the structural files. This assumes that all of the files share a common root - this is generally a safe assumption if the files were created with RHINO.

D = osl_update_inv_dir(D,data_dir);
D.inv{1}.mesh

ans = 

  struct with fields:

       template: 0
           sMRI: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006_CRG.nii'
            def: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/y_3006_CRG.nii'
         Affine: [4×4 double]
          Msize: 2
       tess_mni: [1×1 struct]
       tess_ctx: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006_CRGcortex_8196.surf.gii'
     tess_scalp: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006_CRGscalp_2562.surf.gii'
    tess_oskull: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006_CRGoskull_2562.surf.gii'
    tess_iskull: '/Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/3006_CRGiskull_2562.surf.gii'
            fid: [1×1 struct]
     tess_rhino: [1×1 struct]

Setting chantypes and labels

It's important that you set the channel types and labels correctly in your file, for two reasons

Setting the labels and chantype correctly makes it easier for you and other users to work with the data further down the line, so it is always worth setting them correctly. In this case, suppose we know that the following channel labels correspond to artefact channels:

At this point, we should make sure that these channel types are set correctly, and that the channel labels are informative

D = D.chantype(find(strcmp(D.chanlabels,'EEG060')),'EMG');
D = D.chantype(find(strcmp(D.chanlabels,'EEG059')),'ECG');
D = D.chantype(find(strcmp(D.chanlabels,'EEG057')),'EOG');
D = D.chantype(find(strcmp(D.chanlabels,'EEG058')),'EOG');
D = D.chanlabels(find(strcmp(D.chanlabels,'EEG060')),'EMG');
D = D.chanlabels(find(strcmp(D.chanlabels,'EEG059')),'ECG');
D = D.chanlabels(find(strcmp(D.chanlabels,'EEG057')),'EOG1');
D = D.chanlabels(find(strcmp(D.chanlabels,'EEG058')),'EOG2');

Initial filtering

For filtering, we will use the osl_filter function. This function applies a basic Butterworth filter. You specify both the filter type and key frequencies using two numbers - a lower and upper frequency [f1,f2]. The possibilities are

We would typically apply a high pass filter to remove very slow drift, and then notch filters to remove AC line artefacts and their harmonics (if the sampling frequency of your data is high enough that they may be present)

D = osl_filter(D,[0.1 inf]); % Remove slow drift
D = osl_filter(D,-1*(50+[-2 2])); % Remove 50Hz with notch filter
D = osl_filter(D,-1*(2*100+[-2 2])); % Remove 100Hz with notch filter

SPM12: spm_eeg_filter (v5876)                      16:24:52 - 21/11/2017
========================================================================

SPM12: spm_eeg_filter (v5876)                      16:24:58 - 21/11/2017
========================================================================

SPM12: spm_eeg_filter (v5876)                      16:25:05 - 21/11/2017
========================================================================

Downsampling

The raw files may have an unnecessarily high sampling rate for the analysis you wish to perform, which increases storage requirements and processing time. You can downsample your data at this point if you wish

D = spm_eeg_downsample(struct('D',D,'fsample_new',300));

SPM12: spm_eeg_downsample (v6614)                  16:25:11 - 21/11/2017
========================================================================
Resampling frequency is 300Hz

Remove artefacts

You can remove artefacts using osl_detect_artefacts. Broadly, there are two kinds of artefacts this function identifies and classifies

It can be very important to run osl_detect_artefacts or to otherwise perform artefact rejection before doing ICA, as the presence of known artefacts can significantly degrade the ICA decomposition.

D = osl_detect_artefacts(D);

Detecting artefacts in channel types: ECG,EMG,EOG,MEGGRAD
Bad times - rejected 6.00s (2%) in modality EMG
Bad times - rejected 16.00s (5%) in modality EOG
Bad times - rejected 8.00s (3%) in modality MEGGRAD

By default, osl_detect_artefacts will identify both types of artefacts, although you can optionally specify which ones you want to detect. You can also choose different artefact detection metrics and thresholds, although the function is intended to use sensible defaults. See the artefact detection example for more detailed usage of osl_detect_artefacts. function.

Aside from examining the properties of the MEEG object, it can also be helpful to examine the bad epochs visually. You can do this by opening the MEEG object in oslview

oslview(D)

AFRICA - ICA artefact removal

Most of the information on how to perform ICA artefact removal is provided in the AFRICA example. Here, we are mainly concerned with how to perform automatic artefact rejection. We set artefact_channels to correspond to the chantypes of all the channels whose correlations we want to check. precompute_topos is slow and generates large files, and is only mainly necessary if you are performing manual artefact rejection. So we can save time by skipping that step here. We set the ident_func to the automatic classification function. Lastly, if you examine identify_artefactual_components_auto.m you can see what options are available to control the rejection process. One common setting you might want to change is the correlation threshold for rejection. You might also want to enable or disable other types of ICA component rejection. Finally, call osl_africa

D = osl_africa(D,'artefact_channels',{'EOG1','ECG','EMG','EOG2'},'precompute_topos',false);
has_montage(D)

Number of signals: 274
Number of samples: 87600
Calculating covariance...
Reducing dimension...
Selected [ 150 ] dimensions.
Smallest remaining (non-zero) eigenvalue [ 1.77653 ]
Largest remaining (non-zero) eigenvalue [ 2870.4 ]
Sum of removed eigenvalues [ 173.911 ]
[ 98.6622 ] % of (non-zero) eigenvalues retained.
Whitening...
Check: covariance differs from identity by [ 7.95204e-14 ].
Used approach [ symm ].
Used nonlinearity [ tanh ].
Using stabilized algorithm.
Starting ICA calculation...
Step no. 1
Step no. 2, change in value of estimate: 0.99 
...
...
...
Step no. 585, change in value of estimate: 0.000101 
Convergence after 586 steps
Adding the mean back to the data.
** Saving changes to disk **
Requested to check correlations with chantype EOG1, but no channels have this type in the data
Requested to check correlations with chantype EOG2, but no channels have this type in the data
No ICs correlated with Channel 277_EMG
Rejecting IC 19 due to correlation with Channel 278_ECG (correlation = 0.58)
Rejecting IC 52 due to correlation with Channel 278_ECG (correlation = 0.57)
Rejecting IC 112 due to correlation with Channel 278_ECG (correlation = 0.56)
0 - none (280 channels)
*1 - AFRICA denoised data (280 channels)

As you can see, a new online montage is created corresponding to the ICA-denoised signals.

Forward model

You can check which forward model was used by looking at the contents of D.inv

D.inv{1}.forward.voltype

ans =

    'MEG Local Spheres'

If you want to change the forward model, you can do that now

D = osl_forward_model(D,'forward_meg','Single Shell');
D.inv{1}.forward.voltype

the input is mesh data with 2562 vertices and 5120 triangles

ans =

    'Single Shell'

Beamforming

Finally, you can run the beamformer. You may want to use a filter at this point as well - 1-45Hz is typical.

D = osl_filter(D,[1 45]);
p = parcellation(8);
D = osl_inverse_model(D,p.template_coordinates,struct('pca_order',150));

SPM12: spm_eeg_filter (v5876)                      16:30:25 - 21/11/2017
========================================================================
BF working directory: /Users/romesh/oxford_postdoc/toolboxes/osl/example_data/preproc_simple_example/osl_bf_temp_76c804bd
Initialising batch system... done.


------------------------------------------------------------------------
Running job #1
------------------------------------------------------------------------
Running 'Prepare data'
Done    'Prepare data'
Running 'Define sources'
computing surface normals
Done    'Define sources'
Running 'Covariance features'
Done    'Covariance features'
Running 'Inverse solution'
Done    'Inverse solution'
Running 'Output'
Done    'Output'
Running 'Write'
An online montage is already active, applying beamformer tra to it

SPM12: spm_eeg_montage (v6835)                     16:32:05 - 21/11/2017
========================================================================
An online montage is already active, applying beamformer tra to it

SPM12: spm_eeg_montage (v6835)                     16:32:07 - 21/11/2017
========================================================================
Done    'Write'
Done

Notice that there are now online montages associated with the source reconstruction

has_montage(D)

0 - none (280 channels)
1 - AFRICA denoised data (280 channels)
2 - without weights normalisation, class 1 (3559 channels)
*3 - with weights normalisation, class 1 (3559 channels)

For more information about beamforming, refer to the beamforming example on this website. Now that the MEEG object has a beamformed montage, it is ready to be used in source space analysis e.g. source space connectivity analysis.

Published with MATLAB® R2017a