MoAE#
/demos/MoAE
├── models
└── options
This “Mother of All Experiments” is based on the block design dataset of SPM.
In the options folder has several examples of how to encode the options of your
analysis in a json file.
In the models shows the BIDS statistical model used to run the GLM of this demo.
preprocessing + stats#
- moae_01_bids_app#
# MoAE demo
This script shows how to use the bidspm BIDS app
Download
download the dataset from the FIL for the block design SPM tutorial
Preprocessing
copies the necessary data from the raw to the derivative folder,
runs spatial preprocessing
those are otherwise handled by the workflows:
Stats
This will run the subject level GLM and contrasts on it of the MoaE dataset
GLM specification + estimation
compute contrasts
show results
that are otherwise handled by the workflows
Note
Results might be a bit different from those in the SPM manual as some default options are slightly different in this pipeline (e.g use of FAST instead of AR(1), motion regressors added)
type bidspm help or bidspm(‘action’, ‘help’) or see this page: https://bidspm.readthedocs.io/en/stable/bids_app_api.html for more information on what parameters are obligatory or optional
script
this_dir = fileparts(mfilename('fullpath'));
addpath(fullfile(this_dir, '..', '..'));
bidspm();
% ## Download the dataset
download_data = true;
clean = false;
download_moae_ds(download_data, clean);
% If the `tree` command is installed on your computer, you view it:
system('tree inputs/raw');
% ## Preprocessing
%
% This will run:
%
% - copy the input dataset into a derivative one
% - write a summary description of the data set
% - do slice time correction (if not ignored and if slice timing is specified)
% - realign the functional data (and apply unwarping - if not ignored)
% - coregister the functional to the anatomical one
% - segmentation the anatomical data
% - skullstripping the anatomical data and creation of brain mask in native space
% - normalization to SPM MNI space (IXI549Space)
% - smooth the data
% You can type `bidspm help` to get more info
% on the arguments and parameters needed by the bidspm app.
%
% But it follows the general pattern of any bidsapp:
%
% ```matlab
% bidspm(bids_dir, output_dir, analysis_level, ...)
% ```
bidspm help;
% where the raw bids data is;
bids_dir = fullfile(this_dir, 'inputs', 'raw');
% where we want to output it;
% the data will be saved there in bidspm-preproc subfolder;
output_dir = fullfile(this_dir, 'outputs', 'derivatives');
% the subject we want to analyse;
subject_label = '01';
VERBOSITY = 3;
bidspm(bids_dir, output_dir, 'subject', ...
'participant_label', {subject_label}, ...
'action', 'preprocess', ...
'task', 'auditory', ...
'ignore', {'unwarp', 'slicetiming'}, ...
'space', {'IXI549Space'}, ...
'fwhm', 6, ...
'verbosity', VERBOSITY);
% ## Stats
% for the stats we need to specify where the preprocessed data is;
preproc_dir = fullfile(output_dir, 'bidspm-preproc');
% ### BIDS stats model
%
% The model specification as well as the contrasts to compute
% are defined in a BIDS stats model:
% https://bids-standard.github.io/stats-models/
model_file = fullfile(pwd, 'models', 'model-MoAE_smdl.json');
system('cat models/model-MoAE_smdl.json');
% ### Specify the result to show
%
% Running bidspm for the stats will perform:
%
% - model specification and estimation
% - contrasts computation
% - displaying the results
%
% Hence we need to specify in the options which results
% we want to view and how we want to save it.
%
% The results of a given contrat can be saved as:
% - an png image
% - a SPM montage of slices
% - a thresholded statistical map
% - a binary mask
% - an NIDM results zip file
% - a table of labelled activations
bidspm(bids_dir, output_dir, 'subject', ...
'participant_label', {subject_label}, ...
'action', 'stats', ...
'preproc_dir', preproc_dir, ...
'model_file', model_file, ...
'space', {'IXI549Space'}, ...
'concatenate', false, ...
'fwhm', 6, ...
'verbosity', VERBOSITY);
BIDS stats model
{
"Name": "auditory",
"BIDSModelVersion": "1.0.0",
"Description": "GLM and contrasts to compute for the FIL MoAE dataset",
"Input": {
"task": [
"auditory"
],
"space": [
"IXI549Space"
]
},
"Nodes": [
{
"Level": "Run",
"Name": "run_level",
"GroupBy": [
"run",
"subject"
],
"Model": {
"X": [
"trial_type.listening",
"trans_?",
"rot_?"
],
"HRF": {
"Variables": [
"trial_type.listening"
],
"Model": "spm"
},
"Type": "glm",
"Options": {
"HighPassFilterCutoffHz": 0.0078,
"Mask": {
"suffix": [
"mask"
],
"desc": [
"brain"
]
}
},
"Software": {
"SPM": {
"SerialCorrelation": "AR(1)",
"InclusiveMaskingThreshold": 0.8
},
"bidspm": {
"Results": [
{
"name": [
"listening_1"
],
"p": 0.05,
"MC": "FWE",
"png": true,
"binary": true,
"nidm": true,
"montage": {
"do": true,
"slices": [
-4,
0,
4,
8,
16
],
"background": {
"suffix": "T1w",
"desc": "preproc",
"modality": "anat"
}
}
},
{
"name": [
"listening_inf_baseline"
],
"p": 0.01,
"k": 10,
"MC": "none",
"csv": true,
"atlas": "AAL"
}
]
}
}
},
"DummyContrasts": {
"Test": "t",
"Contrasts": [
"trial_type.listening"
]
},
"Contrasts": [
{
"Name": "listening_inf_baseline",
"ConditionList": [
"trial_type.listening"
],
"Weights": [
-1
],
"Test": "t"
}
]
}
]
}
stats with fmriprep output + default stats model#
- moae_fmriprep#
This script will run the FFX and contrasts on it of the MoAE dataset using the fmriprep preprocessed data
If you want to get the preprocessed data and you have datalad on your computer you can run the following commands to get the necessary data:
datalad install --source git@gin.g-node.org:/SPM_datasets/spm_moae_fmriprep.git \ inputs/fmriprep cd inputs/fmriprep && datalad get *.json \ */*/*tsv \ */*/*json \ */*/*desc-preproc*.nii.gz \ */*/*desc-brain*.nii.gz
Otherwise you also grab the data from OSF: https://osf.io/vufjs/download
script
% (C) Copyright 2019 Remi Gau
clear;
clc;
addpath(fullfile(pwd, '..', '..'));
bidspm();
DOWNLOAD_DATA = true;
SMOOTH = true;
% Getting the raw dataset to get the events.tsv
download_moae_ds(DOWNLOAD_DATA);
%%
subject_label = {'01'};
task = {'auditory'};
space = {'MNI152NLin6Asym'};
WD = fileparts(mfilename('fullpath'));
bids_dir = fullfile(WD, 'inputs', 'raw');
fmriprep_dir = fullfile(WD, 'inputs', 'fmriprep');
% we need to specify where the smoothed data will go
output_dir = fullfile(WD, 'outputs', 'derivatives');
%% Copy and smooth
%
% fmriprep data is not smoothed so we need to do that first
%
if SMOOTH
bidspm(fmriprep_dir, output_dir, 'subject', ...
'action', 'smooth', ...
'participant_label', subject_label, ...
'task', task, ...
'space', space, ...
'fwhm', 8, ...
'verbosity', 3);
end
%% STATS
% create default model
bidspm(bids_dir, output_dir, 'dataset', ...
'action', 'default_model', ...
'participant_label', subject_label, ...
'task', task, ...
'space', space, ...
'ignore', {'contrasts', 'transformations', 'dataset'}, ...
'verbosity', 3);
% run model
model_file = fullfile(output_dir, 'models', 'model-defaultAuditory_smdl.json');
% Specify the result to show
%
% nodeName corresponds to the name of the Node in the BIDS stats model
opt.results(1).nodeName = 'run';
% this results corresponds to the name of the contrast in the BIDS stats model
opt.results(1).name = {'listening_1'};
% Specify how you want your output
% (all the following are on false by default)
opt.results(1).png = true();
opt.results(1).nidm = true();
opt.results(1).csv = true();
opt.results(1).threshSpm = true();
opt.results(1).binary = true();
opt.results(1).montage.do = true();
opt.results(1).montage.background = struct('suffix', 'T1w', ...
'desc', 'preproc', ...
'modality', 'anat');
opt.results(1).montage.slices = -4:2:16;
% the following "bids app" runs:
%
% - GLM specification + estimation,
% - compute contrasts and
% - show results
%
% that are otherwise handled by the bidsFFX.m and bidsResults.m workflows
%
% type bidspm('action', 'help')
% or see this page: https://bidspm.readthedocs.io/en/stable/bids_app_api.html
% for more information on what parameters are obligatory or optional
%
% where the smooth data is
preproc_dir = fullfile(output_dir, 'bidspm-preproc');
bidspm(bids_dir, output_dir, 'subject', ...
'participant_label', subject_label, ...
'action', 'stats', ...
'preproc_dir', preproc_dir, ...
'model_file', model_file, ...
'roi_atlas', 'hcpex', ...
'space', space, ...
'fwhm', 8, ...
'options', opt, ...
'verbosity', 3);
region of interest#
- moae_02_create_roi_extract_data#
This script shows how to create a ROI and extract data from it.
Warning
This is “double dipping” as we use the same data to create the ROI we are going to extract the value from.
script
clc;
addpath(fullfile(pwd, '..', '..'));
bidspm();
subLabel = '01';
opt.dir.derivatives = fullfile(fileparts(mfilename('fullpath')), 'outputs', 'derivatives');
opt.dir.roi = fullfile(opt.dir.derivatives, 'bidspm-roi');
opt.dir.stats = fullfile(opt.dir.derivatives, 'outputs', 'bidspm-stats');
opt.model.file = fullfile(fileparts(mfilename('fullpath')), ...
'models', 'model-MoAE_smdl.json');
opt.pipeline.type = 'stats';
opt = checkOptions(opt);
%% Get the con image to extract data
% we can do this by using the "label-XXXX"
% from the mask created by the contrast estimation in the previous step
ffxDir = getFFXdir(subLabel, opt);
maskImage = spm_select('FPList', ffxDir, '^.*_mask.nii$');
bf = bids.File(spm_file(maskImage, 'filename'));
conImage = spm_select('FPList', ffxDir, ['^con_' bf.entities.label '.nii$']);
%% Create ROI right auditory cortex
saveROI = true;
sphere.radius = 3;
sphere.maxNbVoxels = 200;
sphere.location = [57 -22 11];
specification = struct( ...
'mask1', maskImage, ...
'mask2', sphere);
output_dir = fullfile(opt.dir.roi, ['sub-' subLabel], 'roi');
spm_mkdir(output_dir);
[~, roiFile] = createRoi('expand', specification, conImage, output_dir, saveROI);
bf = bids.File(roiFile);
% rename mask image:
% add a description and remove a whole bunch of entities
% from the original name
newname.entities.hemi = 'R';
newname.entities.desc = 'auditory cortex';
newname.entities.task = '';
newname.entities.label = '';
newname.entities.p = '';
newname.entities.k = '';
newname.entities.MC = '';
newname.entity_order = {'sub', 'hemi', 'space', 'desc'};
rightRoiFile = bf.rename('spec', newname, 'dry_run', false, 'verbose', true);
%% same but with left hemisphere
sphere.location = [-57 -22 11];
specification = struct('mask1', maskImage, ...
'mask2', sphere);
[~, roiFile] = createRoi('expand', specification, conImage, output_dir, saveROI);
bf = bids.File(roiFile);
newname.entities.hemi = 'L';
leftRoiFile = bf.rename('spec', newname, 'dry_run', false, 'verbose', true);
%%
right_data = spm_summarise(conImage, fullfile(output_dir, rightRoiFile.filename));
left_data = spm_summarise(conImage, fullfile(output_dir, leftRoiFile.filename));
BIDS stats model
{
"Name": "auditory",
"BIDSModelVersion": "1.0.0",
"Description": "GLM and contrasts to compute for the FIL MoAE dataset",
"Input": {
"task": [
"auditory"
],
"space": [
"IXI549Space"
]
},
"Nodes": [
{
"Level": "Run",
"Name": "run_level",
"GroupBy": [
"run",
"subject"
],
"Model": {
"X": [
"trial_type.listening",
"trans_?",
"rot_?"
],
"HRF": {
"Variables": [
"trial_type.listening"
],
"Model": "spm"
},
"Type": "glm",
"Options": {
"HighPassFilterCutoffHz": 0.0078,
"Mask": {
"suffix": [
"mask"
],
"desc": [
"brain"
]
}
},
"Software": {
"SPM": {
"SerialCorrelation": "AR(1)",
"InclusiveMaskingThreshold": 0.8
},
"bidspm": {
"Results": [
{
"name": [
"listening_1"
],
"p": 0.05,
"MC": "FWE",
"png": true,
"binary": true,
"nidm": true,
"montage": {
"do": true,
"slices": [
-4,
0,
4,
8,
16
],
"background": {
"suffix": "T1w",
"desc": "preproc",
"modality": "anat"
}
}
},
{
"name": [
"listening_inf_baseline"
],
"p": 0.01,
"k": 10,
"MC": "none",
"csv": true,
"atlas": "AAL"
}
]
}
}
},
"DummyContrasts": {
"Test": "t",
"Contrasts": [
"trial_type.listening"
]
},
"Contrasts": [
{
"Name": "listening_inf_baseline",
"ConditionList": [
"trial_type.listening"
],
"Weights": [
-1
],
"Test": "t"
}
]
}
]
}
slice display#
- moae_03_slice_display#
This script shows how to display the results of a GLM by having on the same image:
the beta estimates
the t statistics
ROI contours
script
clear;
close all;
clc;
addpath(fullfile(pwd, '..', '..'));
bidspm();
this_dir = fileparts(mfilename('fullpath'));
subLabel = '01';
opt.pipeline.type = 'stats';
opt.dir.raw = fullfile(this_dir, 'inputs', 'raw');
opt.dir.derivatives = fullfile(this_dir, 'outputs', 'derivatives');
opt.dir.preproc = fullfile(opt.dir.derivatives, 'bidspm-preproc');
opt.dir.roi = fullfile(opt.dir.derivatives, 'bidspm-roi');
opt.dir.stats = fullfile(opt.dir.derivatives, 'bidspm-stats');
opt.model.file = fullfile(this_dir, 'models', 'model-MoAE_smdl.json');
% Specify the result to compute
opt.results(1).nodeName = 'run_level';
opt.results(1).name = 'listening';
% MONTAGE FIGURE OPTIONS
opt.results(1).montage.do = true();
opt.results(1).montage.slices = -4:2:16; % in mm
% axial is default 'sagittal', 'coronal'
opt.results(1).montage.orientation = 'axial';
% will use the MNI T1 template by default but the underlay image can be changed.
opt.results(1).montage.background = ...
fullfile(spm('dir'), 'canonical', 'avg152T1.nii');
opt = checkOptions(opt);
use_schema = false;
BIDS_ROI = bids.layout(opt.dir.roi, 'use_schema', use_schema);
filter = struct('sub', subLabel, ...
'hemi', 'R', ...
'desc', 'auditoryCortex');
rightRoiFile = bids.query(BIDS_ROI, 'data', filter);
filter.hemi = 'L';
leftRoiFile = bids.query(BIDS_ROI, 'data', filter);
% we get the con image to extract data
ffxDir = getFFXdir(subLabel, opt);
maskImage = spm_select('FPList', ffxDir, '^.*_mask.nii$');
bf = bids.File(spm_file(maskImage, 'filename'));
conImage = spm_select('FPList', ffxDir, ['^con_' bf.entities.label '.nii$']);
%% Layers to put on the figure
layers = sd_config_layers('init', {'truecolor', 'dual', 'contour', 'contour'});
% Layer 1: Anatomical map
[anat_normalized_file, anatRange] = return_normalized_anat_file(opt, subLabel);
layers(1).color.file = anat_normalized_file;
layers(1).color.range = [0 anatRange(2)];
layers(1).color.map = gray(256);
%% Layer 2: Dual-coded layer
%
% - contrast estimates color-coded;
layers(2).color.file = conImage;
color_map_folder = fullfile(fileparts(which('map_luminance')), '..', 'mat_maps');
load(fullfile(color_map_folder, 'diverging_bwr_iso.mat'));
layers(2).color.map = diverging_bwr;
layers(2).color.range = [-4 4];
layers(2).color.label = '\beta_{listening} - \beta_{baseline} (a.u.)';
%% Layer 2: Dual-coded layer
%
% - t-statistics opacity-coded
spmTImage = spm_select('FPList', ffxDir, ['^spmT_' bf.entities.label '.nii$']);
layers(2).opacity.file = spmTImage;
layers(2).opacity.range = [2 3];
layers(2).opacity.label = '| t |';
%% Layer 3 and 4: Contour of ROI
layers(3).color.file = rightRoiFile{1};
layers(3).color.map = [0 0 0];
layers(3).color.line_width = 2;
layers(4).color.file = leftRoiFile{1};
layers(4).color.map = [1 1 1];
layers(4).color.line_width = 2;
%% Settings
settings = sd_config_settings('init');
% we reuse the details for the SPM montage
settings.slice.orientation = opt.results(1).montage.orientation;
settings.slice.disp_slices = -15:3:18;
settings.fig_specs.n.slice_column = 4;
settings.fig_specs.title = opt.results(1).name;
%% Display the layers
[settings, p] = sd_display(layers, settings);
BIDS stats model
{
"Name": "auditory",
"BIDSModelVersion": "1.0.0",
"Description": "GLM and contrasts to compute for the FIL MoAE dataset",
"Input": {
"task": [
"auditory"
],
"space": [
"IXI549Space"
]
},
"Nodes": [
{
"Level": "Run",
"Name": "run_level",
"GroupBy": [
"run",
"subject"
],
"Model": {
"X": [
"trial_type.listening",
"trans_?",
"rot_?"
],
"HRF": {
"Variables": [
"trial_type.listening"
],
"Model": "spm"
},
"Type": "glm",
"Options": {
"HighPassFilterCutoffHz": 0.0078,
"Mask": {
"suffix": [
"mask"
],
"desc": [
"brain"
]
}
},
"Software": {
"SPM": {
"SerialCorrelation": "AR(1)",
"InclusiveMaskingThreshold": 0.8
},
"bidspm": {
"Results": [
{
"name": [
"listening_1"
],
"p": 0.05,
"MC": "FWE",
"png": true,
"binary": true,
"nidm": true,
"montage": {
"do": true,
"slices": [
-4,
0,
4,
8,
16
],
"background": {
"suffix": "T1w",
"desc": "preproc",
"modality": "anat"
}
}
},
{
"name": [
"listening_inf_baseline"
],
"p": 0.01,
"k": 10,
"MC": "none",
"csv": true,
"atlas": "AAL"
}
]
}
}
},
"DummyContrasts": {
"Test": "t",
"Contrasts": [
"trial_type.listening"
]
},
"Contrasts": [
{
"Name": "listening_inf_baseline",
"ConditionList": [
"trial_type.listening"
],
"Weights": [
-1
],
"Test": "t"
}
]
}
]
}