diff --git a/CPAC/anat_preproc/anat_preproc.py b/CPAC/anat_preproc/anat_preproc.py index ac6d8716d7..6ad820cdf1 100644 --- a/CPAC/anat_preproc/anat_preproc.py +++ b/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': diff --git a/CPAC/func_preproc/func_preproc.py b/CPAC/func_preproc/func_preproc.py index aca20c47e8..d49a1206db 100644 --- a/CPAC/func_preproc/func_preproc.py +++ b/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', diff --git a/CPAC/nuisance/bandpass.py b/CPAC/nuisance/bandpass.py index ba707f9982..cd3c0c1ae9 100644 --- a/CPAC/nuisance/bandpass.py +++ b/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 \ No newline at end of file diff --git a/CPAC/nuisance/nuisance.py b/CPAC/nuisance/nuisance.py index 28fec4164b..228d39d6c4 100644 --- a/CPAC/nuisance/nuisance.py +++ b/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 diff --git a/CPAC/pipeline/cpac_pipeline.py b/CPAC/pipeline/cpac_pipeline.py index db3719ba9a..c0c1b27626 100644 --- a/CPAC/pipeline/cpac_pipeline.py +++ b/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', diff --git a/CPAC/pipeline/schema.py b/CPAC/pipeline/schema.py index c82181fa68..881a1d88f7 100644 --- a/CPAC/pipeline/schema.py +++ b/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, diff --git a/CPAC/registration/registration.py b/CPAC/registration/registration.py index d8ef74285e..4cfc879ea3 100644 --- a/CPAC/registration/registration.py +++ b/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') diff --git a/CPAC/resources/configs/pipeline_config_monkey.yml b/CPAC/resources/configs/pipeline_config_monkey.yml index 93416e4a6d..55ee2b79e2 100644 --- a/CPAC/resources/configs/pipeline_config_monkey.yml +++ b/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 diff --git a/CPAC/resources/configs/pipeline_config_rodent.yml b/CPAC/resources/configs/pipeline_config_rodent.yml index db90a1423c..b3fc998565 100644 --- a/CPAC/resources/configs/pipeline_config_rodent.yml +++ b/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 diff --git a/CPAC/utils/configuration.py b/CPAC/utils/configuration.py index 9f6634df54..1a1a763a2f 100644 --- a/CPAC/utils/configuration.py +++ b/CPAC/utils/configuration.py @@ -78,6 +78,7 @@ def check_path(key): template_list = ['template_brain_only_for_anat', 'template_skull_for_anat', + 'template_brain_mask_for_anat', 'ref_mask', 'template_brain_only_for_func', 'template_skull_for_func', diff --git a/dev/docker_data/default_pipeline.yml b/dev/docker_data/default_pipeline.yml index aa3bf50427..08af7284a9 100644 --- a/dev/docker_data/default_pipeline.yml +++ b/dev/docker_data/default_pipeline.yml @@ -139,6 +139,11 @@ acpc_run_preprocessing: after acpc_brainsize: 150 +# 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: robustfov + + # ACPC aligned template acpc_template_skull: /usr/share/fsl/5.0/data/standard/MNI152_T1_1mm.nii.gz acpc_template_brain: None @@ -971,6 +976,10 @@ template_brain_mask_for_func : /usr/share/fsl/5.0/data/standard/MNI152_T1_${res template_epi : s3://fcp-indi/resources/cpac/resources/epi_hbn.nii.gz +# EPI template mask. TODO: Ask for mask, it is a fake mask now. +template_epi_mask : s3://fcp-indi/resources/cpac/resources/epi_hbn.nii.gz + + # Matrix containing all 1's. Used as an identity matrix during registration. # It is not necessary to change this path unless you intend to use non-standard MNI registration. identityMatrix : /usr/share/fsl/5.0/etc/flirtsch/ident.mat @@ -1023,6 +1032,7 @@ Regressors : Bandpass: bottom_frequency: 0.01 top_frequency: 0.1 + method: default - Motion: include_delayed: true @@ -1049,6 +1059,7 @@ Regressors : Bandpass: bottom_frequency: 0.01 top_frequency: 0.1 + method: default # Whether to run frequency filtering before or after nuisance regression. diff --git a/dev/docker_data/required_afni_pkgs.txt b/dev/docker_data/required_afni_pkgs.txt index 2784d7b399..45122ad2aa 100644 --- a/dev/docker_data/required_afni_pkgs.txt +++ b/dev/docker_data/required_afni_pkgs.txt @@ -1,4 +1,5 @@ linux_openmp_64/3dAutomask +linux_openmp_64/1dBandpass linux_openmp_64/3dBandpass linux_openmp_64/3dBlurToFWHM linux_openmp_64/3dBrickStat