From ead41dccd7ab12da572a1ca3da45a8b9ca3c8434 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Wed, 25 May 2022 12:28:40 +0200 Subject: [PATCH 1/4] deal with parametric modulation --- demos/face_repetition/face_rep_02_stats.m | 2 +- .../face_rep_02_stats_parametric.m | 28 +++++ .../model-faceRepetitionParametric_smdl.json | 117 ++++++++++++++++++ src/batches/stats/setBatchEstimateModel.m | 17 ++- .../stats/setBatchSubjectLevelGLMSpec.m | 5 +- src/subject_level/convertOnsetTsvToMat.m | 70 +++++++++-- src/subject_level/createAndReturnOnsetFile.m | 11 +- .../test_convertOnsetTsvToMat.m | 53 +++++++- 8 files changed, 273 insertions(+), 30 deletions(-) create mode 100644 demos/face_repetition/face_rep_02_stats_parametric.m create mode 100644 demos/face_repetition/models/model-faceRepetitionParametric_smdl.json diff --git a/demos/face_repetition/face_rep_02_stats.m b/demos/face_repetition/face_rep_02_stats.m index cd56e3abb..23ebb6628 100644 --- a/demos/face_repetition/face_rep_02_stats.m +++ b/demos/face_repetition/face_rep_02_stats.m @@ -19,7 +19,7 @@ opt = face_rep_get_option_results(); opt.space = 'IXI549Space'; -% bidsFFX('specifyAndEstimate', opt); +bidsFFX('specifyAndEstimate', opt); bidsFFX('contrasts', opt); bidsResults(opt); diff --git a/demos/face_repetition/face_rep_02_stats_parametric.m b/demos/face_repetition/face_rep_02_stats_parametric.m new file mode 100644 index 000000000..c78a941e2 --- /dev/null +++ b/demos/face_repetition/face_rep_02_stats_parametric.m @@ -0,0 +1,28 @@ +% This script will run the FFX and contrasts +% on the the face repetition dataset from SPM. +% +% Results might be a bit different from those in the manual as some +% default options are slightly different in this pipeline (e.g use of FAST +% instead of AR(1)...) +% +% (C) Copyright 2019 Remi Gau + +% TODO +% - create the contrast similar to the tuto + +clear; +clc; + +addpath(fullfile(pwd, '..', '..')); +cpp_spm(); + +opt = face_rep_get_option_results(); +opt.space = 'IXI549Space'; + +opt.model.file = spm_file(opt.model.file, 'basename', 'model-faceRepetitionParametric_smdl'); +opt.model.bm = BidsModel('file', opt.model.file); + +bidsFFX('specifyAndEstimate', opt); +% bidsFFX('contrasts', opt); + +% bidsResults(opt); diff --git a/demos/face_repetition/models/model-faceRepetitionParametric_smdl.json b/demos/face_repetition/models/model-faceRepetitionParametric_smdl.json new file mode 100644 index 000000000..82c42fcb5 --- /dev/null +++ b/demos/face_repetition/models/model-faceRepetitionParametric_smdl.json @@ -0,0 +1,117 @@ +{ + "Name": "parametric modulation", + "BIDSModelVersion": "1.0.0", + "Description": "model for face repetition", + "Input": { + "task": [ + "facerepetition" + ] + }, + "Nodes": [ + { + "Level": "Run", + "Name": "run_level", + "GroupBy": [ + "run", + "subject" + ], + "Transformations": { + "Description": "merge the familiarity and repetition column to create the trial type column", + "Transformer": "cpp_spm", + "Instructions": [ + { + "Name": "Concatenate", + "Input": [ + "face_type", + "repetition_type" + ], + "Output": "trial_type" + }, + { + "Name": "Copy", + "Input": "lag", + "Output": "pmod_lag" + }, + { + "Name": "Power", + "Input": "lag", + "Value": 2, + "Output": "pmod_lag_squared" + } + ] + }, + "Model": { + "X": [ + "trial_type.famous_1", + "trial_type.famous_2", + "trial_type.unfamiliar_1", + "trial_type.unfamiliar_2", + "trans_?", + "rot_?" + ], + "HRF": { + "Variables": [ + "trial_type.famous_1", + "trial_type.famous_2", + "trial_type.unfamiliar_1", + "trial_type.unfamiliar_2" + ], + "Model": "DoubleGamma" + }, + "Type": "glm", + "Options": { + "HighPassFilterCutoffHz": 0.0078 + }, + "Software": { + "SPM": { + "SerialCorrelation": "AR(1)", + "HRFderivatives": "TemporalDispersion" + } + } + }, + "DummyContrasts": { + "Test": "t", + "Contrasts": [ + "trial_type.famous_1", + "trial_type.famous_2", + "trial_type.unfamiliar_1", + "trial_type.unfamiliar_2" + ] + }, + "Contrasts": [ + { + "Name": "faces_gt_baseline", + "ConditionList": [ + "trial_type.famous_1", + "trial_type.famous_2", + "trial_type.unfamiliar_1", + "trial_type.unfamiliar_2" + ], + "Weights": [ + 1, + 1, + 1, + 1 + ], + "Test": "t" + }, + { + "Name": "faces_lt_baseline", + "ConditionList": [ + "trial_type.famous_1", + "trial_type.famous_2", + "trial_type.unfamiliar_1", + "trial_type.unfamiliar_2" + ], + "Weights": [ + -1, + -1, + -1, + -1 + ], + "Test": "t" + } + ] + } + ] +} diff --git a/src/batches/stats/setBatchEstimateModel.m b/src/batches/stats/setBatchEstimateModel.m index 7dc6d4651..a873faa27 100644 --- a/src/batches/stats/setBatchEstimateModel.m +++ b/src/batches/stats/setBatchEstimateModel.m @@ -1,15 +1,24 @@ function matlabbatch = setBatchEstimateModel(matlabbatch, opt, nodeName, contrastsList) % - % Short description of what the function does goes here. + % Set up the estimate model batch for run/subject or group level GLM % % USAGE:: % - % matlabbatch = setBatchEstimateModel(matlabbatch, grpLvlCon) + % matlabbatch = setBatchEstimateModel(matlabbatch, opt) + % matlabbatch = setBatchEstimateModel(matlabbatch, opt, nodeName, contrastsList) % % :param matlabbatch: % :type matlabbatch: structure - % :param grpLvlCon: - % :type grpLvlCon: + % + % :param opt: + % :type opt: structure + % + % :param nodeName: + % :type nodeName: char + % + % :param contrastsList: + % :type contrastsList: cell string + % % % :returns: - :matlabbatch: (structure) % diff --git a/src/batches/stats/setBatchSubjectLevelGLMSpec.m b/src/batches/stats/setBatchSubjectLevelGLMSpec.m index fe20cd95b..da03496b9 100644 --- a/src/batches/stats/setBatchSubjectLevelGLMSpec.m +++ b/src/batches/stats/setBatchSubjectLevelGLMSpec.m @@ -16,7 +16,7 @@ % :type opt: structure % % :param subLabel: - % :type subLabel: string + % :type subLabel: char % % :returns: - :matlabbatch: (structure) % @@ -228,8 +228,7 @@ % get events file from raw data set and convert it to a onsets.mat file % store in the subject level GLM directory - filter = struct( ... - 'sub', subLabel, ... + filter = struct('sub', subLabel, ... 'task', task, ... 'ses', session, ... 'run', run, ... diff --git a/src/subject_level/convertOnsetTsvToMat.m b/src/subject_level/convertOnsetTsvToMat.m index 4682efad8..0c9ca5cdf 100644 --- a/src/subject_level/convertOnsetTsvToMat.m +++ b/src/subject_level/convertOnsetTsvToMat.m @@ -2,11 +2,6 @@ % % Converts an events.tsv file to an onset file suitable for SPM subject level analysis. % - % The function extracts from the events.tsv file the trials (with type, onsets, and durations) - % of the conditions of interest as requested in the model.json. - % It then stores them in the GLM folder in a .mat file - % that can be fed directly in an SPM GLM batch. - % % USAGE:: % % fullpathOnsetFilename = convertOnsetTsvToMat(opt, tsvFile) @@ -17,6 +12,26 @@ % :param tsvFile: % :type tsvFile: string % + % Use a BIDS stats model specified in a JSON file to: + % + % - loads events.tsv and apply the Node.Transformations to its content + % + % - extract the trials (onsets, durations) + % of the conditions that should be convolved + % as requested from Node.HRF.Variables field + % + % It then stores them in in a .mat file + % that can be fed directly in an SPM GLM batch as 'Multiple conditions' + % + % Parametric modulation can be specified via columns in the TSV + % file starting with ``pmod_``. + % These columns can be created via the use of Node.Transformations. + % Only polynomial 1 are supported. + % More complex modulation should be precomputed via the Transformations. + % + % if ``opt.glm.useDummyRegressor`` is set to ``true``, any missing condition + % will be replaced by a DummyRegressor. + % % :returns: :fullpathOnsetFilename: (string) name of the output ``.mat`` file. % % See also: createAndReturnOnsetFile, bids.transformers @@ -52,13 +67,17 @@ names = {}; onsets = {}; durations = {}; + pmod = struct('name', {''}, 'param', {}, 'poly', {}); if ~isfield(opt.model, 'bm') opt.model.bm = BidsModel('file', opt.model.file); end + % TODO get / apply transformers from a specific node transformers = opt.model.bm.getBidsTransformers(); tsvContent = bids.transformers(transformers, tsvContent); + conditionIdx = 1; + for iCond = 1:numel(variablesToConvolve) trialTypeNotFound = false; @@ -74,22 +93,47 @@ trialTypes = tsvContent.(tokens{1}); conditionName = strjoin(tokens(2:end), '.'); - idx = find(strcmp(conditionName, trialTypes)); + rows = find(strcmp(conditionName, trialTypes)); printToScreen(sprintf(' Condition %s: %i trials found.\n', ... conditionName, ... - numel(idx)), ... + numel(rows)), ... opt); - if ~isempty(idx) + if ~isempty(rows) if ~ismember(variablesToConvolve{iCond}, designMatrix) continue end - names{1, end + 1} = conditionName; - onsets{1, end + 1} = tsvContent.onset(idx)'; %#ok<*AGROW,*NASGU> - durations{1, end + 1} = tsvContent.duration(idx)'; + names{1, conditionIdx} = conditionName; + onsets{1, conditionIdx} = tsvContent.onset(rows)'; %#ok<*AGROW,*NASGU> + durations{1, conditionIdx} = tsvContent.duration(rows)'; + + % parametric modulation (pmod) + % + % skipped if parametric modulation is 1 for all onsets + % + % coerces NaNs into 1 + + fields = fieldnames(tsvContent); + pmodIdx = ~cellfun('isempty', regexp(fields, '^pmod_.*', 'match')); + pmodIdx = find(pmodIdx); + + for iMod = 1:numel(pmodIdx) + + thisMod = fields{pmodIdx(iMod)}; + + amplitude = tsvContent.(thisMod)(rows); + amplitude(isnan(amplitude)) = 1; + + if ~all(amplitude == 1) + pmod(1, conditionIdx).name{iMod} = strrep(thisMod, 'pmod_', ''); + pmod(end).param{iMod} = amplitude; + pmod(end).poly{iMod} = 1; + end + + end else @@ -124,6 +168,8 @@ end + conditionIdx = conditionIdx + 1; + end %% save the onsets as a matfile @@ -136,7 +182,7 @@ fullpathOnsetFilename = fullfile(pth, bf.filename); save(fullpathOnsetFilename, ... - 'names', 'onsets', 'durations', ... + 'names', 'onsets', 'durations', 'pmod', ... '-v7'); end diff --git a/src/subject_level/createAndReturnOnsetFile.m b/src/subject_level/createAndReturnOnsetFile.m index 72d607a6d..9b1faca2a 100644 --- a/src/subject_level/createAndReturnOnsetFile.m +++ b/src/subject_level/createAndReturnOnsetFile.m @@ -8,19 +8,22 @@ % % USAGE:: % - % onsetFilename = createAndReturnOnsetFile(opt, subLabel, tsvFile, funcFWHM) + % onsetFilename = createAndReturnOnsetFile(opt, subLabel, tsvFile) % % :param opt: % :type opt: structure + % % :param subLabel: - % :type subLabel: string + % :type subLabel: char + % % :param tsvFile: fullpath name of the tsv file. - % :type tsvFile: string + % :type tsvFile: char % - % :returns: :onsetFilename: (string) fullpath name of the file created. + % :returns: :onsetFilename: (path) fullpath name of the file created. % % See also: convertOnsetTsvToMat % + % % (C) Copyright 2019 CPP_SPM developers if iscell(tsvFile) diff --git a/tests/tests_subject_level/test_convertOnsetTsvToMat.m b/tests/tests_subject_level/test_convertOnsetTsvToMat.m index d0e79b13d..cd420ee85 100644 --- a/tests/tests_subject_level/test_convertOnsetTsvToMat.m +++ b/tests/tests_subject_level/test_convertOnsetTsvToMat.m @@ -11,6 +11,53 @@ end +function test_convertOnsetTsvToMat_parametric_modulation() + + % GIVEN + tsvFile = fullfile(getDummyDataDir(), ... + 'tsv_files', ... + 'sub-01_task-vismotion_events.tsv'); + opt = setOptions('vismotion'); + + opt.model.bm = BidsModel('file', opt.model.file); + + trans = struct('Description', 'add dummy param mod', ... + 'Transformer', 'cpp_spm', ... + 'Instructions', {{ ... + struct('Name', 'Constant', ... + 'Value', 3, ... + 'Output', 'pmod_amp'), ... + struct('Name', 'Power', ... + 'Input', 'pmod_amp', ... + 'Value', 2, ... + 'Output', 'pmod_amp_squared') + }}); + + opt.model.bm.Nodes{1}.Transformations = trans; + + % WHEN + fullpathOnsetFilename = convertOnsetTsvToMat(opt, tsvFile); + + load(fullpathOnsetFilename, 'names', 'onsets', 'durations', 'pmod'); + + assertEqual(names, {'VisMot', 'VisStat'}); + assertEqual(onsets, {2, 4}); + assertEqual(durations, {2, 2}); + + assertEqual(pmod(1).name{1}, 'amp'); + assertEqual(pmod(1).param(1), {3}); + assertEqual(pmod(1).name{2}, 'amp_squared'); + assertEqual(pmod(1).param(2), {9}); + + assertEqual(pmod(2).name{1}, 'amp'); + assertEqual(pmod(2).param(1), {3}); + assertEqual(pmod(2).name{2}, 'amp_squared'); + assertEqual(pmod(2).param(2), {9}); + + cleanUp(fullpathOnsetFilename); + +end + function test_convertOnsetTsvToMat_warning_missing_variable_to_convolve if isOctave @@ -151,12 +198,6 @@ function test_convertOnsetTsvToMat_dummy_regressor() fullpathOnsetFilename = convertOnsetTsvToMat(opt, tsvFile); % THEN - assertEqual(fullfile(getDummyDataDir(), ... - 'tsv_files', ... - 'sub-01_task-vismotion_onsets.mat'), ... - fullpathOnsetFilename); - assertEqual(exist(fullpathOnsetFilename, 'file'), 2); - load(fullpathOnsetFilename); assertEqual(names, {'dummyRegressor'}); From e69a37d918c0833c68dd8e8b5c39c3a369c5f05e Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Wed, 25 May 2022 12:31:51 +0200 Subject: [PATCH 2/4] fix system test --- .github/workflows/system_tests_matlab.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/system_tests_matlab.m b/.github/workflows/system_tests_matlab.m index 98d07aa4c..022e31772 100644 --- a/.github/workflows/system_tests_matlab.m +++ b/.github/workflows/system_tests_matlab.m @@ -5,7 +5,7 @@ addpath(fullfile(root_dir, 'spm12')); -cd demos/MoAE; +cd(fullfile(root_dir, 'demos', 'MoAE')); download_moae_ds(true); cd(fullfile(root_dir, 'manualTests')); From 4b01dbe4e3b589c3d1ff68d2dd472f862cf0d6a5 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Wed, 25 May 2022 12:34:36 +0200 Subject: [PATCH 3/4] refactor --- src/subject_level/convertOnsetTsvToMat.m | 53 +++++++++++++----------- 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/src/subject_level/convertOnsetTsvToMat.m b/src/subject_level/convertOnsetTsvToMat.m index 0c9ca5cdf..ba83a12eb 100644 --- a/src/subject_level/convertOnsetTsvToMat.m +++ b/src/subject_level/convertOnsetTsvToMat.m @@ -109,31 +109,7 @@ names{1, conditionIdx} = conditionName; onsets{1, conditionIdx} = tsvContent.onset(rows)'; %#ok<*AGROW,*NASGU> durations{1, conditionIdx} = tsvContent.duration(rows)'; - - % parametric modulation (pmod) - % - % skipped if parametric modulation is 1 for all onsets - % - % coerces NaNs into 1 - - fields = fieldnames(tsvContent); - pmodIdx = ~cellfun('isempty', regexp(fields, '^pmod_.*', 'match')); - pmodIdx = find(pmodIdx); - - for iMod = 1:numel(pmodIdx) - - thisMod = fields{pmodIdx(iMod)}; - - amplitude = tsvContent.(thisMod)(rows); - amplitude(isnan(amplitude)) = 1; - - if ~all(amplitude == 1) - pmod(1, conditionIdx).name{iMod} = strrep(thisMod, 'pmod_', ''); - pmod(end).param{iMod} = amplitude; - pmod(end).poly{iMod} = 1; - end - - end + pmod = parametricModulation(pmod, tsvContent, rows, conditionIdx); else @@ -187,6 +163,33 @@ end +function pmod = parametricModulation(pmod, tsvContent, rows, conditionIdx) + % parametric modulation (pmod) + % + % skipped if parametric modulation is 1 for all onsets + % + % coerces NaNs into 1 + + fields = fieldnames(tsvContent); + pmodIdx = ~cellfun('isempty', regexp(fields, '^pmod_.*', 'match')); + pmodIdx = find(pmodIdx); + + for iMod = 1:numel(pmodIdx) + + thisMod = fields{pmodIdx(iMod)}; + + amplitude = tsvContent.(thisMod)(rows); + amplitude(isnan(amplitude)) = 1; + + if ~all(amplitude == 1) + pmod(1, conditionIdx).name{iMod} = strrep(thisMod, 'pmod_', ''); + pmod(end).param{iMod} = amplitude; + pmod(end).poly{iMod} = 1; + end + + end +end + function [names, onsets, durations] = addDummyRegressor(names, onsets, durations) names{1, end + 1} = 'dummyRegressor'; From e35e04a87e7e1a10a656c6864bff7286c46e32d3 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Wed, 25 May 2022 14:56:36 +0200 Subject: [PATCH 4/4] update default content onsets.mat --- tests/dummyData/mat_files/onsets.mat | Bin 344 -> 415 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/dummyData/mat_files/onsets.mat b/tests/dummyData/mat_files/onsets.mat index ada01389664b4c7175e58d1608a784c2264d7747..53b59c899105553e6baeb85bd09bce15c4989b62 100644 GIT binary patch delta 115 zcmV-(0F3|G0-pnrG#FQ9WFSppc_1=1ATcyLH8nakH6SuDGBS}-BavVQk#rD~7y&?$ zfS5x+0001Zoa19)U@)Zn9ER*pHm4_4*<-Y2gxK