[FIX] Remove non-steady-state volumes prior to ICA-AROMA

docs/workflows.rst

@@ -588,6 +588,10 @@ be generated, so non-aggressive denoising can be manually performed in the T1w s
         -d sub-<subject_label>_task-<task_id>_bold_MELODICmix.tsv \
         -o sub-<subject_label>_task-<task_id>_bold_space-<space>_AromaNonAggressiveDenoised.nii.gz
 
+*Note*: The non-steady state volumes are removed for the determination of components in melodic.
+Therefore ``*MELODICmix.tsv`` may have zero padded rows to account for the volumes not used in
+melodic's estimation of components.
+
 A visualization of the AROMA component classification is also included in the HTML reports.
 
 .. figure:: _static/aroma.svg

fmriprep/interfaces/confounds.py

@@ -94,6 +94,7 @@ def _run_interface(self, runtime):
 
 class ICAConfoundsInputSpec(BaseInterfaceInputSpec):
     in_directory = Directory(mandatory=True, desc='directory where ICA derivatives are found')
+    skip_vols = traits.Int(desc='number of non steady state volumes identified')
     ignore_aroma_err = traits.Bool(False, usedefault=True, desc='ignore ICA-AROMA errors')
 
 
@@ -111,7 +112,7 @@ class ICAConfounds(SimpleInterface):
 
     def _run_interface(self, runtime):
         aroma_confounds, motion_ics_out, melodic_mix_out = _get_ica_confounds(
-            self.inputs.in_directory, newpath=runtime.cwd)
+            self.inputs.in_directory, self.inputs.skip_vols, newpath=runtime.cwd)
 
         if aroma_confounds is not None:
             self._results['aroma_confounds'] = aroma_confounds
@@ -198,7 +199,7 @@ def _adjust_indices(left_df, right_df):
     return combined_out, confounds_list
 
 
-def _get_ica_confounds(ica_out_dir, newpath=None):
+def _get_ica_confounds(ica_out_dir, skip_vols, newpath=None):
     if newpath is None:
         newpath = os.getcwd()
 
@@ -210,20 +211,21 @@ def _get_ica_confounds(ica_out_dir, newpath=None):
     melodic_mix_out = os.path.join(newpath, 'MELODICmix.tsv')
     motion_ics_out = os.path.join(newpath, 'AROMAnoiseICs.csv')
 
-    # melodic_mix replace spaces with tabs
-    with open(melodic_mix, 'r') as melodic_file:
-        melodic_mix_out_char = melodic_file.read().replace('  ', '\t')
-    # write to output file
-    with open(melodic_mix_out, 'w+') as melodic_file_out:
-        melodic_file_out.write(melodic_mix_out_char)
-
     # copy metion_ics file to derivatives name
     shutil.copyfile(motion_ics, motion_ics_out)
 
     # -1 since python lists start at index 0
     motion_ic_indices = np.loadtxt(motion_ics, dtype=int, delimiter=',', ndmin=1) - 1
     melodic_mix_arr = np.loadtxt(melodic_mix, ndmin=2)
 
+    # pad melodic_mix_arr with rows of zeros corresponding to number non steadystate volumes
+    if skip_vols > 0:
+        zeros = np.zeros([skip_vols, melodic_mix_arr.shape[1]])
+        melodic_mix_arr = np.vstack([zeros, melodic_mix_arr])
+
+    # save melodic_mix_arr
+    np.savetxt(melodic_mix_out, melodic_mix_arr, delimiter='\t')
+
     # Return dummy list of ones if no noise compnents were found
     if motion_ic_indices.size == 0:
         LOGGER.warning('No noise components were classified')

fmriprep/workflows/bold/base.py

@@ -532,6 +532,8 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
             ('outputnode.movpar_file', 'inputnode.movpar_file')]),
         (bold_reg_wf, bold_confounds_wf, [
             ('outputnode.itk_t1_to_bold', 'inputnode.t1_bold_xform')]),
+        (bold_reference_wf, bold_confounds_wf, [
+            ('outputnode.skip_vols', 'inputnode.skip_vols')]),
         (bold_confounds_wf, outputnode, [
             ('outputnode.confounds_file', 'confounds'),
         ]),
@@ -730,6 +732,8 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
                     ('outputnode.bold_mask', 'inputnode.bold_mask')]),
                 (bold_sdc_wf, ica_aroma_wf, [
                     ('outputnode.out_warp', 'inputnode.fieldwarp')]),
+                (bold_reference_wf, ica_aroma_wf, [
+                    ('outputnode.skip_vols', 'inputnode.skip_vols')]),
                 (bold_confounds_wf, join, [
                     ('outputnode.confounds_file', 'in_file')]),
                 (ica_aroma_wf, join,

fmriprep/workflows/bold/confounds.py

@@ -92,6 +92,8 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"):
             BOLD series mask
         movpar_file
             SPM-formatted motion parameters file
+        skip_vols
+            number of non steady state volumes
         t1_mask
             Mask of the skull-stripped template image
         t1_tpms
@@ -136,8 +138,8 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"):
 placed within the corresponding confounds file.
 """
     inputnode = pe.Node(niu.IdentityInterface(
-        fields=['bold', 'bold_mask', 'movpar_file', 't1_mask', 't1_tpms',
-                't1_bold_xform']),
+        fields=['bold', 'bold_mask', 'movpar_file', 'skip_vols',
+                't1_mask', 't1_tpms', 't1_bold_xform']),
         name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(
         fields=['confounds_file']),
@@ -178,7 +180,6 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"):
                     name="fdisp", mem_gb=mem_gb)
 
     # a/t-CompCor
-    non_steady_state = pe.Node(nac.NonSteadyStateDetector(), name='non_steady_state')
     tcompcor = pe.Node(TCompCor(
         components_file='tcompcor.tsv', pre_filter='cosine', save_pre_filter=True,
         percentile_threshold=.05), name="tcompcor", mem_gb=mem_gb)
@@ -250,18 +251,15 @@ def _pick_wm(files):
                             ('bold_mask', 'in_mask')]),
         (inputnode, fdisp, [('movpar_file', 'in_file')]),
 
-        # Calculate nonsteady state
-        (inputnode, non_steady_state, [('bold', 'in_file')]),
-
         # tCompCor
         (inputnode, tcompcor, [('bold', 'realigned_file')]),
-        (non_steady_state, tcompcor, [('n_volumes_to_discard', 'ignore_initial_volumes')]),
+        (inputnode, tcompcor, [('skip_vols', 'ignore_initial_volumes')]),
         (tcc_tfm, tcc_msk, [('output_image', 'roi_file')]),
         (tcc_msk, tcompcor, [('out', 'mask_files')]),
 
         # aCompCor
         (inputnode, acompcor, [('bold', 'realigned_file')]),
-        (non_steady_state, acompcor, [('n_volumes_to_discard', 'ignore_initial_volumes')]),
+        (inputnode, acompcor, [('skip_vols', 'ignore_initial_volumes')]),
         (acc_tfm, acc_msk, [('output_image', 'roi_file')]),
         (acc_msk, acompcor, [('out', 'mask_files')]),
 
@@ -399,6 +397,7 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
 
     The following steps are performed:
 
+    #. Remove non-steady state volumes from the bold series.
     #. Smooth data using FSL `susan`, with a kernel width FWHM=6.0mm.
     #. Run FSL `melodic` outside of ICA-AROMA to generate the report
     #. Run ICA-AROMA
@@ -428,40 +427,53 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
 
     **Parameters**
 
+        aroma_melodic_dim: int or None
+            Set the dimensionality of the Melodic ICA decomposition
+            If None, MELODIC automatically estimates dimensionality.
+        ignore_aroma_err : bool
+            Do not fail on ICA-AROMA errors
+        mem_gb : float
+            Size of BOLD file in GB
+        metadata : dict
+            BIDS metadata for BOLD file
+        name : str
+            Name of workflow (default: ``bold_mni_trans_wf``)
+        omp_nthreads : int
+            Maximum number of threads an individual process may use
         template : str
             Spatial normalization template used as target when that
             registration step was previously calculated with
             :py:func:`~fmriprep.workflows.bold.registration.init_bold_reg_wf`.
             The template must be one of the MNI templates (fMRIPrep uses
             ``MNI152NLin2009cAsym`` by default).
-        metadata : dict
-            BIDS metadata for BOLD file
-        mem_gb : float
-            Size of BOLD file in GB
-        omp_nthreads : int
-            Maximum number of threads an individual process may use
-        name : str
-            Name of workflow (default: ``bold_mni_trans_wf``)
         susan_fwhm : float
             Kernel width (FWHM in mm) for the smoothing step with
             FSL ``susan`` (default: 6.0mm)
         use_fieldwarp : bool
             Include SDC warp in single-shot transform from BOLD to MNI
-        ignore_aroma_err : bool
-            Do not fail on ICA-AROMA errors
-        aroma_melodic_dim: int or None
-            Set the dimensionality of the Melodic ICA decomposition
-            If None, MELODIC automatically estimates dimensionality.
 
 
     **Inputs**
-
-        bold_mni
-            BOLD series, resampled to template space
+        bold_mask
+            BOLD series mask in template space
+        bold_split
+            Individual 3D BOLD volumes, not motion corrected
+        fieldwarp
+            a :abbr:`DFM (displacements field map)` in ITK format
+        hmc_xforms
+            List of affine transforms aligning each volume to ``ref_image`` in ITK format
+        itk_bold_to_t1
+            Affine transform from ``ref_bold_brain`` to T1 space (ITK format)
         movpar_file
             SPM-formatted motion parameters file
-        bold_mask_mni
-            BOLD series mask in template space
+        skip_vols
+            number of non steady state volumes
+        name_source
+            BOLD series NIfTI file
+            Used to recover original information lost during processing
+        t1_2_mni_forward_transform
+            ANTs-compatible affine-and-warp transform file
+
 
     **Outputs**
 
@@ -474,15 +486,15 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
         nonaggr_denoised_file
             BOLD series with non-aggressive ICA-AROMA denoising applied
 
-    .. _ICA-AROMA: https://github.com/rhr-pruim/ICA-AROMA
+    .. _ICA-AROMA: https://github.com/maartenmennes/ICA-AROMA
 
     """
     workflow = Workflow(name=name)
     workflow.__postdesc__ = """\
 Automatic removal of motion artifacts using independent component analysis
 [ICA-AROMA, @aroma] was performed on the *preprocessed BOLD on MNI space*
-time-series after a spatial smoothing with an isotropic, Gaussian kernel
-of 6mm FWHM (full-width half-maximum).
+time-series after removal of non-steady state volumes and spatial smoothing
+with an isotropic, Gaussian kernel of 6mm FWHM (full-width half-maximum).
 Corresponding "non-aggresively" denoised runs were produced after such
 smoothing.
 Additionally, the "aggressive" noise-regressors were collected and placed
@@ -494,6 +506,7 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
             'itk_bold_to_t1',
             't1_2_mni_forward_transform',
             'name_source',
+            'skip_vols',
             'bold_split',
             'bold_mask',
             'hmc_xforms',
@@ -516,6 +529,10 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
     )
     bold_mni_trans_wf.__desc__ = None
 
+    rm_non_steady_state = pe.Node(niu.Function(function=_remove_volumes,
+                                               output_names=['bold_cut']),
+                                  name='rm_nonsteady')
+
     calc_median_val = pe.Node(fsl.ImageStats(op_string='-k %s -p 50'), name='calc_median_val')
     calc_bold_mean = pe.Node(fsl.MeanImage(), name='calc_bold_mean')
 
@@ -539,6 +556,10 @@ def _getusans_func(image, thresh):
         denoise_type='nonaggr', generate_report=True, TR=metadata['RepetitionTime']),
         name='ica_aroma')
 
+    add_non_steady_state = pe.Node(niu.Function(function=_add_volumes,
+                                                output_names=['bold_add']),
+                                   name='add_nonsteady')
+
     # extract the confound ICs from the results
     ica_aroma_confound_extraction = pe.Node(ICAConfounds(ignore_aroma_err=ignore_aroma_err),
                                             name='ica_aroma_confound_extraction')
@@ -562,16 +583,21 @@ def _getbtthresh(medianval):
             ('t1_2_mni_forward_transform', 'inputnode.t1_2_mni_forward_transform'),
             ('fieldwarp', 'inputnode.fieldwarp')]),
         (inputnode, ica_aroma, [('movpar_file', 'motion_parameters')]),
+        (inputnode, rm_non_steady_state, [
+            ('skip_vols', 'skip_vols')]),
+        (bold_mni_trans_wf, rm_non_steady_state, [
+            ('outputnode.bold_mni', 'bold_file')]),
         (bold_mni_trans_wf, calc_median_val, [
-            ('outputnode.bold_mni', 'in_file'),
             ('outputnode.bold_mask_mni', 'mask_file')]),
-        (bold_mni_trans_wf, calc_bold_mean, [
-            ('outputnode.bold_mni', 'in_file')]),
+        (rm_non_steady_state, calc_median_val, [
+            ('bold_cut', 'in_file')]),
+        (rm_non_steady_state, calc_bold_mean, [
+            ('bold_cut', 'in_file')]),
         (calc_bold_mean, getusans, [('out_file', 'image')]),
         (calc_median_val, getusans, [('out_stat', 'thresh')]),
         # Connect input nodes to complete smoothing
-        (bold_mni_trans_wf, smooth, [
-            ('outputnode.bold_mni', 'in_file')]),
+        (rm_non_steady_state, smooth, [
+            ('bold_cut', 'in_file')]),
         (getusans, smooth, [('usans', 'usans')]),
         (calc_median_val, smooth, [(('out_stat', _getbtthresh), 'brightness_threshold')]),
         # connect smooth to melodic
@@ -586,18 +612,63 @@ def _getbtthresh(medianval):
         (melodic, ica_aroma, [('out_dir', 'melodic_dir')]),
         # generate tsvs from ICA-AROMA
         (ica_aroma, ica_aroma_confound_extraction, [('out_dir', 'in_directory')]),
+        (inputnode, ica_aroma_confound_extraction, [
+            ('skip_vols', 'skip_vols')]),
         # output for processing and reporting
         (ica_aroma_confound_extraction, outputnode, [('aroma_confounds', 'aroma_confounds'),
                                                      ('aroma_noise_ics', 'aroma_noise_ics'),
                                                      ('melodic_mix', 'melodic_mix')]),
         # TODO change melodic report to reflect noise and non-noise components
-        (ica_aroma, outputnode, [('nonaggr_denoised_file', 'nonaggr_denoised_file')]),
+        (ica_aroma, add_non_steady_state, [
+            ('nonaggr_denoised_file', 'bold_cut_file')]),
+        (bold_mni_trans_wf, add_non_steady_state, [
+            ('outputnode.bold_mni', 'bold_file')]),
+        (inputnode, add_non_steady_state, [
+            ('skip_vols', 'skip_vols')]),
+        (add_non_steady_state, outputnode, [('bold_add', 'nonaggr_denoised_file')]),
         (ica_aroma, ds_report_ica_aroma, [('out_report', 'in_file')]),
     ])
 
     return workflow
 
 
+def _remove_volumes(bold_file, skip_vols):
+    """remove skip_vols from bold_file"""
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    if skip_vols == 0:
+        return bold_file
+
+    out = fname_presuffix(bold_file, suffix='_cut')
+    bold_img = nb.load(bold_file)
+    bold_img.__class__(bold_img.dataobj[..., skip_vols:],
+                       bold_img.affine, bold_img.header).to_filename(out)
+
+    return out
+
+
+def _add_volumes(bold_file, bold_cut_file, skip_vols):
+    """prepend skip_vols from bold_file onto bold_cut_file"""
+    import nibabel as nb
+    import numpy as np
+    from nipype.utils.filemanip import fname_presuffix
+
+    if skip_vols == 0:
+        return bold_cut_file
+
+    bold_img = nb.load(bold_file)
+    bold_cut_img = nb.load(bold_cut_file)
+
+    bold_data = np.concatenate((bold_img.dataobj[..., :skip_vols],
+                                bold_cut_img.dataobj), axis=3)
+
+    out = fname_presuffix(bold_cut_file, suffix='_addnonsteady')
+    bold_img.__class__(bold_data, bold_img.affine, bold_img.header).to_filename(out)
+
+    return out
+
+
 def _maskroi(in_mask, roi_file):
     import numpy as np
     import nibabel as nb