diff --git a/src/fmripost_aroma/utils/bids.py b/src/fmripost_aroma/utils/bids.py index f0f9560..844fb7e 100644 --- a/src/fmripost_aroma/utils/bids.py +++ b/src/fmripost_aroma/utils/bids.py @@ -7,10 +7,44 @@ from pathlib import Path from bids.layout import BIDSLayout +from bids.utils import listify from fmripost_aroma.data import load as load_data +def extract_entities(file_list): + """Return a dictionary of common entities given a list of files. + + Examples + -------- + >>> extract_entities("sub-01/anat/sub-01_T1w.nii.gz") + {'subject': '01', 'suffix': 'T1w', 'datatype': 'anat', 'extension': '.nii.gz'} + >>> extract_entities(["sub-01/anat/sub-01_T1w.nii.gz"] * 2) + {'subject': '01', 'suffix': 'T1w', 'datatype': 'anat', 'extension': '.nii.gz'} + >>> extract_entities(["sub-01/anat/sub-01_run-1_T1w.nii.gz", + ... "sub-01/anat/sub-01_run-2_T1w.nii.gz"]) + {'subject': '01', 'run': [1, 2], 'suffix': 'T1w', 'datatype': 'anat', 'extension': '.nii.gz'} + + """ + from collections import defaultdict + + from bids.layout import parse_file_entities + + entities = defaultdict(list) + for e, v in [ + ev_pair for f in listify(file_list) for ev_pair in parse_file_entities(f).items() + ]: + entities[e].append(v) + + def _unique(inlist): + inlist = sorted(set(inlist)) + if len(inlist) == 1: + return inlist[0] + return inlist + + return {k: _unique(v) for k, v in entities.items()} + + def collect_derivatives( raw_dir: Path | None, derivatives_dir: Path, diff --git a/src/fmripost_aroma/utils/utils.py b/src/fmripost_aroma/utils/utils.py index cf3ed27..4836a59 100644 --- a/src/fmripost_aroma/utils/utils.py +++ b/src/fmripost_aroma/utils/utils.py @@ -422,3 +422,22 @@ def get_spectrum(data: np.array, tr: float): freqs = np.fft.rfftfreq((power_spectrum.shape[0] * 2) - 1, tr) idx = np.argsort(freqs) return power_spectrum[idx, :], freqs[idx] + + +def _get_wf_name(bold_fname, prefix): + """Derive the workflow name for supplied BOLD file. + + >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_bold.nii.gz", "aroma") + 'aroma_task_nback_wf' + >>> _get_wf_name( + ... "/completely/made/up/path/sub-01_task-nback_run-01_echo-1_bold.nii.gz", + ... "preproc", + ... ) + 'preproc_task_nback_run_01_echo_1_wf' + + """ + from nipype.utils.filemanip import split_filename + + fname = split_filename(bold_fname)[1] + fname_nosub = '_'.join(fname.split('_')[1:-1]) + return f"{prefix}_{fname_nosub.replace('-', '_')}_wf" diff --git a/src/fmripost_aroma/workflows/aroma.py b/src/fmripost_aroma/workflows/aroma.py index 358eeba..8c623bc 100644 --- a/src/fmripost_aroma/workflows/aroma.py +++ b/src/fmripost_aroma/workflows/aroma.py @@ -28,6 +28,7 @@ from fmripost_aroma import config from fmripost_aroma.interfaces.aroma import AROMAClassifier from fmripost_aroma.interfaces.bids import DerivativesDataSink +from fmripost_aroma.utils.utils import _get_wf_name def init_ica_aroma_wf( @@ -47,11 +48,7 @@ def init_ica_aroma_wf( #. Smooth data using FSL `susan`, with a kernel width FWHM=6.0mm. #. Run FSL `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 - #. Calculate ICA-AROMA-identified noise components - (columns named ``AROMAAggrCompXX``) + #. Aggregate components and classifications to TSVs There is a current discussion on whether other confounds should be extracted before or after denoising `here @@ -376,9 +373,10 @@ def init_denoise_wf(bold_file): niu.IdentityInterface( fields=[ 'bold_file', - 'bold_mask_std', + 'bold_mask', 'confounds', 'skip_vols', + 'spatial_reference', ], ), name='inputnode', @@ -403,7 +401,7 @@ def init_denoise_wf(bold_file): workflow.connect([ (inputnode, denoise, [ ('confounds', 'confounds'), - ('bold_mask_std', 'mask_file'), + ('bold_mask', 'mask_file'), ]), (rm_non_steady_state, denoise, [('bold_cut', 'bold_file')]), ]) # fmt:skip @@ -429,7 +427,11 @@ def init_denoise_wf(bold_file): run_without_submitting=True, mem_gb=config.DEFAULT_MEMORY_MIN_GB, ) - workflow.connect([(add_non_steady_state, ds_denoised, [('bold_add', 'denoised_file')])]) + workflow.connect([ + # spatial_reference needs to be parsed into space, cohort, res, den, etc. + (inputnode, ds_denoised, [('spatial_reference', 'space')]), + (add_non_steady_state, ds_denoised, [('bold_add', 'denoised_file')]), + ]) # fmt:skip return workflow @@ -479,26 +481,6 @@ def _add_volumes(bold_file, bold_cut_file, skip_vols): return out -def _get_wf_name(bold_fname, prefix): - """ - Derive the workflow name for supplied BOLD file. - - >>> _get_wf_name("/completely/made/up/path/sub-01_task-nback_bold.nii.gz", "aroma") - 'aroma_task_nback_wf' - >>> _get_wf_name( - ... "/completely/made/up/path/sub-01_task-nback_run-01_echo-1_bold.nii.gz", - ... "preproc", - ... ) - 'preproc_task_nback_run_01_echo_1_wf' - - """ - from nipype.utils.filemanip import split_filename - - fname = split_filename(bold_fname)[1] - fname_nosub = '_'.join(fname.split('_')[1:-1]) - return f"{prefix}_{fname_nosub.replace('-', '_')}_wf" - - def _select_melodic_files(melodic_dir): """Select the mixing and component maps from the Melodic output.""" import os diff --git a/src/fmripost_aroma/workflows/base.py b/src/fmripost_aroma/workflows/base.py index 7b10d13..0e52b95 100644 --- a/src/fmripost_aroma/workflows/base.py +++ b/src/fmripost_aroma/workflows/base.py @@ -29,7 +29,6 @@ """ -import os import sys import warnings from copy import deepcopy @@ -41,17 +40,15 @@ from fmripost_aroma import config from fmripost_aroma.interfaces.bids import DerivativesDataSink from fmripost_aroma.interfaces.reportlets import AboutSummary, SubjectSummary +from fmripost_aroma.utils.utils import _get_wf_name from fmripost_aroma.workflows.resampling import init_resample_volumetric_wf def init_fmripost_aroma_wf(): """Build *fMRIPost-AROMA*'s pipeline. - This workflow organizes the execution of FMRIPREP, with a sub-workflow for - each subject. - - If FreeSurfer's ``recon-all`` is to be run, a corresponding folder is created - and populated with any needed template subjects under the derivatives folder. + This workflow organizes the execution of fMRIPost-AROMA, + with a sub-workflow for each subject. Workflow Graph .. workflow:: @@ -60,33 +57,18 @@ def init_fmripost_aroma_wf(): from fmripost_aroma.workflows.tests import mock_config from fmripost_aroma.workflows.base import init_fmripost_aroma_wf + with mock_config(): wf = init_fmripost_aroma_wf() """ from niworkflows.engine.workflows import LiterateWorkflow as Workflow - from niworkflows.interfaces.bids import BIDSFreeSurferDir ver = Version(config.environment.version) fmripost_aroma_wf = Workflow(name=f'fmripost_aroma_{ver.major}_{ver.minor}_wf') fmripost_aroma_wf.base_dir = config.execution.work_dir - freesurfer = config.workflow.run_reconall - if freesurfer: - fsdir = pe.Node( - BIDSFreeSurferDir( - derivatives=config.execution.output_dir, - freesurfer_home=os.getenv('FREESURFER_HOME'), - spaces=config.workflow.spaces.get_fs_spaces(), - minimum_fs_version='7.0.0', - ), - name=f"fsdir_run_{config.execution.run_uuid.replace('-', '_')}", - run_without_submitting=True, - ) - if config.execution.fs_subjects_dir is not None: - fsdir.inputs.subjects_dir = str(config.execution.fs_subjects_dir.absolute()) - for subject_id in config.execution.participant_label: single_subject_wf = init_single_subject_wf(subject_id) @@ -99,15 +81,7 @@ def init_fmripost_aroma_wf(): for node in single_subject_wf._get_all_nodes(): node.config = deepcopy(single_subject_wf.config) - if freesurfer: - fmripost_aroma_wf.connect( - fsdir, - 'subjects_dir', - single_subject_wf, - 'inputnode.subjects_dir', - ) - else: - fmripost_aroma_wf.add_nodes([single_subject_wf]) + fmripost_aroma_wf.add_nodes([single_subject_wf]) # Dump a copy of the config file into the log directory log_dir = ( @@ -125,12 +99,8 @@ def init_fmripost_aroma_wf(): def init_single_subject_wf(subject_id: str): """Organize the postprocessing pipeline for a single subject. - It collects and reports information about the subject, and prepares - sub-workflows to perform anatomical and functional preprocessing. - Anatomical preprocessing is performed in a single workflow, regardless of - the number of sessions. - Functional preprocessing is performed using a separate workflow for each - individual BOLD series. + It collects and reports information about the subject, + and prepares sub-workflows to postprocess each BOLD series. Workflow Graph .. workflow:: @@ -182,16 +152,15 @@ def init_single_subject_wf(subject_id: str): from niworkflows.utils.misc import fix_multi_T1w_source_name from niworkflows.utils.spaces import Reference - from fmripost_aroma.utils.bids import collect_derivatives + from fmripost_aroma.utils.bids import collect_derivatives, extract_entities from fmripost_aroma.workflows.aroma import init_denoise_wf, init_ica_aroma_wf spaces = config.workflow.spaces workflow = Workflow(name=f'sub_{subject_id}_wf') workflow.__desc__ = f""" -Results included in this manuscript come from preprocessing -performed using *fMRIPost-AROMA* {config.environment.version} -(@fmriprep1; @fmriprep2; RRID:SCR_016216), +Results included in this manuscript come from postprocessing +performed using *fMRIPost-AROMA* {config.environment.version} (@ica_aroma), which is based on *Nipype* {config.environment.nipype_version} (@nipype1; @nipype2; RRID:SCR_002502). @@ -199,8 +168,7 @@ def init_single_subject_wf(subject_id: str): workflow.__postdesc__ = f""" Many internal operations of *fMRIPost-AROMA* use -*Nilearn* {NILEARN_VERSION} [@nilearn, RRID:SCR_001362], -mostly within the functional processing workflow. +*Nilearn* {NILEARN_VERSION} [@nilearn, RRID:SCR_001362]. For more details of the pipeline, see [the section corresponding to workflows in *fMRIPost-AROMA*'s documentation]\ (https://fmripost_aroma.readthedocs.io/en/latest/workflows.html \ @@ -224,9 +192,8 @@ def init_single_subject_wf(subject_id: str): entities=config.execution.bids_filters, ) - anat_only = config.workflow.anat_only # Make sure we always go through these two checks - if not anat_only and not subject_data['bold']: + if not subject_data['bold']: task_id = config.execution.task_id raise RuntimeError( f"No BOLD images found for participant {subject_id} and " @@ -248,7 +215,6 @@ def init_single_subject_wf(subject_id: str): bidssrc = pe.Node( BIDSDataGrabber( subject_data=subject_data, - anat_only=config.workflow.anat_only, subject_id=subject_id, ), name='bidssrc', @@ -321,19 +287,17 @@ def init_single_subject_wf(subject_id: str): ica_aroma_wf = init_ica_aroma_wf(bold_file=bold_file) ica_aroma_wf.__desc__ = func_pre_desc + (ica_aroma_wf.__desc__ or '') + entities = extract_entities(bold_file) + functional_cache = {} if config.execution.derivatives: # Collect native-space derivatives and transforms - from fmripost_aroma.utils.bids import collect_derivatives, extract_entities - - entities = extract_entities(bold_file) - for deriv_dir in config.execution.derivatives: functional_cache.update( collect_derivatives( derivatives_dir=deriv_dir, entities=entities, - ) + ), ) # Resample to MNI152NLin6Asym:res-2, for ICA-AROMA classification @@ -341,7 +305,9 @@ def init_single_subject_wf(subject_id: str): bold_file=bold_file, precomputed=functional_cache, space=Reference.from_string("MNI152NLin6Asym:res-2")[0], + name=_get_wf_name(bold_file, 'resample_raw'), ) + resample_raw_wf.inputs.inputnode.bold_file = bold_file workflow.connect([ (resample_raw_wf, ica_aroma_wf, [ ('outputnode.bold_std', 'inputnode.bold_std'), @@ -351,13 +317,11 @@ def init_single_subject_wf(subject_id: str): else: # Collect MNI152NLin6Asym:res-2 derivatives # Only derivatives dataset was passed in, so we expected standard-space derivatives - from fmripost_aroma.utils.bids import collect_derivatives - functional_cache.update( collect_derivatives( - derivatives_dir=deriv_dir, + derivatives_dir=config.execution.layout, entities=entities, - ) + ), ) ica_aroma_wf.inputs.inputnode.bold_std = functional_cache['bold_std'] ica_aroma_wf.inputs.inputnode.bold_mask_std = functional_cache['bold_mask_std'] @@ -367,53 +331,29 @@ def init_single_subject_wf(subject_id: str): ica_aroma_wf.inputs.inputnode.skip_vols = functional_cache['skip_vols'] ica_aroma_wf.inputs.inputnode.spatial_reference = functional_cache['spatial_reference'] - # Now denoise the native-space BOLD data using ICA-AROMA - denoise_native_wf = init_denoise_wf(bold_file=bold_file) - - # Resample the BOLD series to MNI152NLin6Asym-2mm - - # Run ICA-AROMA using MNI152NLin6Asym-2mm BOLD data - ica_aroma_wf = init_ica_aroma_wf( - bold_file=bold_file, - precomputed=functional_cache, - ) - ica_aroma_wf.__desc__ = func_pre_desc + (ica_aroma_wf.__desc__ or '') - - workflow.connect([ - (ica_aroma_wf, denoise_native_wf, [ - ('outputnode.aroma_noise_ics', 'inputnode.aroma_noise_ics'), - ]), - ]) # fmt:skip - - for space in spaces: - resample_to_space_wf = init_resample_volumetric_wf( - bold_file=bold_file, - functional_cache=functional_cache, - space=space, - ) - workflow.connect([ - (denoise_native_wf, resample_to_space_wf, [ - ('outputnode.denoised_file', 'inputnode.bold_file'), - ]) - ]) - if config.workflow.denoise_method: - # Warp the BOLD series to requested output spaces - # XXX: Probably should just grab the MNI152NLin6Asym-2mm file if that - # space+resolution is requested. + for space in spaces: + # Warp each BOLD run to requested output spaces + resample_to_space_wf = init_resample_volumetric_wf( + bold_file=bold_file, + functional_cache=functional_cache, + space=space, + name=_get_wf_name(bold_file, f'resample_{space}'), + ) - # Run the denoising workflow on each requested BOLD series - denoise_wf = init_denoise_wf(bold_file=bold_file) - workflow.connect([ - (inputnode, denoise_wf, [ - ('bold_std', 'inputnode.bold_std'), - ('bold_mask_std', 'inputnode.bold_mask_std'), - ('spatial_reference', 'inputnode.spatial_reference'), - ]), - (ica_aroma_wf, denoise_wf, [ - ('outputnode.aroma_confounds', 'inputnode.confounds'), - ]), - ]) # fmt:skip + # Now denoise the output-space BOLD data using ICA-AROMA + denoise_wf = init_denoise_wf(bold_file=bold_file) + denoise_wf.inputs.inputnode.skip_vols = functional_cache['skip_vols'] + workflow.connect([ + (resample_to_space_wf, denoise_wf, [ + ('bold_std', 'inputnode.bold_file'), + ('bold_mask_std', 'inputnode.bold_mask'), + ('spatial_reference', 'inputnode.spatial_reference'), + ]), + (ica_aroma_wf, denoise_wf, [ + ('outputnode.aroma_confounds', 'inputnode.confounds'), + ]), + ]) # fmt:skip return clean_datasinks(workflow) diff --git a/src/fmripost_aroma/workflows/resampling.py b/src/fmripost_aroma/workflows/resampling.py index 1c8627f..25e5997 100644 --- a/src/fmripost_aroma/workflows/resampling.py +++ b/src/fmripost_aroma/workflows/resampling.py @@ -4,7 +4,13 @@ from nipype.pipeline import engine as pe -def init_resample_volumetric_wf(bold_file, functional_cache, space, run_stc): +def init_resample_volumetric_wf( + bold_file, + functional_cache, + space, + run_stc, + name='resample_volumetric_wf', +): """Resample raw BOLD data to requested volumetric space space. Parameters @@ -17,13 +23,15 @@ def init_resample_volumetric_wf(bold_file, functional_cache, space, run_stc): Spatial reference. run_stc : bool Whether to run STC. + name : str + Workflow name. """ from fmriprep.workflows.bold.stc import init_bold_stc_wf from niworkflows.engine.workflows import LiterateWorkflow as Workflow from fmripost_aroma.interfaces.resampler import Resampler - workflow = Workflow(name='resample_raw_wf') + workflow = Workflow(name=name) inputnode = pe.Node( niu.IdentityInterface(fields=['bold_file', 'mask_file']),