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.

Skip to end of metadata
Go to start of metadata

You are viewing an old version of this page. View the current version.

Compare with Current View Page History

« Previous Version 5 Next »

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

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.

MEG data

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

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

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

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 will full window
cfg = [];
cfg.toilim = overall_time_window ;
full_data_all = ft_redefinetrial(cfg,trials);

 

Forward operator

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

%% 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 source 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')



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

 

 

  • No labels