[ENH] Enhancements subject / group level GLM and results

miss_hit.cfg

@@ -31,5 +31,5 @@ tab_width: 2
 # metrics limit for the code quality (https://florianschanda.github.io/miss_hit/metrics.html)
 metric "cnest": limit 4
 metric "file_length": limit 400
-metric "cyc": limit 12
+metric "cyc": limit 13
 metric "parameters": limit 7

src/batches/setBatchGroupLevelResults.m

@@ -0,0 +1,35 @@
+function matlabbatch = setBatchGroupLevelResults(varargin)
+  %
+  % Short description of what the function does goes here.
+  %
+  % USAGE::
+  %
+  %   matlabbatch = setBatchGroupLevelResults(matlabbatch, opt, result)
+  %
+  % :param matlabbatch:
+  % :type matlabbatch: structure
+  % :param opt:
+  % :type opt: structure
+  % :param result:
+  % :type result: structure
+  %
+  % :returns: - :matlabbatch: (structure)
+  %
+  % (C) Copyright 2019 CPP_SPM developers
+
+  [matlabbatch, opt, result] = deal(varargin{:});
+
+  result.dir = fullfile(getRFXdir(opt), result.Contrasts.Name);
+
+  load(fullfile(result.dir, 'SPM.mat'));
+  result.nbSubj = SPM.nscan;
+
+  result.contrastNb = 1;
+
+  result.label = 'group';
+
+  result.outputNameStructure = defaultOuputNameStruct(opt, result);
+
+  matlabbatch = setBatchResults(matlabbatch, result);
+
+end

src/batches/setBatchSubjectLevelGLMSpec.m

@@ -1,6 +1,6 @@
 function matlabbatch = setBatchSubjectLevelGLMSpec(varargin)
   %
-  % Short description of what the function does goes here.
+  % Sets up the subject level GLM
   %
   % USAGE::
   %
@@ -15,7 +15,7 @@
   % :param subLabel:
   % :type subLabel: string
   %
-  % :returns: - :argout1: (structure) (matlabbatch)
+  % :returns: - :matlabbatch: (structure)
   %
   % (C) Copyright 2019 CPP_SPM developers
 
@@ -55,13 +55,14 @@
   % Create ffxDir if it doesnt exist
   % If it exists, issue a warning that it has been overwritten
   ffxDir = getFFXdir(subLabel, opt);
-
   overwriteDir(ffxDir, opt);
 
   fmri_spec.dir = {ffxDir};
 
   fmri_spec.fact = struct('name', {}, 'levels', {});
 
+  fmri_spec.mthresh = getInclusiveMaskThreshold(opt.model.file);
+
   fmri_spec.bases.hrf.derivs = getHRFderivatives(opt.model.file);
 
   fmri_spec.cvi = getSerialCorrelationCorrection(opt.model.file);
@@ -175,7 +176,15 @@
     end
     fmriSpec.sess(sesCounter).nscan = numel(hdr);
   else
-    fmriSpec.sess(sesCounter).scans = {fullpathBoldFilename};
+    if opt.glm.maxNbVols == Inf
+      scans = {fullpathBoldFilename};
+    else
+      scans = cellstr(spm_select('ExtFPList', ...
+                                 spm_fileparts(fullpathBoldFilename), ...
+                                 spm_file(fullpathBoldFilename, 'filename'), ...
+                                 1:opt.glm.maxNbVols));
+    end
+    fmriSpec.sess(sesCounter).scans = scans;
   end
 end
 
@@ -216,7 +225,7 @@
   mask = getModelMask(opt.model.file);
 
   if isempty(mask) && ...
-          (strfind(opt.space{1}, 'MNI') || strcmp(opt.space, 'IXI549Space'))
+          (~isempty(strfind(opt.space{1}, 'MNI')) || strcmp(opt.space, 'IXI549Space'))
     mask = spm_select('FPList', fullfile(spm('dir'), 'tpm'), 'mask_ICV.nii');
   end
 

src/batches/setBatchSubjectLevelResults.m

@@ -12,19 +12,14 @@
   % :type opt: structure
   % :param subLabel:
   % :type subLabel: string
-  % :param iNode:
-  % :type iNode: positive integer
-  % :param iCon:
-  % :type iCon: positive integer
   %
   % :returns: - :matlabbatch: (structure)
   %
   % (C) Copyright 2019 CPP_SPM developers
 
-  [matlabbatch, opt, subLabel, iNode, iCon] = deal(varargin{:});
+  [matlabbatch, opt, subLabel, result] = deal(varargin{:});
 
-  result.Contrasts = opt.result.Nodes(iNode).Contrasts(iCon);
-  result.dir = getFFXdir(subLabel, opt);
+  result.nbSubj = 1;
 
   result.contrastNb = getContrastNb(result, opt);
   if isempty(result.contrastNb)
@@ -33,29 +28,9 @@
     return
   end
 
-  if isfield(opt.result.Nodes(iNode), 'Output')
-    result.Output =  opt.result.Nodes(iNode).Output;
-  end
-
-  result.space = opt.space;
-
   result.label = subLabel;
 
-  result.nbSubj = 1;
-
-  result.outputNameStructure = struct( ...
-                                      'suffix', 'spmT', ...
-                                      'ext', '.nii', ...
-                                      'use_schema', 'false', ...
-                                      'entities', struct('sub', '', ...
-                                                         'task', strjoin(opt.taskName, ''), ...
-                                                         'space', result.space, ...
-                                                         'desc', '', ...
-                                                         'label', sprintf('%04.0f', ...
-                                                                          result.contrastNb), ...
-                                                         'p', '', ...
-                                                         'k', '', ...
-                                                         'MC', ''));
+  result.outputNameStructure = defaultOuputNameStruct(opt, result);
 
   matlabbatch = setBatchResults(matlabbatch, result);
 

src/bids_model/getInclusiveMaskThreshold.m

@@ -0,0 +1,39 @@
+function threshold = getInclusiveMaskThreshold(modelFile, nodeType)
+  %
+  % returns the threshold for inclusive masking of subject level GLM
+  % node of a BIDS statistical model
+  %
+  % (C) Copyright 2021 CPP_SPM developers
+
+  if isempty(modelFile)
+    return
+  end
+  if nargin < 2 || isempty(nodeType)
+    nodeType = 'run';
+  elseif ismember(nodeType, 'dataset')
+    warning('No inclusive mask threshold for dataset level GLM.');
+    return
+  end
+
+  model = bids.util.jsondecode(modelFile);
+
+  node = returnModelNode(model, nodeType);
+
+  try
+    threshold = node.Model.Software.SPM.InclusiveMaskingThreshold;
+  catch
+    threshold = '';
+  end
+
+  if isempty(threshold)
+
+    spm_defaults = spm_get_defaults();
+    threshold =  spm_defaults.mask.thresh;
+
+  elseif strcmpi(threshold, '-Inf')
+
+    threshold =  -Inf;
+
+  end
+
+end

src/defaults/checkOptions.m

@@ -83,7 +83,11 @@
   %       if you know the order in which slices were acquired, you should be able to recompute
   %       slice timing and add it to the json files in your BIDS data set.
   %
-  %     - ``opt.glm.roibased.do``
+  %     - ``opt.glm.roibased.do = false`` must be set to ``true`` to use the
+  %       ``bidsRoiBasedGLM`` workflow
+  %     - ``opt.glm.maxNbVols = Inf`` sets the maximum number of volumes to
+  %       include in a run in a subject level GLM. This can be useful if some
+  %       time series have more volumes than necessary.
   %
   %     - ``opt.QA.func.carpetPlot = true`` to plot carpet plot when running ``functionaQA``
   %     - ``opt.QA.func`` contains a lot of options used by ``spmup_first_level_qa``
@@ -215,6 +219,7 @@
   fieldsToSet.QA.func.Basics = 'on';
 
   fieldsToSet.glm.roibased.do = false;
+  fieldsToSet.glm.maxNbVols = Inf;
 
   % specify the results to compute
   fieldsToSet.result.Nodes = returnDefaultResultsStructure();

src/results/defaultOuputNameStruct.m

@@ -0,0 +1,33 @@
+function outputNameStructure = defaultOuputNameStruct(opt, result)
+  %
+  % Short description of what the function does goes here.
+  %
+  % USAGE::
+  %
+  %   matlabbatch = setBatchSubjectLevelResults(matlabbatch, opt, subLabel, funcFWHM, iNode, iCon)
+  %
+  % :param opt:
+  % :type opt: structure
+  % :param result:
+  % :type result: structure
+  %
+  % :returns: - :outputNameStructure: (structure)
+  %
+  % See also: setBatchSubjectLevelResults, bidsResults
+  %
+  % (C) Copyright 2021 CPP_SPM developers
+
+  outputNameStructure = struct( ...
+                               'suffix', 'spmT', ...
+                               'ext', '.nii', ...
+                               'use_schema', 'false', ...
+                               'entities', struct('sub', '', ...
+                                                  'task', strjoin(opt.taskName, ''), ...
+                                                  'space', result.space, ...
+                                                  'desc', '', ...
+                                                  'label', sprintf('%04.0f', result.contrastNb), ...
+                                                  'p', '', ...
+                                                  'k', '', ...
+                                                  'MC', ''));
+
+end

src/utils/regexify.m

@@ -1,8 +1,16 @@
 function string = regexify(string)
   %
-  % Turns a string into a simple regex.
+  % Turns a string into a simple regex. Useful to query bids dataset with
+  % bids.query that by default expects will treat its inputs as regex.
+  %
+  %   Input   -->    Output
+  %
+  %   ``foo`` --> ``^foo$``
+  %
+  % USAGE::
+  %
+  %   string = regexify(string)
   %
-  %   ``foo``  --> ``^foo$``
   %
   % (C) Copyright 2021 CPP_SPM developers
 

src/workflows/bidsFFX.m

@@ -13,7 +13,7 @@
   %
   % USAGE::
   %
-  %  bidsFFX(action, [opt])
+  %  bidsFFX(action, opt)
   %
   % :param action: Action to be conducted:``specifyAndEstimate`` or ``contrasts``.
   % :type action: string
@@ -24,6 +24,8 @@
   % - ``specifyAndEstimate`` for fMRI design + estimate and
   % - ``contrasts`` to estimate contrasts.
   %
+  % See also: setBatchSubjectLevelGLMSpec, setBatchSubjectLevelContrasts
+  %
   % (C) Copyright 2020 CPP_SPM developers
 
   % TODO: get the space to run analysis in from the BIDS stats model input

src/workflows/bidsRFX.m

@@ -4,23 +4,24 @@ function bidsRFX(action, opt)
   %
   % OR
   %
-  % - creates a mean structural image and mean mask over the sample,
+  % - creates a mean structural image and mean mask over the sample
+  %
+  % OR
+  %
   % - specifies and estimates the group level model,
   % - computes the group level contrasts.
   %
   % USAGE::
   %
   %  bidsRFX(action, opt)
   %
-  % :param action: Action to be conducted: ``smoothContrasts`` or ``RFX``
+  % :param action: Action to be conducted: ``'smoothContrasts'`` or ``'RFX'`` or
+  %                ``'meanAnatAndMask'``
   % :type action: string
   % :param opt: structure or json filename containing the options. See
   %             ``checkOptions()`` and ``loadAndCheckOptions()``.
   % :type opt: structure
   %
-  %  - case ``smoothContrasts``: smooth con images
-  %  - case ``RFX``: Mean Struct, MeanMask, Factorial design specification and
-  %    estimation, Contrast estimation
   %
   % (C) Copyright 2020 CPP_SPM developers
 
@@ -36,9 +37,9 @@ function bidsRFX(action, opt)
 
   matlabbatch = {};
 
-  switch action
+  switch lower(action)
 
-    case 'smoothContrasts'
+    case 'smoothcontrasts'
 
       matlabbatch = setBatchSmoothConImages(matlabbatch, opt);
 
@@ -47,7 +48,7 @@ function bidsRFX(action, opt)
                           '_task-', opt.taskName], ...
                          opt);
 
-    case 'MeanAnatAndMask'
+    case 'meananatandmask'
 
       opt.dir.rfx = getRFXdir(opt);
 
@@ -61,10 +62,7 @@ function bidsRFX(action, opt)
                                             fullfile(opt.dir.stats, 'group'));
       saveAndRunWorkflow(matlabbatch, 'create_mean_struc_mask', opt);
 
-    case 'RFX'
-
-      % TODO
-      % saving needs to be improved (maybe??) as the name may vary with FWHM and contrast
+    case 'rfx'
 
       matlabbatch = setBatchFactorialDesign(matlabbatch, opt);
 
@@ -90,7 +88,7 @@ function checks(opt, action)
     errorHandling(mfilename(), 'tooManySpaces', msg, false, opt.verbosity);
   end
 
-  allowedActions = {'smoothContrasts', 'MeanAnatAndMask', 'RFX'};
+  allowedActions = {'smoothContrasts', 'meanAnatAndMask', 'RFX'};
   if ~ismember(action, allowedActions)
     msg = sprintf('action must be: %s.\n%s was given.', createUnorderedList(allowedActions), ...
                   action);