In this tutorial you will find a description of all the steps necessary to obtain the sources localisation of continuous (non averaged) data using a beamformer approach in FieldTrip for a group of subjects in using the MNI template

Prerequisites

Read the previous tutorials (MEG EEG pre-processing, MEG sensors realignment, filter, baseline correction)
Pre process the T1 MRI data or choose a template
Download the latest version of F ieldTrip - Install it in your software folder - Take a look at the FieldTrip tutorials!

Input data

MEG signals (MEG)

FieldTrip can import the fif files coming from our MEG system and those pre processed with our in-house tools. If you acquired simultaneous MEG+EEG data, FieldTrip will import both of them.

MRI

If you have the individual anatomy of each of your subjects (T1 MRI), you should process it with FieldTrip as described below.

If you don't have the individual anatomy, then you have to choose a template in FieldTrip.

Needed software

F ieldTrip

For MEG fif files pre processed outside FieldTrip: You will need two in-house scripts to read the marker files we add to the MEG fif files (put them in your matlab path):

A short script to pre process the MRI in FieldTrip: pre_process_mri_for_ft.m

A simple script to visualize beamformer results: visuBF.m

D.I.C.S. analysis with FieldTrip

Overview of the process:

Suggested data organisation

The mat files generated by FieldTrip can be stored together with the pre-processed data. You may want to create a specific sub folder "FT" to store them in each subject.

Read data

The data should be pre processed (i.e. artefact free) at this point and the FieldTrip folder in your matlab path.

MEG data

Code Block

theme	Confluence
title	Read MEG data

ft_defaults ; %% Load fieldtrip configuration

%% Read two runs, define trials around LEFT_MVT_CLEAN marker (-2s to 2s around the marker), read also BIO005 which is an EMG.
[ trials ] = trial_definition_ft( {'cardiocor_blinkcor_run03_tsss.fif' 'cardiocor_blinkcor_run04_tsss.fif'}, 'LEFT_MVT_CLEAN', 2.0, 2.0, 'BIO005' ) ;

MRI

processing

processin

This step should be done only one time. The results (aligned and segmented MRI, head model and grid) will be loaded when needed.

Code Block

title	MRI processing

mri=ft_read_mri(dicom_file_name);
%% Realign the MRI to the MEG coordinate
cfg=[];
cfg.method = 'interactive' ;
cfg.coordsys = 'neuromag' ;
mri_aligned = ft_volumerealign(cfg,mri) ;

save('mri_aligned.mat', 'mri_aligned');

%% Segmentation
cfg           = [];
cfg.output    = 'brain';
segmentedmri  = ft_volumesegment(cfg, mri_aligned);

save('segmentedmri.mat', 'segmentedmri') ;

 
% construct volume conductor model (i.e. head model) for each subject
cfg = [];
cfg.method = 'singleshell';
vol = ft_prepare_headmodel(cfg, seg);
vol = ft_convert_units(vol, 'cm');
save('vol.mat', 'vol') ;
 
% create the subject specific grid, using the template grid that has just been created
cfg = [];
cfg.grid.warpmni   = 'yes';
cfg.grid.template  = template_grid;
cfg.grid.nonlinear = 'yes'; % use non-linear normalization
cfg.mri            = mri_aligned ;
warped_grid          = ft_prepare_sourcemodel(cfg);
save('warped_grid.m', 'warped_grid')

% make a figure of the single subject headmodel, and grid positions
figure;
ft_plot_vol(vol, 'edgecolor', 'none'); alpha 0.5;
ft_plot_mesh(grid.pos(grid.inside,:));

General configuration

Here we load all the necessary stuff and we set the main parameters of the analysis:

time windows for active state and baseline
frequency of interest
beamformer configuration
MEG channels to be used in the analysis

Code Block

title	Parameters

%% Load the head model created during the MRI pre processing
load vol ; %% load head model located in the subject directory and obtained from the T1 MRI of the subject
 
%% Load the warped MNI grid
load warped_grid ;
%% Load the MRI data for visualisation purpose
load('mri_aligned.mat') ; %% MRI data needed for vizualisation
load('segmentedmri.mat') ;

overall_time_window = [-2 2] ; %% if you modify this modify the trials def above ! Is this used below ???

%% Define time window for active state and baseline
active_time_window = [0.2 0.6] ;
reference_time_window = [-1.4 -1.0] ;

%% Basic configuration of the beamformer
beamformer_lambda_normalization = '10%' ;
lead_field_depth_normalization = 'column' ;

%% Frequency of interest and associated width (we will examine here the 16 Hz to 28 Hz band)
freqofinterest = 22 ;
freqhalfwin = 6 ;

%% Select the channel of interest
channelsofinterest = {'MEG*2', 'MEG*3'}; %% Gradiometers on an Elekta system
%% channelsofinterest = {'MEG*1'} ; %% Magnetometer on an Elekta system

Data preparation

Different sets of trials are prepared for the coming analysis

Code Block

title	Data preparation

% Select time window of interest
cfg = [];
cfg.toilim = active_time_window ;
data_timewindow_active = ft_redefinetrial(cfg,trials);

% Select time window of control
cfg = [];
cfg.toilim = reference_time_window ;
data_timewindow_ref = ft_redefinetrial(cfg,trials);

% Concatenate for common spatial filter computation
data_all = ft_appenddata([], data_timewindow_active, data_timewindow_ref);
 
% Trials willwith full window
cfg = [];
cfg.toilim = overall_time_window ;
full_data_all = ft_redefinetrial(cfg,trials);

Forward operator

We compute the leans field based on the warped MNI grid.

Code Block

title	Forward operator

% Create leadfield grid
cfg                 = [];
cfg.grid        	= warped_grid ;
cfg.vol         	= hdm ;
cfg.channel         = channelsofinterest ;
cfg.grad            = powcsd_active.grad;
cfg.vol             = vol ;
cfg.dics.reducerank = 2; % default for MEG is 2, for EEG is 3
cfg.normalize       =  lead_field_depth_normalization ; %% Depth normalization
[grid] = ft_prepare_leadfield(cfg);

D.I.C.S. Analysis

The spatial filters are computed here. We then contrast the two conditions before saving the results. One very important step necessary to allow an easy visualisation of the results is to change the localisation of the grid points back to the original templates ones.

Code Block

%% Compute DICS filters for the concatenated data
cfg              = [];
cfg.channel      = channelsofinterest ;
cfg.method       = 'dics';
cfg.frequency    = freqofinterest ;
cfg.grid         = grid;
cfg.headmodel    = vol;
cfg.senstype     = 'MEG';
cfg.dics.keepfilter   = 'yes'; % We wish to use the calculated filter later on
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda  = beamformer_lambda_normalization;
source_all = ft_sourceanalysis(cfg, powcsd_all);

%% Apply computed filters to each active and reference windows
cfg              = [];
cfg.channel      = channelsofinterest ;
cfg.method       = 'dics';
cfg.frequency    = freqofinterest ;
cfg.grid         = grid;
cfg.grid.filter  = source_all.avg.filter;
cfg.dics.projectnoise = 'yes';
cfg.headmodel          = vol;
cfg.dics.lambda  = beamformer_lambda_normalization ;
cfg.senstype     ='MEG' ;

%% Source analysis active
source_active = ft_sourceanalysis(cfg, powcsd_active);

%% Source analysis reference
source_reference = ft_sourceanalysis(cfg, powcsd_ref);

%% Compute differences (power) between active and reference
source_diff = source_active ;
source_diff.avg.pow  = (source_active.avg.pow - source_reference.avg.pow) ./ (source_active.avg.pow + source_reference.avg.pow);

%% Display the results
visuBF( source_diff, 'title' )
%% MRI reslicing

mri_resliced = ft_volumereslice([], mri_aligned);

%% Display results superimposed on MRI

cfg            = [];
cfg.parameter = 'avg.pow';
source_active_int  = ft_sourceinterpolate(cfg, source_active, mri_resliced);
source_reference_int  = ft_sourceinterpolate(cfg, source_reference, mri_resliced);
source_diff_int  = source_active_int;
source_diff_int.pow  = (source_active_int.pow - source_reference_int.pow) ./ (source_active_int.pow + source_reference_int.pow);

%% Change sourcegrid locations for display in MNI
sources_exp.pos = template.sourcemodel.pos;
sources_exp.dim = template.sourcemodel.dim;

% 
sources_bsl.pos = template.sourcemodel.pos;
sources_bsl.dim = template.sourcemodel.dim;

%% Save results

save([conf.output '_exp.mat'], 'sources_exp')
save([conf.output '_bsl.mat'], 'sources_bsl')

Image Added

Basic statistical analysis

Code Block




DO THE STATISTICS NOW
list_sujets_HV = {

{{'gamma01_s01'},{}, {'140702'}, {'S01'}, {'sGAMMA_VS_01-0003-00001-000176-01.img'}}
{{'gamma03_s03'},{'140710'}, {'140703'}, {'S03'}, {'sGAMMA_T_02_GI-0003-00001-000176-01.img'}}
{{'gamma09_s09'},{'141020'}, {'141107'}, {'S09'}, {'sGAMMA_BD_S09-0003-00001-000176-01.img'}}
{{'gamma10_s10'},{'141112'}, {'141105'}, {'S10'}, {'sGAMMA_AH_S10-0003-00001-000176-01.img'}}
{{'gamma13_s13'},{'150106'}, {'141209'}, {'S13'}, {'sGAMMA_NJ_S13-0003-00001-000176-01.img'}}
{{'gamma14_s14'},{'150106'}, {'150113'}, {'S14'}, {'sGAMMA_MI_S14-0003-00001-000176-01.img'}}
{{'gamma15_s15'},{'150120'}, {'150108'}, {'S15'}, {'sGAMMA_GC_S15-0003-00001-000176-01.img'}}
{{'gamma19_s19'},{'150922'}, {'150219'}, {'S19'}, {'sGAMMA_CV_S19-0003-00001-000176-01.img'}}
{{'gamma24_s24'},{'150421'}, {'150410'}, {'S24'}, {'sGAMMA_S24_KL-0003-00001-000176-01.img'}}
{{'gamma25_s25'},{'150403'}, {'150413'}, {'S25'}, {'sGAMMA_S25_LP-0003-00001-000176-01.img'}}
{{'gamma28_s28'},{'150519'}, {'150512'}, {'S28'}, {'sGAMMA_S28_GA-0003-00001-000176-01.img'}}
{{'gamma29_s29'},{'150527'}, {'150604'}, {'S29'}, {'sGAMMA_S29_AC-0003-00001-000176-01.img'}}
{{'gamma30_s30'},{'150610'}, {'150528'}, {'S30'}, {'sGAMMA_S20_LM-0003-00001-000176-01.img'}}
{{'gamma32_s32'},{'150626'}, {'150710'}, {'S32'}, {'sGAMMA_S32_GJ-0003-00001-000176-01.img'}}
{{'gamma34_s34'},{'150706'}, {'150715'}, {'S34'}, {'sGAMMA_S34_BF-0003-00001-000176-01.img'}}
{{'gamma35_s35'},{'150709'}, {'150716'}, {'S35'}, {'sGAMMA_S35_BA-0005-00001-000176-01.img'}}
{{'gamma37_s37'},{'150720'}, {'150728'}, {'S37'}, {'sGAMMA_S37_TE-0003-00001-000176-01.img'}}
{{'gamma39_s39'},{'151013'}, {'150929'}, {'S39'}, {'sGAMMA_S39_LP-0003-00001-000176-01.img'}}
{{'gamma41_s41'},{'151012'}, {'151021'}, {'S41'}, {'sGAMMA_S41_DT-0005-00001-000176-01.img'}}
{{'gamma42_s42'},{'151014'}, {'151027'}, {'S42'}, {'sGAMMA_S42_PS-0003-00001-000176-01.img'}}
} ;

%% List des patients -> MEG ID - SHAM tACS - REAL tACS - MRI ID - T1 Name
list_sujets_P = {

{{'gamma04_s04'},{'140722'}, {'151106'}, {'S04'}, {'sGAMMA_WC_P02-0004-00001-000176-01.img'}}
{{'gamma05_s05'},{}, {'141103'}, {'S05'}, {'sGAMMA_OC_S05-0003-00001-000176-01.img'}}
{{'gamma06_s06'},{'141014'}, {'141021'}, {'S06'}, {'sGAMMA_JJ_S06-0003-00001-000176-01.img'}}
{{'gamma07_s07'},{'141016'}, {'141104'}, {'S07'}, {'sGAMMA_CV_S07-0003-00001-000176-01.img'}}
{{'gamma08_s08'},{'141121'}, {'141107'}, {'S08'}, {'sGAMMA_MS_S08-0003-00001-000176-01.img'}}
{{'gamma11_s11'},{'141216'}, {'141125'}, {'S11'}, {'sGAMMA_S11_WP-0003-00001-000176-01.img'}}
{{'gamma12_s12'},{'141218'}, {'150116'}, {'S12'}, {'sGAMMA_CC_S12-0004-00001-000176-01.img'}}
{{'gamma16_s16'},{'150209'}, {'150216'}, {'S16'}, {'sGAMMA_S16-0003-00001-000176-01.img'}}
{{'gamma17_s17'},{'150311'}, {'150218'}, {'S17'}, {'sGAMMA_S17_MM-0003-00001-000176-01.img'}}
{{'gamma18_s18'},{'150217'}, {'150210'}, {'S18'}, {'sGAMMA_S18-0003-00001-000176-01.img'}}
{{'gamma20_s20'},{'150313'}, {'150306'}, {'S20'}, {'sGAMMA_SR_S20-0003-00001-000176-01.img'}}
{{'gamma22_s22'},{'150410'}, {'150326'}, {'S22'}, {'sGAMMA_S22_PS-0003-00001-000176-01.img'}}
{{'gamma26_s26'},{'150430'}, {'150423'}, {'S26'}, {'sGAMMA_S26_SD-0003-00001-000176-01.img'}}
{{'gamma27_s27'},{'150506'}, {'150513'}, {'S27'}, {'sGAMMA_S27_KP-0005-00001-000176-01.img'}}
{{'gamma31_s31'},{'150618'}, {'150703'}, {'S31'}, {'sGAMMA_S31_MA-0003-00001-000176-01.img'}}
{{'gamma33_s33'},{'150623'}, {'150601'}, {'S33'}, {'sGAMMA_S32_GJ-0003-00001-000176-01.img'}}
{{'gamma43_s43'},{'151019'}, {'151026'}, {'S43'}, {'sGAMMA_S43_GA-0003-00001-000176-01.img'}}
} ;

 

%% number of subjects
number_of_subjects_hv = size(list_sujets_HV, 1) ;
number_of_subjects_p = size(list_sujets_P, 1) ;

 

%% number of sessions
number_of_sessions = 2 ; 

%% Construire la grille du beamformer en d?formant le MNI sur le sujet
templatedir  = 'fieldtrip-20160105/external/spm8/templates';
template = ft_read_mri(fullfile(templatedir, 'T1.nii'));

list_labels = {'PREMOV', 'POSTMOV_S', 'POSTMOV_L'} ;
list_markers = {'FIRST_DEV_OK', 'FIRST_NODEV_OK', 'LAST_DEV_OK', 'MID_DEV_OK'} ;
list_sensors = {'MAG', 'GRAD'} ;
list_frequencies = {'THETA', 'ALPHA', 'BETA', 'GAMMA'} ;

number_of_labels = lenght(list_labels) ;
number_of_markers = lenght(list_markers) ;
number_of_sensors = lenght(list_sensors) ;
number_of_frequencies = lenght(list_frequencies) ;

for i_label = 1:number_of_labels
    for i_mk = 1:number_of_mk
        for i_sen = 1:number_of_sen
            for i_freq = 1:number_of_freq  
              
sources_pre_hv = {} ;
sources_pre_p = {} ;

sources_post_hv = {} ;
sources_post_p = {} ;

index_exams = 1 ;

markers = list_markers{i_mk} ;    %% 'FIRST_DEV_OK' ;
label = list_labels{i_label} ;    %% 'PREMOV' ;
sensors = list_sensors{i_sen} ;   %% 'MAG' ; 
freq = list_frequencies{i_freq} ; %%'GAMMA' ;

 

conf = [] ;
conf.markers = markers ;
conf.title_figure = [label '_' conf.markers '_PRE_' sensors '_' freq '_MNI'] ;
conf.template = template ;

 

nb_hv = 0 ;
nb_p = 0 ;

 

%% Main Analysis

for i=1:number_of_subjects_hv
    cd(list_sujets_HV{i}{1}{:})
    list_sujets_HV{i}{1}{:}
    for exam=1
 %% enter exam subdirectory
        if ~isempty([list_sujets_HV{i}{1+exam}{:}])
            cd(list_sujets_HV{i}{1+exam}{:})            
            load([conf.title_figure '_diff.mat']) ;
            sources_post_hv{index_exams} = source_diff ;
            nb_hv = nb_hv + 1 ;
            index_exams = index_exams + 1 ;
            cd ..
        end
    end %% for exam
    cd ..
end %% for subject

 

index_exams = 1 ;
%% Main Analysis
for i=1:number_of_subjects_p
    cd(list_sujets_P{i}{1}{:})
    list_sujets_P{i}{1}{:}
    for exam=1
        %% enter exam subdirectory
        if ~isempty([list_sujets_P{i}{1+exam}{:}])
            cd(list_sujets_P{i}{1+exam}{:})     
            load([conf.title_figure '_diff.mat']) ;
            sources_post_p{index_exams} = source_diff ;
            nb_p = nb_p + 1 ;
            index_exams = index_exams + 1 ;
            cd ..
        end
    end %% for exam
    cd ..
end %% for subject

 

conf = [] ; 
mean_hv = zeros(size(sources_post_hv{1}.avg.pow)) ;
for i=1:nb_hv
    conf.title_figure='test';
    conf.template = template ;
    mean_hv = mean_hv + sources_post_hv{i}.avg.pow ;
end ;

mean_sources_hv = sources_post_hv{i} ;
mean_sources_hv.avg.pow = mean_hv / nb_hv ;

visuBF_MNI( conf, mean_sources_hv )
print -djpeg -r300 ['VS_' conf.title_figure '.jpeg']

conf = [] ; 
mean_p = zeros(size(sources_post_p{1}.avg.pow)) ;
for i=1:nb_p
    conf.title_figure='test';
    conf.template = template ;
    mean_p = mean_p + sources_post_p{i}.avg.pow ;
end ;

mean_sources_p = sources_post_p{i} ;
mean_sources_p.avg.pow = mean_p / nb_p ;
visuBF_MNI( conf, mean_sources_p )
print -djpeg -r300 ['P_' conf.title_figure '.jpeg']

cfg = [];
cfg.dim         = sources_post_hv{1}.dim;
cfg.method      = 'montecarlo';
cfg.parameter   = 'pow';
cfg.correctm    = 'max';
cfg.numrandomization = 2000;
cfg.alpha       = 0.05; % note that this only implies single-sided testing
cfg.tail        = 0;
cfg.statistic   = 'ft_statfun_indepsamplesT';
%% cfg.uvar        = 1; % row of design matrix that contains unit variable (in this case: trials)
cfg.ivar        = 2; % row of design matrix that contains independent variable (the conditions)
cfg.design(1,:) = [ones(1,nb_hv)  ones(1,nb_p)] ;
cfg.design(2,:) = [ones(1,nb_hv)  ones(1,nb_p)*2] ;
stat_hv = ft_sourcestatistics(cfg, sources_post_hv{:}, sources_post_p{:}) ;

conf = [] ;
conf.template = template ;
visuBF_MNI_STAT(conf, stat_hv) ;
print -djpeg -r300 ['STAT_' conf.title_figure '.jpeg']


%% Display the results
visuBF( source_diff_int, 'title' )

Image Removed

Versions Compared

Old Version 5

New Version 6

Key

Prerequisites

Input data

MEG signals (MEG)

MRI

Needed software

D.I.C.S. analysis with FieldTrip

Suggested data organisation

The mat files generated by FieldTrip can be stored together with the pre-processed data. You may want to create a specific sub folder "FT" to store them in each subject.

Read data

MRI

processin

General configuration

Data preparation

Forward operator

Forward operator

D.I.C.S. Analysis

Basic statistical analysis

Page Comparison

Versions Compared

Old Version 5

New Version 6

Key

Prerequisites

Input data

MEG signals (MEG)

MRI

Needed software

D.I.C.S. analysis with FieldTrip

Suggested data organisation

The mat files generated by FieldTrip can be stored together with the pre-processed data. You may want to create a specific sub folder "FT" to store them in each subject.

Read data

MRI

processin

General configuration

Data preparation

Forward operator

Forward operator

D.I.C.S. Analysis

Basic statistical analysis