[ENH] add --return-residuals (#304)

* update FAQ to say old fmriprep outputs no longer work

* add return_residuals to docs

* add --return-residuals option and connect residuals to output

* add test for --return-residuals

* update other tests with return_residuals parameter

Author: jdkent

docs/FAQ.rst

@@ -7,8 +7,7 @@ Frequently Asked Questions
 
 Does NiBetaSeries work with output from old fMRIPrep (< v1.2.0)?
 ----------------------------------------------------------------
-Yes, NiBetaSeries will currently work with older version of fMRIPrep.
-However, this will not continue to be the case.
+No, NiBetaSeries will not work with older version of fMRIPrep.
 To rename your files, you can use renameOldFMRIPREP_.
 
 .. _renameOldFMRIPREP: https://github.com/HBClab/renameOldFMRIPREP

docs/workflows.rst

@@ -26,6 +26,7 @@ Participant Workflow
         norm_betas=False,
         output_dir='.',
         preproc_img_list=[''],
+        return_residuals=False,
         selected_confounds=[''],
         signal_scaling=0,
         smoothing_kernel=None)

src/nibetaseries/cli/run.py

@@ -61,18 +61,22 @@ def get_parser():
     parser.add_argument('-v', '--version', action='version',
                         version=verstr)
 
-    # Atlas Arguments (Required Options)
-    atlas_args = parser.add_argument_group('Required Atlas Arguments')
-    atlas_args.add_argument('-a', '--atlas-img', action='store',
-                            required=('-l' in sys.argv or '--atlas-lut' in sys.argv),
-                            help='input atlas nifti where each voxel within a "region" '
-                                 'is labeled with the same integer and there is a unique '
-                                 'integer associated with each region of interest.')
-    atlas_args.add_argument('-l', '--atlas-lut', action='store',
-                            required=('-a' in sys.argv or '--atlas-img' in sys.argv),
-                            help='atlas look up table (tsv) formatted with the columns: '
-                                  'index, regions which correspond to the regions in the '
-                                  'nifti file specified by --atlas-img.')
+    # Workflow Arguments
+    wf_args = parser.add_argument_group('Workflow Arguments')
+    wf_args.add_argument('-a', '--atlas-img', action='store',
+                         required=('-l' in sys.argv or '--atlas-lut' in sys.argv),
+                         help='input atlas nifti where each voxel within a "region" '
+                              'is labeled with the same integer and there is a unique '
+                              'integer associated with each region of interest.')
+    wf_args.add_argument('-l', '--atlas-lut', action='store',
+                         required=('-a' in sys.argv or '--atlas-img' in sys.argv),
+                         help='atlas look up table (tsv) formatted with the columns: '
+                              'index, regions which correspond to the regions in the '
+                              'nifti file specified by --atlas-img.')
+    wf_args.add_argument('--return-residuals', action='store_true', default=False,
+                         help='setting this option returns the residuals from the model'
+                              'while straightforward for LSA, for any other methods, take'
+                              'the residuals with a grain of salt.')
 
     # preprocessing options
     proc_opts = parser.add_argument_group('Options for processing')
@@ -265,6 +269,7 @@ def main():
             high_pass=opts.high_pass,
             norm_betas=opts.normalize_betas,
             output_dir=output_dir,
+            return_residuals=opts.return_residuals,
             run_label=opts.run_label,
             signal_scaling=signal_scaling,
             selected_confounds=opts.confounds,

src/nibetaseries/cli/tests/test_run.py

@@ -37,17 +37,18 @@ def test_conditional_arguments(monkeypatch):
 
 
 @pytest.mark.parametrize(
-    "use_atlas,estimator,fir_delays,hrf_model,part_label,use_signal_scaling,norm_betas",
+    ("use_atlas,estimator,fir_delays,hrf_model,part_label,"
+     "use_signal_scaling,norm_betas,return_residuals"),
     [
-        (True, 'lsa', None, 'spm', '01', True, True),
-        (False, 'lss', None, 'spm', 'sub-01', False, False),
-        (True, 'lss', [0, 1, 2, 3, 4], 'fir', None, False, True)
+        (True, 'lsa', None, 'spm', '01', True, True, True),
+        (False, 'lss', None, 'spm', 'sub-01', False, False, False),
+        (True, 'lss', [0, 1, 2, 3, 4], 'fir', None, False, True, True),
     ])
 def test_nibs(
         bids_dir, deriv_dir, sub_fmriprep, sub_metadata, bold_file, preproc_file,
         sub_events, confounds_file, brainmask_file, atlas_file, atlas_lut,
         estimator, fir_delays, hrf_model, monkeypatch, part_label, use_atlas,
-        use_signal_scaling, norm_betas):
+        use_signal_scaling, norm_betas, return_residuals):
     import sys
     bids_dir = str(bids_dir)
     out_dir = os.path.join(bids_dir, 'derivatives')
@@ -79,6 +80,8 @@ def test_nibs(
         args.extend([str(d) for d in fir_delays])
     if part_label:
         args.extend(["--participant-label", part_label])
+    if return_residuals:
+        args.append('--return-residuals')
 
     monkeypatch.setattr(sys, 'argv', args)
     assert main() is None

src/nibetaseries/workflows/base.py

@@ -32,8 +32,8 @@
 def init_nibetaseries_participant_wf(
     estimator, atlas_img, atlas_lut, bids_dir,
     database_path, derivatives_pipeline_dir, exclude_description_label,
-    fir_delays, hrf_model, high_pass, norm_betas, output_dir, run_label,
-    selected_confounds, session_label, signal_scaling, smoothing_kernel,
+    fir_delays, hrf_model, high_pass, norm_betas, output_dir, return_residuals,
+    run_label, selected_confounds, session_label, signal_scaling, smoothing_kernel,
     space_label, subject_list, task_label, description_label, work_dir,
         ):
 
@@ -67,6 +67,9 @@ def init_nibetaseries_participant_wf(
             If True, beta estimates are divided by the square root of their variance
         output_dir : str
             Directory where derivatives are saved
+        return_residuals : bool
+            Output the residuals from the betaseries model into the
+            derivatives directory
         run_label : str or None
             Include bold series containing this run label
         selected_confounds : list
@@ -177,6 +180,7 @@ def init_nibetaseries_participant_wf(
             norm_betas=norm_betas,
             output_dir=output_dir,
             preproc_img_list=preproc_img_list,
+            return_residuals=return_residuals,
             selected_confounds=selected_confounds,
             signal_scaling=signal_scaling,
             smoothing_kernel=smoothing_kernel,
@@ -197,8 +201,8 @@ def init_nibetaseries_participant_wf(
 def init_single_subject_wf(
     estimator, atlas_img, atlas_lut, bold_metadata_list, brainmask_list,
     confound_tsv_list, events_tsv_list, fir_delays, hrf_model, high_pass,
-    name, norm_betas, output_dir, preproc_img_list, selected_confounds,
-    signal_scaling, smoothing_kernel,
+    name, norm_betas, output_dir, preproc_img_list, return_residuals,
+    selected_confounds, signal_scaling, smoothing_kernel,
         ):
     """
     This workflow completes the generation of the betaseries files
@@ -257,6 +261,9 @@ def init_single_subject_wf(
             Directory where derivatives are saved
         preproc_img_list : list
             list of preprocessed bold files
+        return_residuals : bool
+            Output the residuals from the betaseries model into the
+            derivatives directory
         selected_confounds : list or None
             the list of confounds to be included in regression
         signal_scaling : False or 0
@@ -311,7 +318,8 @@ def init_single_subject_wf(
 
     output_node = pe.Node(niu.IdentityInterface(fields=['correlation_matrix',
                                                         'correlation_fig',
-                                                        'betaseries_file']),
+                                                        'betaseries_file',
+                                                        'residual_file']),
                           name='output_node')
 
     # initialize the betaseries workflow
@@ -389,6 +397,20 @@ def init_single_subject_wf(
             (output_node, ds_correlation_fig, [('correlation_fig', 'in_file')]),
         ])
 
+    if return_residuals:
+        ds_residual_file = pe.MapNode(DerivativesDataSink(base_directory=output_dir),
+                                      iterfield=['in_file'],
+                                      name='ds_residual_file')
+
+        workflow.connect([
+            (betaseries_wf, output_node,
+                [('output_node.residual_file', 'residual_file')]),
+            (output_node, ds_residual_file,
+                [('residual_file', 'in_file')]),
+            (input_node, ds_residual_file,
+                [('preproc_img', 'source_file')]),
+        ])
+
     return workflow
 
 

src/nibetaseries/workflows/model.py

@@ -85,7 +85,10 @@ def init_betaseries_wf(name="betaseries_wf",
     betaseries_files
         One file per trial type, with each file being
         as long as the number of events for that trial type.
-        (assuming the number of trials for any trial type is above 2)
+    residual_file
+        The residual time series after running beta series.
+        For LSA this is straight forward, but be cautious when
+        interpreting residuals from LSS.
 
     """
     workflow = Workflow(name=name)
@@ -129,7 +132,8 @@ def init_betaseries_wf(name="betaseries_wf",
                 high_pass=high_pass),
             name='betaseries_node')
 
-    output_node = pe.Node(niu.IdentityInterface(fields=['betaseries_files']),
+    output_node = pe.Node(niu.IdentityInterface(fields=['betaseries_files',
+                                                        'residual_file']),
                           name='output_node')
 
     # main workflow
@@ -139,7 +143,8 @@ def init_betaseries_wf(name="betaseries_wf",
                                        ('bold_mask_file', 'mask_file'),
                                        ('bold_metadata', 'bold_metadata'),
                                        ('confounds_file', 'confounds_file')]),
-        (betaseries_node, output_node, [('beta_maps', 'betaseries_files')]),
+        (betaseries_node, output_node, [('beta_maps', 'betaseries_files'),
+                                        ('residual', 'residual_file')]),
     ])
 
     return workflow

src/nibetaseries/workflows/tests/test_base.py

@@ -45,6 +45,7 @@ def test_valid_init_nibetaseries_participant_wf(
         high_pass=0.008,
         norm_betas=norm_betas,
         output_dir=output_dir,
+        return_residuals=False,
         run_label=None,
         selected_confounds=['white_matter', 'csf'],
         session_label=None,
@@ -91,6 +92,7 @@ def test_filters_init_nibetaseries_participant_wf(
             high_pass=0.008,
             norm_betas=False,
             output_dir=output_dir,
+            return_residuals=False,
             run_label=run_label,
             selected_confounds=['white_matter', 'csf'],
             session_label=session_label,