[ENH] add --no-signal-scaling option (#298)

Author: jdkent

docs/workflows.rst

@@ -26,6 +26,7 @@ Participant Workflow
         output_dir='.',
         preproc_img_list=[''],
         selected_confounds=[''],
+        signal_scaling=0,
         smoothing_kernel=None)
 
 The general workflow for a participant models the beta series
@@ -54,6 +55,7 @@ Least Squares- Separate (LSS)
         hrf_model='glover',
         high_pass=0.008,
         smoothing_kernel=0.0,
+        signal_scaling=0,
         selected_confounds=[''])
 
 nistats is used for modeling using the
@@ -75,6 +77,7 @@ Finite BOLD Response- Separate (FS)
         hrf_model='fir',
         high_pass=0.008,
         smoothing_kernel=0.0,
+        signal_scaling=0,
         selected_confounds=[''])
 
 Additionally, NiBetaSeries can be used to perform
@@ -97,6 +100,7 @@ Least Squares- All (LSA)
         hrf_model='glover',
         high_pass=0.008,
         smoothing_kernel=0.0,
+        signal_scaling=0,
         selected_confounds=[''])
 
 For completeness, NiBetaSeries also implements least squares- all (LSA),

src/nibetaseries/cli/run.py

@@ -102,6 +102,10 @@ def get_parser():
                            'are stored (i.e. non-essential files). '
                            'This directory can be deleted once you are reasonably '
                            'certain nibs finished as expected.')
+    proc_opts.add_argument('--no-signal-scaling', action='store_true', default=False,
+                           help='do not scale every voxel with respect to time meaning'
+                                ' beta estimates will be in the same units as the bold'
+                                ' signal')
 
     # Image Selection options
     bids_opts = parser.add_argument_group('Options for selecting images')
@@ -237,6 +241,12 @@ def main():
     else:
         atlas_img = atlas_lut = None
 
+    # check if --no-signal-scaling is set
+    if opts.no_signal_scaling:
+        signal_scaling = False
+    else:
+        signal_scaling = 0
+
     # running participant level
     if opts.analysis_level == "participant":
         nibetaseries_participant_wf = init_nibetaseries_participant_wf(
@@ -252,6 +262,7 @@ def main():
             high_pass=opts.high_pass,
             output_dir=output_dir,
             run_label=opts.run_label,
+            signal_scaling=signal_scaling,
             selected_confounds=opts.confounds,
             session_label=opts.session_label,
             smoothing_kernel=opts.smoothing_kernel,

src/nibetaseries/cli/tests/test_run.py

@@ -36,14 +36,15 @@ def test_conditional_arguments(monkeypatch):
         get_parser().parse_args(no_img)
 
 
-@pytest.mark.parametrize("use_atlas,estimator,fir_delays,hrf_model,part_label",
-                         [(True, 'lsa', None, 'spm', '01'),
-                          (False, 'lss', None, 'spm', 'sub-01'),
-                          (True, 'lss', [0, 1, 2, 3, 4], 'fir', None)])
+@pytest.mark.parametrize("use_atlas,estimator,fir_delays,hrf_model,part_label,use_signal_scaling",
+                         [(True, 'lsa', None, 'spm', '01', True),
+                          (False, 'lss', None, 'spm', 'sub-01', False),
+                          (True, 'lss', [0, 1, 2, 3, 4], 'fir', None, False)])
 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):
+        estimator, fir_delays, hrf_model, monkeypatch, part_label, use_atlas,
+        use_signal_scaling):
     import sys
     bids_dir = str(bids_dir)
     out_dir = os.path.join(bids_dir, 'derivatives')
@@ -64,6 +65,8 @@ def test_nibs(
             out_dir,
             "participant",
             "-c", ".*derivative.*"]
+    if use_signal_scaling:
+        args.extend(["--no-signal-scaling"])
     if use_atlas:
         args.extend(["-a", str(atlas_file), "-l", str(atlas_lut)])
     if fir_delays:

src/nibetaseries/interfaces/nistats.py

@@ -27,6 +27,9 @@ class LSSBetaSeriesInputSpec(BaseInterfaceInputSpec):
                        desc="File that contains all events from the bold run")
     confounds_file = traits.Either(None, File(exists=True),
                                    desc="File that contains all usable confounds")
+    signal_scaling = traits.Enum(False, 0,
+                                 desc="Whether (0) or not (False) to scale each"
+                                      " voxel's timeseries")
     selected_confounds = traits.Either(None, traits.List(),
                                        desc="Column names of the regressors to include")
     hrf_model = traits.String(desc="hemodynamic response model")
@@ -70,7 +73,7 @@ def _run_interface(self, runtime):
             hrf_model=self.inputs.hrf_model,
             mask=self.inputs.mask_file,
             smoothing_fwhm=self.inputs.smoothing_kernel,
-            signal_scaling=0,
+            signal_scaling=self.inputs.signal_scaling,
             high_pass=self.inputs.high_pass,
             drift_model='cosine',
             verbose=1,
@@ -156,6 +159,9 @@ class LSABetaSeriesInputSpec(BaseInterfaceInputSpec):
                        desc="File that contains all events from the bold run")
     confounds_file = traits.Either(None, File(exists=True),
                                    desc="File that contains all usable confounds")
+    signal_scaling = traits.Enum(False, 0,
+                                 desc="Whether (0) or not (False) to scale each"
+                                      " voxel's timeseries")
     selected_confounds = traits.Either(None, traits.List(),
                                        desc="Column names of the regressors to include")
     hrf_model = traits.String(desc="hemodynamic response model")
@@ -196,7 +202,7 @@ def _run_interface(self, runtime):
             hrf_model=self.inputs.hrf_model,
             mask=self.inputs.mask_file,
             smoothing_fwhm=self.inputs.smoothing_kernel,
-            signal_scaling=0,
+            signal_scaling=self.inputs.signal_scaling,
             high_pass=self.inputs.high_pass,
             drift_model='cosine',
             verbose=1,

src/nibetaseries/interfaces/tests/test_nistats.py

@@ -46,6 +46,7 @@ def test_lss_beta_series(sub_metadata, preproc_file, sub_events,
                                 events_file=str(sub_events),
                                 confounds_file=str(confounds_file),
                                 selected_confounds=selected_confounds,
+                                signal_scaling=0,
                                 hrf_model=hrf_model,
                                 smoothing_kernel=None,
                                 high_pass=0.008)
@@ -98,6 +99,7 @@ def test_fs_beta_series(sub_metadata, preproc_file, sub_events,
                                 events_file=str(sub_events),
                                 confounds_file=str(confounds_file),
                                 selected_confounds=selected_confounds,
+                                signal_scaling=0,
                                 hrf_model=hrf_model,
                                 fir_delays=fir_delays,
                                 smoothing_kernel=None,
@@ -149,6 +151,7 @@ def test_lsa_beta_series(sub_metadata, preproc_file, sub_events,
                                 mask_file=mask_file,
                                 events_file=str(sub_events),
                                 confounds_file=str(confounds_file),
+                                signal_scaling=0,
                                 selected_confounds=selected_confounds,
                                 hrf_model=hrf_model,
                                 smoothing_kernel=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,
-    output_dir, run_label, selected_confounds, session_label, smoothing_kernel,
+    fir_delays, hrf_model, high_pass, output_dir, run_label,
+    selected_confounds, session_label, signal_scaling, smoothing_kernel,
     space_label, subject_list, task_label, description_label, work_dir,
         ):
 
@@ -69,6 +69,8 @@ def init_nibetaseries_participant_wf(
             Include bold series containing this run label
         selected_confounds : list
             List of confounds to be included in regression
+        signal_scaling : False or 0
+            Whether (0) or not (False) to scale each voxel's timeseries
         session_label : str or None
             Include bold series containing this session label
         smoothing_kernel : float or None
@@ -173,6 +175,7 @@ def init_nibetaseries_participant_wf(
             output_dir=output_dir,
             preproc_img_list=preproc_img_list,
             selected_confounds=selected_confounds,
+            signal_scaling=signal_scaling,
             smoothing_kernel=smoothing_kernel,
         )
 
@@ -191,8 +194,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, output_dir,
-    preproc_img_list, selected_confounds, smoothing_kernel
+    name, output_dir, preproc_img_list, selected_confounds,
+    signal_scaling, smoothing_kernel,
         ):
     """
     This workflow completes the generation of the betaseries files
@@ -217,6 +220,7 @@ def init_single_subject_wf(
             output_dir='.',
             preproc_img_list=[''],
             selected_confounds=[''],
+            signal_scaling=0,
             smoothing_kernel=0.0)
 
     Parameters
@@ -249,6 +253,8 @@ def init_single_subject_wf(
             list of preprocessed bold files
         selected_confounds : list or None
             the list of confounds to be included in regression
+        signal_scaling : False or 0
+            Whether (0) or not (False) to scale each voxel's timeseries
         smoothing_kernel : float or None
             the size of the smoothing kernel (full width/half max) applied to the bold file (in mm)
 
@@ -308,6 +314,7 @@ def init_single_subject_wf(
                                        hrf_model=hrf_model,
                                        high_pass=high_pass,
                                        selected_confounds=selected_confounds,
+                                       signal_scaling=signal_scaling,
                                        smoothing_kernel=smoothing_kernel)
 
     # initialize the analysis workflow

src/nibetaseries/workflows/model.py

@@ -23,8 +23,9 @@ def init_betaseries_wf(name="betaseries_wf",
                        fir_delays=None,
                        hrf_model='glover',
                        high_pass=0.0078125,
-                       smoothing_kernel=None,
+                       signal_scaling=0,
                        selected_confounds=None,
+                       smoothing_kernel=None,
                        ):
     """Derives Beta Series Maps
     This workflow derives beta series maps from a bold file.
@@ -54,10 +55,12 @@ def init_betaseries_wf(name="betaseries_wf",
     high_pass : float
         high pass filter to apply to bold (in Hertz).
         Reminder - frequencies _lower_ than this number are kept.
-    smoothing_kernel : float or None
-        The size of the smoothing kernel (full width/half max) applied to the bold file (in mm)
     selected_confounds : list or None
         the list of confounds to be included in regression.
+    signal_scaling : False or 0
+        Whether (0) or not (False) to scale each voxel's timeseries
+    smoothing_kernel : float or None
+        The size of the smoothing kernel (full width/half max) applied to the bold file (in mm)
 
     Inputs
     ------
@@ -89,6 +92,7 @@ def init_betaseries_wf(name="betaseries_wf",
         hrf=hrf_model,
         hpf=high_pass,
         selected_confounds=selected_confounds,
+        signal_scaling=signal_scaling,
         estimator=estimator,
         fir_delays=fir_delays,
     )
@@ -105,13 +109,15 @@ def init_betaseries_wf(name="betaseries_wf",
         betaseries_node = pe.Node(LSSBetaSeries(
                 fir_delays=fir_delays,
                 selected_confounds=selected_confounds,
+                signal_scaling=signal_scaling,
                 hrf_model=hrf_model,
                 smoothing_kernel=smoothing_kernel,
                 high_pass=high_pass),
             name='betaseries_node')
     elif estimator == 'lsa':
         betaseries_node = pe.Node(LSABetaSeries(
                 selected_confounds=selected_confounds,
+                signal_scaling=signal_scaling,
                 hrf_model=hrf_model,
                 smoothing_kernel=smoothing_kernel,
                 high_pass=high_pass),
@@ -134,16 +140,17 @@ def init_betaseries_wf(name="betaseries_wf",
 
 
 def gen_wf_description(nistats_ver, fwhm, hrf, hpf,
-                       selected_confounds, estimator,
-                       fir_delays=None):
+                       selected_confounds, signal_scaling,
+                       estimator, fir_delays=None):
     from textwrap import dedent
 
     smooth_str = ('smoothed with a Gaussian kernel with a FWHM of {fwhm} mm,'
                   ' '.format(fwhm=fwhm)
                   if fwhm != 0. else '')
+    signal_scale_str = ', and mean-scaled over time.' if signal_scaling == 0 else '.'
 
-    preproc_str = ('Prior to modeling, preprocessed data were {smooth_str}masked,'
-                   'and mean-scaled over time.'.format(smooth_str=smooth_str))
+    preproc_str = ('Prior to modeling, preprocessed data were {smooth_str}masked{signal_scale_str}'
+                   .format(smooth_str=smooth_str, signal_scale_str=signal_scale_str))
 
     beta_series_tmp = dedent("""
         After fitting {n_models} model, the parameter estimate (i.e., beta) map

src/nibetaseries/workflows/tests/test_base.py

@@ -8,14 +8,14 @@
 from ..base import init_nibetaseries_participant_wf, _check_bs_len
 
 
-@pytest.mark.parametrize("estimator,fir_delays,hrf_model",
-                         [('lsa', None, 'spm'),
-                          ('lss', None, 'spm'),
-                          ('lss', [0, 1, 2, 3, 4], 'fir')])
+@pytest.mark.parametrize("estimator,fir_delays,hrf_model,signal_scaling",
+                         [('lsa', None, 'spm', 0),
+                          ('lss', None, 'spm', False),
+                          ('lss', [0, 1, 2, 3, 4], 'fir', False)])
 def test_valid_init_nibetaseries_participant_wf(
         bids_dir, deriv_dir, sub_fmriprep, sub_top_metadata, bold_file, preproc_file,
         sub_events, confounds_file, brainmask_file, atlas_file, atlas_lut, bids_db_file,
-        estimator, fir_delays, hrf_model):
+        estimator, fir_delays, hrf_model, signal_scaling):
 
     output_dir = op.join(str(bids_dir), 'derivatives', 'atlasCorr')
     work_dir = op.join(str(bids_dir), 'derivatives', 'work')
@@ -43,6 +43,7 @@ def test_valid_init_nibetaseries_participant_wf(
         run_label=None,
         selected_confounds=['white_matter', 'csf'],
         session_label=None,
+        signal_scaling=signal_scaling,
         smoothing_kernel=None,
         space_label=None,
         subject_list=["01"],
@@ -87,6 +88,7 @@ def test_filters_init_nibetaseries_participant_wf(
             run_label=run_label,
             selected_confounds=['white_matter', 'csf'],
             session_label=session_label,
+            signal_scaling=False,
             smoothing_kernel=None,
             space_label=space_label,
             subject_list=["01"],