Merge pull request #465 from cpp-lln-lab/remi-roi_based_glm

[ENH] update roi based glm

.gitignore

@@ -8,10 +8,9 @@
 *.html
 *.ps
 
-
-*.mat
-
-.zenodo.json
+# created by tests
+CHANGES
+README
 
 # Project specific
 options_task-*date*.json

demos/face_repetition/face_rep_02_stats.m

@@ -21,6 +21,7 @@
 end
 
 opt = face_rep_get_option_results();
+opt.space = 'IXI549Space';
 
 bidsFFX('specifyAndEstimate', opt);
 bidsFFX('contrasts', opt);

demos/face_repetition/face_rep_03_roi_analysis.m

@@ -8,6 +8,7 @@
 % (C) Copyright 2019 Remi Gau
 
 clear;
+close all;
 clc;
 
 try
@@ -16,16 +17,14 @@
 end
 
 opt = face_rep_get_option_results();
-opt.taskName = {'facerepetition'};
 
 opt.roi.atlas = 'wang';
 opt.roi.name = {'V1v', 'V1d'};
-opt.roi.space = {'IXI549Space', 'individual'};
+opt.roi.space = {'individual'};
 
-opt.dir.stats = fullfile(opt.dir.raw, '..', 'derivatives', 'cpp_spm-stats');
-
-bidsCreateROI(opt);
+% bidsCreateROI(opt);
 
 opt.glm.roibased.do = true;
+opt.space = 'individual';
 
 bidsRoiBasedGLM(opt);

demos/face_repetition/face_rep_get_option_results.m

@@ -4,6 +4,8 @@
 
   opt = [];
 
+  opt.pipeline.type = 'stats';
+
   opt.verbosity = 1;
 
   opt.dir.raw = fullfile(fileparts(mfilename('fullpath')), 'outputs', 'raw');
@@ -41,7 +43,6 @@
   opt.result.Nodes(1).Output.montage.orientation = 'axial';
 
   %% DO NOT TOUCH
-  opt = checkOptions(opt);
   saveOptions(opt);
 
 end

demos/face_repetition/models/model-faceRepetition_smdl.json

@@ -3,8 +3,7 @@
    "BIDSModelVersion": "1.0.0",
    "Description": "model for face repetition to check resampling effects",
    "Input": {
-      "task": "facerepetition",
-      "space": "IXI549Space"
+      "task": "facerepetition"
    },
    "Nodes": [
       {

docs/source/_static/bidsFFX/Makefile

@@ -6,9 +6,14 @@
 
 # Final output
 # beta_XXXX.nii: mask.nii
-# mask.nii: cpp_spm-stats/sub-01/stats/task-auditory_space-IXI549Space_FWHM-6
+# mask.nii: cpp_spm-stats/sub-01/stats/task-audio_space-IXI549Space_FWHM-6
 
-cpp_spm-stats/sub-01/stats/task-auditory_space-IXI549Space_FWHM-6: SPM.mat sub-01_task-auditory_desc-confounds_regressors.mat sub-01_task-auditory_onsets.mat
+cpp_spm-stats/sub-01/stats/task-audio_space-IXI549Space_FWHM-6: SPM.mat \
+																sub-01_task-audio_regressors.mat \
+																sub-01_task-audio_onsets.mat \
+																sub-01_task-audio_regressors.tsv \
+																sub-01_task-audio_onsets.tsv \
+																sub-01_task-audio_desc-beforeEstimation_designmatrix.png
 
 # Workflows and batches
 setBatchSubjectLevelGLMSpec.m: bidsFFX.m
@@ -23,16 +28,29 @@ convertOnsetTsvToMat.m: createAndReturnOnsetFile.m
 
 # Files
 
-model-auditory_smdl.json: models/
+model-audio_smdl.json: models/
 
-sub-01_task-auditory_desc-confounds_regressors.tsv: preproc/sub-01/func/
+sub-01_task-audio_desc-confounds_regressors.tsv: preproc/sub-01/func/
 
 sub-01_task-audio_space-IXI549Space_desc-smth6_bold.nii: preproc/sub-01/func/
 
 sub-01_task-audio_events.tsv: raw/sub-01/func/
 
-sub-01_task-auditory_desc-confounds_regressors.mat: sub-01_task-auditory_desc-confounds_regressors.tsv model-auditory_smdl.json createAndReturnCounfoundMatFile.m
+sub-01_task-audio_regressors.mat: sub-01_task-audio_desc-confounds_regressors.tsv \
+									model-audio_smdl.json \
+									createAndReturnCounfoundMatFile.m
 
-sub-01_task-auditory_onsets.mat: sub-01_task-audio_events.tsv model-auditory_smdl.json convertOnsetTsvToMat.m
+sub-01_task-audio_regressors.tsv: sub-01_task-audio_regressors.mat
 
-SPM.mat: sub-01_task-audio_space-IXI549Space_desc-smth6_bold.nii sub-01_task-auditory_onsets.mat sub-01_task-auditory_desc-confounds_regressors.mat model-auditory_smdl.json setBatchSubjectLevelGLMSpec.m
+sub-01_task-audio_onsets.mat: sub-01_task-audio_events.tsv \
+									model-audio_smdl.json \
+									convertOnsetTsvToMat.m
+
+sub-01_task-audio_onsets.tsv: sub-01_task-audio_onsets.mat
+
+SPM.mat: sub-01_task-audio_space-IXI549Space_desc-smth6_bold.nii \
+		sub-01_task-audio_onsets.mat \
+		sub-01_task-audio_regressors.mat \
+		model-audio_smdl.json setBatchSubjectLevelGLMSpec.m
+
+sub-01_task-audio_desc-beforeEstimation_designmatrix.png: SPM.mat

docs/source/conf.py

@@ -22,7 +22,7 @@
 author = 'the CPP SPM pipeline dev team'
 
 # The full version, including alpha/beta/rc tags
-release = 'v0.1.0'
+release = 'v1.1.3dev'
 
 
 # -- General configuration ---------------------------------------------------

docs/source/links_and_references.rst

@@ -43,8 +43,11 @@ Design efficiency: https://github.com/Remi-Gau/advanced_fMRI_course/blob/master/
 Content of SPM.mat
 ==================
 
-This is here because SPM has the sad tradition of using variable names that have often attempted to replicate the notation in the papers,
-rather then the ``TypicalLongVariableNames`` that many programmers and new comers would prefer to see.
+This is here because SPM has the sad (and bad) Matlabic tradition of using variable names
+that have often attempted to replicate the notation in the papers to make engineers
+and the generally math enclined happy,
+rather then the ``TypicalLongVariableNames`` that many programmers and new comers
+would prefer to see to help with code readibility.
 
 Adapted from: http://andysbrainblog.blogspot.com/2013/10/whats-in-spmmat-file.html
 
@@ -54,7 +57,7 @@ details on experiment
 
 - ``SPM.xY.RT``                - TR length (RT ="repeat time")
 - ``SPM.xY.P``                 - matrix of file names
-- ``SPM.xY.VY``                - (# of runs x 1) struct array of mapped image volumes (.nii file info)
+- ``SPM.xY.VY``                - (number of runs x 1) struct array of mapped image volumes (.nii file info)
 
 - ``SPM.modality``             - the data you're using (PET, FMRI, EEG)
 
@@ -91,15 +94,15 @@ user-specified covariates/regressors
 
 e.g. motion
 
-- ``SPM.Sess([sesssion]).C.C``    - (n x c) double regressor (c is #covariates, n is #sessions)
+- ``SPM.Sess([sesssion]).C.C``    - (n x c) double regressor (``c`` is number of covariates, ``n`` is number of sessions)
 - ``SPM.Sess([sesssion]).C.name`` - names of covariates
 
 conditions & modulators specified
 +++++++++++++++++++++++++++++++++
 
 i.e. input structure array
 
-- ``SPM.Sess([sesssion]).U(condition).dt``:  - time bin length {seconds}
+- ``SPM.Sess([sesssion]).U(condition).dt``:  - time bin length (seconds)
 - ``SPM.Sess([sesssion]).U(condition).name`` - names of conditions
 - ``SPM.Sess([sesssion]).U(condition).ons``  - onset for condition's trials
 - ``SPM.Sess([sesssion]).U(condition).dur``  - duration for condition's trials
@@ -113,7 +116,7 @@ parameters/modulators specified
 - ``SPM.Sess([sesssion]).U(condition).P.name`` - names of modulators/parameters
 - ``SPM.Sess([sesssion]).U(condition).P.h``    - polynomial order of modulating parameter (order of polynomial expansion where 0 is none)
 - ``SPM.Sess([sesssion]).U(condition).P.P``    - vector of modulating values
-- ``SPM.Sess([sesssion]).U(condition).P.P.i``  - sub-indices of U(i).u for plotting
+- ``SPM.Sess([sesssion]).U(condition).P.P.i``  - sub-indices of ``U(i).u`` for plotting
 
 scan indices for sessions
 +++++++++++++++++++++++++
@@ -156,7 +159,7 @@ design matrix information
 - ``SPM.xX.iB``       - vector of B partition (block effects) indices
 - ``SPM.xX.iG``       - vector of G partition (nuisance variables) indices
 
-- ``SPM.xX.K``        - cell. low frequency confound: high-pass cutoff (secs)
+- ``SPM.xX.K``        - cell. low frequency confound: high-pass cutoff (seconds)
 - ``SPM.xX.K.HParam`` - low frequency cutoff value
 - ``SPM.xX.K.X0``     - cosines (high-pass filter)
 
@@ -327,7 +330,7 @@ added after running contrasts
 
     - If by columns for reduced design matrix then iX0 contains the column indices.
     - Otherwise, it's a string containing the ``spm_FcUtil`` 'Set' action:
-      Usually one of {'c','c+','X0'} defines the indices of the columns that will not be tested. Can be empty.
+      Usually one of ``{'c','c+','X0'}`` defines the indices of the columns that will not be tested. Can be empty.
 
   - ``SPM.xCon.X1o`` - Remaining design space data (X1o is orthogonal to X0)
 

docs/source/statistics.rst

@@ -14,11 +14,44 @@ Subject level
 
 .. autofunction:: bidsFFX
 
+After the specification step an output folder is created.
+To get the fullpath of that folder you can use::
+
+   getFFXdir(subLabel, opt)
+
+A typical folder will contain::
+
+   cpp_spm-stats/sub-01/stats/task-audio_space-IXI549Space_FWHM-6
+   ├── SPM.mat
+   ├── sub-01_task-audio_space-IXI549Space_desc-beforeEstimation_designmatrix.png
+   ├── sub-01_task-audio_run-01_desc-confounds_regressors.mat
+   ├── sub-01_task-audio_run-01_desc-confounds_regressors.tsv
+   ├── sub-01_task-audio_run-01_onsets.mat
+   └── sub-01_task-audio_run-01_onsets.tsv
+
+Each run should have a pair of tsv/mat files:
+
+- One that summarises the onsets used for that design.
+- One that summarises the regressors confounds used for that design.
+
+In most cases those are going to be a subset of the content:
+
+- of the ``_events.tsv`` from the raw BIDS dataset
+- of the ``_regressors.tsv`` from the deriratives BIDS dataset containing
+  the preprocessed data.
+
+What part of the _events.tsv and _regressors.tsv
+gets into the final GLM specification
+depends on the BIDS statistical model used.
+
+The mat files can directly be ingested by SPM:
+the TSV files are there for both logging and interoperability.
+
 .. _fig_FFX-specification:
 .. figure::  _static/bidsFFX/out.png
    :align:   center
 
-   Subject level GLM specification workflow
+   Subject level GLM specification workflow for model specification
 
 .. autofunction:: bidsConcatBetaTmaps
 

src/QA/plotEvents.m

@@ -27,7 +27,7 @@ function plotEvents(eventsFile, modelFile)
     return
   end
 
-  figure('name', figName, 'position', [50 50 1000 200 * numel(conditions)]);
+  spm_figure('Create', 'Tag', figName, 'on');
 
   for iCdt = 1:numel(conditions)
 

src/QA/plotRoiTimeCourse.m

@@ -0,0 +1,82 @@
+function figureFile = plotRoiTimeCourse(varargin)
+  %
+  % Plots the peristimulus histogram from a ROI based GLM
+  %
+  % USAGE::
+  %
+  %   plotRoiTimeCourse(tsvFile)
+  %
+  % :param tsvFile: obligatory argument.
+  % :type tsvFile: path
+  %
+  % See also: bidsRoiBasedGLM(varargin)
+  %
+  % (C) Copyright 2022 CPP_SPM developers
+
+  p = inputParser;
+
+  isFile = @(x) exist(x, 'file') == 2;
+
+  addRequired(p, 'tsvFile', isFile);
+  addOptional(p, 'verbose', true, @islogical);
+
+  parse(p, varargin{:});
+
+  visible = 'off';
+  if p.Results.verbose
+    visible = 'on';
+  end
+
+  tsvFile = p.Results.tsvFile;
+  timeCourse = bids.util.tsvread(tsvFile);
+  conditionNames = fieldnames(timeCourse);
+
+  for i = 1:numel(conditionNames)
+    tmp(:, i) = timeCourse.(conditionNames{i});
+  end
+  timeCourse = tmp;
+
+  jsonFile = bids.internal.file_utils(tsvFile, 'ext', '.json');
+  content = bids.util.jsondecode(jsonFile);
+
+  secs = [0:size(timeCourse, 1) - 1] * content.SamplingFrequency;
+
+  bf = bids.File(tsvFile);
+
+  if isGithubCi()
+    return
+  end
+
+  figName = ['ROI: ' bf.entities.label];
+  if isfield(bf.entities, 'hemi')
+    figName = [figName ' - ' bf.entities.hemi];
+  end
+  spm_figure('Create', 'Tag', figName, visible);
+
+  hold on;
+
+  plot(secs, timeCourse, 'linewidth', 2);
+  plot([0 secs(end)], [0 0], '--k');
+
+  title(['Time courses for ' figName], ...
+        'Interpreter', 'none');
+
+  xlabel('Seconds');
+  ylabel('Percent signal change');
+
+  legend_content = regexprep(conditionNames, '_', ' ');
+
+  axis tight;
+
+  for i = 1:numel(conditionNames)
+    legend_content{i} = sprintf('%s ; max(PSC)= %0.2f', ...
+                                legend_content{i}, ...
+                                content.(conditionNames{i}).percentSignalChange.max);
+  end
+
+  legend(legend_content);
+
+  figureFile = bids.internal.file_utils(tsvFile, 'ext', '.png');
+  print(gcf, figureFile, '-dpng');
+
+end

src/batches/preproc/setBatchCoregistration.m

@@ -8,6 +8,8 @@
   %
   % :param matlabbatch: list of SPM batches
   % :type matlabbatch: structure
+  % :param opt: Other images to apply the coregistration to
+  % :type opt: cell string
   % :param ref: Reference image
   % :type ref: string
   % :param src: Source image
@@ -26,6 +28,7 @@
   p = inputParser;
 
   addRequired(p, 'matlabbatch', @iscell);
+  addRequired(p, 'opt', @isstruct);
   addRequired(p, 'ref', isFile);
   addRequired(p, 'src', isFile);
   addOptional(p, 'other', default_other, @iscell);
@@ -37,7 +40,7 @@
   src = p.Results.src;
   other = p.Results.other;
 
-  printBatchName('coregistration', opt);
+  printBatchName('coregistration', p.Results.opt);
 
   matlabbatch{end + 1}.spm.spatial.coreg.estimate.ref = { ref };
   matlabbatch{end}.spm.spatial.coreg.estimate.source = { src };

src/batches/preproc/setBatchCoregistrationFmap.m

@@ -60,7 +60,10 @@
         filter.suffix = 'magnitude1';
         srcImage = bids.query(BIDS, 'data', filter);
 
-        matlabbatch = setBatchCoregistration(matlabbatch, opt, refImage, srcImage{1}, otherImages);
+        matlabbatch = setBatchCoregistration(matlabbatch, opt, ...
+                                             refImage{1}, ...
+                                             srcImage{1}, ...
+                                             otherImages);
 
       end
 

src/batches/preproc/setBatchCoregistrationFuncToAnat.m

@@ -1,7 +1,7 @@
 function matlabbatch = setBatchCoregistrationFuncToAnat(matlabbatch, BIDS, opt, subLabel)
   %
-  % Set the batch for corregistering the functional images to the
-  % anatomical image
+  % Set the batch for coregistering the functional images
+  % to the anatomical image.
   %
   % USAGE::
   %