FIX: Vastly improve multi-echo handling (#220)

* FIX: Rework BOLD grouping to support multiecho

- Added a new class, `BOLDGrouping`, to better handle binning of BOLD data.
This allows the passing of multiple BOLD file to the functional preprocessing workflow,
which was treating every echo of a multiecho dataset as a unique run.

* FIX: Use more verbose typing for older pythons

* FIX: Correct parameter name

* FIX: Ensure `sessions` are lists, grouping values are selected

* FIX: Update T2star interface, connections

* ENH: Add tedana to dependencies

* STY: formatting

* TST: Add empty echo files for unittest

Author: mgxd

nibabies/data/sub-01_run-01_echo-1_bold.nii.gz

@@ -1,7 +1,25 @@
-#!/usr/bin/env python
-# -*- 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:
+#
+# Copyright 2021 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
 """
 Multi-echo EPI
 ~~~~~~~~~~~~~~
@@ -26,42 +44,51 @@
     traits,
 )
 
-LOGGER = logging.getLogger("nipype.interface")
+LOGGER = logging.getLogger('nipype.interface')
 
 
 class T2SMapInputSpec(CommandLineInputSpec):
     in_files = traits.List(
         File(exists=True),
-        argstr="-d %s",
+        argstr='-d %s',
         position=1,
         mandatory=True,
         minlen=3,
-        desc="multi-echo BOLD EPIs",
+        desc='multi-echo BOLD EPIs',
     )
     echo_times = traits.List(
-        traits.Float, argstr="-e %s", position=2, mandatory=True, minlen=3, desc="echo times"
+        traits.Float,
+        argstr='-e %s',
+        position=2,
+        mandatory=True,
+        minlen=3,
+        desc='echo times',
     )
-    fittype = traits.Enum(
-        "curvefit",
-        "loglin",
-        argstr="--fittype %s",
+    mask_file = File(
+        argstr='--mask %s',
         position=3,
+        desc='mask file',
+        exists=True,
+    )
+    fittype = traits.Enum(
+        'curvefit',
+        'loglin',
+        argstr='--fittype %s',
+        position=4,
         usedefault=True,
-        desc=(
-            "Desired fitting method: "
-            '"loglin" means that a linear model is fit '
-            "to the log of the data. "
-            '"curvefit" means that a more computationally '
-            "demanding monoexponential model is fit "
-            "to the raw data."
-        ),
+        desc='Desired fitting method: '
+        '"loglin" means that a linear model is fit '
+        'to the log of the data. '
+        '"curvefit" means that a more computationally '
+        'demanding monoexponential model is fit '
+        'to the raw data.',
     )
 
 
 class T2SMapOutputSpec(TraitedSpec):
-    t2star_map = File(exists=True, desc="limited T2* map")
-    s0_map = File(exists=True, desc="limited S0 map")
-    optimal_comb = File(exists=True, desc="optimally combined ME-EPI time series")
+    t2star_map = File(exists=True, desc='limited T2* map')
+    s0_map = File(exists=True, desc='limited S0 map')
+    optimal_comb = File(exists=True, desc='optimally combined ME-EPI time series')
 
 
 class T2SMap(CommandLine):
@@ -73,30 +100,29 @@ class T2SMap(CommandLine):
     =======
     >>> from nibabies.interfaces import multiecho
     >>> t2smap = multiecho.T2SMap()
-    >>> t2smap.inputs.in_files = [
-    ...     data_dir / 'sub-01_run-01_echo-1_bold.nii.gz',
-    ...     data_dir / 'sub-01_run-01_echo-2_bold.nii.gz',
-    ...     data_dir / 'sub-01_run-01_echo-3_bold.nii.gz']
+    >>> t2smap.inputs.in_files = ['sub-01_run-01_echo-1_bold.nii.gz', \
+                                  'sub-01_run-01_echo-2_bold.nii.gz', \
+                                  'sub-01_run-01_echo-3_bold.nii.gz']
     >>> t2smap.inputs.echo_times = [0.013, 0.027, 0.043]
-    >>> t2smap.cmdline
-    't2smap -d .../sub-01_run-01_echo-1_bold.nii.gz \
-    .../sub-01_run-01_echo-2_bold.nii.gz \
-    .../sub-01_run-01_echo-3_bold.nii.gz -e 13.0 27.0 43.0 --fittype curvefit'
+    >>> t2smap.cmdline  # doctest: +ELLIPSIS
+    't2smap -d sub-01_run-01_echo-1_bold.nii.gz sub-01_run-01_echo-2_bold.nii.gz \
+sub-01_run-01_echo-3_bold.nii.gz -e 13.0 27.0 43.0 --fittype curvefit'
+
     """
 
-    _cmd = "t2smap"
+    _cmd = 't2smap'
     input_spec = T2SMapInputSpec
     output_spec = T2SMapOutputSpec
 
     def _format_arg(self, name, trait_spec, value):
-        if name == "echo_times":
+        if name == 'echo_times':
             value = [te * 1000 for te in value]
         return super(T2SMap, self)._format_arg(name, trait_spec, value)
 
     def _list_outputs(self):
         outputs = self._outputs().get()
         out_dir = os.getcwd()
-        outputs["t2star_map"] = os.path.join(out_dir, "T2starmap.nii.gz")
-        outputs["s0_map"] = os.path.join(out_dir, "S0map.nii.gz")
-        outputs["optimal_comb"] = os.path.join(out_dir, "desc-optcom_bold.nii.gz")
+        outputs['t2star_map'] = os.path.join(out_dir, 'T2starmap.nii.gz')
+        outputs['s0_map'] = os.path.join(out_dir, 'S0map.nii.gz')
+        outputs['optimal_comb'] = os.path.join(out_dir, 'desc-optcom_bold.nii.gz')
         return outputs

nibabies/data/sub-01_run-01_echo-2_bold.nii.gz

@@ -4,7 +4,27 @@
 import json
 import os
 import sys
+from dataclasses import dataclass, field
 from pathlib import Path
+from typing import IO, List, Union
+
+
+@dataclass
+class BOLDGrouping:
+    """This class is used to facilitate the grouping of BOLD series."""
+
+    session: Union[str, None]
+    pe_dir: str
+    readout: float
+    multiecho_id: str = None
+    files: List[IO] = field(default_factory=list)
+
+    @property
+    def name(self) -> str:
+        return f"{self.session}-{self.pe_dir}-{self.readout}-{self.multiecho_id}"
+
+    def add_file(self, fl) -> None:
+        self.files.append(fl)
 
 
 def write_bidsignore(deriv_dir):
@@ -107,7 +127,7 @@ def _unique(inlist):
     return {k: _unique(v) for k, v in entities.items()}
 
 
-def group_bolds_ref(*, layout, subject, session_id):
+def group_bolds_ref(*, layout, subject, sessions=None):
     """
     Extracts BOLD files from a BIDS dataset and combines them into buckets.
     Files in a bucket share:
@@ -123,8 +143,7 @@ def group_bolds_ref(*, layout, subject, session_id):
         Initialized BIDSLayout
     subject : str
         The subject ID
-    session_id : :obj:`str`, None, or ``False``
-        The session identifier. If ``False``, all sessions will be used (default).
+    sessions : None
 
     Outputs
     -------
@@ -136,8 +155,8 @@ def group_bolds_ref(*, layout, subject, session_id):
     Limitations
     -----------
     Single-band reference (sbref) are excluded.
-    Does not group multi-echo data.
     """
+    import re
     from contextlib import suppress
     from itertools import product
 
@@ -148,19 +167,12 @@ def group_bolds_ref(*, layout, subject, session_id):
         "extension": (".nii", ".nii.gz"),
         "scope": "raw",  # Ensure derivatives are not captured
     }
-    # list of tuples with unique combinations
-    combinations = []
-    # list of lists containing filenames that apply per combination
-    files = []
+    # dictionary containing unique Groupings and files
+    groupings = {}
     # list of all BOLDS encountered
     all_bolds = []
 
-    # TODO: Simplify with pybids.layout.Query.OPTIONAL
-    sessions = (
-        [session_id]
-        if session_id is not False
-        else layout.get_sessions(subject=subject, scope="raw")
-    )
+    sessions = sessions if sessions else layout.get_sessions(subject=subject, scope="raw")
 
     for ses, suffix in sorted(product(sessions or (None,), {"bold"})):
         # bold files same session
@@ -169,7 +181,13 @@ def group_bolds_ref(*, layout, subject, session_id):
         if bolds is None:
             continue
 
-        for bold in bolds:
+        for i, bold in enumerate(bolds):
+            multiecho_id = None
+            # multi-echo should be grouped together
+            if 'echo' in bold.entities:
+                # create unique id by dropping "_echo-{i}"
+                multiecho_id = re.sub(r"_echo-\d+", "", bold.filename)
+
             # session, pe, ro
             meta = bold.get_metadata()
             pe_dir = meta.get("PhaseEncodingDirection")
@@ -180,36 +198,38 @@ def group_bolds_ref(*, layout, subject, session_id):
             if ro is not None:
                 meta.update({"TotalReadoutTime": ro})
 
-            comb = (ses, pe_dir, ro)
+            grouping = BOLDGrouping(
+                session=ses,
+                pe_dir=pe_dir,
+                readout=ro,
+                multiecho_id=multiecho_id,
+            )
 
             if any(v is None for v in (pe_dir, ro)):
                 # cannot be certain so treat as unique
-                combinations.append(comb)
-                files.append([bold.path])
-            elif comb in combinations:
-                # do not add a new entry to the combinations
-                # instead append the file to the existing bucket
-                idx = combinations.index(comb)
-                files[idx].append(bold.path)
+                grouping.add_file(bold.path)
+                groupings[f'unknown{i}'] = grouping
             else:
-                # add a new entry and start a file bucket
-                combinations.append(comb)
-                files.append([bold.path])
+                try:
+                    grouping = groupings[grouping.name]
+                except KeyError:
+                    groupings[grouping.name] = grouping
+
+                grouping.add_file(bold.path)
 
         all_bolds += bolds
 
-        if (len(combinations) != len(files)) or (len(all_bolds) != sum([len(x) for x in files])):
+        if len(all_bolds) != sum([len(g.files) for _, g in groupings.items()]):
             msg = f"""Error encountered when grouping BOLD runs.
-Combinations: {combinations}
-Sorted files: {files}
+Combinations: {groupings}
 BOLD files: {bolds}
 
-Please share this information with the nibabies developers at:
+Please file a bug-report with the nibabies developers at:
 https://github.com/nipreps/nibabies/issues/new/choose
 """
             raise RuntimeError(msg)
 
-    return combinations, files
+    return groupings
 
 
 def validate_input_dir(exec_env, bids_dir, participant_label):

nibabies/data/sub-01_run-01_echo-3_bold.nii.gz

@@ -400,7 +400,7 @@ def init_single_subject_wf(subject_id, session_id=None):
         fmap_estimators = find_estimators(
             layout=config.execution.layout,
             subject=subject_id,
-            sessions=session_id,
+            sessions=[session_id],
             fmapless=False,  # config.workflow.use_syn,
             force_fmapless=False,  # config.workflow.force_syn,
         )
@@ -425,24 +425,25 @@ def init_single_subject_wf(subject_id, session_id=None):
     # 3) total readout time
     from niworkflows.workflows.epi.refmap import init_epi_reference_wf
 
-    _, bold_groupings = group_bolds_ref(
+    bold_groupings = group_bolds_ref(
         layout=config.execution.layout,
         subject=subject_id,
-        session_id=session_id,
+        sessions=[session_id],
     )
-    if any(not x for x in bold_groupings):
-        print("No BOLD files found for one or more reference groupings")
-        return workflow
 
     func_preproc_wfs = []
     has_fieldmap = bool(fmap_estimators)
-    for idx, bold_files in enumerate(bold_groupings):
+    for idx, grouping in enumerate(bold_groupings.values()):
         bold_ref_wf = init_epi_reference_wf(
             auto_bold_nss=True,
             name=f"bold_reference_wf{idx}",
             omp_nthreads=config.nipype.omp_nthreads,
         )
-        bold_ref_wf.inputs.inputnode.in_files = bold_files
+        bold_files = grouping.files
+        bold_ref_wf.inputs.inputnode.in_files = grouping.files
+
+        if grouping.multiecho_id is not None:
+            bold_files = [bold_files]
         for idx, bold_file in enumerate(bold_files):
             func_preproc_wf = init_func_preproc_wf(
                 bold_file,

nibabies/interfaces/multiecho.py

@@ -68,7 +68,7 @@ def init_bold_t2s_wf(echo_times, mem_gb, omp_nthreads, name="bold_t2s_wf"):
 The optimally combined time series was carried forward as the *preprocessed BOLD*.
 """
 
-    inputnode = pe.Node(niu.IdentityInterface(fields=["bold_file"]), name="inputnode")
+    inputnode = pe.Node(niu.IdentityInterface(fields=["bold_file", "bold_mask"]), name="inputnode")
 
     outputnode = pe.Node(niu.IdentityInterface(fields=["bold"]), name="outputnode")
 
@@ -78,7 +78,8 @@ def init_bold_t2s_wf(echo_times, mem_gb, omp_nthreads, name="bold_t2s_wf"):
 
     # fmt: off
     workflow.connect([
-        (inputnode, t2smap_node, [('bold_file', 'in_files')]),
+        (inputnode, t2smap_node, [('bold_file', 'in_files'),
+                                  ('bold_mask', 'mask_file')]),
         (t2smap_node, outputnode, [('optimal_comb', 'bold')]),
     ])
     # fmt: on

nibabies/utils/bids.py

@@ -30,6 +30,7 @@ install_requires =
     scipy ~= 1.6.0; python_version<'3.8'
     sdcflows @ git+https://github.com/nipreps/sdcflows.git@master
     smriprep ~= 0.8.1
+    tedana ~= 0.0.12
     templateflow >= 0.6.0
     toml
 test_requires =