Group analysis in MNE-Python

Group analysis in FieldTrip

Clear Workspace And Set Up Paths

The first part involves setting up the paths such that data can be handled easily and adding the FieldTrip functions to your path. Note that you have to change the path to your FieldTrip and chance the data path to where you have saved the data

%% CLEAN AND SET PATHS

clear variables % clear workspace
restoredefaultpath

% CHANGE THE FieldTrip PATH ACCORDING TO YOUR SYSTEM
addpath /home/lau/matlab/fieldtrip-20180108_workshops_MEG_Nord/
ft_defaults

% CHANGE "data_path" ACCORDING TO YOUR SYSTEM PATH
data_path = '/home/share/workshops_MEG_Nord/mikkel_and_lau/data/';
meg_path = fullfile(data_path, 'sub-01/ses-meg/meg');
mri_path = fullfile(data_path, 'sub-01/ses-mri/anat');

Define And Read In Trials

1. The data set for this subjects contains three files that have undergone some processing already:
   
2. tsss indicates that temporal Signal Space Separation has been applied
   
3. mc indicates that movement correction has been applied
   
4. These two operations have been applied with MaxFilter, which is proprietary software from Elekta
   
5. Three files exist because the fif format is limited to 2GB filesize, thus we need to define the trials from all three files
   
6. Note that the trial duration is 3000 ms, but that we are preparing for a 41 ms shift since we know that there is a delay between trigger and stimulation

%% PREPROCESS RAW FILES

filenames = {
    'oddball_absence-tsss-mc_meg.fif'
    'oddball_absence-tsss-mc_meg-1.fif'
    'oddball_absence-tsss-mc_meg-2.fif'
            };
        
n_filenames = length(filenames);
offset = 41; % the delay in ms between trigger and actual stimulation

header = ft_read_header(fullfile(meg_path, filenames{1})); % recording info
sampling_frequency = header.Fs; % find the sampling frequency (here 1000 Hz)

% process split files separately

split_files = cell(1, n_filenames);

for filename_index = 1:n_filenames
    
    filename = filenames{filename_index};
    full_path = fullfile(meg_path, filename);
    
    % define trials
    
    cfg = [];
    cfg.dataset = full_path;
    cfg.trialfun = 'ft_trialfun_general'; % using the default trial function
    % (s) preparing adjustment of 41 ms; trial duration (3 s)
    cfg.trialdef.prestim = 1.500 - offset / sampling_frequency; 
    cfg.trialdef.poststim = 1.500 +  - offset / sampling_frequency; 
    cfg.trialdef.eventvalue = [1 2 3 4 5]; % all stimulations (800 total)
    cfg.trialdef.eventtype = 'STI101'; %% trigger channel
    
    cfg = ft_definetrial(cfg);
    
    % preprocess defined trials, (no operations applied yet)
    
    cfg.demean = 'no'; % will be done later
    cfg.lpfilter = 'no'; % will be done later
    
    split_files{filename_index} = ft_preprocessing(cfg);
    
end

Append Split Files

After reading in the trials, we can append the split data files for ease of data handling

%% APPEND SPLIT FILES

cfg = [];

preprocessed_data = ft_appenddata(cfg, split_files{:});

Adjust Offset Of Trial

Here, we are moving time zero (stimulation onset) 41 ms back (that’s why we are multiplying with minus 1)

%% ADJUST TRIAL OFFSET

cfg = [];
cfg.offset = -41 * sampling_frequency / 1000; % samples

preprocessed_data = ft_redefinetrial(cfg, preprocessed_data);

Downsample

To cut down on data size and to reduce processing time, we are going to downsample to 200 Hz This will also work as a lowpass filter and will limit the frequency range that we can consider (< 100 Hz)

%% DOWNSAMPLE

cfg = [];
cfg.resamplefs = 200;

preprocessed_downsampled_data = ft_resampledata(cfg, preprocessed_data);

Clean data

Clean the data by removing the high-variance trials by using a summary

%% CLEAN DATA
% remove high-variance trials

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

cleaned_data = ft_rejectvisual(cfg, preprocessed_downsampled_data);

cfg.channel = 'MEGGRAD';

cleaned_data = ft_rejectvisual(cfg, cleaned_data);

Save Preprocessed Data

This will be a good time to save. (Users with memory problems can start from here)

%% SAVE PREPROCESSED DATA

save(fullfile(meg_path, 'cleaned_data'), 'cleaned_data');

Inspect Data

We can inspect the data by using ft_databrowser

%% INSPECT DATA

cfg = [];
cfg.continuous = 'no';
cfg.viewmode = 'vertical'; % also try 'butterfly'
cfg.channel = 'MEG';

ft_databrowser(cfg, cleaned_data);

Data browser

Data browser; single channel

Timelocked Analysis

Now, we will do a timelocked analysis. (Averaging all responses to a given timepoint)
(Also called Evoked Response or Event-Related Field (ERF))
Here, we are applying two kinds of preprocessing before taking the average
1. The data is demeaned using a “baseline” period before the onset of stimulation
2. A low-pass filter of 70 Hz is applied, aiming at filtering away any contributions higher than 70 Hz

We are not going to process these any further, plots are included below for illustration and for sanity checks

%% TIMELOCKED ANALYSIS
% for the timelocked analysis, we are cropping the responses, because we
% there are not many timelocked responses after 600 ms
cfg = [];
cfg.toilim = [-0.200 0.600];

cropped_data = ft_redefinetrial(cfg, cleaned_data);
      
cfg = [];
cfg.preproc.demean = 'yes'; % demeaning ...
cfg.preproc.baselinewindow = [-0.200 0.000]; % ... based on this interval
cfg.preproc.lpfilter = 'yes'; % applying low-pass filter ...
cfg.preproc.lpfreq = 70; % ... with this frequency (Hz)

timelocked = ft_timelockanalysis(cfg, cropped_data);

Save Timelocked

%% SAVE TIMELOCKED

save(fullfile(meg_path, 'timelocked'), 'timelocked');

Plot Magnetometers

Draw boxes around sensors to investigate single sensors (or means of several)
From single sensors, draw boxes around responses to investigate topographies

%% PLOT TIMELOCKED MAGNETOMETERS

figure

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

ft_multiplotER(cfg, timelocked);

Time-Frequency Representation

We can represent the fluctuations in power in different frequency bands over time 1. To remove any specific trial offsets in power, we are first demeaning the data by using the whole time window time as the “baseline” data. 2. After that we apply a wavelet analysis to the frequencies in the range of 1 to 40 Hz
3. We consider a time range of 3 s, estimating power for every 5 milliseconds.

%% TFR ANALYSIS

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

demeaned_data = ft_preprocessing(cfg, cleaned_data);

cfg = [];
cfg.method = 'wavelet';
cfg.foilim = [1 40];
cfg.toi = -1.500:0.005:1.500;
cfg.width = 7;
cfg.pad = 'nextpow2';

tfr = ft_freqanalysis(cfg, demeaned_data);

Save TFR

%% SAVE TFR

save(fullfile(meg_path, 'tfr'), 'tfr');

Combine Planar Gradiometers

1. Our MEG system contains 102 sensor triplets, each containing a magnetometer (plotted above) and two planar gradiometers orthogonal to one another
   
2. Due to their orthogonality, they see differently directed gradients
   
3. To summarize what the gradiometers on a sensor triplet “see”, it is useful to combine the absolute responses from each gradiometer. (Note that any information about polarity is lost in this step)

%% COMBINE PLANAR GRADIOMETERS

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

combined_tfr = ft_combineplanar(cfg, tfr);

Save Combined TFR

%% SAVE COMBINED TFR

save(fullfile(meg_path, 'combined_tfr'), 'combined_tfr')

Plot Multiplot TFR

1. Since absolute numbers of power are hard to interpret, we are going to demean every estimated frequency by the mean power in that frequency over the whole trial
   
2. Furthermore, since there is more power in low frequencies than in high frequencies, we are using a relative baseline type
3. To zoom in on the beta rebound, we are reducing the xlim and the ylim
   
4. Try for yourself to see what it looks like, if you use absolute as baseline type, or if you plot the whole time and frequency range

%% PLOT TFR (GRADIOMETERS)

figure

cfg = [];
cfg.layout = 'neuromag306cmb.lay';
cfg.baseline = [-Inf Inf];
cfg.baselinetype = 'relative';
cfg.xlim = [0.300 1.500]; % s
cfg.ylim = [10 40]; % Hz
cfg.colorbar = 'yes';

ft_multiplotTFR(cfg, combined_tfr);

Multiplot TFR

z-axis (colour coding) shows relative chance to the mean power in that frequency band

Single plot TFR

Steps From Here

Now, we have identified the beta rebound. To estimate its neural origin, we need to preprocess the MRI data, which is the subject of the next tutorial