Set up general paths

Change these to appropriate paths for your operating system and setup

clear variables
restoredefaultpath
addpath /home/lau/matlab/fieldtrip-20180129_workshop_teaching/
ft_defaults

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 if you have more subjects than one)

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 = 8; % right index finger
% to do all five fingers set events as below:
% events = [1 2 4 8 16]; %% right little, ring, middle, index, 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 code is read in (8), but note that all can be 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
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(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{:});

Inspect raw data

Try to inspect the data and see if you can identify some bad electrodes (and MEG sensors)


cfg = [];
cfg.continuous = 'yes';
cfg.viewmode = 'vertical'; % also try 'butterfly'

ft_databrowser(cfg, preprocessed_data_MEG);

ft_databrowser(cfg, preprocessed_data_EEG);

Remove bad electrodes

These would mess up the source reconstructions if left in
Please note that you need not find them now, as they are specified in the next step (interpolation)


% 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);

Plot neighbours

Define the neighbours of each electrode, so bad ones can be interpolated

Do the actual interpolation

cfg = [];
cfg.method                    = 'spline';
cfg.neighbours                = neighbours_EEG;
cfg.badchannel                = badchannels;
cfg.senstype                  = 'EEG';
interpolated_data = ft_channelrepair(cfg, preprocessed_data_EEG);
Plot the interpolated and non-interplolated

Can you spot the differences?

cfg = [];
cfg.viewmode = 'butterfly';

ft_databrowser(cfg, preprocessed_data_EEG); % not interpolated
ft_databrowser(cfg, interpolated_data);

NOT interpolated

Interpolated

Re-reference to common average

This is necessary as electric potentials can only compared by reference to a common reference

cfg = [];
cfg.reref                  = 'yes';
cfg.refchannel             = 'all';
 
rereferenced_interpolated_data = ft_preprocessing(cfg, interpolated_data);
Plot the re-referenced data

Append MEG and EEG data

Now the data can be handled together

cfg = [];

preprocessed_data = ft_appenddata(cfg, preprocessed_data_MEG, ...
                                        rereferenced_interpolated_data);

Clean the data by trial

  1. This is visually guided way to reject trials by removing those showing high variance
  2. Note, if you have several experimental conditions, always collapse all conditions before summarising and rejecting to avoid subjective bias
  3. You can also remove specific artefacts by using the the ft_artifact_xxx and ft_rejectartifact functions
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 for the offset between trigger value (event code) and the actual delivery of the stimulation

This moves the offset (0 msec) 41 msec back in time. This has been measured with an accelerometer

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

cleaned_data = ft_redefinetrial(cfg, cleaned_data);

Downsample the data (speeds up processing time, but decreases the number of frequencies we can look at)

Downsample from 1000 Hz to 200 Hz

cfg = [];
cfg.resamplefs = 200;

cleaned_downsampled_data = ft_resampledata(cfg, cleaned_data);

% save(fullfile(output_path, 'cleaned_downsampled_data'), ... % uncomment to save the data in your current path
  %                      'cleaned_downsampled_data');

Average data to obtain evoked responses (timelockeds)

  1. Loop through events to obtain the evoked response for each event (we only have one event, but this shows how you can make loops for several events)
  2. We have included noise covariance estimation, which is necessary if you want to combine magnetometers, gradiometers and electrodes when doing source analysis
  3. The responses are also demeaned, where a baseline period is used to adjust the whole time period. This is necessary for source reconstruction
  4. A lowpass-filter of 70 Hz is also included

Q: Why is it necessary to include noise covariance estimation if you want to combine magnetometers, gradiometers and electrodes when doing source analsis?

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

cfg = [];
cfg.toilim = [-0.200 0.600];
    
cropped_data = ft_redefinetrial(cfg, cleaned_downsampled_data);

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'); %% uncomment to
%save on output_path

Compare with non-demeaned

Just for making a comparison

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

cfg = [];
cfg.toilim = [-0.200 0.600];
    
cropped_data = ft_redefinetrial(cfg, cleaned_downsampled_data);

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

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

Plot demeaned and non-demeaned

See the difference that demeaning does

Q: Which time course is demeaned and which is not? Which topography is demeaned and which not?

Combine gradiometers and plot the three kinds of channels

Planar gradiometers are combined by taking the absolute sum of each pair of channels
After that EEG, gradiometers and magnetometers are plotted separately on a 2d plot
Then they are plotted on 3d surfaces. Use the rotation tools to answer the questions below

Q: Describe the three topographies in as much detail as possible answering:
a) do the topographies point towards a single dipolar source or several distrubuted sources?
b) is/are the underlying source(s) most likely to be radial or tangential dipoles?
c) what is are the position(s) of the underlying source(s)? (Indicate on an image with LibreOffice/Microsoft Office/etc.)
d) what is/are the orientation(s) of the source(s)? (Indicate on an image with LibreOffice/Microsoft Office/etc.)
e) on one of the plots it is not possible to get any information about the orientation of the source. Which one is that?
f) what are the channels showing the maximum/minimum responses?
g) what are the maximum/minimum values approximately for each of the plots?

%% combine gradiometers and subsequently plot

close all 

cfg = [];
cfg.method = 'sum';

combined_planar = ft_combineplanar(cfg, timelockeds{1});

% plot eeg

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306eeg1005_natmeg.lay';
cfg.channel = 'EEG';
cfg.xlim = [0.057 0.057];

ft_topoplotER(cfg, timelockeds{1});

% plot gradiometers

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306cmb.lay';
cfg.channel = 'MEGGRAD';
cfg.xlim = [0.057 0.057];

ft_topoplotER(cfg, combined_planar);

% plot magnetometers

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306mag.lay';
cfg.channel = 'MEGMAG';
cfg.xlim = [0.057 0.057];

ft_topoplotER(cfg, timelockeds{1});
%% 3d topo plots

close all

elec_pos = combined_planar.elec.chanpos;
sensors = combined_planar.grad;
mag_indices = 103:204;
grad_indices = 1:102;
mag_pos = sensors.chanpos(mag_indices, :);
grad_pos = sensors.chanpos(grad_indices, :);

time_index = 52;

figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(elec_pos, combined_planar.avg(205:end, time_index));
ft_plot_sens(combined_planar.elec, 'elecsize', 10, 'label', 'label')
c = colorbar;
c.Label.String = 'Electric Potential (V)';
title('EEG 57 msec', 'fontsize', 40);


figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(mag_pos, combined_planar.avg(mag_indices, time_index));
ft_plot_sens(combined_planar.grad, 'label', 'label', 'chantype', 'megmag')
c = colorbar;
c.Label.String = 'Magnetic Field Strength (T)';
title('Magnetometers 57 msec', 'fontsize', 40);


figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(grad_pos, combined_planar.avg(grad_indices, time_index));
ft_plot_sens(combined_planar.grad, 'label', 'label', 'chantype', 'unknown')
c = colorbar;
c.Label.String = 'Magnetic Field Gradient (T/m)';
title('Gradiometers 57 msec', 'fontsize', 40);