Source code for fmriprep.workflows.bold.confounds

# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
Calculate BOLD confounds

.. autofunction:: init_bold_confs_wf
.. autofunction:: init_ica_aroma_wf

from niworkflows.nipype.pipeline import engine as pe
from niworkflows.nipype.interfaces import utility as niu, fsl
from niworkflows.nipype.interfaces.nilearn import SignalExtraction
from niworkflows.nipype.algorithms import confounds as nac

from niworkflows.interfaces.segmentation import ICA_AROMARPT
from niworkflows.interfaces.masks import ROIsPlot
from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms

from ...interfaces import (
    TPM2ROI, AddTPMs, AddTSVHeader, GatherConfounds, ICAConfounds
from ...interfaces.patches import RobustACompCor as ACompCor

[docs]def init_bold_confs_wf(mem_gb, use_aroma, ignore_aroma_err, metadata, 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 regressors in a :abbr:`GLM (general linear model)`. The following confounds are calculated, with column headings in parentheses: #. Region-wise average signal (``CSF``, ``WhiteMatter``, ``GlobalSignal``) #. DVARS - standard, nonstandard, and voxel-wise standard variants (``stdDVARS``, ``non-stdDVARS``, ``vx-wisestdDVARS``) #. Framewise displacement, based on MCFLIRT motion parameters (``FramewiseDisplacement``) #. Temporal CompCor (``tCompCorXX``) #. Anatomical CompCor (``aCompCorXX``) #. Cosine basis set for high-pass filtering w/ 0.008 Hz cut-off (``CosineXX``) #. Non-steady-state volumes (``NonSteadyStateXX``) #. Estimated head-motion parameters, in mm and rad (``X``, ``Y``, ``Z``, ``RotX``, ``RotY``, ``RotZ``) #. ICA-AROMA-identified noise components, if enabled (``AROMAAggrCompXX``) Prior to estimating aCompCor and tCompCor, non-steady-state volumes are censored and high-pass filtered using a :abbr:`DCT (discrete cosine transform)` basis. The cosine basis, as well as one regressor per censored volume, are included for convenience. .. workflow:: :graph2use: orig :simple_form: yes from fmriprep.workflows.bold.confounds import init_bold_confs_wf wf = init_bold_confs_wf( mem_gb=1, use_aroma=True, ignore_aroma_err=True, metadata={}) **Parameters** mem_gb : float Size of BOLD file in GB - please note that this size should be calculated after resamplings that may extend the FoV use_aroma : bool Perform ICA-AROMA on MNI-resampled functional series ignore_aroma_err : bool Do not fail on ICA-AROMA errors metadata : dict BIDS metadata for BOLD file **Inputs** bold BOLD image, after the prescribed corrections (STC, HMC and SDC) when available. bold_mask BOLD series mask movpar_file SPM-formatted motion parameters file t1_mask Mask of the skull-stripped template image t1_tpms List of tissue probability maps in T1w space t1_bold_xform Affine matrix that maps the T1w space into alignment with the native BOLD space bold_mni BOLD image resampled in MNI space (only if ``use_aroma`` enabled) bold_mask_mni Brain mask corresponding to the BOLD image resampled in MNI space (only if ``use_aroma`` enabled) **Outputs** confounds_file TSV of all aggregated confounds confounds_list List of calculated confounds for reporting acompcor_report Reportlet visualizing white-matter/CSF mask used for aCompCor tcompcor_report Reportlet visualizing ROI identified in tCompCor ica_aroma_report Reportlet visualizing MELODIC ICs, with ICA-AROMA signal/noise labels aroma_noise_ics CSV of noise components identified by ICA-AROMA melodic_mix FSL MELODIC mixing matrix nonaggr_denoised_file BOLD series with non-aggressive ICA-AROMA denoising applied **Subworkflows** * :py:func:`~fmriprep.workflows.bold.confounds.init_ica_aroma_wf` """ inputnode = pe.Node(niu.IdentityInterface( fields=['bold', 'bold_mask', 'movpar_file', 't1_mask', 't1_tpms', 't1_bold_xform', 'bold_mni', 'bold_mask_mni']), name='inputnode') outputnode = pe.Node(niu.IdentityInterface( fields=['confounds_file', 'confounds_list', 'rois_report', 'ica_aroma_report', 'aroma_noise_ics', 'melodic_mix', 'nonaggr_denoised_file']), name='outputnode') # Get masks ready in T1w space acc_tpm = pe.Node(AddTPMs(indices=[0, 2]), name='tpms_add_csf_wm') # acc stands for aCompCor csf_roi = pe.Node(TPM2ROI(erode_mm=0, mask_erode_mm=30), name='csf_roi') wm_roi = pe.Node(TPM2ROI( erode_prop=0.6, mask_erode_prop=0.6**3), # 0.6 = radius; 0.6^3 = volume name='wm_roi') acc_roi = pe.Node(TPM2ROI( erode_prop=0.6, mask_erode_prop=0.6**3), # 0.6 = radius; 0.6^3 = volume name='acc_roi') # Map ROIs in T1w space into BOLD space csf_tfm = pe.Node(ApplyTransforms(interpolation='NearestNeighbor', float=True), name='csf_tfm', mem_gb=0.1) wm_tfm = pe.Node(ApplyTransforms(interpolation='NearestNeighbor', float=True), name='wm_tfm', mem_gb=0.1) acc_tfm = pe.Node(ApplyTransforms(interpolation='NearestNeighbor', float=True), name='acc_tfm', mem_gb=0.1) tcc_tfm = pe.Node(ApplyTransforms(interpolation='NearestNeighbor', float=True), name='tcc_tfm', mem_gb=0.1) # Ensure ROIs don't go off-limits (reduced FoV) csf_msk = pe.Node(niu.Function(function=_maskroi), name='csf_msk') wm_msk = pe.Node(niu.Function(function=_maskroi), name='wm_msk') acc_msk = pe.Node(niu.Function(function=_maskroi), name='acc_msk') tcc_msk = pe.Node(niu.Function(function=_maskroi), name='tcc_msk') # DVARS dvars = pe.Node(nac.ComputeDVARS(save_all=True, remove_zerovariance=True), name="dvars", mem_gb=mem_gb) # Frame displacement fdisp = pe.Node(nac.FramewiseDisplacement(parameter_source="SPM"), name="fdisp", mem_gb=mem_gb) # a/t-CompCor non_steady_state = pe.Node(nac.NonSteadyStateDetector(), name='non_steady_state') tcompcor = pe.Node(nac.TCompCor( components_file='tcompcor.tsv', pre_filter='cosine', save_pre_filter=True, percentile_threshold=.05), name="tcompcor", mem_gb=mem_gb) acompcor = pe.Node(ACompCor( components_file='acompcor.tsv', pre_filter='cosine', save_pre_filter=True), name="acompcor", mem_gb=mem_gb) # Set TR if present if 'RepetitionTime' in metadata: tcompcor.inputs.repetition_time = metadata['RepetitionTime'] acompcor.inputs.repetition_time = metadata['RepetitionTime'] # Global and segment regressors mrg_lbl = pe.Node(niu.Merge(3), name='merge_rois', run_without_submitting=True) signals = pe.Node(SignalExtraction( detrend=True, class_labels=["CSF", "WhiteMatter", "GlobalSignal"]), name="signals", mem_gb=mem_gb) # Arrange confounds add_header = pe.Node(AddTSVHeader(columns=["X", "Y", "Z", "RotX", "RotY", "RotZ"]), name="add_header", mem_gb=0.01, run_without_submitting=True) concat = pe.Node(GatherConfounds(), name="concat", mem_gb=0.01, run_without_submitting=True) # Generate reportlet mrg_compcor = pe.Node(niu.Merge(2), name='merge_compcor', run_without_submitting=True) rois_plot = pe.Node(ROIsPlot(compress_report=True, colors=['r', 'b', 'magenta'], generate_report=True), name='rois_plot') def _pick_csf(files): return files[0] def _pick_wm(files): return files[-1] workflow = pe.Workflow(name=name) workflow.connect([ # Massage ROIs (in T1w space) (inputnode, acc_tpm, [('t1_tpms', 'in_files')]), (inputnode, csf_roi, [(('t1_tpms', _pick_csf), 'in_tpm'), ('t1_mask', 'in_mask')]), (inputnode, wm_roi, [(('t1_tpms', _pick_wm), 'in_tpm'), ('t1_mask', 'in_mask')]), (inputnode, acc_roi, [('t1_mask', 'in_mask')]), (acc_tpm, acc_roi, [('out_file', 'in_tpm')]), # Map ROIs to BOLD (inputnode, csf_tfm, [('bold_mask', 'reference_image'), ('t1_bold_xform', 'transforms')]), (csf_roi, csf_tfm, [('roi_file', 'input_image')]), (inputnode, wm_tfm, [('bold_mask', 'reference_image'), ('t1_bold_xform', 'transforms')]), (wm_roi, wm_tfm, [('roi_file', 'input_image')]), (inputnode, acc_tfm, [('bold_mask', 'reference_image'), ('t1_bold_xform', 'transforms')]), (acc_roi, acc_tfm, [('roi_file', 'input_image')]), (inputnode, tcc_tfm, [('bold_mask', 'reference_image'), ('t1_bold_xform', 'transforms')]), (csf_roi, tcc_tfm, [('eroded_mask', 'input_image')]), # Mask ROIs with bold_mask (inputnode, csf_msk, [('bold_mask', 'in_mask')]), (inputnode, wm_msk, [('bold_mask', 'in_mask')]), (inputnode, acc_msk, [('bold_mask', 'in_mask')]), (inputnode, tcc_msk, [('bold_mask', 'in_mask')]), # connect inputnode to each non-anatomical confound node (inputnode, dvars, [('bold', 'in_file'), ('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')]), (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')]), (acc_tfm, acc_msk, [('output_image', 'roi_file')]), (acc_msk, acompcor, [('out', 'mask_files')]), # Global signals extraction (constrained by anatomy) (inputnode, signals, [('bold', 'in_file')]), (csf_tfm, csf_msk, [('output_image', 'roi_file')]), (csf_msk, mrg_lbl, [('out', 'in1')]), (wm_tfm, wm_msk, [('output_image', 'roi_file')]), (wm_msk, mrg_lbl, [('out', 'in2')]), (inputnode, mrg_lbl, [('bold_mask', 'in3')]), (mrg_lbl, signals, [('out', 'label_files')]), # Collate computed confounds together (inputnode, add_header, [('movpar_file', 'in_file')]), (signals, concat, [('out_file', 'signals')]), (dvars, concat, [('out_all', 'dvars')]), (fdisp, concat, [('out_file', 'fd')]), (tcompcor, concat, [('components_file', 'tcompcor'), ('pre_filter_file', 'cos_basis')]), (acompcor, concat, [('components_file', 'acompcor')]), (add_header, concat, [('out_file', 'motion')]), # Set outputs (concat, outputnode, [('confounds_file', 'confounds_file'), ('confounds_list', 'confounds_list')]), (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, outputnode, [('out_report', 'rois_report')]), ]) if use_aroma: # ICA-AROMA ica_aroma_wf = init_ica_aroma_wf(name='ica_aroma_wf', ignore_aroma_err=ignore_aroma_err) workflow.connect([ (inputnode, ica_aroma_wf, [('bold_mni', 'inputnode.bold_mni'), ('bold_mask_mni', 'inputnode.bold_mask_mni'), ('movpar_file', 'inputnode.movpar_file')]), (ica_aroma_wf, concat, [('outputnode.aroma_confounds', 'aroma')]), (ica_aroma_wf, outputnode, [('outputnode.out_report', 'ica_aroma_report'), ('outputnode.aroma_noise_ics', 'aroma_noise_ics'), ('outputnode.melodic_mix', 'melodic_mix'), ('outputnode.nonaggr_denoised_file', 'nonaggr_denoised_file')]) ]) return workflow
[docs]def init_ica_aroma_wf(name='ica_aroma_wf', ignore_aroma_err=False): ''' This workflow wraps `ICA-AROMA`_ to identify and remove motion-related independent components from a BOLD time series. The following steps are performed: #. Smooth data using SUSAN #. Run MELODIC outside of ICA-AROMA to generate the report #. Run ICA-AROMA #. Aggregate identified motion components (aggressive) to TSV #. Return classified_motion_ICs and melodic_mix for user to complete non-aggressive denoising in T1w space Additionally, non-aggressive denoising is performed on the BOLD series resampled into MNI space. .. workflow:: :graph2use: orig :simple_form: yes from fmriprep.workflows.bold.confounds import init_ica_aroma_wf wf = init_ica_aroma_wf() **Parameters** ignore_aroma_err : bool Do not fail on ICA-AROMA errors **Inputs** bold_mni BOLD series, resampled to template space movpar_file SPM-formatted motion parameters file bold_mask_mni BOLD series mask in template space **Outputs** aroma_confounds TSV of confounds identified as noise by ICA-AROMA aroma_noise_ics CSV of noise components identified by ICA-AROMA melodic_mix FSL MELODIC mixing matrix nonaggr_denoised_file BOLD series with non-aggressive ICA-AROMA denoising applied out_report Reportlet visualizing MELODIC ICs, with ICA-AROMA signal/noise labels .. _ICA-AROMA: ''' workflow = pe.Workflow(name=name) inputnode = pe.Node(niu.IdentityInterface( fields=['bold_mni', 'movpar_file', 'bold_mask_mni']), name='inputnode') outputnode = pe.Node(niu.IdentityInterface( fields=['aroma_confounds', 'out_report', 'aroma_noise_ics', 'melodic_mix', 'nonaggr_denoised_file']), name='outputnode') 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') def getusans_func(image, thresh): return [tuple([image, thresh])] getusans = pe.Node(niu.Function(function=getusans_func, output_names=['usans']), name='getusans', mem_gb=0.01) smooth = pe.Node(fsl.SUSAN(fwhm=6.0), name='smooth') # melodic node melodic = pe.Node(fsl.MELODIC(no_bet=True, no_mm=True), name="melodic") # ica_aroma node ica_aroma = pe.Node(ICA_AROMARPT(denoise_type='nonaggr', generate_report=True), name='ica_aroma') # 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') def _getbtthresh(medianval): return 0.75 * medianval # connect the nodes workflow.connect([ # Connect input nodes to complete smoothing (inputnode, calc_median_val, [('bold_mni', 'in_file'), ('bold_mask_mni', 'mask_file')]), (inputnode, calc_bold_mean, [('bold_mni', 'in_file')]), (calc_bold_mean, getusans, [('out_file', 'image')]), (calc_median_val, getusans, [('out_stat', 'thresh')]), (inputnode, smooth, [('bold_mni', 'in_file')]), (getusans, smooth, [('usans', 'usans')]), (calc_median_val, smooth, [(('out_stat', _getbtthresh), 'brightness_threshold')]), # connect smooth to melodic (smooth, melodic, [('smoothed_file', 'in_files')]), (inputnode, melodic, [('bold_mask_mni', 'mask')]), # connect nodes to ICA-AROMA (smooth, ica_aroma, [('smoothed_file', 'in_file')]), (inputnode, ica_aroma, [('bold_mask_mni', 'report_mask'), ('movpar_file', 'motion_parameters')]), (melodic, ica_aroma, [('out_dir', 'melodic_dir')]), # generate tsvs from ICA-AROMA (ica_aroma, ica_aroma_confound_extraction, [('out_dir', 'in_directory')]), # 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, [('out_report', 'out_report'), ('nonaggr_denoised_file', 'nonaggr_denoised_file')]), ]) return workflow
def _maskroi(in_mask, roi_file): import numpy as np import nibabel as nb from niworkflows.nipype.utils.filemanip import fname_presuffix roi = nb.load(roi_file) roidata = roi.get_data().astype(np.uint8) msk = nb.load(in_mask).get_data().astype(bool) roidata[~msk] = 0 roi.set_data_dtype(np.uint8) out = fname_presuffix(roi_file, suffix='_boldmsk') roi.__class__(roidata, roi.affine, roi.header).to_filename(out) return out