[skip tests][skip docs][skip ds210][skip ds054] wip: remove non steady state volumes from aroma calculation

Author: jdkent

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')
+    n_volumes = 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.n_volumes, 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, n_volumes, 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 n_volumes > 0:
+        zeros = np.zeros([n_volumes, 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

@@ -18,6 +18,7 @@
 from nipype.interfaces.fsl import Split as FSLSplit
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
+from nipype.algorithms import confounds as nac
 
 from ...interfaces import (
     DerivativesDataSink,
@@ -424,6 +425,9 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
                                              omp_nthreads=omp_nthreads,
                                              use_compression=False)
 
+    # get non steady-state volumes
+    non_steady_state = pe.Node(nac.NonSteadyStateDetector(), name='non_steady_state')
+
     # get confounds
     bold_confounds_wf = init_bold_confs_wf(
         mem_gb=mem_gb['largemem'],
@@ -525,13 +529,18 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
             ('outputnode.out_warp', 'inputnode.fieldwarp'),
             ('outputnode.bold_mask', 'inputnode.bold_mask')]),
         (bold_sdc_wf, summary, [('outputnode.method', 'distortion_correction')]),
+        # Calculate nonsteady state
+        (bold_bold_trans_wf, non_steady_state, [
+            ('outputnode.bold', 'in_file')]),
         # Connect bold_confounds_wf
         (inputnode, bold_confounds_wf, [('t1_tpms', 'inputnode.t1_tpms'),
                                         ('t1_mask', 'inputnode.t1_mask')]),
         (bold_hmc_wf, bold_confounds_wf, [
             ('outputnode.movpar_file', 'inputnode.movpar_file')]),
         (bold_reg_wf, bold_confounds_wf, [
             ('outputnode.itk_t1_to_bold', 'inputnode.t1_bold_xform')]),
+        (non_steady_state, bold_confounds_wf, [
+            ('n_volumes_to_discard', 'inputnode.n_volumes_to_discard')]),
         (bold_confounds_wf, outputnode, [
             ('outputnode.confounds_file', 'confounds'),
         ]),
@@ -730,6 +739,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')]),
+                (non_steady_state, ica_aroma_wf, [
+                    ('n_volumes_to_discard', 'inputnode.n_volumes_to_discard')]),
                 (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
+        n_volumes_to_discard
+            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', 'n_volumes_to_discard', 
+                '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, [('n_volumes_to_discard', '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, [('n_volumes_to_discard', 'ignore_initial_volumes')]),
         (acc_tfm, acc_msk, [('output_image', 'roi_file')]),
         (acc_msk, acompcor, [('out', 'mask_files')]),
 
@@ -494,6 +492,7 @@ def init_ica_aroma_wf(template, metadata, mem_gb, omp_nthreads,
             'itk_bold_to_t1',
             't1_2_mni_forward_transform',
             'name_source',
+            'n_volumes_to_discard',
             'bold_split',
             'bold_mask',
             'hmc_xforms',
@@ -516,6 +515,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 +542,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 +569,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, [
+            ('n_volumes_to_discard', 'n_volumes')]),
+        (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
@@ -585,17 +597,65 @@ 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, [
+            ('n_volumes_to_discard', 'n_volumes')]),
         # 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, [
+            ('n_volumes_to_discard', 'n_volumes')]),
+        (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, n_volumes):
+        import nibabel as nb
+        from nipype.utils.filemanip import fname_presuffix
+
+        # load the bold file and get the 4d matrix
+        bold_img = nb.load(bold_file)
+        bold_data = bold_img.get_data()
+
+        # cut off the beginning volumes
+        bold_data_cut = bold_data[..., n_volumes:]
+
+        # modify header with new shape (fewer volumes)
+        data_shape = list(bold_img.header.get_data_shape())
+        data_shape[-1] -= n_volumes
+        bold_img.header.set_data_shape(tuple(data_shape))
+
+        # save the resulting bold file
+        out = fname_presuffix(bold_file, suffix='_cut')
+        bold_img.__class__(bold_data_cut, bold_img.affine, bold_img.header).to_filename(out)
+        return out
+
+
+def _add_volumes(bold_file, bold_cut_file, n_volumes):
+    """prepend n_volumes from bold_file onto bold_cut_file"""
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    # load the data
+    bold_img = nb.load(bold_file)
+    bold_data = bold_img.get_data()
+    bold_cut_img = nb.load(bold_cut_file)
+    bold_cut_data = bold_cut_img.get_data()
+
+    # assign everything from n_volumes foward to bold_cut_data
+    bold_data[..., n_volumes:] = bold_cut_data
+
+    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