This is a collaborative space. In order to contribute, send an email to maximilien.chaumon@icm-institute.org
On any page, type the letter L on your keyboard to add a "Label" to the page, which will make search easier.

Single trials - Sources analysis - Group Level

Owned by Former user (Deleted)
Last updated: Sept 26, 2017

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

Prerequisites

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

FieldTrip

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.

FieldTrip initialisation.

Read MEG data
ft_defaults ; %% Load fieldtrip configuration

 

 

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.

MRI processing
mri=ft_read_mri(dicom_file_name);
 
%% Construire la 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') ;
 
% 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,:));

 

 

 

Data preparation

 

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
  • Different sets of trials are prepared for the coming analysis

Data preparation
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 ;

% Define the list of subjects (local structure with for each selected subject, meg label, date of first meg session, date of second meg session, mri label, t1 name)
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) ;
 
%% Charger la grille du beamformer en d?formant le MNI sur le sujet
templatedir  = 'fieldtrip-20160105/template/sourcemodel';

template = load(fullfile(templatedir, 'standard_sourcemodel3d8mm'));

 

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. 

%% 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
            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.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);
            for freq_index=1:number_of_fois

				foi = fois{freq_index}{1} ;
                wof = fois{freq_index}{2} ;
 
				%% Analysis on MAG
				%% Select the channel of interest
				channelsofinterest = {'MEG*1'} ; %% Magnetometer 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}]) ;

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

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