diff --git a/src/rbc/bids/functional.py b/src/rbc/bids/functional.py index 1991fd18..69b8a460 100644 --- a/src/rbc/bids/functional.py +++ b/src/rbc/bids/functional.py @@ -158,6 +158,12 @@ def export_functional( desc=bids_safe_label(reg), extension=".1D", ) + func.save( + outputs.bpf_regressor_file[reg], + suffix="regressors+filtered", + desc=bids_safe_label(reg), + extension=".1D", + ) mni = func.derive(space=TemplateSpace.MNI152NLIN6ASYM) for reg in regressors: diff --git a/src/rbc/bids/longitudinal.py b/src/rbc/bids/longitudinal.py index b0f52032..b46fa9be 100644 --- a/src/rbc/bids/longitudinal.py +++ b/src/rbc/bids/longitudinal.py @@ -2,11 +2,13 @@ from __future__ import annotations -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, TypedDict -from rbc.bids import Suffix, TemplateSpace +from rbc.bids import Suffix, TemplateSpace, bids_safe_label +from rbc.bids.metrics import export_metrics if TYPE_CHECKING: + from collections.abc import Sequence from pathlib import Path import polars as pl @@ -14,6 +16,7 @@ from rbc.bids import Bids from rbc.workflows.anatomical import AnatomicalLongOutputs from rbc.workflows.functional import FunctionalLongOutputs + from rbc.workflows.metrics import MetricsOutputs def _require_file(path: Path | None, field: str) -> Path: @@ -107,6 +110,18 @@ def export_longitudinal_anat(aex: Bids, outputs: AnatomicalLongOutputs) -> None: ) +class LongitudinalFuncInputs(TypedDict): + """Resolved inputs for the longitudinal functional workflow.""" + + template: Path + anat_to_template_xfm: Path + bold_to_anat_itk: Path + sbref: Path + bold: Path + bold_mask: Path + regressor_files: dict[str, Path] + + def resolve_longitudinal_func( func_q: Bids, tpl_q: Bids, @@ -114,7 +129,8 @@ def resolve_longitudinal_func( tpl_df: pl.DataFrame, *, ses: str, -) -> dict[str, Path | None]: + regressors: Sequence[str], +) -> LongitudinalFuncInputs: """Resolve inputs for longitudinal functional processing. Args: @@ -123,6 +139,7 @@ def resolve_longitudinal_func( func_df: DataFrame of functional derivatives. tpl_df: DataFrame of longitudinal template files. ses: Session label (used for template xfm lookup). + regressors: Regressor names (e.g. ``["36-parameter"]``). Returns: Dict with keys matching ``longitudinal_process`` parameters. @@ -146,9 +163,15 @@ def resolve_longitudinal_func( "bold": func_q.expect( func_df, suffix=Suffix.BOLD, desc="preproc", without=["space"] ), - "bold_mask": func_q.find( + "bold_mask": func_q.expect( func_df, suffix=Suffix.MASK, desc="brain", without=["space"] ), + "regressor_files": { + reg: func_q.expect( + func_df, suffix="regressors", desc=bids_safe_label(reg), extension=".1D" + ) + for reg in regressors + }, } @@ -167,5 +190,83 @@ def export_longitudinal_func(fex: Bids, outputs: FunctionalLongOutputs) -> None: desc="composite", extra={"from": "bold", "to": "longitudinal", "mode": "image"}, ) + fex.save(outputs.bold_mask, suffix=Suffix.MASK, desc="brain") + for reg, path in outputs.regressed_bold.items(): + fex.save( + path, + suffix=Suffix.BOLD, + desc="regressed", + extra={"reg": bids_safe_label(reg)}, + ) + for reg, path in outputs.cleaned_bold.items(): + fex.save( + path, + suffix=Suffix.BOLD, + desc="preproc", + extra={"reg": bids_safe_label(reg)}, + ) if outputs.bold_mask: fex.save(outputs.bold_mask, suffix=Suffix.MASK, desc="brain") + + +class LongitudinalMetricsInputs(TypedDict): + """Resolved inputs for longitudinal metrics.""" + + regressed_bold: Path + cleaned_bold: Path + template_brain_mask: Path + + +def resolve_longitudinal_metrics( + long_q: Bids, + func_df: pl.DataFrame, + *, + regressor: str, +) -> LongitudinalMetricsInputs: + """Resolve inputs for longitudinal metrics. + + Args: + long_q: Bids builder with ``space="longitudinal"`` and identity entities. + func_df: DataFrame of longitudinal functional derivatives. + regressor: Regressor name (e.g. ``"36-parameter"``). + + Returns: + Dict with keys matching ``single_session_metrics`` parameters. + """ + return { + "regressed_bold": long_q.expect( + func_df, + suffix=Suffix.BOLD, + desc="regressed", + extra={"reg": bids_safe_label(regressor)}, + ), + "cleaned_bold": long_q.expect( + func_df, + suffix=Suffix.BOLD, + desc="preproc", + extra={"reg": bids_safe_label(regressor)}, + ), + "template_brain_mask": long_q.expect( + func_df, + suffix=Suffix.MASK, + desc="brain", + ), + } + + +def export_longitudinal_metrics( + long_q: Bids, + outputs: MetricsOutputs, + *, + regressor: str, + atlases: list[str], +) -> None: + """Export longitudinal metrics outputs. + + Args: + long_q: Bids builder with ``space="longitudinal"`` and identity entities. + outputs: Results from the metrics workflow. + regressor: Regressor name. + atlases: Atlas labels. + """ + export_metrics(long_q, outputs, regressor=regressor, atlases=atlases) diff --git a/src/rbc/cli/longitudinal.py b/src/rbc/cli/longitudinal.py index f6630157..2f730d6d 100644 --- a/src/rbc/cli/longitudinal.py +++ b/src/rbc/cli/longitudinal.py @@ -3,12 +3,13 @@ from __future__ import annotations from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal from rbc.cli.base import BaseArgs, _or_default, _validate_nifti_path +from rbc.cli.metrics import _resolve_atlas_args from rbc.orchestration import Filters, RunnerConfig from rbc.orchestration.longitudinal import run -from rbc_resources import REGISTRATION_TEMPLATES +from rbc_resources import ATLAS_REGISTRY, REGISTRATION_TEMPLATES if TYPE_CHECKING: import argparse @@ -22,22 +23,29 @@ class LongitudinalArgs(BaseArgs): anatomical: bool functional: bool + metrics: bool registration_template: Path + regressor: Sequence[Literal["36-parameter", "aCompCor"]] + atlas_files: dict[str, Path] @classmethod def validate_namespace(cls, ns: argparse.Namespace) -> LongitudinalArgs: """Validation of longitudinal workflow specific arguments to NamedTuple.""" - if not ns.functional and not ns.anatomical: + if not ns.functional and not ns.anatomical and not ns.metrics: raise ValueError( - "At least one of '--anatomical' or '--functional' is required." + "At least one of '--anatomical', '--functional', " + "or '--metrics' is required." ) return cls( **BaseArgs.validate_namespace(ns).__dict__, anatomical=ns.anatomical, functional=ns.functional, + metrics=ns.metrics, registration_template=_or_default( ns.anat_template, REGISTRATION_TEMPLATES.brain_1mm ), + regressor=ns.regressor, + atlas_files=_resolve_atlas_args(ns.atlas), ) @@ -52,7 +60,10 @@ def main(args: LongitudinalArgs) -> int: ), anatomical=args.anatomical, functional=args.functional, + metrics=args.metrics, registration_template=args.registration_template, + regressors=args.regressor, + atlas_files=args.atlas_files, runner_config=RunnerConfig( runner=args.runner, verbose=bool(args.verbose), @@ -86,7 +97,30 @@ def register_command( action="store_true", help="Use functional longitudinal pipeline for processing", ) - + parser.add_argument( + "--regressor", + nargs="+", + choices=["36-parameter", "aCompCor"], + default=["36-parameter"], + help="Space-delimited nuisance regression method(s) to apply.", + ) + parser.add_argument( + "--metrics", + default=False, + action="store_true", + help="Compute longitudinal metrics (ALFF, ReHo, timeseries).", + ) + parser.add_argument( + "--atlas", + nargs="+", + default=["schaefer_200"], + metavar="ATLAS", + help=( + "Atlas(es) for timeseries extraction. Accepts registry names " + f"({', '.join(sorted(ATLAS_REGISTRY))}) or paths to custom NIfTI " + "atlas files." + ), + ) templates = parser.add_argument_group("template overrides") templates.add_argument( "--anat-template", diff --git a/src/rbc/orchestration/longitudinal.py b/src/rbc/orchestration/longitudinal.py index ffadf2b7..84c17eac 100644 --- a/src/rbc/orchestration/longitudinal.py +++ b/src/rbc/orchestration/longitudinal.py @@ -17,20 +17,25 @@ load_table, ) from rbc.bids.longitudinal import ( + LongitudinalFuncInputs, export_longitudinal_anat, export_longitudinal_func, + export_longitudinal_metrics, resolve_longitudinal_anat, resolve_longitudinal_func, + resolve_longitudinal_metrics, ) from rbc.bids.session import iter_session_files, load_session from rbc.context import RunContext from rbc.orchestration import Filters, RunnerConfig, init_runner +from rbc.orchestration.metrics import _read_header_tr from rbc.workflows.anatomical import longitudinal_process as anatomical_longitudinal from rbc.workflows.functional import longitudinal_process as functional_longitudinal +from rbc.workflows.metrics import single_session_metrics as metrics_longitudinal from rbc_resources import REGISTRATION_TEMPLATES if TYPE_CHECKING: - from collections.abc import Sequence + from collections.abc import Mapping, Sequence from pathlib import Path _logger = logging.getLogger(__name__) @@ -75,6 +80,7 @@ def process_func( pipe_ctx: RunContext, func_df: pl.DataFrame, tpl_df: pl.DataFrame, + regressors: Sequence[str], ) -> None: """Handle functional longitudinal processing for one BOLD run. @@ -82,6 +88,7 @@ def process_func( pipe_ctx: RunContext bound to this subject/session. func_df: Functional derivative DataFrame for this run. tpl_df: Longitudinal template DataFrame. + regressors: Regressor names (e.g. ``["36-parameter"]``). """ row = func_df.filter(suffix=Suffix.BOLD).row(0, named=True) ents = extract_entities(row, ["task", "run"]) @@ -89,18 +96,55 @@ def process_func( func_q = pipe_ctx.bids(datatype=Datatype.FUNC, entities=ents) tpl_q = pipe_ctx.bids(datatype=Datatype.ANAT).derive(ses="longitudinal") - resolved = resolve_longitudinal_func( + resolved: LongitudinalFuncInputs = resolve_longitudinal_func( func_q, tpl_q, func_df, tpl_df, ses=pipe_ctx.ses, # type: ignore[arg-type] + regressors=regressors, ) func_outputs = functional_longitudinal(**resolved) # type: ignore[arg-type] fex = func_q.derive(space="longitudinal") export_longitudinal_func(fex, func_outputs) +def process_metrics( + pipe_ctx: RunContext, + func_df: pl.DataFrame, + *, + regressors: Sequence[str], + atlas_files: Mapping[str, Path], + fwhm: float = 6.0, +) -> None: + """Handle longitudinal metrics for one BOLD run. + + Args: + pipe_ctx: RunContext bound to this subject/session. + func_df: Functional derivative DataFrame for this run. + regressors: Regressor names. + atlas_files: Mapping of atlas labels to NIfTI file paths. + fwhm: Smoothing kernel FWHM in mm. + """ + row = func_df.filter(suffix=Suffix.BOLD).row(0, named=True) + ents = extract_entities(row, ["task", "run"]) + + long_q = pipe_ctx.bids(datatype=Datatype.FUNC, entities=ents, space="longitudinal") + + for regressor in regressors: + resolved = resolve_longitudinal_metrics(long_q, func_df, regressor=regressor) + run_tr = _read_header_tr(resolved["regressed_bold"]) + outputs = metrics_longitudinal( + **resolved, + tr=run_tr, + atlas_files=atlas_files, + fwhm=fwhm, + ) + export_longitudinal_metrics( + long_q, outputs, regressor=regressor, atlases=list(atlas_files) + ) + + def run( input_dirs: Sequence[Path], output_dir: Path, @@ -108,7 +152,10 @@ def run( filters: Filters, anatomical: bool = True, functional: bool = True, + metrics: bool = False, registration_template: Path = REGISTRATION_TEMPLATES.brain_1mm, + regressors: Sequence[str] = ("36-parameter",), + atlas_files: Mapping[str, Path] | None = None, runner_config: RunnerConfig | None = None, ) -> None: """Run the longitudinal pipeline for all matching subjects/sessions. @@ -119,7 +166,10 @@ def run( filters: Participant/session/task filters. anatomical: Run anatomical longitudinal processing. functional: Run functional longitudinal processing. + metrics: Run longitudinal metrics (ALFF, ReHo, timeseries). + atlas_files: Mapping of atlas labels to NIfTI file paths. registration_template: Brain template for ANTs registration. + regressors: Nuisance regressor strategies to apply (e.g. ``["36-parameter"]``). runner_config: Execution backend configuration. """ config = runner_config or RunnerConfig() @@ -173,7 +223,21 @@ def run( if functional: for func_df, _ in iter_session_files(session, groupby=FUNC_GROUP_ENTITIES): - process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) + process_func( + pipe_ctx=pipe_ctx, + func_df=func_df, + tpl_df=tpl_df, + regressors=regressors, + ) + + if metrics: + for func_df, _ in iter_session_files(session, groupby=FUNC_GROUP_ENTITIES): + process_metrics( + pipe_ctx=pipe_ctx, + func_df=func_df, + regressors=regressors, + atlas_files=atlas_files or {}, + ) pipe_ctx.ensure_dataset_description() diff --git a/src/rbc/workflows/functional.py b/src/rbc/workflows/functional.py index ce854aef..c1ca61a2 100644 --- a/src/rbc/workflows/functional.py +++ b/src/rbc/workflows/functional.py @@ -81,7 +81,8 @@ class FunctionalOutputs(NamedTuple): template_bold: BOLD resampled to template space. regressed_bold: Nuisance-regressed & non-bandpassed BOLD. cleaned_bold: Nuisance-regressed & bandpass-filtered BOLD. - regressor_file: Bandpass-filtered nuisance regressor ``.1D`` file. + regressor_file: Nuisance regressor ``.1D`` files. + bpf_regressor_file: Bandpass-filtered nuisance regressor ``.1D`` files. template_brain_mask: Brain mask warped to template space. """ @@ -105,6 +106,7 @@ class FunctionalOutputs(NamedTuple): regressed_bold: dict[str, Path] cleaned_bold: dict[str, Path] regressor_file: dict[str, Path] + bpf_regressor_file: dict[str, Path] template_brain_mask: Path @@ -378,7 +380,8 @@ def single_session_preprocess( template_bold=template_bold, regressed_bold={r: regression[r].regressed_bold for r in regressor_set}, cleaned_bold={r: cleaned[r].regressed_bold for r in regressor_set}, - regressor_file=filtered_regressors, + regressor_file={r: regressors[r].regressor_file for r in regressor_set}, + bpf_regressor_file=filtered_regressors, template_brain_mask=tmpl_brain, ) @@ -390,14 +393,19 @@ class FunctionalLongOutputs(NamedTuple): bold_to_long_xfm: BOLD-to-longitudinal-template composite warp. sbref: Motion reference volume warped to longitudinal template space. bold: Preprocessed BOLD warped to longitudinal template space. - bold_mask: Brain mask warped to longitudinal template space, - or *None* if no mask was provided. + bold_mask: Brain mask warped to longitudinal template space. + regressed_bold: Nuisance-regressed (non-bandpassed) BOLD in longitudinal + template space. Suitable for ALFF/fALFF. + cleaned_bold: Nuisance-regressed + bandpass-filtered BOLD in longitudinal + template space (Hallquist 2013). """ bold_to_long_xfm: Path sbref: Path bold: Path - bold_mask: Path | None = None + bold_mask: Path + regressed_bold: dict[str, Path] + cleaned_bold: dict[str, Path] def longitudinal_process( @@ -407,7 +415,8 @@ def longitudinal_process( bold_to_anat_itk: Path, sbref: Path, bold: Path, - bold_mask: Path | None, + bold_mask: Path, + regressor_files: dict[str, Path], ) -> FunctionalLongOutputs: """Transform preprocessed functional outputs to longitudinal template space. @@ -415,33 +424,74 @@ def longitudinal_process( composite warp is available, and anatomical data has already been processed to longitudinal template space. + Regressors are computed once during cross-sectional preprocessing and passed + in via ``regressor_files``. Only the regression steps are re-run against the + longitudinal space BOLD. + + Steps: + 1. Compose BOLD-to-anatomical + anatomical-to-longitudinal-template transforms. + 2. Warp sbref (3D) and preproc BOLD (4D) to longitudinal template space. + 3. Warp brain mask to longitudinal template space. + 4. Nuisance regression without bandpass on longitudinal-space BOLD + (for ALFF/fALFF). + 5. Nuisance regression with simultaneous bandpass filtering on longitudinal-space + BOLD (Hallquist 2013). + Args: template: Longitudinal template image. anat_to_template_xfm: T1w-to-longitudinal-template composite warp. bold_to_anat_itk: BOLD-to-T1w affine in ITK format. sbref: Motion reference (single-band reference) volume. bold: Preprocessed bold image. - bold_mask: Bold brain mask, if available. + bold_mask: Bold brain mask in native space. + regressor_files: Per-regressor nuisance regressor ``.1D`` files. + Returns: :class:`FunctionalLongOutputs` with all non-null inputs transformed to template space. """ + # 1. Compose full BOLD -> longitudinal template transform bold_to_tpl_xfm = compose_transform( ref=template, bold_to_anat_itk=bold_to_anat_itk, anat_to_tpl_xfm=anat_to_template_xfm, ) + # 2. Warp sbref & bold to longitudinal space + warped_sbref = func_transform( + in_file=sbref, template=template, xfm=bold_to_tpl_xfm, strategy="single" + ) + warped_bold = func_transform( + in_file=bold, template=template, xfm=bold_to_tpl_xfm, strategy="chunked" + ) + + # 3. Warp bold mask to longitudinal space + warped_mask = mask_transform(mask=bold_mask, template=template, xfm=bold_to_tpl_xfm) + + regression: dict[str, ApplyRegressionOutputs] = {} + cleaned: dict[str, ApplyRegressionOutputs] = {} + for reg, reg_file in regressor_files.items(): + # 4. Nuisance regression without bandpass + _logger.info("%s nuisance regression (no bandpass)", reg) + regression[reg] = apply_regression( + bold_file=warped_bold, + brain_mask_file=warped_mask, + regressor_file=reg_file, + ) + # 5. Simultaneous regression + bandpass filtering (Hallquist 2013) + _logger.info("%s nuisance regression + bandpass filtering", reg) + cleaned[reg] = apply_regression_bandpass( + bold_file=warped_bold, + brain_mask_file=warped_mask, + regressor_file=reg_file, + ) + return FunctionalLongOutputs( - sbref=func_transform( # 3D volume - in_file=sbref, template=template, xfm=bold_to_tpl_xfm, strategy="single" - ), - bold=func_transform( - in_file=bold, template=template, xfm=bold_to_tpl_xfm, strategy="chunked" - ), - bold_mask=mask_transform(mask=bold_mask, template=template, xfm=bold_to_tpl_xfm) - if bold_mask - else None, bold_to_long_xfm=bold_to_tpl_xfm, + sbref=warped_sbref, + bold=warped_bold, + bold_mask=warped_mask, + regressed_bold={r: regression[r].regressed_bold for r in regressor_files}, + cleaned_bold={r: cleaned[r].regressed_bold for r in regressor_files}, ) diff --git a/tests/unit/bids/test_exports.py b/tests/unit/bids/test_exports.py index 1c96faab..53b65d1a 100644 --- a/tests/unit/bids/test_exports.py +++ b/tests/unit/bids/test_exports.py @@ -74,6 +74,7 @@ def _make_func_outputs(w: Path, regressors: list[str]) -> FunctionalOutputs: regressed_bold={r: _dummy(w, f"regressed_{r}.nii.gz") for r in regressors}, cleaned_bold={r: _dummy(w, f"cleaned_{r}.nii.gz") for r in regressors}, regressor_file={r: _dummy(w, f"regressors_{r}.1D") for r in regressors}, + bpf_regressor_file={r: _dummy(w, f"bpf_regressors_{r}.1D") for r in regressors}, template_brain_mask=_dummy(w, "template_mask.nii.gz"), ) @@ -201,7 +202,7 @@ def test_file_count_single_regressor( outputs = _make_func_outputs(workdir, ["36-parameter"]) export_functional(func_bids, outputs, regressors=["36-parameter"]) saved = list(pipe_ctx.output_dir.rglob("*.*")) - assert len(saved) == 13 + assert len(saved) == 14 def test_file_count_two_regressors( self, func_bids: Bids, workdir: Path, pipe_ctx: RunContext @@ -214,7 +215,7 @@ def test_file_count_two_regressors( outputs = _make_func_outputs(workdir, regs) export_functional(func_bids, outputs, regressors=regs) saved = list(pipe_ctx.output_dir.rglob("*.*")) - assert len(saved) == 16 + assert len(saved) == 18 # --------------------------------------------------------------------------- diff --git a/tests/unit/cli/test_longitudinal.py b/tests/unit/cli/test_longitudinal.py index da966521..00c26c36 100644 --- a/tests/unit/cli/test_longitudinal.py +++ b/tests/unit/cli/test_longitudinal.py @@ -27,9 +27,12 @@ def base_args(tmp_path: Path) -> argparse.Namespace: session_label=[], anatomical=True, functional=False, + metrics=False, tmp_dir=None, anat_template=None, ants_threads=1, + atlas=["schaefer_200"], + regressor=["36-parameter"], ) diff --git a/tests/unit/orchestration/test_longitudinal.py b/tests/unit/orchestration/test_longitudinal.py index ca28e510..0e436e9a 100644 --- a/tests/unit/orchestration/test_longitudinal.py +++ b/tests/unit/orchestration/test_longitudinal.py @@ -68,6 +68,8 @@ def _mock_func_outputs(*, with_bold_mask: bool = True) -> Mock: m.bold = fake / "bold.nii.gz" m.bold_to_long_xfm = fake / "bold_to_long_xfm.nii.gz" m.bold_mask = (fake / "bold_mask.nii.gz") if with_bold_mask else None + m.regressed_bold = {"36-parameter": fake / "regressed.nii.gz"} + m.cleaned_bold = {"36-parameter": fake / "cleaned.nii.gz"} return m @@ -319,7 +321,12 @@ def test_calls_functional_longitudinal( ), patch("rbc.bids.builder.shutil.copy2"), ): - process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) + process_func( + pipe_ctx=pipe_ctx, + func_df=func_df, + tpl_df=tpl_df, + regressors=["36-parameter"], + ) assert mock_func.call_count == 1 @pytest.mark.parametrize( @@ -359,26 +366,12 @@ def test_missing_required_file_raises( ), pytest.raises(FileNotFoundError), ): - process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) - - def test_optional_bold_mask_file_not_found( - self, func_df: pl.DataFrame, tpl_df: pl.DataFrame, tmp_path: Path - ) -> None: - """Optional bold_mask not found is caught; 3 exports emitted.""" - pipe_ctx = RunContext(sub="01", ses="baseline", output_dir=tmp_path) - with ( - patch( - "rbc.orchestration.longitudinal.functional_longitudinal", - return_value=_mock_func_outputs(with_bold_mask=False), - ), - patch( - "rbc.bids.query.find_file", - side_effect=_none_for(suffix="mask", desc="brain"), - ), - patch("rbc.bids.builder.shutil.copy2") as mock_copy, - ): - process_func(pipe_ctx=pipe_ctx, func_df=func_df, tpl_df=tpl_df) - assert mock_copy.call_count == 3 + process_func( + pipe_ctx=pipe_ctx, + func_df=func_df, + tpl_df=tpl_df, + regressors=["36-parameter"], + ) class TestLongitudinalDispatch: