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 FieldTrip - Install it in your software folder - Take a look at the FieldTrip tutorials!
Input data
EEG signals (EEG)Fieldtrip can import the data coming from our BrainAmps system. If a 3D digitalisation of the EEG electrodes coordinates was performed you should insert this information in the file using multiconv (see the appropriate tutorial). It is also possible to use an EEG template (you can do that after the EEG importation in FieldTrip).
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
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.
MRI processing and grid normalisation to the MNI space
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 | ||
|---|---|---|
| ||
mri=ft_read_mri(dicom_file_name); %% ReadConstruire 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
| Code Block | ||
|---|---|---|
| ||
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('segmentedmrila grille du beamformer en d?formant le MNI sur le sujet templatedir = 'fieldtrip-20160105/template/sourcemodel'; template = load(fullfile(templatedir, 'standard_sourcemodel3d8mm')); %% 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') ; |
General configuration
| Code Block | ||
|---|---|---|
| ||
%% Load % create the headsubject modelspecific createdgrid, duringusing the MRItemplate pregrid processingthat loadhas voljust ;been %%created loadcfg head= model located in the subject directory and obtained from the T1 MRI of the subject %% 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
| Code Block | ||
|---|---|---|
| ||
% 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 will full window
cfg = [];
cfg.toilim = overall_time_window ;
full_data_all = ft_redefinetrial(cfg,trials);
|
Time-Frequency analysis
| Code Block | ||
|---|---|---|
| ||
cfg = []; cfg.output = 'pow'; cfg.channel = 'MEG'; cfg.method = 'mtmconvol'; cfg.taper = 'hanning'; cfg.foi = 7:2:40; % analysis 2 to 30 Hz in steps of 2 Hz cfg.t_ftimwin = ones(length(cfg.foi),1).*0.5; % length of time window = 0.5 sec cfg.toi = -2.0:0.05:2.0; % time window "slides" from -2.0 to 2.0 sec in steps of 0.05 sec (50 ms) cfg.trials = 'all'; TFRhann = ft_freqanalysis(cfg, full_data_all); %% Visualization cfg = []; cfg.xlim = [0.1 0.2]; cfg.ylim = [20 20]; cfg.zlim = [-1e-28 1e-28]; cfg.baseline = [-0.2 -0.0]; cfg.baselinetype = 'absolute'; cfg.layout = 'neuromag306mag.lay'; figure; ft_topoplotTFR(cfg,TFRhann); % for the multiple plots also: cfg = []; cfg.ylim = [20 20]; cfg.baseline = [-0.2 -0.0]; cfg.baselinetype = 'absolute'; cfg.layout = 'neuromag306mag.lay'; cfg.xlim = [-0.4:0.1:0.4]; cfg.comment = 'xlim'; cfg.commentpos = 'title'; figure; ft_topoplotTFR(cfg,TFRhann); cfg = []; cfg.baseline = [-1.0 -0.8]; cfg.baselinetype = 'absolute'; %% cfg.zlim = [-1e-27 1e-27]; []; 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,:)); |
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 |
|---|
%% Load the MRI data for visualisation purpose load('mri_aligned.mat') ; %% MRI data needed for vizualisation load('segmentedmri.mat') ; index = 1 ; conf = [] ; %% marker of interest markers = 'Move_OK' ; %% output label conf.label = 'WIKI_' ; %% Main Analysis for i=1:number_of_subjects %% Enter subject dir cd(list_sujets{i}{1}{:}) %% 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 ; for exam=1:2 %% Is the subject/session really here? if ~isempty([list_sujets{i}{1+exam}{:}]) %% enter exam subdirectory cfg.showlabels = 'yes'; cfg.layout = 'neuromag306mag.lay'; figure; ft_multiplotTFR(cfg, TFRhann); |
CSD Computation
| Code Block | ||
|---|---|---|
| ||
%% Compute cross-sepctral density matrices
powcsd_all = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_all) ;
powcsd_active = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_timewindow_active) ;
powcsd_ref = compute_crosprect(channelsofinterest, freqofinterest, freqhalfwin, data_timewindow_ref) ; |
Forward operator
| Code Block | ||
|---|---|---|
| ||
% Create leadfield grid
cfg = [];
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.grid.resolution = 1.0; % use a 3-D grid with a 0.5 cm resolution
cfg.grid.unit = 'cm';
cfg.grid.tight = 'yes';
cfg.normalize = lead_field_depth_normalization ; %% Depth normalization
[grid] = ft_prepare_leadfield(cfg); |
D.I.C.S. Analysis
| Code Block |
|---|
%% Compute DICS filters for the concatenated data cfg cd(list_sujets{i}{1+exam}{:}) %% Datasets to be scanned for the marker of interest datasets = {'B_cardiocor_blinkcor_run01_pre_tsss.fif' 'B_cardiocor_blinkcor_run02_pre_tsss.fif' 'B_cardiocor_blinkcor_run03_pre_tsss.fif'} ; %% 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( datasets, markers, 2.0, 2.0, 'BIO005' ) ; % 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 with full window cfg = []; cfg.toilim = overall_time_window ; full_data_all = ft_redefinetrial(cfg,trials); % Create leadfield grid cfg = []; cfg.grid = warped_grid ; cfg.vol = =hdm []; cfg.channel = channelsofinterest ; cfg.methodgrad = 'dics'; cfg.frequency = freqofinterest ; cfg.grid powcsd_active.grad; cfg.vol = grid; cfg.headmodel = vol; cfg.senstype = 'MEG'; cfg.dics.keepfilter reducerank = 'yes'2; % Wedefault wishfor toMEG useis the2, calculatedfor filterEEG lateris on3 cfg.dics.projectnoise = 'yes'; cfg.dics.lambdanormalize = beamformer_lambda lead_field_depth_normalization ; source_all %% Depth normalization [grid] = ft_prepare_sourceanalysisleadfield(cfg, powcsd_all); %% Apply computed filters to each active and reference windows cfg for freq_index=1:number_of_fois foi = fois{freq_index}{1} ; = []; cfg.channel = channelsofinterest wof = fois{freq_index}{2} ; cfg.method = 'dics'; cfg.frequency = freqofinterest ; cfg.grid = grid; cfg.grid.filter = source_all.avg.filter; cfg.dics.projectnoise = 'yes'; cfg.headmodel %% Analysis on MAG %% Select the channel of interest channelsofinterest = {'MEG*1'} ; %% Magnetometer on an Elekta system %% Compute DICS filters for the concatenated data cfg = vol[]; cfg.dics.lambdachannel = beamformer_lambda_normalizationchannelsofinterest ; cfg.senstypemethod = 'MEGdics' ; cfg.frequency %% Source analysis active= source_activefreqofinterest = 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.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.parameter = 'avg.pow'; source_active_int cfg.channel = 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); %% Display the results visuBF( source_diff_int, 'title' ) |
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);
%% Change grid locations for display in MNI
source_diff.pos = template.sourcemodel.pos;
source_diff.dim = template.sourcemodel.dim;
save('source_diff', [conf.label '_' conf.markers '_PRE_MAG_' fois{freq_index}{3}]) ;
%% Analysis on MAG
%% Select the channel of interest
channelsofinterest = {'MEG*2', 'MEG*3'}; %% Gradiometers on an Elekta system
%% 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);
%% Change grid locations for display in MNI
source_diff.pos = template.sourcemodel.pos;
source_diff.dim = template.sourcemodel.dim;
save('source_diff', [conf.label '_' conf.markers '_PRE_MAG_' fois{freq_index}{3}]) ;
end
cd ..
end
end %% for exam
cd ..
end %% for subject |
Basic statistical analysis
| Code Block |
|---|
ft_defaults ;
list_sujets = {
{{'gamma03_s03'},{'140703'}, {'140710'}, {'S03'}, {'sGAMMA_T_02_GI-0003-00001-000176-01.img'}}
{{'gamma07_s07'},{'141016'}, {'141104'}, {'S07'}, {'sGAMMA_CV_S07-0003-00001-000176-01.img'}}
{{'gamma08_s08'},{'141107'}, {'141121'}, {'S08'}, {'sGAMMA_MS_S08-0003-00001-000176-01.img'}}
{{'gamma10_s10'},{'141105'}, {'141112'}, {'S10'}, {'sGAMMA_AH_S10-0003-00001-000176-01.img'}}
{{'gamma18_s18'},{'150210'}, {'150217'}, {'S18'}, {'sGAMMA_S18-0003-00001-000176-01.img'}}
{{'gamma19_s19'},{'150219'}, {'150922'}, {'S19'}, {'sGAMMA_CV_S19-0003-00001-000176-01.img'}}
{{'gamma20_s20'},{'150306'}, {'150313'}, {'S20'}, {'sGAMMA_SR_S20-0003-00001-000176-01.img'}}
{{'gamma22_s22'},{'150410'}, {'150326'}, {'S22'}, {'sGAMMA_S22_PS-0003-00001-000176-01.img'}}
{{'gamma24_s24'},{'150410'}, {'150421'}, {'S24'}, {'sGAMMA_S24_KL-0003-00001-000176-01.img'}}
{{'gamma25_s25'},{'150403'}, {'150413'}, {'S25'}, {'sGAMMA_S25_LP-0003-00001-000176-01.img'}}
{{'gamma26_s26'},{'150423'}, {'150430'}, {'S26'}, {'sGAMMA_S26_SD-0003-00001-000176-01.img'}}
{{'gamma27_s27'},{'150506'}, {'150513'}, {'S27'}, {'sGAMMA_S27_KP-0005-00001-000176-01.img'}}
{{'gamma29_s29'},{'150527'}, {'150604'}, {'S29'}, {'sGAMMA_S29_AC-0003-00001-000176-01.img'}}
{{'gamma31_s31'},{'150618'}, {'150703'}, {'S31'}, {'sGAMMA_S31_MA-0003-00001-000176-01.img'}}
{{'gamma39_s39'},{'150929'}, {'151013'}, {'S39'}, {'sGAMMA_S39_LP-0003-00001-000176-01.img'}}
} ;
%% number of subjects
number_of_subjects = size(list_sujets, 1) ;
%% Frequencies of interest (here only one, we can define several if needed)
fois={{22,6, 'BETA_MNI'}} ;
number_of_fois = length(fois) ;
%% number of subjects
number_of_subjects = size(list_sujets, 1) ;
%% number of sessions
number_of_sessions = 2 ;
list_labels = {'WIKI_'} ;
list_markers = {'Move_OK'} ;
list_frequencies = {'BETA'} ;
number_of_labels = length(list_labels) ;
number_of_markers = length(list_markers) ;
number_of_frequencies = length(list_frequencies) ;
for i_label = 1:number_of_labels
for i_mk = 1:number_of_markers
for i_freq = 1:number_of_frequencies
sources_mag = {} ;
sources_grad = {} ;
index_exams = 1 ;
markers = list_markers{i_mk} ; %%
label = list_labels{i_label} ; %%
freq = list_frequencies{i_freq} ; %%
%% Main Analysis
%% Read MAG Results
nb_hv = 0 ;
conf.title_figure = [label '_' markers '_PRE_MAG_' freq '_MNI'] ;
for i=1:number_of_subjects
cd(list_sujets{i}{1}{:})
list_sujets{i}{1}{:}
for exam=1
%% enter exam subdirectory
if ~isempty([list_sujets{i}{1+exam}{:}])
cd(list_sujets{i}{1+exam}{:})
load([conf.title_figure '_diff.mat']) ;
nb_hv = nb_hv + 1 ;
sources_mag{nb_hv} = source_diff ;
index_exams = index_exams + 1 ;
cd ..
end
end %% for exam
cd ..
end %% for subject
conf = [] ;
mean_hv = zeros(size(sources{1}.avg.pow)) ;
for i=1:nb_hv
mean_hv = mean_hv + sources_mag{i}.avg.pow ;
end ;
mean_sources_mag = sources{i} ;
mean_sources_mag.avg.pow = mean_hv / nb_hv ;
visuBF_MNI(mean_sources_mag, 'test')
%% Read GRAD Results
nb_hv = 0 ;
conf.title_figure = [label '_' markers '_PRE_GRAD_' freq '_MNI'] ;
for i=1:number_of_subjects
cd(list_sujets{i}{1}{:})
list_sujets{i}{1}{:}
for exam=1
%% enter exam subdirectory
if ~isempty([list_sujets{i}{1+exam}{:}])
cd(list_sujets{i}{1+exam}{:})
load([conf.title_figure '_diff.mat']) ;
nb_hv = nb_hv + 1 ;
sources_grad{nb_hv} = source_diff ;
index_exams = index_exams + 1 ;
cd ..
end
end %% for exam
cd ..
end %% for subject
conf = [] ;
mean_hv = zeros(size(sources{1}.avg.pow)) ;
for i=1:nb_hv
mean_hv = mean_hv + sources{i}.avg.pow ;
end ;
mean_sources_grad = sources{i} ;
mean_sources_grad.avg.pow = mean_hv / nb_hv ;
visuBF_MNI(mean_sources_grad, 'test')
%% print -djpeg -r300 ['VS_' conf.title_figure '.jpeg']
end %% Frequency
end %% Markers
end %% LABEL
cfg = [];
cfg.dim = sources_mag{1}.dim;
cfg.method = 'montecarlo';
cfg.parameter = 'pow';
cfg.correctm = 'max';
cfg.numrandomization = 2000;
cfg.alpha = 0.01; % note that this only implies single-sided testing
cfg.tail = 0;
cfg.statistic = 'ft_statfun_depsamplesT';
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_hv)] ;
cfg.design(2,:) = [ones(1,nb_hv) ones(1,nb_hv)*2] ;
stat_hv = ft_sourcestatistics(cfg, sources_mag{:}, sources_grad{:}) ;
visuBF_MNI_STAT(stat_hv, 'test') ;
|
