Please read this about paths and files before you proceed!

Since the workshop was not intended to cover the preprocessing steps, the raw files were not actually provided.
Therefore, the preprocessing tutorials rather serve the purpose of showing what was done to the preprocessed data.
They of course also serve as examples for you to model your own analysis if deemed useful

Set up general paths

Change these to appropriate paths for your operating system and setup

addpath /home/lau/matlab/fieldtrip-20170613_workshop_natmeg/
ft_defaults

raw_meg_path = '/archive/20067_workshop_source_reconstruction/MEG/';
meg_path = '/home/share/workshop_source_reconstruction/data/MEG/';
mri_path = '/home/share/workshop_source_reconstruction/data/MRI/';

Set up subject specific paths

Make subject and recording specific paths (the cell array “subjects_and_dates” can be expanded)

%% subjects and dates

subjects_and_dates = ...
                    {
                        'NatMEG_0177/170424/'
                    };
                
filenames = {...
        'tactile_stim_raw_tsss_mc.fif'
        'tactile_stim_raw_tsss_mc-1.fif'
        'tactile_stim_raw_tsss_mc-2.fif'
            };
        
n_filenames = length(filenames);  

output_path = fullfile(meg_path, subjects_and_dates{1});
events = [1 2 4 8 16]; %% right little, right ring, right middle, right index, right thumb

Set preprocessing options

  1. Elekta data files (.fif) are split into several files due to a limitation of 2 GB for each files
  2. Event codes are read in (1, 2, 4, 8, 16)
  3. We anticipate the measured delay of 41 ms between the trigger and the actual stimulation
  4. MEG data and EEG data are handled separately since EEG data need to be referenced to a reference “channel” (here the average of all channels)
  5. Not much preprocessing is done to the data (we’ll do that later for the respective analyses (timelocked and time frequency representation (TFR)))
  6. The split files are appended to one another to ease data handling
%% read in data

split_files_MEG = cell(1, n_filenames);
split_files_EEG = split_files_MEG;

for filename_index = 1:n_filenames

    filename = filenames{filename_index};
    full_path = fullfile(raw_meg_path, subjects_and_dates{1}, filename);
 
    % define trials
    
    cfg = [];
    cfg.dataset = full_path;
    cfg.trialdef.prestim = 1.609; % seconds % adjusted with 41 ms
    cfg.trialdef.poststim = 1.691; % seconds
    cfg.trialdef.eventvalue = events;
    cfg.trialdef.eventtype = 'STI101';
    
    cfg = ft_definetrial(cfg);
    
    % preprocess
 
    cfg.demean = 'no';
    cfg.lpfilter = 'no';
    cfg.hpfilter = 'no';
    cfg.dftfilter = 'no';
    cfg.channel = 'MEG';

    split_files_MEG{filename_index} = ft_preprocessing(cfg);
    
    cfg.channel = 'EEG';
    cfg.reref = 'no';
    cfg.refchannel = 'all';
    
    split_files_EEG{filename_index} = ft_preprocessing(cfg);
         
end

% append split files

cfg = [];

preprocessed_data_MEG = ft_appenddata(cfg, split_files_MEG{:});
preprocessed_data_EEG = ft_appenddata(cfg, split_files_EEG{:});

Remove bad electrodes

These would mess up the source reconstructions if left in

%% remove bad electrodes

% Find bad channels (#27+34+40+79+78+96) % numbers of bad channels found
cfg = [];
cfg.method      = 'summary';
cfg.keepchannel = 'no';
cfg.channel     = 'EEG';
cfg.layout      = 'neuromag306eeg1005_natmeg.lay';

temp = ft_rejectvisual(cfg, preprocessed_data_EEG);

Interpolate bad channels and re-reference data

  1. The data for the removed channels is interpolated back in based on their neighbours
  2. The data is rereferenced to the average of all the channels
%% Interpolate bad channels
% Get the names of the selected channels
% badchannels = setdiff(preprocessed_data_EEG.label, temp.label);
badchannels = union(preprocessed_data_MEG.label, {'EEG027' 'EEG034' 'EEG040' 'EEG079' 'EEG078' 'EEG096'});

% Interpolate bad channels: Find neighbours
cfg             = [];
cfg.method      = 'triangulation';
cfg.senstype   = 'EEG'; % 
neighbours_EEG = ft_prepare_neighbours(cfg, preprocessed_data_EEG);
 
% plotting neighbours for inspection
cfg            = [];
cfg.neighbours = neighbours_EEG;
cfg.senstype   = 'EEG';
ft_neighbourplot(cfg, preprocessed_data_EEG);

% Interpolate channels
cfg = [];
cfg.method                    = 'spline';
cfg.neighbours                = neighbours_EEG;
cfg.badchannel                = badchannels;
cfg.senstype                  = 'EEG';
interpolated_data = ft_channelrepair(cfg, preprocessed_data_EEG);
 
% Rereference to common average
cfg = [];
cfg.reref                  = 'yes';
cfg.refchannel             = 'all';
 
rereferenced_interpolated_data = ft_preprocessing(cfg, interpolated_data);

Append MEG and EEG data so that it can all be handled in one data object

%% append MEG and EEG

cfg = [];

preprocessed_data = ft_appenddata(cfg, preprocessed_data_MEG, rereferenced_interpolated_data);

Clean data by removing high-variance trials based on summary plots

  1. The three kinds of sensors are handled separately since the units are of different orders
    ‘MEGMAG’ = magnetometers
    ‘MEGGRAD’ = gradiometers
    ‘EEG’ = electrodes
  2. Note that keepchannel is necessary, otherwise the channels not chosen will be discarded (and you will thus end up with an empty dataset if it is set to false)
%% clean data by trial

cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
cfg.channel = 'MEGMAG';
cfg.layout = 'neuromag306eeg_1005_natmeg.lay';

cleaned_data = ft_rejectvisual(cfg, preprocessed_data);

cfg.channel = 'MEGGRAD';

cleaned_data = ft_rejectvisual(cfg, cleaned_data);

cfg.channel = 'EEG';

cleaned_data = ft_rejectvisual(cfg, cleaned_data);

Adjust the timeline

There was a mismatch in timing between the onset of the trigger and the actual stimulation of 41 msec

cfg = [];
cfg.offset = -41; % ms

cleaned_data = ft_redefinetrial(cfg, cleaned_data);

Downsample data to a sampling frequency to 200 Hz

This is done to optimise processing speed and data handling
Data is also saved at this point

cfg = [];
cfg.resamplefs = 200;

cleaned_downsampled_data = ft_resampledata(cfg, cleaned_data);

save(fullfile(output_path, 'cleaned_downsampled_data'), 'cleaned_downsampled_data');

Plot and identify components

%% plot components

% MEG
cfg = [];
cfg.component = 1:60;
cfg.layout = 'neuromag306mag.lay';
cfg.comment = 'no';

ft_topoplotIC(cfg, ica);

cfg = [];
cfg.layout = 'neuromag306mag.lay';
cfg.viewmode = 'component';

ft_databrowser(cfg, ica);

Remove the appropriate comments

%% remove appropriate components

cfg = [];
cfg.component = [1 2]; % NB! these should only be the eye related components
cfg.demean = 'no';

ica_removed = ft_rejectcomponent(cfg, ica, cleaned_downsampled_data);

Do a timelocked analysis (average all trials to increase signal-to-noise ratio)

Preprocessing options etc. are applied

1, Data is cropped, so we only look at a relevant period for timelocked components
2. Data is demeaned with the baseline period specified
3. Data is low-pass filtered
4. The covariance is calculated (necessary for MNE source reconstructions)
5. Timelockeds are also saved

%% timelocked

n_events = length(events);
timelockeds = cell(1, n_events);

cfg = [];
cfg.toilim = [-0.200 1.000];
    
cropped_data = ft_redefinetrial(cfg, ica_removed);

for event_index = 1:n_events
    
    event = events(event_index);
    
    cfg = [];
    cfg.covariance = 'yes';
    cfg.covariancewindow = 'prestim';
    cfg.preproc.demean = 'yes';
    cfg.preproc.baselinewindow = [-0.200 0];
    cfg.preproc.lpfilter = 'yes';
    cfg.preproc.lpfreq = 70;
    cfg.trials = cropped_data.trialinfo == event;

    timelockeds{event_index} = ft_timelockanalysis(cfg, cropped_data);
    
end

save(fullfile(output_path, 'timelockeds'), 'timelockeds');  

Calculate the time-frequency representation (TFR)

  1. First, we baseline the data that we are going to use for the TFRs
  2. We are calculating power for the frequencies in the range from 1 to 40 Hz (foilim (frequencies of interest limits))
  3. We are using the wavelet method
  4. We are calculating the power for each time sample (toi (times of interest))
  5. Gradiometers are combined for ease of interpretation
  6. TFRs are saved
%% tfrs

n_events = length(events);
tfrs = cell(1, n_events);

cfg = [];
cfg.demean = 'yes';
cfg.baselinewindow = 'all';

baseline_data = ft_preprocessing(cfg, ica_removed);

save(fullfile(output_path, 'baseline_data'), 'baseline_data');

for event_index = 1:n_events
    
    event = events(event_index);
    
    cfg = [];
    cfg.method = 'wavelet';
    cfg.foilim = [1 40];
    cfg.toi = -1.650:0.005:1.650;
    cfg.width = 7;
    cfg.pad = 'nextpow2';
    cfg.trials = baseline_data.trialinfo == event;

    tfrs{event_index} = ft_freqanalysis(cfg, baseline_data);
    
end

save(fullfile(output_path, 'tfrs'), 'tfrs', '-v7.3');

% combine tfrs

combined_tfrs = cell(1, n_events);

for event_index = 1:n_events
        
    cfg = [];
    
    combined_tfrs{event_index} = ft_combineplanar(cfg, tfrs{event_index});

end

save(fullfile(output_path, 'combined_tfrs'), 'combined_tfrs', '-v7.3');