Differences between revisions 94 and 105 (spanning 11 versions)

MEG visual tutorial: Group analysis

Authors: Francois Tadel, Elizabeth Bock.

The aim of this tutorial is to reproduce in the Brainstorm environment the analysis described in the SPM tutorial "Multimodal, Multisubject data fusion". The data processed here consists in simultaneous MEG/EEG recordings of 19 subjects performing a simple visual task on a large number of famous, unfamiliar and scrambled faces. This tutorial follows another page that explains how to process one single subject in details.

Contents

License
Download and installation
Overview of the analysis
1. Coregistration of the acquisition runs
2. Objectives
Subject averages: Famous, Unfamiliar, Scrambled
Subject averages: Faces
Subject averages: Within-subject differences
1. Sources: Faces - Scrambled
2. Sources: Famous - Unfamiliar
Subject averages: Filter and normalize
Subject averages: Screen captures
Group analysis: MEG/EEG
Group analysis: Sources
Scripting

License

This dataset was obtained from the OpenfMRI project (http://www.openfmri.org), accession #ds117. It is made available under the Creative Commons Attribution 4.0 International Public License. Please cite the following reference if you use these data:
Wakeman DG, Henson RN, A multi-subject, multi-modal human neuroimaging dataset, Scientific Data (2015)

Any questions, please contact: rik.henson@mrc-cbu.cam.ac.uk

Download and installation

First, make sure you have enough space on your hard drive, at least 30Gb:

The Brainstorm database with all the data imported, downloaded from this website: 20Gb
The same database after this tutorial: 30Gb

You can follow this tutorial after processing the recordings for the 16 good subjects (6 runs per subject) as illustrated in the single subject tutorial. Otherwise, we provide a Brainstorm protocol that includes all the imported data, downsampled at 275Hz:

Go to the Download page, download the file TutorialGroup.zip (20Gb).
Unzip this file in your Brainstorm database folder (brainstorm_db).
In Brainstorm, menu File > Load protocol > Load from folder > Select brainstorm_db/TutorialGroup

The database you need in order to follow this tutorial should contain the following:

The individual anatomy imported from FreeSurfer or BrainSuite for each subject (16 subjects).
The sensor level averages (MEG+EEG) for each run (downsampled to 275Hz or not).
The forward model, noise covariance and inverse models for each subject and each run.
This protocol TutorialGroup is produced from the single subject protocol TutorialVisual with the script: brainstorm3/toolbox/script/tutorial_visual_copy.m

1 function tutorial_visual_copy(ProtocolNameSingle, ProtocolNameGroup, reports_dir) 2 % TUTORIAL_VISUAL_COPY: Copy the subject averages for the Brainstorm/SPM group tutorial into a new protocol (BIDS VERSION) 3 % 4 % ONLINE TUTORIALS: 5 % - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle 6 % - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup 7 % 8 % INPUTS: 9 % - ProtocolNameSingle : Name of the protocol created with all the data imported (TutorialVisual) 10 % - ProtocolNameGroup : Name of the protocol with just the averages, downsampled to 275Hz (TutorialGroup) 11 % - reports_dir : If defined, exports all the reports as HTML to this folder 12 13 % @============================================================================= 14 % This function is part of the Brainstorm software: 15 % https://neuroimage.usc.edu/brainstorm 16 % 17 % Copyright (c) University of Southern California & McGill University 18 % This software is distributed under the terms of the GNU General Public License 19 % as published by the Free Software Foundation. Further details on the GPLv3 20 % license can be found at http://www.gnu.org/copyleft/gpl.html. 21 % 22 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 23 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 24 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 25 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 26 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 27 % 28 % For more information type "brainstorm license" at command prompt. 29 % =============================================================================@ 30 % 31 % Author: Francois Tadel, Elizabeth Bock, 2016-2018 32 33 % ===== CHECK PROTOCOLS ===== 34 % Start brainstorm without the GUI 35 if ~brainstorm('status') 36 brainstorm nogui 37 end 38 % Output folder for reports 39 if (nargin < 3) || isempty(reports_dir) || ~isdir(reports_dir) 40 reports_dir = []; 41 end 42 % You have to specify the folder in which the tutorial dataset is unzipped 43 if (nargin < 2) || isempty(ProtocolNameSingle) || isempty(ProtocolNameGroup) 44 ProtocolNameSingle = 'TutorialVisual'; 45 ProtocolNameGroup = 'TutorialGroup'; 46 end 47 48 % Output protocol: Delete existing protocol 49 gui_brainstorm('DeleteProtocol', ProtocolNameGroup); 50 % Output protocol: Create new protocol 51 iProtocolGroup = gui_brainstorm('CreateProtocol', ProtocolNameGroup, 0, 0); 52 % Output protocol: Get protocol information 53 ProtocolInfoGroup = bst_get('ProtocolInfo'); 54 55 % Input protocol: Check that it exists 56 iProtocolSingle = bst_get('Protocol', ProtocolNameSingle); 57 if isempty(iProtocolSingle) 58 error(['Unknown protocol: ' ProtocolNameSingle]); 59 end 60 % Select input protocol 61 gui_brainstorm('SetCurrentProtocol', iProtocolSingle); 62 % Input protocol: Get protocol information 63 ProtocolInfoSingle = bst_get('ProtocolInfo'); 64 65 66 % ===== COPY ONLY GOOD SUBJECTS ===== 67 % Start a new report (one report per subject) 68 bst_report('Start'); 69 % Loop on subjects 70 for iSubj = 1:16 71 % Subject folders 72 SubjectName = sprintf('sub-%02d', iSubj); 73 AnatSrc = bst_fullfile(ProtocolInfoSingle.SUBJECTS, SubjectName); 74 DataSrc = bst_fullfile(ProtocolInfoSingle.STUDIES, SubjectName); 75 AnatDest = bst_fullfile(ProtocolInfoGroup.SUBJECTS, SubjectName); 76 DataDest = bst_fullfile(ProtocolInfoGroup.STUDIES, SubjectName); 77 % If subject folder doesn't exist: skip 78 if ~file_exist(AnatSrc) || ~file_exist(DataSrc) 79 disp(['Subject "' SubjectName '" does not exist or is incomplete.']); 80 continue; 81 end 82 % Copy anatomy files 83 mkdir(AnatDest); 84 disp(['Copying: ' AnatSrc ' to ' AnatDest '...']); 85 file_copy(bst_fullfile(AnatSrc, '*.mat'), AnatDest); 86 % Copy analysis folders 87 mkdir(bst_fullfile(DataDest, '@default_study')); 88 mkdir(bst_fullfile(DataDest, '@intra')); 89 disp(['Copying: ' DataSrc ' to ' DataDest '...']); 90 file_copy(bst_fullfile(DataSrc, '@default_study', '*.mat'), bst_fullfile(DataDest, '@default_study')); 91 file_copy(bst_fullfile(DataSrc, '@intra', '*.mat'), bst_fullfile(DataDest, '@intra')); 92 % Loop on runs 93 for iRun = 1:6 94 % Run folders 95 RunName = sprintf('sub-%02d_ses-meg_task-facerecognition_run-%02d_proc-sss_meg_notch', iSubj, iRun); 96 RunSrc = bst_fullfile(DataSrc, RunName); 97 RunDest = bst_fullfile(DataDest, RunName); 98 % If run folder doesn't exist: skip 99 if ~file_exist(RunSrc) 100 disp(['Run "' SubjectName '/' RunName '" does not exist or is incomplete.']); 101 continue; 102 end 103 % Copy files 104 mkdir(RunDest); 105 disp(['Copying: ' RunSrc ' to ' RunDest '...']); 106 file_copy(bst_fullfile(RunSrc, 'brainstormstudy.mat'), RunDest); 107 file_copy(bst_fullfile(RunSrc, 'channel_*.mat'), RunDest); 108 file_copy(bst_fullfile(RunSrc, '*_average_*.mat'), RunDest); 109 if ~isempty(dir(bst_fullfile(RunSrc, 'headmodel_*.mat'))) 110 file_copy(bst_fullfile(RunSrc, 'headmodel_*.mat'), RunDest); 111 end 112 if ~isempty(dir(bst_fullfile(RunSrc, 'noisecov_full.mat'))) 113 file_copy(bst_fullfile(RunSrc, 'noisecov_full.mat'), RunDest); 114 end 115 if ~isempty(dir(bst_fullfile(RunSrc, 'results_*.mat'))) 116 file_copy(bst_fullfile(RunSrc, 'results_*.mat'), RunDest); 117 end 118 if ~isempty(dir(bst_fullfile(RunSrc, 'timefreq_*.mat'))) 119 file_copy(bst_fullfile(RunSrc, 'timefreq_*.mat'), RunDest); 120 end 121 end 122 end 123 124 % ===== RELOAD ===== 125 % Reload output protocol 126 db_reload_database(iProtocolGroup); 127 % Select output protocol 128 gui_brainstorm('SetCurrentProtocol', iProtocolGroup); 129 130 % ===== DOWNSAMPLE TO 275HZ ===== 131 % Process: Select data files in: */* 132 sDataAll = bst_process('CallProcess', 'process_select_files_data', [], []); 133 % Process: Resample: 275Hz 134 sDataAll = bst_process('CallProcess', 'process_resample', sDataAll, [], ... 135 'freq', 275, ... 136 'overwrite', 1); 137 % Process: Select time-frequency files in: */* 138 sTimefreqAll = bst_process('CallProcess', 'process_select_files_timefreq', [], []); 139 % Process: Resample: 275Hz 140 if ~isempty(sTimefreqAll) 141 sTimefreqAll = bst_process('CallProcess', 'process_resample', sTimefreqAll, [], ... 142 'freq', 275, ... 143 'overwrite', 1); 144 end 145 146 % ===== RENAME: DATA ===== 147 % Process: Select data files in: */* 148 sDataAll = bst_process('CallProcess', 'process_select_files_data', [], []); 149 % Rename data files 150 for i = 1:length(sDataAll) 151 % Remove all the processing tags 152 iTag = strfind(sDataAll(i).Comment, ' |'); 153 if isempty(iTag) 154 continue; 155 end 156 newComment = sDataAll(i).Comment(1:iTag-1); 157 % Process: Set comment: AA 158 bst_process('CallProcess', 'process_set_comment', sDataAll(i), [], ... 159 'tag', newComment, ... 160 'isindex', 0); 161 end 162 163 % ===== RENAME: TIME-FREQ ===== 164 % Process: Select time-frequency files in: */* 165 sTimefreqAll = bst_process('CallProcess', 'process_select_files_timefreq', [], []); 166 % Rename timefreq files 167 %AllConditions = {'Famous', 'Scrambled', 'Unfamiliar'}; 168 for i = 1:length(sTimefreqAll) 169 % Remove all the processing tags 170 iTag = strfind(sTimefreqAll(i).Comment, ' |'); 171 if isempty(iTag) 172 continue; 173 end 174 newComment = sTimefreqAll(i).Comment(1:iTag-1); 175 %newComment = ['Avg: ', AllConditions{sTimefreqAll(i).iItem}, ', Power, 6-60Hz']; 176 % Process: Set comment 177 bst_process('CallProcess', 'process_set_comment', sTimefreqAll(i), [], ... 178 'tag', newComment, ... 179 'isindex', 0); 180 end 181 182 % Save report 183 ReportFile = bst_report('Save', []); 184 if ~isempty(reports_dir) && ~isempty(ReportFile) 185 bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolNameGroup '_copy.html'])); 186 end 187 188 189 190 191
If you are not using the package TutorialGroup.zip, make sure you exclude the bad subjects from the analysis: sub001, sub005, sub016. You would also need to add some selection processes to select only the averages (and not the continuous recordings or the individual trials).

Overview of the analysis

Coregistration of the acquisition runs

For each subject, all the runs have been registered to a common head position with MaxFilter. To verify this, select all the channel files within one subject, right-click > Display sensors > MEG (all). The surfaces representing the MEG helmet are perfectly overlapping for all the runs. When the runs are not aligned, it looks like this.

This means we can safely average or compare the MEG sensor values across runs within one subject. However, it is not reliable to average MEG recordings across subjects, because of the anatomical differences between subjects.

This doesn't mean we can estimate the sources only once per subject. We have computed different SSP projectors and selected different bad channels for each acquisition run. To be able to use this information efficiently we should estimate the sources for each run separately, then average the sources across runs.

The forward model is the same for all the runs within one subject, therefore it can be computed for the first run and copied to all the other runs.

Objectives

The objectives for this tutorial are to reproduce the analysis presented in the following documents:

Wakeman DG, Henson RN, A multi-subject, multi-modal human neuroimaging dataset
Scientific Data (2015)
Wakeman DG, Henson RN, Functional and structural connrctivity in face processing: MEG, EEG, fMRI, MRI and DWI data, Biomag 2016
SPM12 manual, chapter 41: Multimodal, Multisubject data fusion

Summary of the results we will compute in this tutorial:

Subject and grand averages for each condition (Famous, Unfamiliar, Scrambled).
Normalization of these averages (Z-score for the sources, or ERS/D for the time-frequency maps).
Projection of the sources results on a template and spatial smoothing of the source maps.
Contrast between faces (famous+unfamiliar) and non-faces (scrambled): difference of averages and significance test.
Contrast between famous faces and unfamiliar faces: difference of averages and significance test.
Sensors of interest: EEG070 (or EEG060 or EEG065 if EEG070 is marked as bad)

Expected effects:

OFA: Bilateral occipital face area
FFA: Bilateral fusiform face area
STS: Right superior temporal sulcus
Around 170ms (N170): Difference between Faces and Scrambled.
After 250ms: Difference between Famous and Unfamiliar.

The methodology that we will follow for computing the averages and the other statistics is described in the tutorial "Workflows".

Subject averages: Famous, Unfamiliar, Scrambled

We will start by computing the subject-level averages for all the data types we have: sensor-level recordings, source maps and time-frequency results. We will use a weighted average to group the results we have for each run (weighted by the number of good trials). We want to compute the averages for each experimental condition separately (famous, unfamiliar, scrambled).

Since the files are already selected in the interface, we will also compute a grand average across subject. In this context, we will consider that each subject has the same weight in the grand average (option "weighted" not selected). Note that it is not accurate to average MEG recordings across subjects, but we can do it just to get a general idea of the group effects (more information).

MEG/EEG

Drag and drop all the subjects in Process1. Select button [Process recordings].
Select process Average > Average files: By trial group (subject average)
Arithmetic average, Weighted
Add process Average > Average files: By trial group (grand average)
Arithmetic average, Not weighted
In output, you get three grand averages in the "Group analysis/Inter-subject" folder (top), and three subject averages per subject in the "intra-subject" folders (bottom).

Sources

Keep all the subjects selected in Process1. Select button [Process sources].
If you run the averaging process directly with this file selection, it would average together the MEG and the EEG sources (36 files selected per subject instead of 18). To average separately the two modalities, we can use the Filter box at the bottom-right corner of the Brainstorm window. Enter "EEG" and it would select only the files with a comment including the tag "EEG".
Select process Average > Average files: By trial group (subject average)
EEG source files only, Arithmetic average, Weighted
Add process File > Add tag: EEG.
Repeat the same steps with the MEG source files (select all subjects, filter, average, tag).
This time, we cannot compute a grand average directly because the source spaces used to estimate the brain activity do not match between subjects. We have to go through an additional step of projection on a template anatomy. We will also wait to have normalized the source maps before displaying them.

Time-frequency

Keep all the subjects selected in Process1. Select button [Process time-freq].
Make sure you reset the Filter search box at the bottom-right of the Brainstorm window.
Run process Average > Average files: By trial group (subject average)
Arithmetic average, Weighted, Match signals
Like for the sources, we will first normalize the time-frequency maps with respect with a baseline before computing grand averages. Like the sensor data: it is valid to average EEG recordings across subjects, but it is not accurate to average MEG sensors across subjects.
We can average the time-frequency files across runs (it is as valid as averaging the MEG or EEG signals), wut we need to select this option "match signals between files using their names" to account for the variable number of bad channels we may have across runs or subjects.

Subject averages: Faces

One of the contrast we want to study is faces (Famous and Unfamiliar) vs non-faces (Scrambled). We need to re-average the Famous and Unfamiliar averages together.

MEG/EEG

In Process2: FilesA = all the Famous subject averages, FilesB = all the Unfamiliar subject averages. The subjects must be in the same order in the two lists because the process will process the files by pairs of files (A,B).
Select process Other > Average A&B: Weighted.
Add process File > Set comment: Comment = "WAvg: Avg: Faces"
Add process Average > Average files: By trial group (subject average): Not weighted
Selecting that many files from the database explorer can be long when done manually, but can be very efficient from a script. Check out the script equivalent at the end of this page.
It creates one new average for each pair (A,B) and computes the grand average for condition Faces.

Sources

In Process2: FilesA = all the Famous subject averages (EEG only) [Process sources]
In Process2: FilesB = all the Unfamiliar subject averages (EEG only) [Process sources]
Make sure the order of the files matches in the two lists.
Select process Other > Average A&B: Weighted.
Add process File > Set comment: Comment = "WAvg: Avg: Faces | EEG"
Repeat the same steps for the MEG sources.

Time-frequency

In Process2: FilesA = all the Famous subject averages (EEG only) [Process timefreq]
In Process2: FilesB = all the Unfamiliar subject averages (EEG only) [Process timefreq]
Select process Other > Average A&B: Weighted.
Add process File > Set comment: Comment = "WAvg: Avg: Faces"
In output: the average for the Face condition, for each subject.

Subject averages: Within-subject differences

To detect correctly the differences between two conditions at the source level, we need to estimate the differences of the conditions for each subject before normalizing them (see tutorial Workflows).

Sources: Faces - Scrambled

In Process2: FilesA = all the Faces subject averages (EEG only) [Process sources]
In Process2: FilesB = all the Scrambled subject averages (EEG only) [Process sources]
Select process Difference > Difference A-B: Do not use absolute values.
Add process File > Set comment: Comment=Faces - Scrambled | EEG
Repeat the same steps for the MEG sources.

Sources: Famous - Unfamiliar

In Process2: FilesA = all the Famous subject averages (EEG only) [Process sources]
In Process2: FilesB = all the Unfamiliar subject averages (EEG only) [Process sources]
Select process Difference > Difference A-B: Do not use absolute values.
Add process File > Set comment: Comment=Famous - Unfamiliar | EEG
Repeat the same steps for the MEG sources.

Subject averages: Filter and normalize

Before comparing the averages across subjects we are going to low-pass filter the signals below 32Hz (to smooth possible latency differences between subjects) and normalize the source and time-frequency values with respect with a baseline (see tutorial Workflows).

MEG/EEG

In Process1, select all the Intra-subject folders from all the subjects, select [Process recordings]. For a faster selection, you can use the view "Functional data (sorted by conditions)".
Select process Pre-process > Band-pass filter: 0Hz-32Hz, MEG,EEG, Mirror, Overwrite
Add process Extract > Extract time: Time window=[-200,900]ms
Two tags are added at the end of the comments of the averaged recordings.

Sources

In Process1, select all the Intra-subject folders from all the subjects, select [Process sources].
Select process Pre-process > Band-pass filter: 0Hz-32Hz, MEG,EEG, Mirror, Overwrite
Add process Extract > Extract time: Time window=[-200,900]ms
Add process Standardize > Baseline normalization: Baseline=[-200,-5]ms, Z-score, Overwrite
Three tags are added at the end of the comments of the averaged sources.

Time-frequency

In Process1, select all the Intra-subject folders from all the subjects, select [Process time-freq].
Run process Standardize > Baseline normalization: Baseline=[-200,-5]ms, ERS/ERD, Overwrite
One tag is added at the end of the comments of the averaged time-frequency maps.

Subject averages: Screen captures

Now we have all the measures ready to be compared across subjects: MEG/EEG, sources, time-frequency. Let's take a few screen captures to make sure the primary visual response look good for all the subjects. The screen captures below represent from left to right:

Registration sensors/MRI: Head surface, MEG (inside of the helmet), EEG+headpoints (green points).
Faces: EEG recordings
Faces: MEG sources at 107ms
Faces: Time-frequency for EEG070

sub002

sub003

sub004

sub006

sub007

sub008

sub009

sub010

sub011

sub012

sub013

sub014

sub015

sub017

sub018

sub019

Group analysis: MEG/EEG

Grand averages

We have already computed the group averages for all the conditions. We will not look much further at these results as we are more interested in the contrasts between conditions. Below are screen captures for all group averages (top=MEG MAG, bottom=EEG). Topography at [50,100,150,200,250,300]ms.

Faces: Famous / Unfamiliar

Scrambled

Faces - Scrambled: Differences of averages

We could compute the contrasts directly from the grand averages, but we will do it from the subject averages because it will be the same file selection as for the statistics.

In Process2: FilesA = all the Faces subject averages (from the Intra-subject folders).
In Process2: FilesB = all the Scrambled subject averages (from the Intra-subject folders).
Run process: Test > Difference of means: Arithmetic average, Not weighted.
Rename the file: Faces - Scrambled (top=MEG MAG, bottom=EEG).

Faces - Scrambled: Significance testing

We have computed the amplitude of the different between the two conditions, and we will try now to estimate a p-value indicating the significance of each of these values. In all the screen captures below: top=MEG MAG, bottom=EEG.

In Process2: Keep the same file selection.
Run process: Test > Parametric test: Paired: All file, All sensors, No average, two-tailed.
Rename the file: Faces - Scrambled: Parametric t-test. Display with p<0.05, FDR-corrected.
You can run other tests in a similar way, with almost identical results.
Process: Test > Permutation test: Paired: All file, All sensors, Paired t-test, 1000 randomizations. Display with p<0.05, FDR-corrected.
Process: Test > FieldTrip: ft_timelockstatistics: All file, EEG, Paired t-test, 1000 randomizations, correction=cluster, cluster alpha=0.05. The cluster-based statistics must be executed on one type of sensors at a time (EEG, MEG MAG, MEG GRAD), because it tries to identify spatio-temporal clusters which groups adjacent sensors.

Famous - Unfamiliar: Differences of averages

In Process2: FilesA = all the Famous subject averages (from the Intra-subject folders).
In Process2: FilesB = all the Unfamiliar subject averages (from the Intra-subject folders).
Run process: Test > Difference of means: Arithmetic average, Not weighted.
Rename the file: Famous - Unfamiliar.

Famous - Unfamiliar: Significance testing

In Process2: Keep the same file selection.
Run process: Test > Parametric test: Paired: All file, All sensors, No average, two-tailed.
Rename the file: Faces - Scrambled: Parametric t-test. Display with p<0.05, FDR-corrected.
Run process: Test > FieldTrip: ft_timelockstatistics: All file, EEG, Paired t-test, 1000 randomizations, correction=cluster, cluster alpha=0.05.

Group analysis: Sources

Projet sources on template

The sources were estimated on the individual anatomy of each subject, the resulting cortical source maps cannot be averaged directly. We need first to re-interpolate all the individual results on a common template (the ICBM152 brain, available in the "default anatomy" folder of the protocol). We also need to extract the absolute values for these source maps: the sign of the minimum norm maps are relative to the orientation of the current with respect to the surface normal, which can vary between subjects.

In Process1, select all the Intra-subject folders from all the subjects, select [Process recordings]. For a faster selection, you can use the view "Functional data (sorted by conditions)".
Select process Pre-process > Absolute values: Sensor types=MEG,EEG, Overwrite.
Add process Sources > Project on default anatomy.
All the 192 subject averages are projected in the same folder "Group analysis / Intra-subject". Re-organize all the projected source files in new folders, one folder per condition.
To create a folder, right-click on Group analysis > New folder. To move a file to the new folder: drag and drop it in the database explorer or use the keyboard shortcuts Ctrl+X/Ctrl+V.

Spatial smoothing

The source maps estimated with constrained orientations can show very focal activity: two adjacent vertices may have very different normals, and therefore very different current values. When averaging multiple subjects, the peaks of activity may not align very well across subjects. Smoothing spatially the source maps may help obtaining better group results.

In Process1, select all the source maps in Group analysis.
Run process Sources > Spatial smoothing:

MEG: mean(|Faces-Scrambled|)

In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG.
Run process Average > Average files: By folder (grand average), Not weighted.
Display with an arbitrary threshold: Z-score>3.
Regions of interest: OFA (Occipital Face Area), FFA (Fusiform Face Area), V1.
This approach identifies correctly the effect size at the subject level. But because the sign is ambiguous across subjects, we need to average the absolute values of the differences. Therefore we are able to identify when and where the two conditions are different, but not which condition leads to a stronger brain response.

MEG: Chi2-test |Faces-Scrambled|=0

In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG.
Run process Test > Parametric test against zero: All file, One-sample Chi2-test two-tailed.
Screen capture: p<0.05 FDR-corrected
This tests the hypotesis H0:[mean_subject(|Faces-Scrambled|) = 0]
This test is appropriate to detect the differences between the two conditions, but cannot determine which condition is stronger.

MEG: Chi2-test log(|Faces-Scrambled|)=0

Right-click on the folder Group_analysis/Faces-Scrambled_MEG > File > Duplicate folder.
Rename the new folder: Faces-Scrambled_MEG_log.
In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG_log.
Select process Pre-process > Run Matlab command > "Data = log(Data)";
Add process Test > Parametric test against zero: All file, One-sample Chi2-test two-tailed.
Screen capture: p<0.05, FDR-corrected

MEG: mean(|Faces|)-mean(|Scrambled|)

In Process2: FilesA = all the source files in Group_analysis/Faces_MEG.
In Process2: FilesB = all the source files in Group_analysis/Scrambled_MEG.
Run process Test > Difference of means: Arithmetic average, Not weighted.
Display with an arbitrary threshold: Z-score>3
Regions of interest: OFA (Occipital Face Area), FFA (Fusiform Face Area), V1.
This approach is not able to detect correctly all the differences between the two conditions. At the subject level,we discarded the sign of each condition, therefore we are not able to

MEG: Student t-test |Faces|=|Scrambled|

In Process2: FilesA = all the source files in Group_analysis/Faces_MEG.
In Process2: FilesB = all the source files in Group_analysis/Scrambled_MEG.
Run process Test > Parametric test: Paired: All file, No average, two-tailed.
Screen capture: p<0.05 FDR-corrected
This test the hypothesis H0:[mean(abs(Faces))=mean(abs(Scrambled))].
This measure is not appropriate to properly identify the differences between the two conditions because it is insensitive to the sign of the sources within each subject.

EEG: Faces-Scrambled

Repeat the same operations with the EEG sources.
Difference of averages: |Faces-Scrambled| (Z-score>5):
Paramtric Chi2 t-test: mean_subject(|Faces-Scrambled|) = 0, p<0.05 FDR-corrected
Paramtric Chi2 t-test: log(mean_subject(|Faces-Scrambled|)) = 0, p<0.05 FDR-corrected
Difference of averages: |Faces|-|Scrambled| (Z-score>3)
Parametric Student t-test: mean(abs(Faces))=mean(abs(Scrambled)), p<0.05 FDR-corrected

Scripting

The following script from the Brainstorm distribution reproduces the analysis presented in this tutorial page: brainstorm3/toolbox/script/tutorial_visual_group.m

function tutorial_visual_group(ProtocolName, reports_dir)
% TUTORIAL_VISUAL_GROUP: Runs the Brainstorm/SPM group analysis pipeline (group analysis, BIDS VERSION)
%
% ONLINE TUTORIALS: 
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup
%
% INPUTS:
%    - ProtocolName : Name of the protocol in which the recordings for the 19 subjects have been imported (TutorialVisual or TutorialGroup)
%                     This folder must include: the run-level averages (recordings and time-frequency) and the source kernels for each run
%    - reports_dir  : If defined, exports all the reports as HTML to this folder

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, Elizabeth Bock, 2016-2018

% ===== CHECK PROTOCOL =====
% Start Brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(ProtocolName)
    ProtocolName = 'TutorialGroup';
end
% Select current protocol
iProtocol = bst_get('Protocol', ProtocolName);
if isempty(iProtocol)
    error(['Unknown protocol: ' ProtocolName]);
end
if (iProtocol ~= bst_get('iProtocol'))
    gui_brainstorm('SetCurrentProtocol', iProtocol);
end
% Process: Select results files in: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch
sFiles = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname', 'sub-01', ...
    'condition',   'sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch');
if isempty(sFiles)
    error(['No source files available in folder: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch.' 10 ...
           'You should run tutorial_visual_single first, or download the ' 10 ...
           'protocol TutorialGroup.zip from the Brainstorm website.']);
end
% Process: Select file comments with tag: Unconstr
sFiles = bst_process('CallProcess', 'process_select_tag', sFiles, [], 'tag', 'Unconstr');
if ~isempty(sFiles)
    isUnconstrained = 1;
else
    isUnconstrained = 0;
end
% Start a new report
bst_report('Start');


% ===== DELETE ALL THE EXISTING FILES =====
GroupSubjectName = bst_get('NormalizedSubjectName');
DirIntra = bst_get('DirAnalysisIntra');
% Delete group analysis
[sSubject, iSubject] = bst_get('Subject', GroupSubjectName);
if ~isempty(sSubject)
    db_delete_subjects(iSubject);
end
% Get all the files in the intra subjects
sAvgSrc = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% If there are existing files: delete the contents of the folders
if ~isempty(sAvgSrc)
    % Get the base folder of the protocol
    ProtocolInfo = bst_get('ProtocolInfo', iProtocol);
    % Get all the folders
    AllIntraDir = unique(cellfun(@bst_fileparts, {sAvgSrc.FileName}, 'UniformOutput', 0));
    % Loop over the intra folders
    for iSubj = 1:length(AllIntraDir)
        % Get all the files
        subjIntraDir = dir(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, '*.mat'));
        subjIntraFiles = setdiff({subjIntraDir.name}, 'brainstormstudy.mat');
        % Delete all the files
        for iFile = 1:length(subjIntraFiles)
            delete(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, subjIntraFiles{iFile}));
        end
    end
    % Reload protocol
    db_reload_database(iProtocol);
end


%% ===== SUBJECT AVERAGES: FAMOUS, UNFAMILIAR, SCRAMBLED ==============================
%  ====================================================================================

% ===== SUBJECT AVERAGE: MEG+EEG =====
% Process: Select data files in: */*/Avg
sAvgRun = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', 'All', ...
    'tag',         'Avg');
% Process: Weighted Average: By trial group (subject average)
sAvgSubj = bst_process('CallProcess', 'process_average', sAvgRun, [], ...
    'avgtype',    6, ...  % By trial group (subject average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   1, ...
    'keepevents', 1);
% Process: Average: By trial group (grand average)
sAvgGroup = bst_process('CallProcess', 'process_average', sAvgSubj, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== SUBJECT AVERAGES: SOURCES (EEG) =====
% Process: Select source files in: */*/EEG
sAvgRunSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'EEG');
% Process: Weighted Average: By trial group (subject average) - EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_average', sAvgRunSrcEeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 0);
% Process: Add tag: EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcEeg, [], ...
    'tag',    'EEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: SOURCES (MEG) =====
% Process: Select source files in: */*/MEG
sAvgRunSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'MEG');
% Process: Weighted Average: By trial group (subject average) - MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_average', sAvgRunSrcMeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcMeg, [], ...
    'tag',    'MEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: TIMEFREQ =====
% Process: Select time-frequency files in: */*
sAvgRunTf = bst_process('CallProcess', 'process_select_files_timefreq', [], []);
% Process: Weighted Average: By trial group (subject average)
sAvgSubjTf = bst_process('CallProcess', 'process_average', sAvgRunTf, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);



%% ===== SUBJECT AVERAGES: FACES ======================================================
%  ====================================================================================

% ===== FACES AVERAGE: MEG+EEG ======
% Process: Select data files in: */*
sAllData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous (6 files)"
sAvgSubjFamous = bst_process('CallProcess', 'process_select_tag', sAllData, [], 'tag', 'WAvg: Avg: Famous (6 files)');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar (6 files)"
sAvgSubjUnfamiliar = bst_process('CallProcess', 'process_select_tag', sAllData, [],  'tag', 'WAvg: Avg: Unfamiliar (6 files)');
% Process: Weighted Average A&B
sAvgSubjFaces = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamous, sAvgSubjUnfamiliar, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFaces = bst_process('CallProcess', 'process_add_tag', sAvgSubjFaces, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFaces = bst_process('CallProcess', 'process_set_comment', sAvgSubjFaces, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);
% Process: Average: By trial group (grand average)
sAvgGroupFaces = bst_process('CallProcess', 'process_average', sAvgSubjFaces, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== FACES AVERAGE: SOURCES (EEG) ======
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcEeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcEeg, [], ...
    'tag',     'WAvg: Avg: Faces | EEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: SOURCES (MEG) ======
% Process: Select source files in: */*
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: Famous
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Famous');
% Process: Select file comments with tag: Unfamiliar
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcMeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcMeg, [], ...
    'tag',     'WAvg: Avg: Faces | MEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: TIMEFREQ ======
% Process: Select time-frequency files in: */*
sAllTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesTf = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousTf, sAvgSubjUnfamiliarTf, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFacesTf = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesTf, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesTf = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesTf, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);



%% ===== SUBJECT AVERAGES: WITHIN-SUBJECT DIFFERENCES =================================
%  ====================================================================================

% ===== FACES - SCRAMBLED: SOURCES EEG =====
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcEeg, sAvgSubjScrambledSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcEeg, [], ...
    'tag',     'Faces - Scrambled | EEG', ...
    'isindex', 0);

% ===== FACES - SCRAMBLED: SOURCES MEG =====
% Process: Select source files in: */*/MEG
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcMeg, sAvgSubjScrambledSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcMeg, [], ...
    'tag',     'Faces - Scrambled | MEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE EEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcEeg, [], ...
    'tag',     'Famous - Unfamiliar | EEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE MEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcMeg, [], ...
    'tag',     'Famous - Unfamiliar | MEG', ...
    'isindex', 0);



%% ===== SUBJECT AVERAGES: FILTER AND NORMALIZE =======================================
%  ====================================================================================

% ===== FILTER: MEG+EEG =====
% Process: Select data files in: */Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgData = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'sensortypes', 'MEG, EEG', ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgData = bst_process('CallProcess', 'process_extract_time', sAvgData, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);

% ===== FILTER AND NORMALIZE: SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgSrc = bst_process('CallProcess', 'process_bandpass', sAvgSrc, [], ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgSrc = bst_process('CallProcess', 'process_extract_time', sAvgSrc, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);
% Process: Z-score transformation: [-202ms,-5ms]
sAvgSrc = bst_process('CallProcess', 'process_baseline_norm', sAvgSrc, [], ...
    'baseline',   [-0.2018, -0.005], ...
    'source_abs', 0, ...
    'method',     'zscore', ...  % Z-score transformation:    x_std = (x - μ) / σ
    'overwrite',  1);

% ===== NORMALIZE: TIMEFREQ =====
% Process: Select time-frequency files in: */Intra
sAvgTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Event-related perturbation (ERS/ERD): [-200ms,-4ms]
sAvgTf = bst_process('CallProcess', 'process_baseline_norm', sAvgTf, [], ...
    'baseline',  [-0.2, -0.004], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - μ) / μ * 100
    'overwrite', 1);

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_1prepare.html']));
end



%% ===== SUBJECT AVERAGES: SCREEN CAPTURES ============================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: local color scale
bst_colormaps('SetMaxMode', 'stat2', 'local');

% ===== FACES: RECORDINGS ======
% Process: Select data files in: */*/WAvg: Avg: Faces |
sAvgFaces = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces |', ...
    'includeintra',  1);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         1, ...  % Sensor/MRI registration
    'modality',       1);     % MEG
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         5, ...  % Recordings time series
    'modality',       4, ...  % EEG
    'time',           0.1055);

% ===== FACES: SOURCES MEG ======
% Process: Select source files in: */*/WAvg: Avg: Faces | MEG
sAvgFacesSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Faces | MEG', ...
    'includeintra',  1);
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sAvgFacesSrcMeg, [], ...
    'target',         8, ...  % Sources (one time)
    'orient',         6, ...  % back
    'time',           0.1055, ...
    'threshold',      10);

% ===== FACES: TIME-FREQ EEG ======
% Process: Select time-frequency files in: */*/WAvg: Avg: Faces
sAvgFacesTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces', ...
    'includeintra',  1);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sAvgFacesTf, [], ...
    'target',         14, ...  % Time-frequency maps
    'time',           0.1055, ...
    'rowname',        'EEG070');

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_2snapshots.html']));
end



%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: global color scale
bst_colormaps('SetMaxMode', 'meg', 'global');
bst_colormaps('SetMaxMode', 'eeg', 'global');
% Set display properties for statistics: p<0.05, FDR-corrected
StatThreshOptions = bst_get('StatThreshOptions');
StatThreshOptions.pThreshold = 0.05;
StatThreshOptions.Correction = 'fdr';
StatThreshOptions.Control    = [1 2 3];
bst_set('StatThreshOptions', StatThreshOptions);

% ===== GRAND AVERAGES =====
% Process: Select data files in: Group_analysis/Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   GroupSubjectName, ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Take screen captures
DataScreenCapture(sAvgData(1), 1, 1, '');
DataScreenCapture(sAvgData(2), 1, 1, '');
DataScreenCapture(sAvgData(3), 1, 1, '');
DataScreenCapture(sAvgData(4), 1, 1, '');
% Define function for screen captures
function DataScreenCapture(sFiles, isMEG, isEEG, Comment)
    if isMEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       3, ...  % MEG (Magnetometers)
            'time',           0.1055, ...
            'Comment',        ['MEG MAG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       3, ...  % MEG (Magnetometers)
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['MEG MAG: ' Comment]);
    end
    if isEEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       4, ...  % EEG
            'time',           0.1055, ...
            'Comment',        ['EEG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       4, ...  % EEG
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['EEG: ' Comment]);
    end
end

% ===== FACES-SCRAMBLED: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFacesData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFacesData = bst_process('CallProcess', 'process_set_comment', sDiffFacesData, [], ...
    'tag',     'Faces - Scrambled', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFacesData, 1, 1, sDiffFacesData.Comment);

% ===== FACES-SCRAMBLED: PARAMETRIC T-TEST =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: t-test paired [ALL]          H0:(A=B), H1:(A<>B)
sStatParamFacesData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFacesData = bst_process('CallProcess', 'process_set_comment', sStatParamFacesData, [], ...
    'tag',     'Faces - Scrambled: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFacesData, 1, 1, sStatParamFacesData.Comment);

% ===== FACES-SCRAMBLED: NON-PARAMETRIC T-TEST =====
% Process: Perm t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatPermFacesData = bst_process('CallProcess', 'process_test_permutation2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',     [], ...
    'sensortypes',    '', ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgrow',         0, ...
    'iszerobad',      0, ...
    'Comment',        '', ...
    'test_type',      'ttest_paired', ...  % Paired Student's t-test T = mean(A-B) / std(A-B) * sqrt(n)
    'randomizations', 1000, ...
    'tail',           'two');  % Two-tailed
% Process: Set comment
sStatPermFacesData = bst_process('CallProcess', 'process_set_comment', sStatPermFacesData, [], ...
    'tag',     'Faces - Scrambled: Permutation t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatPermFacesData, 1, 1, sStatPermFacesData.Comment);

% ===== FACES-SCRAMBLED: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFacesData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFacesData = bst_process('CallProcess', 'process_set_comment', sStatClustFacesData, [], ...
        'tag',     'Faces - Scrambled: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFacesData, 0, 1, sStatClustFacesData.Comment);
end


%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% ===== FAMOUS-UNFAMILIAR: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Famous | "
sSubjAvgFamousData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Famous (6 files)', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Unfamiliar (6 files)"
sSubjAvgUnfamiliarData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Unfamiliar (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFamousData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFamousData = bst_process('CallProcess', 'process_set_comment', sDiffFamousData, [], ...
    'tag',     'Famous - Unfamiliar', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFamousData, 1, 1, sDiffFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: PARAMETRIC T-TEST =====
% Process: t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatParamFamousData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFamousData = bst_process('CallProcess', 'process_set_comment', sStatParamFamousData, [], ...
    'tag',     'Famous - Unfamiliar: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFamousData, 1, 1, sStatParamFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFamousData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFamousData = bst_process('CallProcess', 'process_set_comment', sStatClustFamousData, [], ...
        'tag',     'Famous - Unfamiliar: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFamousData, 0, 1, sStatClustFamousData.Comment);
end

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_3meeg.html']));
end



%% ===== FACES-SCRAMBLED: SOURCES (MEG) ===============================================
%  ====================================================================================
% Start a new report
bst_report('Start');

% ===== PROJECT SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Absolute values / Norm for unconstrained
if isUnconstrained
    % Process: Unconstrained to flat map
    sAvgSrc = bst_process('CallProcess', 'process_source_flat', sAvgSrc, [], ...
        'method', 1);  % Norm: sqrt(x^2+y^2+z^2)
else
    sAvgSrc = bst_process('CallProcess', 'process_absolute', sAvgSrc, [], ...
        'overwrite', 1);
end
% Process: Project on default anatomy
sProjSrc = bst_process('CallProcess', 'process_project_sources', sAvgSrc, [], ...
    'headmodeltype', 'surface');  % Cortex surface
% Process: Spatial smoothing (3.00,abs)
sProjSrc = bst_process('CallProcess', 'process_ssmooth_surfstat', sProjSrc, [], ...
    'fwhm',       3, ...
    'overwrite',  1, ...
    'source_abs', 1);

% ===== MOVE FILES =====
% File tags / destination folder
ListTags = {'Avg: Faces | EEG',                 'Faces_EEG'; ...
            'Avg: Faces | MEG',                 'Faces_MEG'; ...
            'Avg: Famous (6 files) | EEG',      'Famous_EEG'; ...
            'Avg: Famous (6 files) | MEG',      'Famous_MEG'; ...
            'Avg: Unfamiliar (6 files) | EEG',  'Unfamiliar_EEG'; ...
            'Avg: Unfamiliar (6 files) | MEG',  'Unfamiliar_MEG'; ...
            'Avg: Scrambled (6 files) | EEG',   'Scrambled_EEG'; ...
            'Avg: Scrambled (6 files) | MEG',   'Scrambled_MEG'; ...
            'Faces - Scrambled | EEG',          'Faces-Scrambled_EEG'; ...
            'Faces - Scrambled | MEG',          'Faces-Scrambled_MEG'; ...
            'Famous - Unfamiliar | EEG',        'Famous-Unfamiliar_EEG'; ...
            'Famous - Unfamiliar | MEG',        'Famous-Unfamiliar_MEG'};
% Move group by group
for i = 1:size(ListTags,1)
    % Find names of files to move
    iTag = find(~cellfun(@(c)isempty(strfind(c, ListTags{i,1})), {sProjSrc.Comment}));
    if isempty(iTag)
        continue;
    end
    % Process: Move files
    bst_process('CallProcess', 'process_movefile', {sProjSrc(iTag).FileName}, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      ListTags{i,2});
end


% ===== LOOP ON MODALITY =====
AllModalities = {'MEG','EEG'};
for iMod = 1:length(AllModalities)
    % Selected modality for this loop
    Mod = AllModalities{iMod};
    
    % ===== |FACES-SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_', Mod]);
    % Process: Weighted Average: By folder (grand average)
    sAbsDiffFaces = bst_process('CallProcess', 'process_average', sDiffFaces, [], ...
        'avgtype',         4, ...  % By folder (grand average)
        'avg_func',        1, ...  % Arithmetic average:  mean(x)
        'weighted',        0, ...
        'scalenormalized', 0);
    % Process: Set comment
    sAbsDiffFaces = bst_process('CallProcess', 'process_set_comment', sAbsDiffFaces, [], ...
        'tag',     ['mean(|Faces-Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sAbsDiffFaces = bst_process('CallProcess', 'process_movefile', sAbsDiffFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);

    % ===== |FACES-SCRAMBLED|: PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFaces = bst_process('CallProcess', 'process_test_parametric1', sDiffFaces, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFaces = bst_process('CallProcess', 'process_set_comment', sChiParamFaces, [], ...
        'tag',     ['|Faces-Scrambled|=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFaces = bst_process('CallProcess', 'process_movefile', sChiParamFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);

    % ===== LOG(|FACES-SCRAMBLED|): PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Duplicate folders: Add tag "_log"
    sDiffFacesLog = bst_process('CallProcess', 'process_duplicate', sDiffFaces, [], ...
        'target', 2, ...  % Duplicate folders
        'tag',    '_log');
    % Process: Run Matlab command
    sDiffFacesLog = bst_process('CallProcess', 'process_matlab_eval', sDiffFacesLog, [], ...
        'matlab',    'Data = log(Data);', ...
        'overwrite', 1);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFacesLog = bst_process('CallProcess', 'process_test_parametric1', sDiffFacesLog, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFacesLog = bst_process('CallProcess', 'process_set_comment', sChiParamFacesLog, [], ...
        'tag',     ['log(|Faces-Scrambled|)=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFacesLog = bst_process('CallProcess', 'process_movefile', sChiParamFacesLog, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);


    % ===== |FACES|-|SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: Difference of means [abs(mean)]
    sDiffAbsFaces = bst_process('CallProcess', 'process_diff_mean', sAvgFaces, sAvgScrambled, ...
        'avg_func',   2, ...  % Absolute value of average abs(mean(A)) - abs(mean(B))
        'weighted',   0);
    % Process: Set comment
    sDiffAbsFaces = bst_process('CallProcess', 'process_set_comment', sDiffAbsFaces, [], ...
        'tag',     ['mean(|Faces|)-mean(|Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);

    % ===== |FACES|-|SCRAMBLED|: PARAMETRIC T-TEST ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: t-test paired [-202ms,900ms]          H0:(A=B), H1:(A<>B)
    sTtestParamFaces = bst_process('CallProcess', 'process_test_parametric2p', sAvgFaces, sAvgScrambled, ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sTtestParamFaces = bst_process('CallProcess', 'process_set_comment', sTtestParamFaces, [], ...
        'tag',     ['|Faces|=|Scrambled|: Parametric t-test | ' Mod], ...
        'isindex', 0);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);    
end



% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_4sources.html']));
end




end


-  ⇤ ← Revision 94 as of 2016-07-15 19:07:32 → 
  Size: 27440
  Editor: FrancoisTadel
  Comment:
+   ← Revision 105 as of 2016-07-15 21:31:17 → ⇥
  Size: 29260
  Editor: FrancoisTadel
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 4:
-The aim of this tutorial is to reproduce in the Brainstorm environment the analysis described in the SPM tutorial "[[ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.henson/wakemandg_hensonrn/Publications/SPM12_manual_chapter.pdf|Multimodal, Multisubject data fusion]]".  The data processed here consists in simultaneous MEG/EEG recordings of  19 subjects performing simple visual task on a large number of famous,  unfamiliar and scrambled faces. This tutorial follows another page that explains how to process [[Tutorials/VisualSingle|one single subject]] in details.
+The aim of this tutorial is to reproduce in the Brainstorm environment the analysis described in the SPM tutorial "[[ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.henson/wakemandg_hensonrn/Publications/SPM12_manual_chapter.pdf|Multimodal, Multisubject data fusion]]".  The data processed here consists in simultaneous MEG/EEG recordings of  19 subjects performing a simple visual task on a large number of famous,  unfamiliar and scrambled faces. This tutorial follows another page that explains how to process [[Tutorials/VisualSingle|one single subject]] in details.
 Line 19:
-You can follow this tutorial after processing the recordings for the 19 subjects (6 runs per subject) as illustrated in the [[Tutorials/VisualSingle|single subject tutorial]]. Otherwise, we provide a Brainstorm protocol that includes all the imported data, downsampled at 275Hz:
+You can follow this tutorial after processing the recordings for the 16 good subjects (6 runs per subject) as illustrated in the [[Tutorials/VisualSingle|single subject tutorial]]. Otherwise, we provide a Brainstorm protocol that includes all the imported data, downsampled at 275Hz:
 Line 27:
- * The individual anatomy imported from FreeSurfer or BrainSuite for each subject (19 subjects).
+ * The individual anatomy imported from FreeSurfer or BrainSuite for each subject (16 subjects).
 Line 31:
- * If you are not using the package TutorialGroup.zip, make sure you exclude the bad subjects from the analysis: sub001, sub005, sub016. You would also need to add some selection processes to select only the averages (and not the continuous recordings or the trials).
+ * If you are not using the package TutorialGroup.zip, make sure you exclude the bad subjects from the analysis: sub001, sub005, sub016. You would also need to add some selection processes to select only the averages (and not the continuous recordings or the individual trials).
 Line 39:
-This means that we can safely average or compare the MEG sensor values across runs within one subject. However, it is not reliable to average MEG recordings across subjects, because of the [[http://neuroimage.usc.edu/brainstorm/Tutorials/Workflows#MEG_recordings|anatomical differences between subjects]].

This doesn't mean that we can estimate the sources only once per subject. We have computed different SSP projectors and different bad channels for the each acquisition run. To be able to use this information efficiently we should estimate the sources for the trial averages for each run separately, then average the sources across runs.
+This means we can safely average or compare the MEG sensor values across runs within one subject. However, it is not reliable to average MEG recordings across subjects, because of the [[http://neuroimage.usc.edu/brainstorm/Tutorials/Workflows#MEG_recordings|anatomical differences between subjects]].

This doesn't mean we can estimate the sources only once per subject. We have computed different SSP projectors and selected different bad channels for each acquisition run. To be able to use this information efficiently we should estimate the sources for each run separately, then average the sources across runs.
 Line 45:
-=== Objectives [TODO] ===
+=== Objectives ===
 Line 56:
- * Contrast between faces (Famous+Unfamiliar) and non-faces (Scrambled).
 * Contrast between famous faces and unfamiliar faces.
 * Two contrast measures: first a difference of grand averages, then a t-test across subjects.
+ * Projection of the sources results on a template and spatial smoothing of the source maps.
 * Contrast between faces (famous+unfamiliar) and non-faces (scrambled): difference of averages and significance test.
 * Contrast between famous faces and unfamiliar faces: difference of averages and significance test.
 Line 69:
-The methodology that we will follow for computing the averages and the other statistics are described in the tutorial "[[Tutorials/Workflows|Workflows]]".
+The methodology that we will follow for computing the averages and the other statistics is described in the tutorial "[[Tutorials/Workflows|Workflows]]".
 Line 74:
-Since the files are already selected in the interface, we will also compute a grand average across subject. In this context, we will consider that each subject has the same weight in the gran average (option "weighted" not selected). Note that it is not accurate to average MEG recordings across subjects, but we can do it just to get a general idea ([[http://neuroimage.usc.edu/brainstorm/Tutorials/Workflows#MEG_recordings|more information]]).
+Since the files are already selected in the interface, we will also compute a grand average across subject. In this context, we will consider that each subject has the same weight in the grand average (option "weighted" not selected). Note that it is not accurate to average MEG recordings across subjects, but we can do it just to get a general idea of the group effects ([[http://neuroimage.usc.edu/brainstorm/Tutorials/Workflows#MEG_recordings|more information]]).
 Line 86:
- * Add process '''File > Set tag''': '''EEG'''.<<BR>><<BR>> {{attachment:average_process_sources.gif||height="504",width="603"}}
+ * Add process '''File > Add tag''': '''EEG'''.<<BR>><<BR>> {{attachment:average_process_sources.gif||height="504",width="603"}}
 Line 88:
- * This time, we cannot compute a grand average directly because the source spaces used to estimate the brain activity do not match between subjects. We will have to go through an additional step of projection on a template anatomy. We will also wait to have normalized the source maps before displaying them.<<BR>><<BR>> {{attachment:average_files_sources.gif||height="208",width="243"}}
+ * This time, we cannot compute a grand average directly because the source spaces used to estimate the brain activity do not match between subjects. We have to go through an additional step of projection on a template anatomy. We will also wait to have normalized the source maps before displaying them.<<BR>><<BR>> {{attachment:average_files_sources.gif||height="208",width="243"}}
 Line 93:
- * Like for the sources, we will first normalize the time-frequency maps with respect with a baseline before computing grand averages.
 * We can average the time-frequency files across runs (it is as valid as averaging the MEG or EEG signals across subjects). But we need to select this option "match signals between files using their names" to account for the variable number of bad channels we may have across runs or subjects.  <<BR>><<BR>> {{attachment:average_files_tf.gif||height="250",width="633"}}
+ * Like for the sources, we will first normalize the time-frequency maps with respect with a baseline before computing grand averages. Like the sensor data: it is valid to average EEG recordings across subjects, but it is not accurate to average MEG sensors across subjects.
 * We can average the time-frequency files across runs (it is as valid as averaging the MEG or EEG signals), wut we need to select this option "match signals between files using their names" to account for the variable number of bad channels we may have across runs or subjects.  <<BR>><<BR>> {{attachment:average_files_tf.gif||height="250",width="633"}}
 Line 102:
- * Selecting that many files from the database explorer can be long when done manually, but can be done very efficiently from a script. Check out the script equivalent at the end of this page.
 * It creates one new average for each pair (A,B) and compute the grand average for condition Faces. <<BR>><<BR>> {{attachment:average_faces_data2.gif||height="234",width="299"}}
+ * Selecting that many files from the database explorer can be long when done manually, but can be very efficient from a script. Check out the script equivalent at the end of this page.
 * It creates one new average for each pair (A,B) and computes the grand average for condition Faces. <<BR>><<BR>> {{attachment:average_faces_data2.gif||height="234",width="299"}}
 Line 214:
-=== Faces - Scrambled: Significance testing [TODO] ===
+=== Faces - Scrambled: Significance testing ===
 Line 228:
- * Rename the file: '''Famous - Unfamiliar'''. <<BR>><<BR>> {{attachment:diff_data_famous.gif||height="303",width="450"}}

=== Famous - Unfamiliar: Significance testing [TODO] ===
+ * Rename the file: '''Famous - Unfamiliar'''. <<BR>><<BR>> {{attachment:diff_data_famous.gif||width="450"}}

=== Famous - Unfamiliar: Significance testing ===
 Line 235:
- * Run process: '''Test > FieldTrip: ft_timelockstatistics''': All file, EEG, Paired t-test, 1000 randomizations, correction=cluster, cluster alpha=0.05. <<BR>><<BR>> {{attachment:cluster_data_famous.gif}}
+ * Run process: '''Test > FieldTrip: ft_timelockstatistics''': All file, EEG, Paired t-test, 1000 randomizations, correction=cluster, cluster alpha=0.05. <<BR>><<BR>> {{attachment:cluster_data_famous.gif||width="450"}}
 Line 247:
-=== Spatial smoothing [TODO] ===
+=== Spatial smoothing ===
 Line 251:
- * Run process Sources > Spatial smoothing: <<BR>><<BR>> {{attachment:sources_smooth.gif}}

=== MEG: mean(|Faces-Scrambled|) [TODO] ===
+ * Run process Sources > Spatial smoothing: <<BR>><<BR>> {{attachment:sources_smooth.gif||height="348",width="504"}}

=== MEG: mean(|Faces-Scrambled|) ===
 Line 255:
- * Run process '''Average > Average files''': By folder (grand average), Not weighted <<BR>><<BR>> {{attachment:sources_absdiff_meg.gif}}
+ * Run process '''Average > Average files''': By folder (grand average), Not weighted. <<BR>>Display with an arbitrary threshold: Z-score>3. <<BR>><<BR>> {{attachment:sources_absdiff_meg.gif||width="300"}}  {{attachment:sources_absdiff_meg_right.gif||width="300"}}
 Line 259:
-=== MEG: mean(|Faces|)-mean(|Scrambled|) [TODO] ===
+=== MEG: Chi2-test |Faces-Scrambled|=0 ===
 * In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG.
 * Run process '''Test > Parametric test against zero''': All file, One-sample Chi2-test two-tailed.
 * Screen capture: p<0.05 FDR-corrected<<BR>><<BR>> {{attachment:sources_chi2_meg.gif||width="300"}} {{attachment:sources_chi2_meg_right.gif||width="300"}}
 * This tests the hypotesis H0:[mean_subject(|Faces-Scrambled|) = 0]<<BR>>This  test is appropriate to detect the differences between the two  conditions, but cannot determine which condition is stronger.

=== MEG: Chi2-test log(|Faces-Scrambled|)=0 ===
 * Right-click on the folder Group_analysis/Faces-Scrambled_MEG > '''File > Duplicate folder'''. <<BR>><<BR>> {{attachment:duplicate_folder.gif||height="235",width="351"}}
 * Rename the new folder: '''Faces-Scrambled_MEG_log'''.
 * In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG_log.
 * Select process '''Pre-process > Run Matlab command > "Data = log(Data)";''' <<BR>> Add process '''Test > Parametric test against zero''': All file, One-sample Chi2-test two-tailed. <<BR>><<BR>> {{attachment:log_process.gif||height="378",width="654"}}
 * Screen capture: p<0.05, FDR-corrected<<BR>><<BR>>  {{attachment:sources_chi2log_meg.gif||width="300"}} {{attachment:sources_chi2log_meg_right.gif||width="300"}}

=== MEG: mean(|Faces|)-mean(|Scrambled|) ===
-Line 262:
+Line 275:
- * Run process '''Test > Difference of means''': Arithmetic average, __Not__ weighted. <<BR>><<BR>> {{attachment:sources_diffmean_meg.gif}}
+ * Run process '''Test > Difference of means''': Arithmetic average, __Not__ weighted. <<BR>>Display with an arbitrary threshold: Z-score>3<<BR>><<BR>> {{attachment:sources_diffmean_meg.gif||width="300"}}   {{attachment:sources_diffmean_meg_right.gif||width="300"}}
-Line 266:
+Line 279:
-=== MEG: Student t-test |Faces|=|Scrambled| [TODO] ===
+=== MEG: Student t-test |Faces|=|Scrambled| ===
-Line 270:
+Line 283:
- * Screen capture: p<0.05 FDR-corrected <<BR>><<BR>> {{attachment:sources_ttest_meg.gif}}
+ * Screen capture: p<0.05 FDR-corrected <<BR>><<BR>> {{attachment:sources_ttest_meg.gif||width="300"}} {{attachment:sources_ttest_meg_right.gif||width="300"}}
-Line 273:
+Line 286:
-=== MEG: Chi2-test |Faces-Scrambled|=0 ===
 * In Process1, select all the source files in Group_analysis/Faces-Scrambled_MEG.
 * Run process '''Test > Parametric test against zero''': All file, One-sample Chi2-test two-tailed.
 * Screen capture: p<0.05 FDR-corrected<<BR>><<BR>> {{attachment:sources_chi2_meg.gif}}
 * This tests the hypotesis H0:[mean_subject(|Faces-Scrambled|) = 0]<<BR>>This test is appropriate to detect the differences between the two conditions, but cannot determine which condition is stronger.
 * '''ERROR''': There is something wrong here, because everything is significant even with a Bonferroni correction for multiple comparisons. Either the test is not adapted, or it's not coded properly.

=== Faces-Scrambled (EEG) [TODO] ===
+=== EEG: Faces-Scrambled ===
-Line 282:
+Line 288:
- * Difference of averages: <<BR>><<BR>> {{attachment:sources_diffmean_eeg.gif}}
 * Parametric Student t-test: [mean(abs(Faces))=mean(abs(Scrambled))], p<0.05 FDR-corrected <<BR>><<BR>> {{attachment:sources_ttest_eeg.gif}}
 * Paramtric Chi2 t-test: [mean_subject(|Faces-Scrambled|) = 0], p<0.05 FDR-corrected <<BR>><<BR>> {{attachment:sources_chi2_eeg.gif}}
+ * Difference of averages: |Faces-Scrambled| (Z-score>5): <<BR>>  {{attachment:sources_diffmean_eeg.gif||width="300"}} {{attachment:sources_diffmean_eeg_right.gif||width="300"}}
 * Paramtric Chi2 t-test: mean_subject(|Faces-Scrambled|) = 0, p<0.05 FDR-corrected <<BR>> {{attachment:sources_chi2_eeg.gif||width="300"}} {{attachment:sources_chi2_eeg_right.gif||width="300"}}
 * Paramtric Chi2 t-test: '''log'''(mean_subject(|Faces-Scrambled|)) = 0, p<0.05 FDR-corrected <<BR>> {{attachment:sources_chi2log_eeg.gif||width="300"}} {{attachment:sources_chi2log_eeg_right.gif||width="300"}}
 * Difference of averages: |Faces|-|Scrambled| (Z-score>3) <<BR>> {{attachment:sources_absdiff_eeg.gif|sources_ttest_eeg.gif|width="300"}} {{attachment:sources_absdiff_eeg_right.gif|sources_ttest_eeg_right.gif|width="300"}}
 * Parametric Student t-test: mean(abs(Faces))=mean(abs(Scrambled)), p<0.05 FDR-corrected <<BR>> {{attachment:sources_ttest_eeg.gif||width="300"}} {{attachment:sources_ttest_eeg_right.gif||width="300"}}

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):