Remove sdcflows and model ASLPrep on fMRIPrep/next

[tsalo]

Closes #219, closes #297, closes #346, closes #340, closes #341, and closes #199.

Changes proposed in this pull request

• Model ASLPrep on fMRIPrep 23.2.0a2.
• Outsource everything I can to fMRIPrep.
• New parameters: --bids-database-dir, --level (which is mostly untested), --medial-surface-nan, --project-goodvoxels, --cifti-output, --no-msm, --no-submm-recon, --fs-subjects-dir, --fs-no-reconall, --config-file.

Breaking changes

• Users must point to the actual output directory. aslprep will no longer be appended automatically.
• The parameter --anat-derivatives is replaced with --derivatives, which takes one or more paths to anatomical or ASLPrep derivatives.
• The parameter --dummy-vols is now --dummy-scans to match fMRIPrep.
• Native-space aslref image split into desc-hmc and desc-coreg versions.
• Add coverage files for mean CBF and basil CBF parcellations.
• desc-pvGM_cbf renamed to desc-basilGM_cbf
• desc-pvWM_cbf renamed to desc-basilWM_cbf
• desc-confounds_regressors.tsv renamed to desc-confounds_timeseries.tsv
• from-T1w_to-scanner_mode-image_xfm.txt removed. The inverse is invertible, so it's fine.
• from-scanner_to-T1w_mode-image_xfm.txt split into from-orig_to-aslref_mode-image_xfm.txt and from-aslref_to-T1w_mode-image_xfm.txt

[tsalo]

I'm getting the following error now that it's installing a version of niworkflows from GitHub:

Collecting niworkflows@ git+https://github.com/nipreps/niworkflows.git@master
  Cloning https://github.com/nipreps/niworkflows.git (to revision master) to /tmp/pip-install-enby527a/niworkflows_56e388d470f74c29b63b04d7b4f3bdab
  Running command git clone --filter=blob:none --quiet https://github.com/nipreps/niworkflows.git /tmp/pip-install-enby527a/niworkflows_56e388d470f74c29b63b04d7b4f3bdab
The authenticity of host 'github.com (140.82.112.3)' can't be established.

My solution: Use the most recent releases from each of the dependencies and remove anything that requires the pinned GitHub versions (e.g., the MSMSulc elements).

[tsalo]

Next problem:

/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/workflows.py:628: in run
    execgraph = generate_expanded_graph(deepcopy(flatgraph))
/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/utils.py:967: in generate_expanded_graph
    graph_in = _remove_nonjoin_identity_nodes(graph_in, keep_iterables=True)
/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/utils.py:847: in _remove_nonjoin_identity_nodes
    _remove_identity_node(graph, node)
/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/utils.py:873: in _remove_identity_node
    _propagate_internal_output(graph, node, field, connections, portinputs)
/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/utils.py:947: in _propagate_internal_output
    value = evaluate_connect_function(src[1], src[2], value)
/usr/local/miniconda/lib/python3.9/site-packages/nipype/pipeline/engine/utils.py:697: in evaluate_connect_function
    output_value = func(first_arg, *list(args))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

lst = <undefined>

>   ???
E   TypeError: '_Undefined' object is not subscriptable

<string>:3: TypeError

So far my debugging process is to mount my local clone of Nipype within run_local_tests.py. I then added some print/raise statements into Nipype to improve the error messages, which let me pin down the problem to

https://github.com/PennLINC/aslprep/blob/3a34405b90488eeb89f502bc7e43c1a3cee1906c/aslprep/workflows/asl/qc.py#L182

Basically, I had accidentally removed the connection of asl_mask_std from asl_std_trans_wf to compute_cbf_qc_wf, so that input was _Undefined, which the connection function _select_last_in_list could not handle. This would have been much easier to debug if Nipype printed out the connection where the error occurred, or even the connection function's name.

[tsalo]

The current problem is that fMRIPrep's init_bold_reg_wf (which I'm using now because it's basically the same as ASLPrep's copy) writes out a report file, but the fMRIPrep DerivativesDataSink config file doesn't allow for the asl suffix. I think I have a few options:

1. Rename the input file(s) so they're called "bold" instead of "asl" to trick the fMRIPrep DerivativesDataSink.
2. Open an issue in niworkflows about adding asl to the list of valid suffixes for figures.
3. Set write_report in the init_bold_reg_wf call to False and write out the report within init_asl_preproc_wf instead.
   • This is the one I'll probably go with for this PR, but I think modifying niworkflows' config file is a good idea in the long run.
4. Open an issue in fMRIPrep about moving some of the more generic workflows (e.g., init_bold_reg_wf, init_bold_t1_trans_wf, init_bold_preproc_trans_wf) from fMRIPrep to niworkflows and making them even more generic.

[tsalo]

fMRIPrep's init_carpetplot_wf uses the fMRIPrep DerivativesDataSink, which has the same suffix problem I mentioned in https://github.com/PennLINC/aslprep/pull/338#issuecomment-1806836077. However, unlike init_bold_reg_wf, there is no write_report parameter. It seems like I should open an issue to add asl to the fMRIPrep/niworkflows datasink config file.

[tsalo]

The carpetplots from fMRIPrep look bad for ASL data because detrend is automatically enabled in the fMRIPrep fMRISummary (which is what is used in init_carpetplot_wf. A few options:

1. Try using a context manager to swap out the fMRISummary version using by init_carpetplot_wf.
   • Given the simplicity of the workflow, this seems like overkill.
2. Open an issue in fMRIPrep requesting that init_carpetplot_wf make detrend a parameter.
   • The fMRIPrep devs may not want to do this.
3. Use a copy of init_carpetplot_wf in ASLPrep. I don't love this idea, but it is a very simple workflow...
   • I'll need an updated clone of fMRISummary.

[tsalo]

Here's a trick ChatGPT gave me to patch in ASLPrep's DerivativesDataSink into an fMRIPrep workflow with a context manager. I'm not using it for init_carpetplot_wf because another parameter (detrend) is buried too deeply into the workflow to make it useful to call directly.

class OverrideDerivativesDataSink:
    """A context manager for temporarily overriding the definition of SomeClass.

    Parameters
    ----------
    None

    Attributes
    ----------
    original_class (type): The original class that is replaced during the override.

    Methods
    -------
    __init__()
        Initialize the context manager.
    __enter__()
        Enters the context manager and performs the class override.
    __exit__(exc_type, exc_value, traceback)
        Exits the context manager and restores the original class definition.
    """

    def __init__(self, module):
        """Initialize the context manager with the target module.

        Parameters
        -----------
        module
            The module where SomeClass should be overridden.
        """
        self.module = module

    def __enter__(self):
        """Enter the context manager and perform the class override.

        Returns
        -------
        OverrideConfoundsDerivativesDataSink
            The instance of the context manager.
        """
        # Save the original class
        self.original_class = self.module.DerivativesDataSink
        # Replace SomeClass with YourOwnClass
        self.module.DerivativesDataSink = DerivativesDataSink
        return self

    def __exit__(self, exc_type, exc_value, traceback):  # noqa: U100
        """Exit the context manager and restore the original class definition.

        Parameters
        ----------
        exc_type : type
            The type of the exception (if an exception occurred).
        exc_value : Exception
            The exception instance (if an exception occurred).
        traceback : traceback
            The traceback information (if an exception occurred).

        Returns
        -------
        None
        """
        # Restore the original class
        self.module.DerivativesDataSink = self.original_class

Usage:

with OverrideDerivativesDataSink(package.module):
    # do some stuff

[tsalo]

It looks like fMRIPrep is currently undergoing a massive refactor, and that's going to have upstream impacts on SDCFlows, (probably) niworkflows, and sMRIPrep. Basically everything I'm doing is based on the version prior to this refactor, which unfortunately doesn't have any releases related to it.

SDCFlows has a bunch of backwards-incompatible changes, and I also need to make a bunch of changes to support ASL data because SDCFlows relies on hardcoded suffixes and datatypes.

[tsalo]

I think I might be stuck.

1. I can't pin to an older version of SDCFlows because I need to make changes to the codebase to support ASL and M0scan files for fieldmap-less and PEPOLAR distortion correction.
2. I can't pin to a modified branch on my fork because CircleCI refuses to install dependencies from GitHub repos.
3. I can update SDCFlows and wait for a new release, but (1) that will delay this PR substantially and (2) I don't have good example code to use because the version of fMRIPrep that uses the new version of SDCFlows doesn't appear to be compatible with ASLPrep (e.g., the fit/apply structure won't work with ASL since we need to calculate CBF and the code was completely rewritten).
   • Actually, since fMRIPrep does include one boldref-space non-BOLD derivatives (i.e., the T2* map for multi-echo data), its new structure may work reasonably well with ASL data.
   • However, other than SDCFlows, fMRIPrep's next branch relies on a bunch of unreleased versions (e.g., niworkflows and smriprep). I can't pin to these versions.

UPDATE:

1. I've been mounting my SDCFlows branch in my local tests. I'll need to open a PR to add ASL to SDCFlows and wait for a release before this refactor is fully buttoned-up.
2. I just gave up on this.
3. Local testing has been a useful workaround. I also ended up fully adopting the fit/apply structure.

[tsalo]

Currently, at least some calls to antsApplyTransforms are failing with the following error:

Inverse does not exist for /aslprep/aslprep/tests/test_data/test_001/derivatives/smriprep/sub-01/anat/sub-01_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5

It looks like feeding in h5 transforms with the invert transforms flag doesn't work? I'll have to ask Matt about it.

UPDATE: Per Matt, h5s don't work with the invert transforms flag. I had to pass in ght eanat-to-standard transform, which fixed the problem.

[tsalo]

Some of the jobs are being killed:

INFO     nipype.workflow:nodes.py:453 [Node] Setting-up "aslprep_0_6_wf.sub_A00086748_wf.asl_preproc_ses_BAS1_wf.asl_anat_wf.resample" in "/src/aslprep/.circleci/work/test_003/aslprep_0_6_wf/sub_A00086748_wf/asl_preproc_ses_BAS1_wf/asl_anat_wf/resample".
INFO     nipype.workflow:nodes.py:716 [Node] Executing "resample" <fmriprep.interfaces.resampling.ResampleSeries>
WARNING  py.warnings:_warnings.py:37 WARNING: Reference space not set
WARNING  py.warnings:_warnings.py:37 WARNING: Reference space not set
/bin/bash: line 1:   434 Killed                  pytest -rP -o log_cli=true -m "test_003" --cov-append --cov-report term-missing --cov=aslprep --data_dir=/src/aslprep/.circleci/data --output_dir=/src/aslprep/.circleci/out --working_dir=/src/aslprep/.circleci/work aslprep

[tsalo]

I'm noticing that, at least in our two GE test datasets, the deltam volumes exhibit a weird background noise pattern. The changes in this PR default to using the deltam volume(s) for the ASL reference image, which is causing some problems. I need to check if the same problem is present in non-GE data (i.e., that deltam, control, or other volume types perform poorly as reference images). If it isn't, then I need to add a GE-specific case for reference image estimation, to use the M0 scan instead.

EDIT: The reference looks fine for test_001, test_003, pcasl_singlepld_siemens, pcasl_singlepld_philips, and pcasl_multipld

UPDATE: It appears GE-specific, and having a step that selects the M0 scan to use as the reference for GE data seems to have fixed the problem.

[tsalo]

The Philips results look terrible (see below), but most of the others are looking pretty good. Actually, I'm seeing a little bottom-row weirdness in the Siemens PVC BASIL CBF results too.

[codecov-commenter]

Codecov Report

Attention: 203 lines in your changes are missing coverage. Please review.

Comparison is base (10fc6ba) 62.80% compared to head (d3b027c) 75.99%.

Files	Patch %	Lines
aslprep/cli/run.py	0.00%	40 Missing :warning:
aslprep/workflows/asl/fit.py	87.07%	14 Missing and 5 partials :warning:
aslprep/workflows/asl/base.py	83.92%	14 Missing and 4 partials :warning:
aslprep/workflows/asl/confounds.py	82.17%	17 Missing and 1 partial :warning:
aslprep/utils/plotting.py	68.62%	8 Missing and 8 partials :warning:
aslprep/workflows/asl/apply.py	25.00%	15 Missing :warning:
aslprep/config.py	86.66%	7 Missing and 7 partials :warning:
aslprep/utils/bids.py	70.58%	5 Missing and 5 partials :warning:
aslprep/interfaces/plotting.py	74.28%	8 Missing and 1 partial :warning:
aslprep/interfaces/reference.py	81.81%	4 Missing and 4 partials :warning:
... and 7 more

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #338 +/- ## =========================================== + Coverage 62.80% 75.99% +13.19% =========================================== Files 56 39 -17 Lines 5186 4120 -1066 Branches 722 595 -127 =========================================== - Hits 3257 3131 -126 + Misses 1676 791 -885 + Partials 253 198 -55 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[tsalo]

I think the problem stems from a mismatch in the coverage between the M0 scan and the ASL scan. The high values are from voxels outside the M0 scan (i.e., those values are like 0.0001) but inside the ASL scan.

[tsalo]

Yup, everything's looking much better now that I take the intersection of the T1w, ASL, and M0 masks.

[tsalo]

I'm testing locally on a PNC subject with three sessions, using a mounted fork of SDCFlows:

1. Field map specified through IntendedFor.
2. Field map specified through B0FieldSource/B0FieldIdentifier.
3. No field map.

2 and 3 work as expected, so far, but 1 isn't finding the field map. I still need to try out SyN, though I'm not optimistic with these versions given https://github.com/nipreps/fmriprep/issues/3158.

Also, the aslref-to-fmap transform was written out to the anat folder, so it looks like the OverrideDerivativesDataSink patch in ds_fmapreg_wf didn't work. And some steps in the surface projection were killed on my laptop. Need to see if I can just run with more memory.

[tsalo]

Error with CIFTI output:

nipype.pipeline.engine.nodes.NodeExecutionError: Exception raised while executing Node ds_asl_cifti.

Traceback:
    Traceback (most recent call last):
      File "/usr/local/miniconda/lib/python3.10/site-packages/nipype/interfaces/base/core.py", line 397, in run
        runtime = self._run_interface(runtime)
      File "/usr/local/miniconda/lib/python3.10/site-packages/niworkflows/interfaces/bids.py", line 630, in _run_interface
        raise ValueError(f"Could not build path with entities {out_entities}.")
    ValueError: Could not build path with entities {'subject': '01', 'session': '2', 'acquisition': 'se', 'suffix': 'asl', 'datatype': 'perf', 'extension': 'dtseries.nii', 'density': '91k', 'space': 'fsLR'}.

[tsalo]

It looks like the only thing I need to get working is the resampling workflow!

EDIT: Also, there aren't workflow plots for the workflows in aslprep.workflows.asl.outputs or aslprep.workflows.asl.apply.init_asl_cifti_resample_wf.

[tsalo]

😭 I got an error with the CIFTI resampling:

While running:
        wb_command -volume-to-surface-mapping /src/aslprep/.circleci/work/test_003_full/aslprep_0_6_wf/sub_01_wf/asl_preproc_ses_1_wf/ds_asl_cifti_wf/warp_mean_cbf_to_anat/sub-01_ses-1_asl_DeltaMOrCBF_meancbf_trans.nii.gz /src/aslprep/.circleci/data/anatomical/smriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii sub-01_hemi-L_midthickness.surf_mapped.func.gii -ribbon-constrained /src/aslprep/.circleci/data/anatomical/smriprep/sub-01/anat/sub-01_hemi-L_white.surf.gii /src/aslprep/.circleci/data/anatomical/smriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii -volume-roi /src/aslprep/.circleci/work/test_003_full/aslprep_0_6_wf/sub_01_wf/asl_preproc_ses_1_wf/asl_cifti_resample_wf/goodvoxels_bold_mask_wf/goodvoxels_mask/sub-01_ses-1_asl_reduced_std_maths_maths_maths_thresh_maths.nii.gz

        ERROR: roi volume is not in the same volume space as input volume

I think I need the working directory to figure out what's going wrong, but that's not really feasible with local testing or using CircleCI, so I'm leaning toward merging this once I have CI passing, building the Singularity image on CUBIC, and testing there.