diff --git a/.circleci/ds005_outputs.txt b/.circleci/ds005_outputs.txt index 18d51e6c7..d3f07df52 100644 --- a/.circleci/ds005_outputs.txt +++ b/.circleci/ds005_outputs.txt @@ -53,6 +53,7 @@ fmriprep/sub-01/func fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_AROMAnoiseICs.csv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.dtseries.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_bold.dtseries.nii +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_desc-confounds_regressors.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_desc-confounds_regressors.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_desc-MELODIC_mixing.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_space-fsaverage5_hemi-L.func.gii @@ -85,6 +86,7 @@ fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-01_space-T1w_desc-preproc_ fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_AROMAnoiseICs.csv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.dtseries.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.dtseries.nii +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-confounds_regressors.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-confounds_regressors.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-MELODIC_mixing.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-fsaverage5_hemi-L.func.gii diff --git a/.circleci/ds005_partial_outputs.txt b/.circleci/ds005_partial_outputs.txt index 740930cbf..1cc0abed2 100644 --- a/.circleci/ds005_partial_outputs.txt +++ b/.circleci/ds005_partial_outputs.txt @@ -53,6 +53,7 @@ fmriprep/sub-01/func fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_AROMAnoiseICs.csv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.dtseries.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_bold.dtseries.nii +fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-confounds_regressors.json fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-confounds_regressors.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_desc-MELODIC_mixing.tsv fmriprep/sub-01/func/sub-01_task-mixedgamblestask_run-02_space-fsaverage5_hemi-L.func.gii diff --git a/.circleci/ds054_outputs.txt b/.circleci/ds054_outputs.txt index 4939bb6db..17b81eb92 100644 --- a/.circleci/ds054_outputs.txt +++ b/.circleci/ds054_outputs.txt @@ -27,6 +27,7 @@ fmriprep/sub-100185/anat/sub-100185_space-MNI152NLin2009cAsym_label-CSF_probseg. fmriprep/sub-100185/anat/sub-100185_space-MNI152NLin2009cAsym_label-GM_probseg.nii.gz fmriprep/sub-100185/anat/sub-100185_space-MNI152NLin2009cAsym_label-WM_probseg.nii.gz fmriprep/sub-100185/func +fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_desc-confounds_regressors.json fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_desc-confounds_regressors.tsv fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-MNI152NLin2009cAsym_boldref.nii.gz fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-MNI152NLin2009cAsym_desc-brain_mask.json @@ -38,6 +39,7 @@ fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-T1w_desc-brain fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-T1w_desc-brain_mask.nii.gz fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-T1w_desc-preproc_bold.json fmriprep/sub-100185/func/sub-100185_task-machinegame_run-01_space-T1w_desc-preproc_bold.nii.gz +fmriprep/sub-100185/func/sub-100185_task-machinegame_run-02_desc-confounds_regressors.json fmriprep/sub-100185/func/sub-100185_task-machinegame_run-02_desc-confounds_regressors.tsv fmriprep/sub-100185/func/sub-100185_task-machinegame_run-02_space-MNI152NLin2009cAsym_boldref.nii.gz fmriprep/sub-100185/func/sub-100185_task-machinegame_run-02_space-MNI152NLin2009cAsym_desc-brain_mask.json diff --git a/.circleci/ds210_outputs.txt b/.circleci/ds210_outputs.txt index 0f5019ae3..a6ab90994 100644 --- a/.circleci/ds210_outputs.txt +++ b/.circleci/ds210_outputs.txt @@ -27,6 +27,7 @@ fmriprep/sub-02/anat/sub-02_space-MNI152NLin2009cAsym_label-CSF_probseg.nii.gz fmriprep/sub-02/anat/sub-02_space-MNI152NLin2009cAsym_label-GM_probseg.nii.gz fmriprep/sub-02/anat/sub-02_space-MNI152NLin2009cAsym_label-WM_probseg.nii.gz fmriprep/sub-02/func +fmriprep/sub-02/func/sub-02_task-cuedSGT_run-01_desc-confounds_regressors.json fmriprep/sub-02/func/sub-02_task-cuedSGT_run-01_desc-confounds_regressors.tsv fmriprep/sub-02/func/sub-02_task-cuedSGT_run-01_space-MNI152NLin2009cAsym_boldref.nii.gz fmriprep/sub-02/func/sub-02_task-cuedSGT_run-01_space-MNI152NLin2009cAsym_desc-brain_mask.json diff --git a/.zenodo.json b/.zenodo.json index 9b9d987a7..0df51fd7a 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -116,6 +116,11 @@ "affiliation": "Neuroscience Program, University of Iowa", "orcid": "0000-0002-3339-4857" }, + { + "name": "Ciric, Rastko", + "affiliation": "Stanford University", + "orcid": "0000-0001-6347-7939" + }, { "name": "Poldrack, Russell A.", "affiliation": "Department of Psychology, Stanford University", diff --git a/docs/outputs.rst b/docs/outputs.rst index 5dc791c6b..3a02e6347 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -173,7 +173,7 @@ motion-correction parameters estimated by fMRIPrep; if present, ``non_steady_state_outlier_XX`` columns indicate non-steady state volumes with a single ``1`` value and ``0`` elsewhere (*i.e.*, there is one ``non_steady_state_outlier_XX`` column per outlier/volume); -six noise components are calculated using :abbr:`CompCor (Component Based Noise Correction)`, +additional noise components are calculated using :abbr:`CompCor (Component Based Noise Correction)`, according to both the anatomical (``a_comp_cor_XX``) and temporal (``t_comp_cor_XX``) variants; and the motion-related components identified by :abbr:`ICA (independent components analysis)`-:abbr:`AROMA (Automatic Removal Of Motion Artifacts)` diff --git a/docs/workflows.rst b/docs/workflows.rst index a65176cf2..ecfcdcef6 100644 --- a/docs/workflows.rst +++ b/docs/workflows.rst @@ -46,6 +46,9 @@ is presented below: ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}), ('T1w', {}), ('fsnative', {})]), reportlets_dir='.', + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, skull_strip_fixed_seed=False, skull_strip_template='OASIS30ANTs', subject_id='test', @@ -299,6 +302,9 @@ BOLD preprocessing reportlets_dir='.', t2s_coreg=False, use_aroma=False, + regressors_all_comps=False, + regressors_fd_th=0.5, + regressors_dvars_th=1.5, use_bbr=True, use_syn=True, layout=BIDSLayout('.'), @@ -543,7 +549,11 @@ Confounds estimation name="discover_wf", mem_gb=1, metadata={"RepetitionTime": 2.0, - "SliceTiming": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]}) + "SliceTiming": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]}, + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, + ) Given a motion-corrected fMRI, a brain mask, ``mcflirt`` movement parameters and a segmentation, the `discover_wf` sub-workflow calculates potential diff --git a/fmriprep/__about__.py b/fmriprep/__about__.py index 32a42a93d..5e6e829ee 100644 --- a/fmriprep/__about__.py +++ b/fmriprep/__about__.py @@ -86,7 +86,8 @@ 'indexed_gzip>=0.8.8', 'nibabel>=2.2.1', 'nilearn!=0.5.0,!=0.5.1', - 'nipype>=1.1.6', + 'nipype @ git+https://github.com/nipy/nipype.git@' + 'd353f0d879826031334b09d33e9443b8c9b3e7fe', 'nitime', 'niworkflows<0.10.0a0,>=0.9.1.post1', 'numpy', @@ -103,6 +104,8 @@ LINKS_REQUIRES = [ + 'git+https://github.com/nipy/nipype.git@' + 'd353f0d879826031334b09d33e9443b8c9b3e7fe#egg=nipype', ] TESTS_REQUIRES = [ diff --git a/fmriprep/cli/run.py b/fmriprep/cli/run.py index 71c77535f..db9dd44de 100755 --- a/fmriprep/cli/run.py +++ b/fmriprep/cli/run.py @@ -182,6 +182,22 @@ def get_parser(): help='Exact or maximum number of MELODIC components to estimate ' '(positive = exact, negative = maximum)') + # Confounds options + g_confounds = parser.add_argument_group('Specific options for estimating confounds') + g_confounds.add_argument( + '--return-all-components', required=False, action='store_true', default=False, + help='Include all components estimated in CompCor decomposition in the confounds ' + 'file instead of only the components sufficient to explain 50 percent of ' + 'BOLD variance in each CompCor mask') + g_confounds.add_argument( + '--fd-spike-threshold', required=False, action='store', default=0.5, type=float, + help='Threshold for flagging a frame as an outlier on the basis of framewise ' + 'displacement') + g_confounds.add_argument( + '--dvars-spike-threshold', required=False, action='store', default=1.5, type=float, + help='Threshold for flagging a frame as an outlier on the basis of standardised ' + 'DVARS') + # ANTs options g_ants = parser.add_argument_group('Specific options for ANTs registrations') g_ants.add_argument('--skull-strip-template', action='store', default='OASIS30ANTs', @@ -574,6 +590,9 @@ def build_workflow(opts, retval): output_dir=str(output_dir), output_spaces=output_spaces, run_uuid=run_uuid, + regressors_all_comps=opts.return_all_components, + regressors_fd_th=opts.fd_spike_threshold, + regressors_dvars_th=opts.dvars_spike_threshold, skull_strip_fixed_seed=opts.skull_strip_fixed_seed, skull_strip_template=opts.skull_strip_template, subject_list=subject_list, diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib index bcf0d1107..8a391aee4 100644 --- a/fmriprep/data/boilerplate.bib +++ b/fmriprep/data/boilerplate.bib @@ -241,6 +241,20 @@ @article{power_fd_dvars year = 2014 } +@article{confounds_satterthwaite_2013, + author = {Satterthwaite, Theodore D. and Elliott, Mark A. and Gerraty, Raphael T. and Ruparel, Kosha and Loughead, James and Calkins, Monica E. and Eickhoff, Simon B. and Hakonarson, Hakon and Gur, Ruben C. and Gur, Raquel E. and Wolf, Daniel H.}, + doi = {10.1016/j.neuroimage.2012.08.052}, + issn = {10538119}, + journal = {NeuroImage}, + number = 1, + pages = {240--256}, + title = {{An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data}}, + url = {http://linkinghub.elsevier.com/retrieve/pii/S1053811912008609}, + volume = 64, + year = 2013 +} + + @article{nilearn, author = {Abraham, Alexandre and Pedregosa, Fabian and Eickenberg, Michael and Gervais, Philippe and Mueller, Andreas and Kossaifi, Jean and Gramfort, Alexandre and Thirion, Bertrand and Varoquaux, Gael}, doi = {10.3389/fninf.2014.00014}, diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 07bffbdb7..797043673 100755 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -55,6 +55,9 @@ def init_fmriprep_wf( omp_nthreads, output_dir, output_spaces, + regressors_all_comps, + regressors_dvars_th, + regressors_fd_th, run_uuid, skull_strip_fixed_seed, skull_strip_template, @@ -105,6 +108,9 @@ def init_fmriprep_wf( output_spaces=OrderedDict([ ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}), ('T1w', {}), ('fsnative', {})]), + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, run_uuid='X', skull_strip_fixed_seed=False, skull_strip_template='OASIS30ANTs', @@ -120,43 +126,44 @@ def init_fmriprep_wf( Parameters - layout : BIDSLayout object - BIDS dataset layout - subject_list : list - List of subject labels - task_id : str or None - Task ID of BOLD series to preprocess, or ``None`` to preprocess all + anat_only : bool + Disable functional workflows + bold2t1w_dof : 6, 9 or 12 + Degrees-of-freedom for BOLD-T1w registration + cifti_output : bool + Generate bold CIFTI file in output spaces + debug : bool + Enable debugging outputs echo_idx : int or None Index of echo to preprocess in multiecho BOLD series, or ``None`` to preprocess all - run_uuid : str - Unique identifier for execution instance - work_dir : str - Directory in which to store workflow execution state and temporary files - output_dir : str - Directory in which to save derivatives + err_on_aroma_warn : bool + Do not fail on ICA-AROMA errors + fmap_bspline : bool + **Experimental**: Fit B-Spline field using least-squares + fmap_demean : bool + Demean voxel-shift map during unwarp + force_syn : bool + **Temporary**: Always run SyN-based SDC + freesurfer : bool + Enable FreeSurfer surface reconstruction (may increase runtime) + hires : bool + Enable sub-millimeter preprocessing in FreeSurfer ignore : list Preprocessing steps to skip (may include "slicetiming", "fieldmaps") - debug : bool - Enable debugging outputs - low_mem : bool - Write uncompressed .nii files in some cases to reduce memory usage - anat_only : bool - Disable functional workflows + layout : BIDSLayout object + BIDS dataset layout longitudinal : bool Treat multiple sessions as longitudinal (may increase runtime) See sub-workflows for specific differences - t2s_coreg : bool - For multi-echo EPI, use the calculated T2*-map for T2*-driven coregistration + low_mem : bool + Write uncompressed .nii files in some cases to reduce memory usage + medial_surface_nan : bool + Replace medial wall values with NaNs on functional GIFTI files omp_nthreads : int Maximum number of threads an individual process may use - skull_strip_template : str - Name of ANTs skull-stripping template ('OASIS30ANTs' or 'NKI') - skull_strip_fixed_seed : bool - Do not use a random seed for skull-stripping - will ensure - run-to-run replicability when used with --omp-nthreads 1 - freesurfer : bool - Enable FreeSurfer surface reconstruction (may increase runtime) + output_dir : str + Directory in which to save derivatives output_spaces : OrderedDict Ordered dictionary where keys are TemplateFlow ID strings (e.g. ``MNI152Lin``, ``MNI152NLin6Asym``, ``MNI152NLin2009cAsym``, or ``fsLR``) strings designating @@ -165,30 +172,35 @@ def init_fmriprep_wf( Values of the dictionary aggregate modifiers (e.g. the value for the key ``MNI152Lin`` could be ``{'resolution': 2}`` if one wants the resampling to be done on the 2mm resolution version of the selected template). - medial_surface_nan : bool - Replace medial wall values with NaNs on functional GIFTI files - cifti_output : bool - Generate bold CIFTI file in output spaces - hires : bool - Enable sub-millimeter preprocessing in FreeSurfer + regressors_all_comps + Return all CompCor component time series instead of the top fraction + regressors_dvars_th + Criterion for flagging DVARS outliers + regressors_fd_th + Criterion for flagging framewise displacement outliers + run_uuid : str + Unique identifier for execution instance + skull_strip_template : str + Name of ANTs skull-stripping template ('OASIS30ANTs' or 'NKI') + skull_strip_fixed_seed : bool + Do not use a random seed for skull-stripping - will ensure + run-to-run replicability when used with --omp-nthreads 1 + subject_list : list + List of subject labels + t2s_coreg : bool + For multi-echo EPI, use the calculated T2*-map for T2*-driven coregistration + task_id : str or None + Task ID of BOLD series to preprocess, or ``None`` to preprocess all + use_aroma : bool + Perform ICA-AROMA on MNI-resampled functional series use_bbr : bool or None Enable/disable boundary-based registration refinement. If ``None``, test BBR result for distortion before accepting. - bold2t1w_dof : 6, 9 or 12 - Degrees-of-freedom for BOLD-T1w registration - fmap_bspline : bool - **Experimental**: Fit B-Spline field using least-squares - fmap_demean : bool - Demean voxel-shift map during unwarp use_syn : bool **Experimental**: Enable ANTs SyN-based susceptibility distortion correction (SDC). If fieldmaps are present and enabled, this is not run, by default. - force_syn : bool - **Temporary**: Always run SyN-based SDC - use_aroma : bool - Perform ICA-AROMA on MNI-resampled functional series - err_on_aroma_warn : bool - Do not fail on ICA-AROMA errors + work_dir : str + Directory in which to store workflow execution state and temporary files """ fmriprep_wf = Workflow(name='fmriprep_wf') @@ -206,36 +218,39 @@ def init_fmriprep_wf( reportlets_dir = os.path.join(work_dir, 'reportlets') for subject_id in subject_list: single_subject_wf = init_single_subject_wf( - layout=layout, - subject_id=subject_id, - task_id=task_id, + anat_only=anat_only, + aroma_melodic_dim=aroma_melodic_dim, + bold2t1w_dof=bold2t1w_dof, + cifti_output=cifti_output, + debug=debug, echo_idx=echo_idx, - name="single_subject_" + subject_id + "_wf", - reportlets_dir=reportlets_dir, - output_dir=output_dir, + err_on_aroma_warn=err_on_aroma_warn, + fmap_bspline=fmap_bspline, + fmap_demean=fmap_demean, + force_syn=force_syn, + freesurfer=freesurfer, + hires=hires, ignore=ignore, - debug=debug, - low_mem=low_mem, - anat_only=anat_only, + layout=layout, longitudinal=longitudinal, - t2s_coreg=t2s_coreg, + low_mem=low_mem, + medial_surface_nan=medial_surface_nan, + name="single_subject_" + subject_id + "_wf", omp_nthreads=omp_nthreads, - skull_strip_template=skull_strip_template, - skull_strip_fixed_seed=skull_strip_fixed_seed, - freesurfer=freesurfer, + output_dir=output_dir, output_spaces=output_spaces, - medial_surface_nan=medial_surface_nan, - cifti_output=cifti_output, - hires=hires, + regressors_all_comps=regressors_all_comps, + regressors_dvars_th=regressors_dvars_th, + regressors_fd_th=regressors_fd_th, + reportlets_dir=reportlets_dir, + skull_strip_fixed_seed=skull_strip_fixed_seed, + skull_strip_template=skull_strip_template, + subject_id=subject_id, + t2s_coreg=t2s_coreg, + task_id=task_id, + use_aroma=use_aroma, use_bbr=use_bbr, - bold2t1w_dof=bold2t1w_dof, - fmap_bspline=fmap_bspline, - fmap_demean=fmap_demean, use_syn=use_syn, - force_syn=force_syn, - use_aroma=use_aroma, - aroma_melodic_dim=aroma_melodic_dim, - err_on_aroma_warn=err_on_aroma_warn, ) single_subject_wf.config['execution']['crashdump_dir'] = ( @@ -275,6 +290,9 @@ def init_single_subject_wf( output_dir, output_spaces, reportlets_dir, + regressors_all_comps, + regressors_dvars_th, + regressors_fd_th, skull_strip_fixed_seed, skull_strip_template, subject_id, @@ -326,6 +344,9 @@ def init_single_subject_wf( ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}), ('T1w', {}), ('fsnative', {})]), reportlets_dir='.', + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, skull_strip_fixed_seed=False, skull_strip_template='OASIS30ANTs', subject_id='test', @@ -392,6 +413,12 @@ def init_single_subject_wf( resolution version of the selected template). reportlets_dir : str Directory in which to save reportlets + regressors_all_comps + Return all CompCor component time series instead of the top fraction + regressors_fd_th + Criterion for flagging framewise displacement outliers + regressors_dvars_th + Criterion for flagging DVARS outliers skull_strip_fixed_seed : bool Do not use a random seed for skull-stripping - will ensure run-to-run replicability when used with --omp-nthreads 1 @@ -412,6 +439,7 @@ def init_single_subject_wf( **Experimental**: Enable ANTs SyN-based susceptibility distortion correction (SDC). If fieldmaps are present and enabled, this is not run, by default. + Inputs subjects_dir @@ -560,6 +588,9 @@ def init_single_subject_wf( output_dir=output_dir, output_spaces=output_spaces, reportlets_dir=reportlets_dir, + regressors_all_comps=regressors_all_comps, + regressors_fd_th=regressors_fd_th, + regressors_dvars_th=regressors_dvars_th, t2s_coreg=t2s_coreg, use_aroma=use_aroma, use_bbr=use_bbr, diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index 716e64077..fe7f3aafa 100755 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -69,6 +69,9 @@ def init_func_preproc_wf( omp_nthreads, output_dir, output_spaces, + regressors_all_comps, + regressors_dvars_th, + regressors_fd_th, reportlets_dir, t2s_coreg, use_aroma, @@ -106,6 +109,9 @@ def init_func_preproc_wf( output_spaces=OrderedDict([ ('MNI152Lin', {}), ('fsaverage', {'density': '10k'}), ('T1w', {}), ('fsnative', {})]), + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, reportlets_dir='.', t2s_coreg=False, use_aroma=False, @@ -157,6 +163,12 @@ def init_func_preproc_wf( Values of the dictionary aggregate modifiers (e.g. the value for the key ``MNI152Lin`` could be ``{'resolution': 2}`` if one wants the resampling to be done on the 2mm resolution version of the selected template). + regressors_all_comps + Return all CompCor component time series instead of the top fraction + regressors_dvars_th + Criterion for flagging DVARS outliers + regressors_fd_th + Criterion for flagging framewise displacement outliers reportlets_dir : str Absolute path of a directory in which reportlets will be temporarily stored t2s_coreg : bool @@ -384,7 +396,8 @@ def init_func_preproc_wf( fields=['bold_t1', 'bold_t1_ref', 'bold_mask_t1', 'bold_aseg_t1', 'bold_aparc_t1', 'bold_std', 'bold_std_ref' 'bold_mask_std', 'bold_aseg_std', 'bold_aparc_std', 'bold_cifti', 'cifti_variant', 'cifti_variant_key', 'confounds', 'surfaces', - 'aroma_noise_ics', 'melodic_mix', 'nonaggr_denoised_file']), + 'aroma_noise_ics', 'melodic_mix', 'nonaggr_denoised_file', + 'confounds_metadata']), name='outputnode') # BOLD buffer: an identity used as a pointer to either the original BOLD @@ -433,7 +446,8 @@ def init_func_preproc_wf( ('nonaggr_denoised_file', 'inputnode.nonaggr_denoised_file'), ('bold_cifti', 'inputnode.bold_cifti'), ('cifti_variant', 'inputnode.cifti_variant'), - ('cifti_variant_key', 'inputnode.cifti_variant_key') + ('cifti_variant_key', 'inputnode.cifti_variant_key'), + ('confounds_metadata', 'inputnode.confounds_metadata'), ]), ]) @@ -475,6 +489,9 @@ def init_func_preproc_wf( bold_confounds_wf = init_bold_confs_wf( mem_gb=mem_gb['largemem'], metadata=metadata, + regressors_all_comps=regressors_all_comps, + regressors_fd_th=regressors_fd_th, + regressors_dvars_th=regressors_dvars_th, name='bold_confounds_wf') bold_confounds_wf.get_node('inputnode').inputs.t1_transform_flags = [False] @@ -625,6 +642,9 @@ def init_func_preproc_wf( (bold_confounds_wf, outputnode, [ ('outputnode.confounds_file', 'confounds'), ]), + (bold_confounds_wf, outputnode, [ + ('outputnode.confounds_metadata', 'confounds_metadata'), + ]), # Connect bold_bold_trans_wf (bold_split, bold_bold_trans_wf, [ ('out_files', 'inputnode.bold_file')]), diff --git a/fmriprep/workflows/bold/confounds.py b/fmriprep/workflows/bold/confounds.py index e0a740574..bfd61da74 100755 --- a/fmriprep/workflows/bold/confounds.py +++ b/fmriprep/workflows/bold/confounds.py @@ -15,6 +15,7 @@ from templateflow.api import get as get_template from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from niworkflows.interfaces.confounds import ExpandModel, SpikeRegressors from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms from niworkflows.interfaces.images import SignalExtraction from niworkflows.interfaces.masks import ROIsPlot @@ -23,9 +24,12 @@ RobustACompCor as ACompCor, RobustTCompCor as TCompCor, ) +from niworkflows.interfaces.plotting import ( + CompCorVariancePlot, ConfoundsCorrelationPlot +) from niworkflows.interfaces.segmentation import ICA_AROMARPT from niworkflows.interfaces.utils import ( - TPM2ROI, AddTPMs, AddTSVHeader + TPM2ROI, AddTPMs, AddTSVHeader, TSV2JSON, DictMerge ) from ...interfaces import ( @@ -36,7 +40,14 @@ DEFAULT_MEMORY_MIN_GB = 0.01 -def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): +def init_bold_confs_wf( + mem_gb, + metadata, + regressors_all_comps, + regressors_dvars_th, + regressors_fd_th, + name="bold_confs_wf", +): """ This workflow calculates confounds for a BOLD series, and aggregates them into a :abbr:`TSV (tab-separated value)` file, for use as nuisance @@ -70,7 +81,11 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): from fmriprep.workflows.bold.confounds import init_bold_confs_wf wf = init_bold_confs_wf( mem_gb=1, - metadata={}) + metadata={}, + regressors_all_comps=False, + regressors_dvars_th=1.5, + regressors_fd_th=0.5, + ) **Parameters** @@ -82,6 +97,15 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): BIDS metadata for BOLD file name : str Name of workflow (default: ``bold_confs_wf``) + regressors_all_comps: bool + Indicates whether CompCor decompositions should return all + components instead of the minimal number of components necessary + to explain 50 percent of the variance in the decomposition mask. + regressors_dvars_th + Criterion for flagging DVARS outliers + regressors_fd_th + Criterion for flagging framewise displacement outliers + **Inputs** @@ -109,6 +133,8 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): rois_report Reportlet visualizing white-matter/CSF mask used for aCompCor, the ROI for tCompCor and the BOLD brain mask. + confounds_metadata + Confounds metadata dictionary. """ workflow = Workflow(name=name) @@ -126,23 +152,29 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): *preprocessed BOLD* time-series (using a discrete cosine filter with 128s cut-off) for the two *CompCor* variants: temporal (tCompCor) and anatomical (aCompCor). -Six tCompCor components are then calculated from the top 5% variable +tCompCor components are then calculated from the top 5% variable voxels within a mask covering the subcortical regions. This subcortical mask is obtained by heavily eroding the brain mask, which ensures it does not include cortical GM regions. -For aCompCor, six components are calculated within the intersection of +For aCompCor, components are calculated within the intersection of the aforementioned mask and the union of CSF and WM masks calculated in T1w space, after their projection to the native space of each -functional run (using the inverse BOLD-to-T1w transformation). +functional run (using the inverse BOLD-to-T1w transformation). Components +are also calculated separately within the WM and CSF masks. The head-motion estimates calculated in the correction step were also placed within the corresponding confounds file. -""" +The confound time series derived from head motion estimates and global +signals were expanded with the inclusion of temporal derivatives and +quadratic terms for each [@confounds_satterthwaite_2013]. +Frames that exceeded a threshold of {fd} mm FD or {dv} standardised DVARS +were classified as motion outliers. +""".format(fd=regressors_fd_th, dv=regressors_dvars_th) inputnode = pe.Node(niu.IdentityInterface( 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']), + fields=['confounds_file', 'confounds_metadata']), name='outputnode') # Get masks ready in T1w space @@ -180,16 +212,28 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): name="fdisp", mem_gb=mem_gb) # a/t-CompCor + mrg_lbl_cc = pe.Node(niu.Merge(3), name='merge_rois_cc', run_without_submitting=True) + tcompcor = pe.Node( TCompCor(components_file='tcompcor.tsv', header_prefix='t_comp_cor_', pre_filter='cosine', - save_pre_filter=True, percentile_threshold=.05), + save_pre_filter=True, save_metadata=True, percentile_threshold=.05, + failure_mode='NaN'), name="tcompcor", mem_gb=mem_gb) acompcor = pe.Node( ACompCor(components_file='acompcor.tsv', header_prefix='a_comp_cor_', pre_filter='cosine', - save_pre_filter=True), + save_pre_filter=True, save_metadata=True, mask_names=['combined', 'CSF', 'WM'], + merge_method='none', failure_mode='NaN'), name="acompcor", mem_gb=mem_gb) + # Set number of components + if regressors_all_comps: + acompcor.inputs.num_components = 'all' + tcompcor.inputs.num_components = 'all' + else: + acompcor.inputs.variance_threshold = 0.5 + tcompcor.inputs.variance_threshold = 0.5 + # Set TR if present if 'RepetitionTime' in metadata: tcompcor.inputs.repetition_time = metadata['RepetitionTime'] @@ -212,7 +256,32 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): name="add_motion_headers", mem_gb=0.01, run_without_submitting=True) concat = pe.Node(GatherConfounds(), name="concat", mem_gb=0.01, run_without_submitting=True) - # Generate reportlet + # CompCor metadata + tcc_metadata_fmt = pe.Node( + TSV2JSON(index_column='component', drop_columns=['mask'], output=None, + additional_metadata={'Method': 'tCompCor'}, enforce_case=True), + name='tcc_metadata_fmt') + acc_metadata_fmt = pe.Node( + TSV2JSON(index_column='component', output=None, + additional_metadata={'Method': 'aCompCor'}, enforce_case=True), + name='acc_metadata_fmt') + mrg_conf_metadata = pe.Node(niu.Merge(2), name='merge_confound_metadata', + run_without_submitting=True) + mrg_conf_metadata2 = pe.Node(DictMerge(), name='merge_confound_metadata2', + run_without_submitting=True) + + # Expand model to include derivatives and quadratics + model_expand = pe.Node(ExpandModel( + model_formula='(dd1(rps + wm + csf + gsr))^^2 + others'), + name='model_expansion') + + # Add spike regressors + spike_regress = pe.Node(SpikeRegressors( + fd_thresh=regressors_fd_th, + dvars_thresh=regressors_dvars_th), + name='spike_regressors') + + # Generate reportlet (ROIs) mrg_compcor = pe.Node(niu.Merge(2), name='merge_compcor', run_without_submitting=True) rois_plot = pe.Node(ROIsPlot(colors=['b', 'magenta'], generate_report=True), name='rois_plot', mem_gb=mem_gb) @@ -222,6 +291,27 @@ def init_bold_confs_wf(mem_gb, metadata, name="bold_confs_wf"): name='ds_report_bold_rois', run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB) + # Generate reportlet (CompCor) + mrg_cc_metadata = pe.Node(niu.Merge(2), name='merge_compcor_metadata', + run_without_submitting=True) + compcor_plot = pe.Node( + CompCorVariancePlot(variance_thresholds=(0.5, 0.7, 0.9), + metadata_sources=['tCompCor', 'aCompCor']), + name='compcor_plot') + ds_report_compcor = pe.Node( + DerivativesDataSink(desc='compcorvar', keep_dtype=True), + name='ds_report_compcor', run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB) + + # Generate reportlet (Confound correlation) + conf_corr_plot = pe.Node( + ConfoundsCorrelationPlot(reference_column='global_signal', max_dim=70), + name='conf_corr_plot') + ds_report_conf_corr = pe.Node( + DerivativesDataSink(desc='confoundcorr', keep_dtype=True), + name='ds_report_conf_corr', run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB) + def _pick_csf(files): return files[0] @@ -270,7 +360,10 @@ def _pick_wm(files): (inputnode, acompcor, [('bold', 'realigned_file')]), (inputnode, acompcor, [('skip_vols', 'ignore_initial_volumes')]), (acc_tfm, acc_msk, [('output_image', 'roi_file')]), - (acc_msk, acompcor, [('out', 'mask_files')]), + (acc_msk, mrg_lbl_cc, [('out', 'in1')]), + (csf_msk, mrg_lbl_cc, [('out', 'in2')]), + (wm_msk, mrg_lbl_cc, [('out', 'in3')]), + (mrg_lbl_cc, acompcor, [('out', 'mask_files')]), # Global signals extraction (constrained by anatomy) (inputnode, signals, [('bold', 'in_file')]), @@ -294,14 +387,32 @@ def _pick_wm(files): (add_dvars_header, concat, [('out_file', 'dvars')]), (add_std_dvars_header, concat, [('out_file', 'std_dvars')]), + # Confounds metadata + (tcompcor, tcc_metadata_fmt, [('metadata_file', 'in_file')]), + (acompcor, acc_metadata_fmt, [('metadata_file', 'in_file')]), + (tcc_metadata_fmt, mrg_conf_metadata, [('output', 'in1')]), + (acc_metadata_fmt, mrg_conf_metadata, [('output', 'in2')]), + (mrg_conf_metadata, mrg_conf_metadata2, [('out', 'in_dicts')]), + + # Expand the model with derivatives, quadratics, and spikes + (concat, model_expand, [('confounds_file', 'confounds_file')]), + (model_expand, spike_regress, [('confounds_file', 'confounds_file')]), + # Set outputs - (concat, outputnode, [('confounds_file', 'confounds_file')]), + (spike_regress, outputnode, [('confounds_file', 'confounds_file')]), + (mrg_conf_metadata2, outputnode, [('out_dict', 'confounds_metadata')]), (inputnode, rois_plot, [('bold', 'in_file'), ('bold_mask', 'in_mask')]), (tcompcor, mrg_compcor, [('high_variance_masks', 'in1')]), (acc_msk, mrg_compcor, [('out', 'in2')]), (mrg_compcor, rois_plot, [('out', 'in_rois')]), (rois_plot, ds_report_bold_rois, [('out_report', 'in_file')]), + (tcompcor, mrg_cc_metadata, [('metadata_file', 'in1')]), + (acompcor, mrg_cc_metadata, [('metadata_file', 'in2')]), + (mrg_cc_metadata, compcor_plot, [('out', 'metadata_files')]), + (compcor_plot, ds_report_compcor, [('out_file', 'in_file')]), + (concat, conf_corr_plot, [('confounds_file', 'confounds_file')]), + (conf_corr_plot, ds_report_conf_corr, [('out_file', 'in_file')]), ]) return workflow diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 546587c0a..19b3d3835 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -41,15 +41,12 @@ def init_func_derivatives_wf( from smriprep.workflows.outputs import _bids_relative workflow = Workflow(name=name) - inputnode = pe.Node( - niu.IdentityInterface( - fields=['template', 'source_file', - 'bold_t1', 'bold_t1_ref', 'bold_mask_t1', - 'bold_std', 'bold_std_ref', 'bold_mask_std', - 'bold_aseg_t1', 'bold_aparc_t1', 'bold_aseg_std', - 'bold_aparc_std', 'cifti_variant_key', - 'confounds', 'surfaces', 'aroma_noise_ics', 'melodic_mix', - 'nonaggr_denoised_file', 'bold_cifti', 'cifti_variant']), + inputnode = pe.Node(niu.IdentityInterface(fields=[ + 'aroma_noise_ics', 'bold_aparc_std', 'bold_aparc_t1', 'bold_aseg_std', + 'bold_aseg_t1', 'bold_cifti', 'bold_mask_std', 'bold_mask_t1', 'bold_std', + 'bold_std_ref', 'bold_t1', 'bold_t1_ref', 'cifti_variant', 'cifti_variant_key', + 'confounds', 'confounds_metadata', 'melodic_mix', 'nonaggr_denoised_file', + 'source_file', 'surfaces', 'template']), name='inputnode') raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources') @@ -62,7 +59,8 @@ def init_func_derivatives_wf( workflow.connect([ (inputnode, raw_sources, [('source_file', 'in_files')]), (inputnode, ds_confounds, [('source_file', 'source_file'), - ('confounds', 'in_file')]), + ('confounds', 'in_file'), + ('confounds_metadata', 'meta_dict')]), ]) # Resample to T1w space diff --git a/requirements.txt b/requirements.txt index 5b2ef799e..05bdbc2d3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,5 @@ nilearn!=0.5.0,!=0.5.1 +git+https://github.com/nipy/nipype.git@d353f0d879826031334b09d33e9443b8c9b3e7fe#egg=nipype niworkflows<0.10.0a0,>=0.9.1.post1 smriprep<0.3.0a0,>=0.2.0 templateflow<0.2.0a0,>=0.1.3