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
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/';
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
%% 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{:});
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
% 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
cfg = [];
preprocessed_data = ft_appenddata(cfg, preprocessed_data_MEG, rereferenced_interpolated_data);
%% 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);
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);
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');
Eye blinks and movements cause big artefacts that mask the signal from the brain
%% find ica components
cfg = [];
cfg.method = 'runica';
cfg.numcomponent = 60;
cfg.demean = 'yes';
cfg.channel = {'MEG' 'EEG'};
ica = ft_componentanalysis(cfg, cleaned_downsampled_data);
save(fullfile(output_path, 'ica'), 'ica');
%% 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 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);
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');
%% 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');