optional fov tool for acpc

CPAC/anat_preproc/anat_preproc.py

@@ -55,14 +55,36 @@ def acpc_alignment(skullstrip_tool='afni', config=None, acpc_target='whole-head'
                                                          'acpc_brain_mask']),
                          name='outputspec')
 
-    robust_fov = pe.Node(interface=fsl_utils.RobustFOV(),
-                            name='anat_acpc_1_robustfov')
-    robust_fov.inputs.brainsize = config.acpc_brainsize
-    robust_fov.inputs.out_transform = 'fov_xfm.mat'
+    if config.acpc_FOV_crop == 'robustfov':
+        robust_fov = pe.Node(interface=fsl_utils.RobustFOV(),
+                                name='anat_acpc_1_fov')
+        robust_fov.inputs.brainsize = config.acpc_brainsize
+        robust_fov.inputs.out_transform = 'fov_xfm.mat'
+
+        fov, in_file = (robust_fov, 'in_file')
+        fov, fov_mtx = (robust_fov, 'out_transform')
+        fov, fov_outfile = (robust_fov, 'out_roi')
+
+    elif config.acpc_FOV_crop == 'flirt':
+        # robustfov doesn't work on some monkey data. prefer using flirt.
+        # ${FSLDIR}/bin/flirt -in "${Input}" -applyxfm -ref "${Input}" -omat "$WD"/roi2full.mat -out "$WD"/robustroi.nii.gz
+        # adopted from DCAN NHP https://github.com/DCAN-Labs/dcan-macaque-pipeline/blob/master/PreFreeSurfer/scripts/ACPCAlignment.sh#L80-L81
+        flirt_fov = pe.Node(interface=fsl.FLIRT(),
+                                name='anat_acpc_1_fov')
+        flirt_fov.inputs.args = '-applyxfm'
+
+        fov, in_file = (flirt_fov, 'in_file')
+        fov, ref_file = (flirt_fov, 'reference')
+        fov, fov_mtx = (flirt_fov, 'out_matrix_file')
+        fov, fov_outfile = (flirt_fov, 'out_file')
+    else: 
+        raise 'ACPC alignment FOV crop tool should be either robustfov or flirt. Please check your pipeline config "acpc_FOV_crop" field.'
 
     # align head-to-head to get acpc.mat (for human)
     if acpc_target == 'whole-head':
-        preproc.connect(inputnode, 'anat_leaf', robust_fov, 'in_file')
+        preproc.connect(inputnode, 'anat_leaf', fov, in_file)
+        if config.acpc_FOV_crop == 'flirt':
+            preproc.connect(inputnode, 'anat_leaf', fov, ref_file)
 
     # align brain-to-brain to get acpc.mat (for monkey)
     if acpc_target == 'brain':
@@ -79,13 +101,15 @@ def acpc_alignment(skullstrip_tool='afni', config=None, acpc_target='whole-head'
             preproc.connect(inputnode, 'template_skull_for_anat',
                             initial_skullstrip, 'inputspec.template_skull_for_anat') 
         preproc.connect(initial_skullstrip, 'outputspec.brain', 
-                        robust_fov, 'in_file')
+                        fov, in_file)
+        if config.acpc_FOV_crop == 'flirt':
+            preproc.connect(initial_skullstrip, 'outputspec.brain', fov, ref_file)
 
     convert_fov_xfm = pe.Node(interface=fsl_utils.ConvertXFM(),
                                 name='anat_acpc_2_fov_convertxfm')
     convert_fov_xfm.inputs.invert_xfm = True
 
-    preproc.connect(robust_fov, 'out_transform',
+    preproc.connect(fov, fov_mtx,
                     convert_fov_xfm, 'in_file')
 
     align = pe.Node(interface=fsl.FLIRT(),
@@ -95,7 +119,7 @@ def acpc_alignment(skullstrip_tool='afni', config=None, acpc_target='whole-head'
     align.inputs.searchr_y = [30, 30]
     align.inputs.searchr_z = [30, 30]
 
-    preproc.connect(robust_fov, 'out_roi', align, 'in_file')
+    preproc.connect(fov, fov_outfile, align, 'in_file')
 
     # align head-to-head to get acpc.mat (for human)
     if acpc_target == 'whole-head':

CPAC/func_preproc/func_preproc.py

@@ -904,7 +904,7 @@ def create_func_preproc(skullstrip_tool, motion_correct_tool,
                                                         'TR']),
                          name='inputspec')
 
-    output_node = pe.Node(util.IdentityInterface(fields=['refit',
+    output_node = pe.Node(util.IdentityInterface(fields=['refit', 
                                                          'reorient',
                                                          'reorient_mean',
                                                          'motion_correct',
@@ -915,6 +915,7 @@ def create_func_preproc(skullstrip_tool, motion_correct_tool,
                                                          'mask',
                                                          'skullstrip',
                                                          'func_mean',
+                                                         'get_func_volume',
                                                          'preprocessed',
                                                          'preprocessed_median',
                                                          'preprocessed_mask',

CPAC/nuisance/bandpass.py

@@ -133,4 +133,35 @@ def bandpass_voxels(realigned_file, regressor_file, bandpass_freqs,
                         delimiter='\t')
 
     return bandpassed_file, regressor_bandpassed_file
+
+
+def afni_1dBandpass(in_file, highpass, lowpass, tr=1):
+
+    '''
+    Perform AFNI 1dBandpass
+
+    Parameters
+    ----------
+    in_file : string
+        Path of an input 1D file
+    highpass : float
+        LowCutoff/HighPass
+    lowpass : float
+        HighCutoff/LowPass
 
+    Returns
+    -------
+    out_file : string
+        Path of an output 1D file
+    '''
+
+    import os
+
+    basename = os.path.basename(in_file)
+    filename, file_extension = os.path.splitext(basename)
+    out_file = os.path.join(os.getcwd(), filename + '_bp' + file_extension)
+
+    cmd = '1dBandpass -dt %f %f %f %s > %s' % (tr, highpass, lowpass, in_file, out_file)
+    os.system(cmd)
+
+    return out_file

CPAC/nuisance/nuisance.py

@@ -25,7 +25,7 @@
     cosine_filter,
     TR_string_to_float)
 from CPAC.utils.datasource import check_for_s3
-from .bandpass import bandpass_voxels
+from .bandpass import (bandpass_voxels, afni_1dBandpass)
 
 
 def gather_nuisance(functional_file_path,
@@ -1508,7 +1508,9 @@ def filtering_bold_and_regressors(nuisance_selectors,
 
     inputspec = pe.Node(util.IdentityInterface(fields=[
         'functional_file_path',
-        'regressors_file_path'
+        'regressors_file_path',
+        'functional_brain_mask_file_path',
+        'tr'
     ]), name='inputspec')
 
     outputspec = pe.Node(util.IdentityInterface(fields=['residual_file_path',
@@ -1518,33 +1520,89 @@ def filtering_bold_and_regressors(nuisance_selectors,
     filtering_wf = pe.Workflow(name=name)
     bandpass_selector = nuisance_selectors.get('Bandpass')
 
-    frequency_filter = pe.Node(
-                Function(input_names=['realigned_file',
-                                      'regressor_file',
-                                      'bandpass_freqs',
-                                      'sample_period'],
-                         output_names=['bandpassed_file',
-                                       'regressor_file'],
-                         function=bandpass_voxels,
-                         as_module=True),
-                name='frequency_filter'
-            )
+    if bandpass_selector.get('method'):
+        bandpass_method = bandpass_selector.get('method')
+    else:
+        bandpass_method = 'default'
+
+    if bandpass_method == 'default':
+
+        frequency_filter = pe.Node(
+                    Function(input_names=['realigned_file',
+                                        'regressor_file',
+                                        'bandpass_freqs',
+                                        'sample_period'],
+                            output_names=['bandpassed_file',
+                                        'regressor_file'],
+                            function=bandpass_voxels,
+                            as_module=True),
+                    name='frequency_filter'
+                )
 
-    frequency_filter.inputs.bandpass_freqs = [
-                bandpass_selector.get('bottom_frequency'),
-                bandpass_selector.get('top_frequency')
-            ]
+        frequency_filter.inputs.bandpass_freqs = [
+                    bandpass_selector.get('bottom_frequency'),
+                    bandpass_selector.get('top_frequency')
+                ]
+
+        filtering_wf.connect(inputspec, 'functional_file_path',
+                            frequency_filter, 'realigned_file')
+
+        filtering_wf.connect(inputspec, 'regressors_file_path',
+                            frequency_filter, 'regressor_file')
+
+        filtering_wf.connect(frequency_filter, 'bandpassed_file',
+                            outputspec, 'residual_file_path')
+
+        filtering_wf.connect(frequency_filter, 'regressor_file',
+                            outputspec, 'residual_regressor')
+
+    elif bandpass_method == 'AFNI':
+
+        bandpass_ts = pe.Node(interface=afni.Bandpass(),
+                                    name='bandpass_ts')
+
+        bandpass_ts.inputs.highpass = bandpass_selector.get('bottom_frequency')
+        bandpass_ts.inputs.lowpass = bandpass_selector.get('top_frequency')
+        bandpass_ts.inputs.outputtype = 'NIFTI_GZ'
+
+        tr_string2float_node = pe.Node(util.Function(input_names=['tr'],
+                                                     output_names=['tr_float'],
+                                                     function=TR_string_to_float),
+                                        name='tr_string2float')
+
+        filtering_wf.connect(inputspec, 'tr',
+                            tr_string2float_node, 'tr')
+
+        filtering_wf.connect(tr_string2float_node, 'tr_float',
+                            bandpass_ts, 'tr')
+
+        filtering_wf.connect(inputspec, 'functional_file_path',
+                            bandpass_ts, 'in_file')
+
+        filtering_wf.connect(inputspec, 'functional_brain_mask_file_path',
+                            bandpass_ts, 'mask')
+
+        filtering_wf.connect(bandpass_ts, 'out_file',
+                            outputspec, 'residual_file_path')
+
+        bandpass_regressor = pe.Node(Function(input_names=['in_file',
+                                                           'highpass',
+                                                           'lowpass',
+                                                           'tr'],
+                                              output_names=['out_file'],
+                                              function=afni_1dBandpass),
+                                     name='bandpass_regressor')
 
-    filtering_wf.connect(inputspec, 'functional_file_path',
-                         frequency_filter, 'realigned_file')
+        bandpass_regressor.inputs.highpass = bandpass_selector.get('bottom_frequency')
+        bandpass_regressor.inputs.lowpass = bandpass_selector.get('top_frequency')
 
-    filtering_wf.connect(inputspec, 'regressors_file_path',
-                         frequency_filter, 'regressor_file')
+        filtering_wf.connect(inputspec, 'regressors_file_path',
+                            bandpass_regressor, 'in_file')
 
-    filtering_wf.connect(frequency_filter, 'bandpassed_file',
-                         outputspec, 'residual_file_path')
+        filtering_wf.connect(tr_string2float_node, 'tr_float',
+                            bandpass_regressor, 'tr')
 
-    filtering_wf.connect(frequency_filter, 'regressor_file',
-                         outputspec, 'residual_regressor')
+        filtering_wf.connect(bandpass_regressor, 'out_file',
+                            outputspec, 'residual_regressor')
 
     return filtering_wf

CPAC/pipeline/cpac_pipeline.py

@@ -819,7 +819,9 @@ def build_workflow(subject_id, sub_dict, c, pipeline_name=None, num_ants_cores=1
         (c.resolution_for_func_preproc, c.template_skull_for_func, 'template_skull_for_func_preproc', 'resolution_for_func_preproc'),
         (c.resolution_for_func_preproc, c.template_brain_mask_for_func, 'template_brain_mask_for_func_preproc', 'resolution_for_func_preproc'),
         (c.resolution_for_func_preproc, c.template_epi, 'template_epi', 'resolution_for_func_preproc'),  # no difference of skull and only brain
+        (c.resolution_for_func_preproc, c.template_epi_mask, 'template_epi_mask', 'resolution_for_func_preproc'),  
         (c.resolution_for_func_derivative, c.template_epi, 'template_epi_derivative', 'resolution_for_func_derivative'),  # no difference of skull and only brain
+        (c.resolution_for_func_derivative, c.template_epi_mask, 'template_epi_mask_derivative', 'resolution_for_func_derivative'),  
         (c.resolution_for_func_derivative, c.template_brain_only_for_func, 'template_brain_for_func_derivative', 'resolution_for_func_preproc'),
         (c.resolution_for_func_derivative, c.template_skull_for_func, 'template_skull_for_func_derivative', 'resolution_for_func_preproc'),
     ]
@@ -1398,7 +1400,7 @@ def build_workflow(subject_id, sub_dict, c, pipeline_name=None, num_ants_cores=1
                     ants_reg_anat_mni, 'inputspec.reference_mask'
                     )
 
-                # pass the reference mask file
+                # pass the moving mask file
                 node, out_file = strat['anatomical_brain_mask']
                 workflow.connect(
                     node, out_file,
@@ -2376,11 +2378,41 @@ def build_workflow(subject_id, sub_dict, c, pipeline_name=None, num_ants_cores=1
                     )
 
                     if 'Bandpass' in regressors_selector:
+
+                        bandpass_selector = regressors_selector.get('Bandpass')
+                        if bandpass_selector.get('method'):
+                            bandpass_method = bandpass_selector.get('method')
+                        else:
+                            bandpass_method = 'default'
+
+                        if not bandpass_method in ['default', 'AFNI']:
+                            err = "\n\n[!] C-PAC says: You must choose 'default' or 'AFNI' as filtering method, " \
+                                    "but you provided:\n{0}\n\n".format(bandpass_method)
+                            raise Exception(err)
+
                         if 'Before' in c.filtering_order:
                             nuisance_regression_after_workflow = create_nuisance_regression_workflow(
                                 regressors_selector,
                                 name='nuisance_regression_after-filt_{0}_'
                                      '{1}'.format(regressors_selector_i, num_strat))
+
+                            if bandpass_method == 'AFNI':
+
+                                node, out_file = new_strat['functional_brain_mask']
+
+                                workflow.connect(
+                                    node, out_file,
+                                    filtering,
+                                    'inputspec.functional_brain_mask_file_path'
+                                )
+
+                                node, node_out = new_strat['tr']
+
+                                workflow.connect(
+                                    node, node_out,
+                                    filtering,
+                                    'inputspec.tr'
+                                )
 
                             workflow.connect(
                                 filtering,
@@ -2453,6 +2485,25 @@ def build_workflow(subject_id, sub_dict, c, pipeline_name=None, num_ants_cores=1
                             new_strat.append_name(nuisance_regression_after_workflow.name)
 
                         elif 'After' in c.filtering_order:
+
+                            if bandpass_method == 'AFNI':
+
+                                node, out_file = new_strat['functional_brain_mask']
+
+                                workflow.connect(
+                                    node, out_file,
+                                    filtering,
+                                    'inputspec.functional_brain_mask_file_path'
+                                )
+
+                                node, node_out = new_strat['tr']
+
+                                workflow.connect(
+                                    node, node_out,
+                                    filtering,
+                                    'inputspec.tr'
+                                )
+
                             workflow.connect(
                                 nuisance_regression_before_workflow,
                                 'outputspec.residual_file_path',

CPAC/pipeline/schema.py

@@ -33,6 +33,7 @@
     'resolution_for_anat': All(str, Match(r'^[0-9]+mm$')),
     'template_brain_only_for_anat': str,
     'template_skull_for_anat': str,
+    'template_brain_mask_for_anat': str,
     'template_symmetric_brain_only': str,
     'template_symmetric_skull': str,
     'dilated_symmetric_brain_mask': str,

CPAC/registration/registration.py

@@ -552,6 +552,7 @@ def create_register_func_to_epi(name='register_func_to_epi', reg_option='ANTS',
                                                        'func_3d',
                                                        'func_3d_mask',
                                                        'epi',
+                                                       'epi_mask',
                                                        'interp',
                                                        'ants_para']),
                         name='inputspec')
@@ -582,6 +583,8 @@ def create_register_func_to_epi(name='register_func_to_epi', reg_option='ANTS',
                 ('epi', 'inputspec.reference_brain'),
                 ('func_3d', 'inputspec.moving_skull'),
                 ('epi', 'inputspec.reference_skull'),
+                ('epi_mask', 'inputspec.reference_mask'),
+                ('func_3d_mask', 'inputspec.moving_mask'),
                 ('interp', 'inputspec.interp'),
                 ('ants_para', 'inputspec.ants_para')
             ]),
@@ -1230,6 +1233,9 @@ def connect_func_to_template_reg(workflow, strat_list, c):
                 node, out_file = strat['template_epi']
                 workflow.connect(node, out_file, func_to_epi, 'inputspec.epi')
 
+                node, out_file = strat['template_epi_mask']
+                workflow.connect(node, out_file, func_to_epi, 'inputspec.epi_mask')
+
                 node, out_file = strat['functional_brain_mask']
                 workflow.connect(node, out_file, func_to_epi, 'inputspec.func_3d_mask')
 

CPAC/resources/configs/pipeline_config_monkey.yml

@@ -126,6 +126,11 @@ acpc_align: True
 acpc_brainsize: 70
 
 
+# Choose a tool to crop the FOV in ACPC alignment. 
+# Using FSL's robustfov or flirt command. Default: robustfov for human data, flirt for monkey data. 
+acpc_FOV_crop: flirt
+
+
 # ACPC aligned template
 acpc_template_skull: s3://fcp-indi/resources/cpac/resources/MacaqueYerkes19_T1w_0.5mm.nii.gz
 acpc_template_brain: s3://fcp-indi/resources/cpac/resources/MacaqueYerkes19_T1w_0.5mm_brain.nii.gz

CPAC/resources/configs/pipeline_config_rodent.yml

@@ -803,6 +803,7 @@ Regressors:
 - Bandpass:
     bottom_frequency: 0.01
     top_frequency: 0.1
+    method: AFNI
   CerebrospinalFluid:
       erode_mask: false
       summary: Mean