diff --git a/docs/builtin.rst b/docs/builtin.rst index ef45e3d54..6ad2eed67 100644 --- a/docs/builtin.rst +++ b/docs/builtin.rst @@ -209,22 +209,30 @@ Available - Apply parcellation and perform aggregation function - Done - 0.0.1 - * - :class:`.FunctionalConnectivityParcels` - - Compute functional connectivity over parcellation - - Done - - 0.0.1 - * - :class:`.CrossParcellationFC` - - Compute functional connectivity across two parcellations - - Done - - 0.0.1 * - :class:`.SphereAggregation` - Spherical aggregation using mean - Done - 0.0.1 + * - :class:`.MapsAggregation` + - Apply maps (probabilistic atlas) and extract mean time series per region + - Done + - 0.0.7 + * - :class:`.FunctionalConnectivityParcels` + - Compute functional connectivity over parcellation + - Done + - 0.0.1 * - :class:`.FunctionalConnectivitySpheres` - Compute functional connectivity over spheres placed on coordinates - Done - 0.0.1 + * - :class:`.FunctionalConnectivityMaps` + - Compute functional connectivity over maps (probabilistic atlas) + - Done + - 0.0.7 + * - :class:`.CrossParcellationFC` + - Compute functional connectivity across two parcellations + - Done + - 0.0.1 * - :class:`.RSSETSMarker` - Compute root sum of squares of edgewise timeseries - Done @@ -237,6 +245,10 @@ Available - Calculate regional homogeneity over spheres placed on coordinates - Done - 0.0.1 + * - :class:`.ReHoMaps` + - Calculate regional homogeneity over maps (probabilistic atlas) + - Done + - 0.0.7 * - :class:`.ALFFParcels` - Calculate (f)ALFF and aggregate using parcellations - Done @@ -245,6 +257,10 @@ Available - Calculate (f)ALFF and aggregate using spheres placed on coordinates - Done - 0.0.1 + * - :class:`.ALFFMaps` + - Calculate (f)ALFF and aggregate using maps (probabilistic atlas) + - Done + - 0.0.7 * - :class:`.EdgeCentricFCParcels` - | Calculate edge-centric functional connectivity over parcellation, as | found in @@ -257,6 +273,12 @@ Available | `Jo et al. (2021) `_ - Done - 0.0.2 + * - :class:`.EdgeCentricFCMaps` + - | Calculate edge-centric functional connectivity over maps (probabilistic atlas), + | as found in + | `Jo et al. (2021) `_ + - Done + - 0.0.7 * - :class:`.TemporalSNRParcels` - Calculate temporal signal-to-noise ratio using parcellations - Done @@ -266,6 +288,10 @@ Available | coordinates - Done - 0.0.2 + * - :class:`.TemporalSNRMaps` + - Calculate temporal signal-to-noise ratio using maps (probabilistic atlas) + - Done + - 0.0.7 * - :class:`.HurstExponent` - | Calculate Hurst exponent of a time series as found in | `Peng et al. (1995) `_ diff --git a/docs/changes/newsfragments/460.feature b/docs/changes/newsfragments/460.feature new file mode 100644 index 000000000..db09db463 --- /dev/null +++ b/docs/changes/newsfragments/460.feature @@ -0,0 +1 @@ +Introduce :class:`.MapsAggregation`, :class:`.FunctionalConnectivityMaps`, :class:`.EdgeCentricFCMaps`, :class:`.ReHoMaps`, :class:`.ALFFMaps` and :class:`.TemporalSNRMaps` by `Synchon Mandal`_ diff --git a/docs/extending/maps.rst b/docs/extending/maps.rst index 74516adfa..154c27bb7 100644 --- a/docs/extending/maps.rst +++ b/docs/extending/maps.rst @@ -116,7 +116,7 @@ file, we can save the above code in a Python file, say Afterwards continue configuring the rest of the pipeline in this YAML file, and you will be able to use this maps using the name you gave the maps when registering it. For example, we can add a -``MapsAggregation`` Marker to demonstrate how this can be done: +:class:`.MapsAggregation` Marker to demonstrate how this can be done: .. code-block:: yaml @@ -124,7 +124,6 @@ maps when registering it. For example, we can add a - name: CustomMaps_mean kind: MapsAggregation maps: my_custom_maps - method: mean Now, you can simply use this YAML file to run your pipeline. diff --git a/junifer/conftest.py b/junifer/conftest.py index 1a2e7ba80..28e9dcec1 100644 --- a/junifer/conftest.py +++ b/junifer/conftest.py @@ -1,10 +1,14 @@ """Provide conftest for pytest.""" # Authors: Federico Raimondo +# Synchon Mandal # License: AGPL +from pathlib import Path + import pytest +from junifer.datagrabber import PatternDataladDataGrabber from junifer.utils.singleton import Singleton @@ -23,3 +27,31 @@ def reset_singletons() -> None: # Force deleting the singletons for elem in to_remove: del elem + + +@pytest.fixture +def maps_datagrabber(tmp_path: Path) -> PatternDataladDataGrabber: + """Return a PatternDataladDataGrabber for maps testing. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + dg = PatternDataladDataGrabber( + uri="https://github.com/OpenNeuroDatasets/ds005226.git", + types=["BOLD"], + patterns={ + "BOLD": { + "pattern": ( + "derivatives/pre-processed_data/space-MNI/{subject}/" + "{subject-padded}_task-{task}_run-{run}_space-MNI152NLin6Asym" + "_res-2_desc-preproc_bold.nii.gz" + ), + "space": "MNI152NLin6Asym", + }, + }, + replacements=["subject", "subject-padded", "task", "run"], + ) + return dg diff --git a/junifer/markers/__init__.pyi b/junifer/markers/__init__.pyi index 73462565f..34c35c7d0 100644 --- a/junifer/markers/__init__.pyi +++ b/junifer/markers/__init__.pyi @@ -1,17 +1,23 @@ __all__ = [ "BaseMarker", "RSSETSMarker", + "MapsAggregation", "ParcelAggregation", "SphereAggregation", + "FunctionalConnectivityMaps", "FunctionalConnectivityParcels", "FunctionalConnectivitySpheres", "CrossParcellationFC", + "EdgeCentricFCMaps", "EdgeCentricFCParcels", "EdgeCentricFCSpheres", + "ReHoMaps", "ReHoParcels", "ReHoSpheres", + "ALFFMaps", "ALFFParcels", "ALFFSpheres", + "TemporalSNRMaps", "TemporalSNRParcels", "TemporalSNRSpheres", "BrainPrint", @@ -19,18 +25,22 @@ __all__ = [ from .base import BaseMarker from .ets_rss import RSSETSMarker +from .maps_aggregation import MapsAggregation from .parcel_aggregation import ParcelAggregation from .sphere_aggregation import SphereAggregation from .functional_connectivity import ( + FunctionalConnectivityMaps, FunctionalConnectivityParcels, FunctionalConnectivitySpheres, CrossParcellationFC, + EdgeCentricFCMaps, EdgeCentricFCParcels, EdgeCentricFCSpheres, ) -from .reho import ReHoParcels, ReHoSpheres -from .falff import ALFFParcels, ALFFSpheres +from .reho import ReHoMaps, ReHoParcels, ReHoSpheres +from .falff import ALFFMaps, ALFFParcels, ALFFSpheres from .temporal_snr import ( + TemporalSNRMaps, TemporalSNRParcels, TemporalSNRSpheres, ) diff --git a/junifer/markers/falff/__init__.pyi b/junifer/markers/falff/__init__.pyi index e70e972cf..1b003e374 100644 --- a/junifer/markers/falff/__init__.pyi +++ b/junifer/markers/falff/__init__.pyi @@ -1,4 +1,5 @@ -__all__ = ["ALFFParcels", "ALFFSpheres"] +__all__ = ["ALFFMaps", "ALFFParcels", "ALFFSpheres"] +from .falff_maps import ALFFMaps from .falff_parcels import ALFFParcels from .falff_spheres import ALFFSpheres diff --git a/junifer/markers/falff/falff_maps.py b/junifer/markers/falff/falff_maps.py new file mode 100644 index 000000000..1d6434044 --- /dev/null +++ b/junifer/markers/falff/falff_maps.py @@ -0,0 +1,151 @@ +"""Provide class for ALFF / fALFF on maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, Optional, Union + +from ...api.decorators import register_marker +from ...utils import logger +from ..maps_aggregation import MapsAggregation +from .falff_base import ALFFBase + + +__all__ = ["ALFFMaps"] + + +@register_marker +class ALFFMaps(ALFFBase): + """Class for ALFF / fALFF on maps. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + using : {"junifer", "afni"} + Implementation to use for computing ALFF: + + * "junifer" : Use ``junifer``'s own ALFF implementation + * "afni" : Use AFNI's ``3dRSFC`` + + highpass : positive float, optional + The highpass cutoff frequency for the bandpass filter. If 0, + it will not apply a highpass filter (default 0.01). + lowpass : positive float, optional + The lowpass cutoff frequency for the bandpass filter (default 0.1). + tr : positive float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIfTI header (default None). + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + Notes + ----- + The ``tr`` parameter is crucial for the correctness of fALFF/ALFF + computation. If a dataset is correctly preprocessed, the ``tr`` should be + extracted from the NIfTI without any issue. However, it has been + reported that some preprocessed data might not have the correct ``tr`` in + the NIfTI header. + + ALFF/fALFF are computed using a bandpass butterworth filter. See + :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more + details. + + """ + + def __init__( + self, + maps: str, + using: str, + highpass: float = 0.01, + lowpass: float = 0.1, + tr: Optional[float] = None, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + name: Optional[str] = None, + ) -> None: + # Superclass init first to validate `using` parameter + super().__init__( + highpass=highpass, + lowpass=lowpass, + using=using, + tr=tr, + name=name, + ) + self.maps = maps + self.masks = masks + + def compute( + self, + input: dict[str, Any], + extra_input: Optional[dict[str, Any]] = None, + ) -> dict[str, Any]: + """Compute. + + Parameters + ---------- + input : dict + The BOLD data as dictionary. + extra_input : dict, optional + The other fields in the pipeline data object (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``alff`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + * ``falff`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + """ + logger.info("Calculating ALFF / fALFF for maps") + + # Compute ALFF + fALFF + alff_output, falff_output, alff_output_path, falff_output_path = ( + self._compute(input_data=input) + ) + + # Perform aggregation on ALFF + fALFF + aggregation_alff_input = dict(input.items()) + aggregation_falff_input = dict(input.items()) + aggregation_alff_input["data"] = alff_output + aggregation_falff_input["data"] = falff_output + aggregation_alff_input["path"] = alff_output_path + aggregation_falff_input["path"] = falff_output_path + + return { + "alff": { + **MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute( + input=aggregation_alff_input, + extra_input=extra_input, + )["aggregation"], + }, + "falff": { + **MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute( + input=aggregation_falff_input, + extra_input=extra_input, + )["aggregation"], + }, + } diff --git a/junifer/markers/falff/tests/test_falff_maps.py b/junifer/markers/falff/tests/test_falff_maps.py new file mode 100644 index 000000000..aee7a92b2 --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_maps.py @@ -0,0 +1,167 @@ +"""Provide tests for ALFFMaps.""" + +# Authors: Synchon Mandal +# License: AGPL + +import logging +from pathlib import Path + +import pytest +import scipy as sp + +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import ALFFMaps +from junifer.pipeline import WorkDirManager +from junifer.pipeline.utils import _check_afni +from junifer.storage import HDF5FeatureStorage + + +MAPS = "Smith_rsn_10" + + +@pytest.mark.parametrize( + "feature", + [ + "alff", + "falff", + ], +) +def test_ALFFMaps_get_output_type(feature: str) -> None: + """Test ALFFMaps get_output_type(). + + Parameters + ---------- + feature : str + The parametrized feature name. + + """ + assert "vector" == ALFFMaps( + maps=MAPS, + using="junifer", + ).get_output_type(input_type="BOLD", output_feature=feature) + + +def test_ALFFMaps( + caplog: pytest.LogCaptureFixture, + tmp_path: Path, + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test ALFFMaps. + + Parameters + ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path + + with caplog.at_level(logging.DEBUG): + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + + # Initialize marker + marker = ALFFMaps( + maps=MAPS, + using="junifer", + ) + # Check correct output + for name in ["alff", "falff"]: + assert "vector" == marker.get_output_type( + input_type="BOLD", output_feature=name + ) + + # Fit transform marker on data + output = marker.fit_transform(element_data) + + assert "Creating cache" in caplog.text + + # Get BOLD output + assert "BOLD" in output + for feature in output["BOLD"].keys(): + output_bold = output["BOLD"][feature] + # Assert BOLD output keys + assert "data" in output_bold + assert "col_names" in output_bold + + output_bold_data = output_bold["data"] + # Assert BOLD output data dimension + assert output_bold_data.ndim == 2 + assert output_bold_data.shape == (1, 10) + + # Reset log capture + caplog.clear() + # Initialize storage + storage = HDF5FeatureStorage(tmp_path / "falff_maps.hdf5") + # Fit transform marker on data with storage + marker.fit_transform( + input=element_data, + storage=storage, + ) + # Check stored feature name + features = storage.list_features() + assert any( + x["name"] in ["BOLD_ALFFMaps_alff", "BOLD_ALFFMaps_falff"] + for x in features.values() + ) + # Cache working correctly + assert "Creating cache" not in caplog.text + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires AFNI to be in PATH" +) +def test_ALFFMaps_comparison( + tmp_path: Path, maps_datagrabber: PatternDataladDataGrabber +) -> None: + """Test ALFFMaps implementation comparison. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path + + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + + # Initialize marker + junifer_marker = ALFFMaps( + maps=MAPS, + using="junifer", + ) + # Fit transform marker on data + junifer_output = junifer_marker.fit_transform(element_data) + # Get BOLD output + junifer_output_bold = junifer_output["BOLD"] + + # Initialize marker + afni_marker = ALFFMaps( + maps=MAPS, + using="afni", + ) + # Fit transform marker on data + afni_output = afni_marker.fit_transform(element_data) + # Get BOLD output + afni_output_bold = afni_output["BOLD"] + + for feature in afni_output_bold.keys(): + # Check for Pearson correlation coefficient + r, _ = sp.stats.pearsonr( + junifer_output_bold[feature]["data"][0], + afni_output_bold[feature]["data"][0], + ) + assert r > 0.97 diff --git a/junifer/markers/functional_connectivity/__init__.pyi b/junifer/markers/functional_connectivity/__init__.pyi index dd7790d35..6df28438f 100644 --- a/junifer/markers/functional_connectivity/__init__.pyi +++ b/junifer/markers/functional_connectivity/__init__.pyi @@ -1,13 +1,17 @@ __all__ = [ + "FunctionalConnectivityMaps", "FunctionalConnectivityParcels", "FunctionalConnectivitySpheres", "CrossParcellationFC", + "EdgeCentricFCMaps", "EdgeCentricFCParcels", "EdgeCentricFCSpheres", ] +from .functional_connectivity_maps import FunctionalConnectivityMaps from .functional_connectivity_parcels import FunctionalConnectivityParcels from .functional_connectivity_spheres import FunctionalConnectivitySpheres from .crossparcellation_functional_connectivity import CrossParcellationFC +from .edge_functional_connectivity_maps import EdgeCentricFCMaps from .edge_functional_connectivity_parcels import EdgeCentricFCParcels from .edge_functional_connectivity_spheres import EdgeCentricFCSpheres diff --git a/junifer/markers/functional_connectivity/edge_functional_connectivity_maps.py b/junifer/markers/functional_connectivity/edge_functional_connectivity_maps.py new file mode 100644 index 000000000..9f833a203 --- /dev/null +++ b/junifer/markers/functional_connectivity/edge_functional_connectivity_maps.py @@ -0,0 +1,115 @@ +"""Provide class for edge-centric functional connectivity using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, Optional, Union + +from ...api.decorators import register_marker +from ..maps_aggregation import MapsAggregation +from ..utils import _ets +from .functional_connectivity_base import FunctionalConnectivityBase + + +__all__ = ["EdgeCentricFCMaps"] + + +@register_marker +class EdgeCentricFCMaps(FunctionalConnectivityBase): + """Class for edge-centric FC using maps. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (default None). + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + name : str, optional + The name of the marker. If None, will use + ``BOLD_EdgeCentricFCParcels`` (default None). + + References + ---------- + .. [1] Jo et al. (2021) + Subject identification using edge-centric functional connectivity. + https://doi.org/10.1016/j.neuroimage.2021.118204 + + """ + + def __init__( + self, + maps: str, + conn_method: str = "correlation", + conn_method_params: Optional[dict] = None, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + name: Optional[str] = None, + ) -> None: + self.maps = maps + super().__init__( + conn_method=conn_method, + conn_method_params=conn_method_params, + masks=masks, + name=name, + ) + + def aggregate( + self, input: dict[str, Any], extra_input: Optional[dict] = None + ) -> dict: + """Perform maps aggregation and ETS computation. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``aggregation`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + """ + # Perform aggregation + aggregation = MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute(input, extra_input=extra_input) + # Compute edgewise timeseries + ets, edge_names = _ets( + bold_ts=aggregation["aggregation"]["data"], + roi_names=aggregation["aggregation"]["col_names"], + ) + + return { + "aggregation": { + "data": ets, + "col_names": edge_names, + }, + } diff --git a/junifer/markers/functional_connectivity/functional_connectivity_maps.py b/junifer/markers/functional_connectivity/functional_connectivity_maps.py new file mode 100644 index 000000000..90889869b --- /dev/null +++ b/junifer/markers/functional_connectivity/functional_connectivity_maps.py @@ -0,0 +1,95 @@ +"""Provide class for functional connectivity using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, Optional, Union + +from ...api.decorators import register_marker +from ..maps_aggregation import MapsAggregation +from .functional_connectivity_base import FunctionalConnectivityBase + + +__all__ = ["FunctionalConnectivityMaps"] + + +@register_marker +class FunctionalConnectivityMaps(FunctionalConnectivityBase): + """Class for functional connectivity using maps. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (default None). + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + name : str, optional + The name of the marker. If None, will use + ``BOLD_FunctionalConnectivityMaps`` (default None). + + """ + + def __init__( + self, + maps: str, + conn_method: str = "correlation", + conn_method_params: Optional[dict] = None, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + name: Optional[str] = None, + ) -> None: + self.maps = maps + super().__init__( + conn_method=conn_method, + conn_method_params=conn_method_params, + masks=masks, + name=name, + ) + + def aggregate( + self, input: dict[str, Any], extra_input: Optional[dict] = None + ) -> dict: + """Perform maps aggregation. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``aggregation`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + """ + return MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute(input=input, extra_input=extra_input) diff --git a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_maps.py b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_maps.py new file mode 100644 index 000000000..ae44c9d07 --- /dev/null +++ b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_maps.py @@ -0,0 +1,74 @@ +"""Provide tests for edge-centric functional connectivity using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from pathlib import Path + +import pytest + +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import EdgeCentricFCMaps +from junifer.storage import HDF5FeatureStorage + + +@pytest.mark.parametrize( + "conn_method_params", + [ + {"empirical": False}, + {"empirical": True}, + ], +) +def test_EdgeCentricFCMaps( + tmp_path: Path, + maps_datagrabber: PatternDataladDataGrabber, + conn_method_params: dict[str, bool], +) -> None: + """Test EdgeCentricFCMaps. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + conn_method_params : dict + The parametrized parameters to connectivity measure method. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Setup marker + marker = EdgeCentricFCMaps( + maps="Smith_rsn_10", + conn_method="correlation", + conn_method_params=conn_method_params, + ) + # Check correct output + assert "matrix" == marker.get_output_type( + input_type="BOLD", output_feature="functional_connectivity" + ) + + # Fit-transform the data + edge_fc = marker.fit_transform(element_data) + edge_fc_bold = edge_fc["BOLD"]["functional_connectivity"] + + # For 10 ROIs we should get (10 * (10 -1) / 2) edges in the ETS + n_edges = int(10 * (10 - 1) / 2) + assert "data" in edge_fc_bold + assert "row_names" in edge_fc_bold + assert "col_names" in edge_fc_bold + assert edge_fc_bold["data"].shape == (n_edges, n_edges) + assert len(set(edge_fc_bold["row_names"])) == n_edges + assert len(set(edge_fc_bold["col_names"])) == n_edges + + # Store + storage = HDF5FeatureStorage(uri=tmp_path / "test_edge_fc_maps.hdf5") + marker.fit_transform(input=element_data, storage=storage) + features = storage.list_features() + assert any( + x["name"] == "BOLD_EdgeCentricFCMaps_functional_connectivity" + for x in features.values() + ) diff --git a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_maps.py b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_maps.py new file mode 100644 index 000000000..c4f92ccc3 --- /dev/null +++ b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_maps.py @@ -0,0 +1,110 @@ +"""Provide tests for functional connectivity using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from pathlib import Path +from typing import TYPE_CHECKING + +import pytest +from nilearn.connectome import ConnectivityMeasure +from nilearn.maskers import NiftiMapsMasker +from numpy.testing import assert_array_almost_equal +from sklearn.covariance import EmpiricalCovariance, LedoitWolf + +from junifer.data import MapsRegistry +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import FunctionalConnectivityMaps +from junifer.storage import HDF5FeatureStorage + + +if TYPE_CHECKING: + from sklearn.base import BaseEstimator + + +@pytest.mark.parametrize( + "conn_method_params, cov_estimator", + [ + ({"empirical": False}, LedoitWolf(store_precision=False)), + ({"empirical": True}, EmpiricalCovariance(store_precision=False)), + ], +) +def test_FunctionalConnectivityMaps( + tmp_path: Path, + maps_datagrabber: PatternDataladDataGrabber, + conn_method_params: dict[str, bool], + cov_estimator: type["BaseEstimator"], +) -> None: + """Test FunctionalConnectivityParcels. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + conn_method_params : dict + The parametrized parameters to connectivity measure method. + cov_estimator : estimator object + The parametrized covariance estimator. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Setup marker + marker = FunctionalConnectivityMaps( + maps="Smith_rsn_10", + conn_method="correlation", + conn_method_params=conn_method_params, + ) + # Check correct output + assert "matrix" == marker.get_output_type( + input_type="BOLD", output_feature="functional_connectivity" + ) + + # Fit-transform the data + fc = marker.fit_transform(element_data) + fc_bold = fc["BOLD"]["functional_connectivity"] + + assert "data" in fc_bold + assert "row_names" in fc_bold + assert "col_names" in fc_bold + assert fc_bold["data"].shape == (10, 10) + assert len(set(fc_bold["row_names"])) == 10 + assert len(set(fc_bold["col_names"])) == 10 + + # Compare with nilearn + # Load testing maps for the target data + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Extract timeseries + masker = NiftiMapsMasker(maps_img=maps, standardize=False) + extracted_timeseries = masker.fit_transform( + element_data["BOLD"]["data"] + ) + # Compute the connectivity measure + connectivity_measure = ConnectivityMeasure( + cov_estimator=cov_estimator, + kind="correlation", # type: ignore + ).fit_transform([extracted_timeseries])[0] + + # Check that FC are almost equal + assert_array_almost_equal( + connectivity_measure, fc_bold["data"], decimal=3 + ) + + # Store + storage = HDF5FeatureStorage( + uri=tmp_path / "test_fc_maps.hdf5", + ) + marker.fit_transform(input=element_data, storage=storage) + features = storage.list_features() + assert any( + x["name"] + == "BOLD_FunctionalConnectivityMaps_functional_connectivity" + for x in features.values() + ) diff --git a/junifer/markers/maps_aggregation.py b/junifer/markers/maps_aggregation.py new file mode 100644 index 000000000..2d231f67a --- /dev/null +++ b/junifer/markers/maps_aggregation.py @@ -0,0 +1,201 @@ +"""Provide class for maps aggregation.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, ClassVar, Optional, Union + +from nilearn.maskers import NiftiMapsMasker + +from ..api.decorators import register_marker +from ..data import get_data +from ..stats import get_aggfunc_by_name +from ..typing import Dependencies, MarkerInOutMappings +from ..utils import logger, raise_error, warn_with_log +from .base import BaseMarker + + +__all__ = ["MapsAggregation"] + + +@register_marker +class MapsAggregation(BaseMarker): + """Class for maps aggregation. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + time_method : str, optional + The method to use to aggregate the time series over the time points, + after aggregation (only applicable to BOLD data). If None, + it will not operate on the time dimension (default None). + time_method_params : dict, optional + The parameters to pass to the time aggregation method (default None). + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + on : {"T1w", "T2w", "BOLD", "VBM_GM", "VBM_WM", "VBM_CSF", "fALFF", \ + "GCOR", "LCOR"} or list of the options, optional + The data types to apply the marker to. If None, will work on all + available data (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + Raises + ------ + ValueError + If ``time_method`` is specified for non-BOLD data or if + ``time_method_params`` is not None when ``time_method`` is None. + + """ + + _DEPENDENCIES: ClassVar[Dependencies] = {"nilearn", "numpy"} + + _MARKER_INOUT_MAPPINGS: ClassVar[MarkerInOutMappings] = { + "T1w": { + "aggregation": "vector", + }, + "T2w": { + "aggregation": "vector", + }, + "BOLD": { + "aggregation": "timeseries", + }, + "VBM_GM": { + "aggregation": "vector", + }, + "VBM_WM": { + "aggregation": "vector", + }, + "VBM_CSF": { + "aggregation": "vector", + }, + "fALFF": { + "aggregation": "vector", + }, + "GCOR": { + "aggregation": "vector", + }, + "LCOR": { + "aggregation": "vector", + }, + } + + def __init__( + self, + maps: str, + time_method: Optional[str] = None, + time_method_params: Optional[dict[str, Any]] = None, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + on: Union[list[str], str, None] = None, + name: Optional[str] = None, + ) -> None: + self.maps = maps + self.masks = masks + super().__init__(on=on, name=name) + + # Verify after super init so self._on is set + if "BOLD" not in self._on and time_method is not None: + raise_error( + "`time_method` can only be used with BOLD data. " + "Please remove `time_method` parameter." + ) + if time_method is None and time_method_params is not None: + raise_error( + "`time_method_params` can only be used with `time_method`. " + "Please remove `time_method_params` parameter." + ) + self.time_method = time_method + self.time_method_params = time_method_params or {} + + def compute( + self, input: dict[str, Any], extra_input: Optional[dict] = None + ) -> dict: + """Compute. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``aggregation`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + Warns + ----- + RuntimeWarning + If time aggregation is required but only time point is available. + + """ + t_input_img = input["data"] + logger.debug("Maps aggregation") + + # Get maps tailored to target image + maps_img, labels = get_data( + kind="maps", + names=self.maps, + target_data=input, + extra_input=extra_input, + ) + + # Load mask + mask_img = None + if self.masks is not None: + logger.debug(f"Masking with {self.masks}") + # Get tailored mask + mask_img = get_data( + kind="mask", + names=self.masks, + target_data=input, + extra_input=extra_input, + ) + + # Initialize masker + logger.debug("Masking") + masker = NiftiMapsMasker( + maps_img=maps_img, + mask_img=mask_img, + target_affine=t_input_img.affine, + ) + # Mask the input data and extract data + data = masker.fit_transform(t_input_img) + + # Apply time dimension aggregation if required + if self.time_method is not None: + if data.shape[0] > 1: + logger.debug("Aggregating time dimension") + time_agg_func = get_aggfunc_by_name( + self.time_method, func_params=self.time_method_params + ) + data = time_agg_func(data, axis=0) + else: + warn_with_log( + "No time dimension to aggregate as only one time point is " + "available." + ) + # Format the output + return { + "aggregation": { + "data": data, + "col_names": labels, + }, + } diff --git a/junifer/markers/reho/__init__.pyi b/junifer/markers/reho/__init__.pyi index 6f1ad0f7c..82d880733 100644 --- a/junifer/markers/reho/__init__.pyi +++ b/junifer/markers/reho/__init__.pyi @@ -1,4 +1,5 @@ -__all__ = ["ReHoParcels", "ReHoSpheres"] +__all__ = ["ReHoMaps", "ReHoParcels", "ReHoSpheres"] +from .reho_maps import ReHoMaps from .reho_parcels import ReHoParcels from .reho_spheres import ReHoSpheres diff --git a/junifer/markers/reho/reho_maps.py b/junifer/markers/reho/reho_maps.py new file mode 100644 index 000000000..50ed38f7e --- /dev/null +++ b/junifer/markers/reho/reho_maps.py @@ -0,0 +1,161 @@ +"""Provide class for regional homogeneity (ReHo) on maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, Optional, Union + +import numpy as np + +from ...api.decorators import register_marker +from ...utils import logger +from ..maps_aggregation import MapsAggregation +from .reho_base import ReHoBase + + +__all__ = ["ReHoMaps"] + + +@register_marker +class ReHoMaps(ReHoBase): + """Class for regional homogeneity on maps. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + using : {"junifer", "afni"} + Implementation to use for computing ReHo: + + * "junifer" : Use ``junifer``'s own ReHo implementation + * "afni" : Use AFNI's ``3dReHo`` + + reho_params : dict, optional + Extra parameters for computing ReHo map as a dictionary (default None). + If ``using="afni"``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + - 7 : for facewise neighbours only + - 19 : for face- and edge-wise neighbours + - 27 : for face-, edge-, and node-wise neighbors + + * ``neigh_rad`` : positive float, optional + The radius of a desired neighbourhood (default None). + * ``neigh_x`` : positive float, optional + The semi-radius for x-axis of ellipsoidal volumes (default None). + * ``neigh_y`` : positive float, optional + The semi-radius for y-axis of ellipsoidal volumes (default None). + * ``neigh_z`` : positive float, optional + The semi-radius for z-axis of ellipsoidal volumes (default None). + * ``box_rad`` : positive int, optional + The number of voxels outward in a given cardinal direction for a + cubic box centered on a given voxel (default None). + * ``box_x`` : positive int, optional + The number of voxels for +/- x-axis of cuboidal volumes + (default None). + * ``box_y`` : positive int, optional + The number of voxels for +/- y-axis of cuboidal volumes + (default None). + * ``box_z`` : positive int, optional + The number of voxels for +/- z-axis of cuboidal volumes + (default None). + + else if ``using="junifer"``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27, 125}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + * 7 : for facewise neighbours only + * 19 : for face- and edge-wise neighbours + * 27 : for face-, edge-, and node-wise neighbors + * 125 : for 5x5 cuboidal volume + + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + name : str, optional + The name of the marker. If None, it will use the class name + (default None). + + """ + + def __init__( + self, + maps: str, + using: str, + reho_params: Optional[dict] = None, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + name: Optional[str] = None, + ) -> None: + # Superclass init first to validate `using` parameter + super().__init__(using=using, name=name) + self.maps = maps + self.reho_params = reho_params + self.masks = masks + + def compute( + self, + input: dict[str, Any], + extra_input: Optional[dict[str, Any]] = None, + ) -> dict[str, Any]: + """Compute. + + Parameters + ---------- + input : dict + The BOLD data as dictionary. + extra_input : dict, optional + The other fields in the pipeline data object (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``reho`` : dictionary with the following keys: + + - ``data`` : ROI values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + """ + logger.info("Calculating ReHo for maps") + + # Compute voxelwise reho + # If the input data space is "native", then reho_file_path points to + # the input data path as it might be required for parcellation + # transformation to native space. + if self.reho_params is not None: + reho_map, reho_file_path = self._compute( + input_data=input, **self.reho_params + ) + else: + reho_map, reho_file_path = self._compute(input_data=input) + + # Perform aggregation on reho map + aggregation_input = dict(input.items()) + aggregation_input["data"] = reho_map + aggregation_input["path"] = reho_file_path + maps_aggregation = MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute( + input=aggregation_input, + extra_input=extra_input, + ) + + return { + "reho": { + # Only use the first row and expand row dimension + "data": maps_aggregation["aggregation"]["data"][0][ + np.newaxis, : + ], + "col_names": maps_aggregation["aggregation"]["col_names"], + } + } diff --git a/junifer/markers/reho/tests/test_reho_maps.py b/junifer/markers/reho/tests/test_reho_maps.py new file mode 100644 index 000000000..42db2d987 --- /dev/null +++ b/junifer/markers/reho/tests/test_reho_maps.py @@ -0,0 +1,135 @@ +"""Provide tests for ReHoMaps.""" + +# Authors: Synchon Mandal +# License: AGPL + +import logging +from pathlib import Path + +import pytest +import scipy as sp + +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import ReHoMaps +from junifer.pipeline import WorkDirManager +from junifer.pipeline.utils import _check_afni +from junifer.storage import HDF5FeatureStorage + + +def test_ReHoMaps( + caplog: pytest.LogCaptureFixture, + tmp_path: Path, + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test ReHoMaps. + + Parameters + ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path + + with caplog.at_level(logging.DEBUG): + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + + # Initialize marker + marker = ReHoMaps( + maps="Smith_rsn_10", + using="junifer", + ) + # Check correct output + assert "vector" == marker.get_output_type( + input_type="BOLD", output_feature="reho" + ) + + # Fit transform marker on data + output = marker.fit_transform(element_data) + + assert "Creating cache" in caplog.text + + # Get BOLD output + assert "BOLD" in output + output_bold = output["BOLD"]["reho"] + # Assert BOLD output keys + assert "data" in output_bold + assert "col_names" in output_bold + + output_bold_data = output_bold["data"] + # Assert BOLD output data dimension + assert output_bold_data.ndim == 2 + # Assert BOLD output data is normalized + assert (output_bold_data > 0).all() and ( + output_bold_data < 1 + ).all() + + # Reset log capture + caplog.clear() + # Initialize storage + storage = HDF5FeatureStorage(tmp_path / "reho_maps.hdf5") + # Fit transform marker on data with storage + marker.fit_transform( + input=element_data, + storage=storage, + ) + # Check stored feature name + features = storage.list_features() + assert any( + x["name"] == "BOLD_ReHoMaps_reho" for x in features.values() + ) + # Cache working correctly + assert "Creating cache" not in caplog.text + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires AFNI to be in PATH" +) +def test_ReHoMaps_comparison( + tmp_path: Path, maps_datagrabber: PatternDataladDataGrabber +) -> None: + """Test ReHoMaps implementation comparison. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path + + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + + # Initialize marker + junifer_marker = ReHoMaps(maps="Smith_rsn_10", using="junifer") + # Fit transform marker on data + junifer_output = junifer_marker.fit_transform(element_data) + # Get BOLD output + junifer_output_bold = junifer_output["BOLD"]["reho"] + + # Initialize marker + afni_marker = ReHoMaps(maps="Smith_rsn_10", using="afni") + # Fit transform marker on data + afni_output = afni_marker.fit_transform(element_data) + # Get BOLD output + afni_output_bold = afni_output["BOLD"]["reho"] + + # Check for Pearson correlation coefficient + r, _ = sp.stats.pearsonr( + junifer_output_bold["data"].flatten(), + afni_output_bold["data"].flatten(), + ) + assert r >= 0.2 # this is very bad, but they differ... diff --git a/junifer/markers/temporal_snr/__init__.pyi b/junifer/markers/temporal_snr/__init__.pyi index 74022fbd1..4b542648e 100644 --- a/junifer/markers/temporal_snr/__init__.pyi +++ b/junifer/markers/temporal_snr/__init__.pyi @@ -1,4 +1,5 @@ -__all__ = ["TemporalSNRParcels", "TemporalSNRSpheres"] +__all__ = ["TemporalSNRMaps", "TemporalSNRParcels", "TemporalSNRSpheres"] +from .temporal_snr_maps import TemporalSNRMaps from .temporal_snr_parcels import TemporalSNRParcels from .temporal_snr_spheres import TemporalSNRSpheres diff --git a/junifer/markers/temporal_snr/temporal_snr_maps.py b/junifer/markers/temporal_snr/temporal_snr_maps.py new file mode 100644 index 000000000..e5408c50e --- /dev/null +++ b/junifer/markers/temporal_snr/temporal_snr_maps.py @@ -0,0 +1,80 @@ +"""Provide class for temporal SNR using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Any, Optional, Union + +from ...api.decorators import register_marker +from ..maps_aggregation import MapsAggregation +from .temporal_snr_base import TemporalSNRBase + + +__all__ = ["TemporalSNRMaps"] + + +@register_marker +class TemporalSNRMaps(TemporalSNRBase): + """Class for temporal signal-to-noise ratio using maps. + + Parameters + ---------- + maps : str + The name of the map(s) to use. + See :func:`.list_data` for options. + masks : str, dict or list of dict or str, optional + The specification of the masks to apply to regions before extracting + signals. Check :ref:`Using Masks ` for more details. + If None, will not apply any mask (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + """ + + def __init__( + self, + maps: str, + masks: Union[str, dict, list[Union[dict, str]], None] = None, + name: Optional[str] = None, + ) -> None: + self.maps = maps + super().__init__( + masks=masks, + name=name, + ) + + def aggregate( + self, input: dict[str, Any], extra_input: Optional[dict] = None + ) -> dict: + """Perform maps aggregation. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which the data + is the voxelwise temporal SNR map. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + dict + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``aggregation`` : dictionary with the following keys: + + - ``data`` : ROI-wise tSNR values as ``numpy.ndarray`` + - ``col_names`` : ROI labels as list of str + + """ + return MapsAggregation( + maps=self.maps, + masks=self.masks, + on="BOLD", + ).compute(input=input, extra_input=extra_input) diff --git a/junifer/markers/temporal_snr/tests/test_temporal_snr_maps.py b/junifer/markers/temporal_snr/tests/test_temporal_snr_maps.py new file mode 100644 index 000000000..173509fa2 --- /dev/null +++ b/junifer/markers/temporal_snr/tests/test_temporal_snr_maps.py @@ -0,0 +1,67 @@ +"""Provide tests for temporal signal-to-noise ratio using maps.""" + +# Authors: Synchon Mandal +# License: AGPL + +from pathlib import Path + +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import TemporalSNRMaps +from junifer.storage import HDF5FeatureStorage + + +def test_TemporalSNRMaps_computation( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test TemporalSNRMaps fit-transform. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + marker = TemporalSNRMaps(maps="Smith_rsn_10") + # Check correct output + assert "vector" == marker.get_output_type( + input_type="BOLD", output_feature="tsnr" + ) + + # Fit-transform the data + tsnr_parcels = marker.fit_transform(element_data) + tsnr_parcels_bold = tsnr_parcels["BOLD"]["tsnr"] + + assert "data" in tsnr_parcels_bold + assert "col_names" in tsnr_parcels_bold + assert tsnr_parcels_bold["data"].shape == (1, 10) + assert len(set(tsnr_parcels_bold["col_names"])) == 10 + + +def test_TemporalSNRMaps_storage( + tmp_path: Path, maps_datagrabber: PatternDataladDataGrabber +) -> None: + """Test TemporalSNRMaps store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + marker = TemporalSNRMaps(maps="Smith_rsn_10") + # Store + storage = HDF5FeatureStorage(tmp_path / "test_tsnr_maps.hdf5") + marker.fit_transform(input=element_data, storage=storage) + features = storage.list_features() + assert any( + x["name"] == "BOLD_TemporalSNRMaps_tsnr" for x in features.values() + ) diff --git a/junifer/markers/tests/test_maps_aggregation.py b/junifer/markers/tests/test_maps_aggregation.py new file mode 100644 index 000000000..9c2b3e8bb --- /dev/null +++ b/junifer/markers/tests/test_maps_aggregation.py @@ -0,0 +1,430 @@ +"""Provide test for maps aggregation.""" + +# Authors: Synchon Mandal +# License: AGPL + +from copy import deepcopy +from pathlib import Path + +import pytest +from nilearn.maskers import NiftiMapsMasker +from numpy.testing import assert_array_almost_equal, assert_array_equal + +from junifer.data import MapsRegistry, MaskRegistry +from junifer.data.masks import compute_brain_mask +from junifer.datagrabber import PatternDataladDataGrabber +from junifer.datareader import DefaultDataReader +from junifer.markers import MapsAggregation +from junifer.storage import HDF5FeatureStorage + + +@pytest.mark.parametrize( + "input_type, storage_type", + [ + ( + "T1w", + "vector", + ), + ( + "T2w", + "vector", + ), + ( + "BOLD", + "timeseries", + ), + ( + "VBM_GM", + "vector", + ), + ( + "VBM_WM", + "vector", + ), + ( + "VBM_CSF", + "vector", + ), + ( + "fALFF", + "vector", + ), + ( + "GCOR", + "vector", + ), + ( + "LCOR", + "vector", + ), + ], +) +def test_MapsAggregation_input_output( + input_type: str, storage_type: str +) -> None: + """Test MapsAggregation input and output types. + + Parameters + ---------- + input_type : str + The parametrized input type. + storage_type : str + The parametrized storage type. + + """ + assert storage_type == MapsAggregation( + maps="Smith_rsn_10", on=input_type + ).get_output_type(input_type=input_type, output_feature="aggregation") + + +def test_MapsAggregation_3D( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test MapsAggregation on 3D data. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Deepcopy data for later use + input_data = deepcopy(element_data) + # Modify for junifer + element_data["BOLD"]["data"] = element_data["BOLD"]["data"].slicer[ + ..., + 0:1, + ] + marker = MapsAggregation(maps="Smith_rsn_10") + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + assert maps_agg_bold_data.ndim == 2 + + # Compare with nilearn + # Load testing maps + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Extract data + masker = NiftiMapsMasker(maps_img=maps) + nifti_maps_masked_bold = masker.fit_transform( + input_data["BOLD"]["data"].slicer[..., 0:1] + ) + + assert_array_equal( + nifti_maps_masked_bold.shape, maps_agg_bold_data.shape + ) + assert_array_equal(nifti_maps_masked_bold, maps_agg_bold_data) + + +def test_MapsAggregation_4D( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test MapsAggregation on 4D data. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Deepcopy data for later use + input_data = deepcopy(element_data) + marker = MapsAggregation(maps="Smith_rsn_10") + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + assert maps_agg_bold_data.ndim == 2 + + # Compare with nilearn + # Load testing maps + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Extract data + masker = NiftiMapsMasker(maps_img=maps) + nifti_maps_masked_bold = masker.fit_transform( + input_data["BOLD"]["data"] + ) + + assert_array_equal( + nifti_maps_masked_bold.shape, maps_agg_bold_data.shape + ) + assert_array_equal(nifti_maps_masked_bold, maps_agg_bold_data) + + +def test_MapsAggregation_storage( + maps_datagrabber: PatternDataladDataGrabber, tmp_path: Path +) -> None: + """Test MapsAggregation storage. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + tmp_path : pathlib.Path + The path to the test directory. + + """ + # Store 3D + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + storage = HDF5FeatureStorage( + uri=tmp_path / "test_maps_storage_3D.hdf5" + ) + marker = MapsAggregation( + maps="Smith_rsn_10", + on="BOLD", + ) + element_data["BOLD"]["data"] = element_data["BOLD"]["data"].slicer[ + ..., 0:1 + ] + marker.fit_transform(input=element_data, storage=storage) + features = storage.list_features() + assert any( + x["name"] == "BOLD_MapsAggregation_aggregation" + for x in features.values() + ) + + # Store 4D + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + storage = HDF5FeatureStorage( + uri=tmp_path / "test_maps_storage_4D.sqlite" + ) + marker = MapsAggregation( + maps="Smith_rsn_10", + on="BOLD", + ) + marker.fit_transform(input=element_data, storage=storage) + features = storage.list_features() + assert any( + x["name"] == "BOLD_MapsAggregation_aggregation" + for x in features.values() + ) + + +def test_MapsAggregation_3D_mask( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test MapsAggregation on 3D data with mask. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Deepcopy data for later use + input_data = deepcopy(element_data) + # Modify for junifer + element_data["BOLD"]["data"] = element_data["BOLD"]["data"].slicer[ + ..., + 0:1, + ] + marker = MapsAggregation( + maps="Smith_rsn_10", + on="BOLD", + masks="compute_brain_mask", + ) + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + assert maps_agg_bold_data.ndim == 2 + + # Compare with nilearn + # Load testing maps + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Load mask + mask_img = MaskRegistry().get( + "compute_brain_mask", target_data=element_data["BOLD"] + ) + # Extract data + masker = NiftiMapsMasker(maps_img=maps, mask_img=mask_img) + nifti_maps_masked_bold = masker.fit_transform( + input_data["BOLD"]["data"].slicer[..., 0:1] + ) + + assert_array_equal( + nifti_maps_masked_bold.shape, maps_agg_bold_data.shape + ) + assert_array_equal(nifti_maps_masked_bold, maps_agg_bold_data) + + +def test_MapsAggregation_3D_mask_computed( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test MapsAggregation on 3D data with computed masks. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + + # Compare with nilearn + # Load testing maps + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Get a mask + mask_img = compute_brain_mask(element_data["BOLD"], threshold=0.2) + # Create NiftiMapsMasker + masker = NiftiMapsMasker(maps_img=maps, mask_img=mask_img) + nifti_maps_masked_bold_good = masker.fit_transform( + element_data["BOLD"]["data"] + ) + + # Get another mask + mask_img = compute_brain_mask(element_data["BOLD"], threshold=0.5) + # Create NiftiMapsMasker + masker = NiftiMapsMasker(maps_img=maps, mask_img=mask_img) + nifti_maps_masked_bold_bad = masker.fit_transform(mask_img) + + # Use the MapsAggregation object + marker = MapsAggregation( + maps="Smith_rsn_10", + masks={"compute_brain_mask": {"threshold": 0.2}}, + on="BOLD", + ) + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + + assert maps_agg_bold_data.ndim == 2 + assert_array_almost_equal( + nifti_maps_masked_bold_good, maps_agg_bold_data + ) + + with pytest.raises(AssertionError): + assert_array_almost_equal( + maps_agg_bold_data, nifti_maps_masked_bold_bad + ) + + +def test_MapsAggregation_4D_agg_time( + maps_datagrabber: PatternDataladDataGrabber, +): + """Test MapsAggregation on 4D data, aggregating time. + + Parameters + ---------- + maps_datagrabber : PatternDataladDataGrabber + The testing PatternDataladDataGrabber, as fixture. + + """ + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + # Create MapsAggregation object + marker = MapsAggregation( + maps="Smith_rsn_10", + time_method="mean", + on="BOLD", + ) + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + + # Compare with nilearn + # Loading testing maps + maps, _ = MapsRegistry().get( + maps="Smith_rsn_10", + target_data=element_data["BOLD"], + ) + # Extract data + masker = NiftiMapsMasker(maps_img=maps) + nifti_maps_masked_bold = masker.fit_transform( + element_data["BOLD"]["data"] + ) + nifti_maps_masked_bold_mean = nifti_maps_masked_bold.mean(axis=0) + + assert maps_agg_bold_data.ndim == 1 + assert_array_equal( + nifti_maps_masked_bold_mean.shape, maps_agg_bold_data.shape + ) + assert_array_almost_equal( + nifti_maps_masked_bold_mean, maps_agg_bold_data, decimal=2 + ) + + # Test picking first time point + nifti_maps_masked_bold_pick_0 = nifti_maps_masked_bold[:1, :] + marker = MapsAggregation( + maps="Smith_rsn_10", + time_method="select", + time_method_params={"pick": [0]}, + on="BOLD", + ) + maps_agg_bold_data = marker.fit_transform(element_data)["BOLD"][ + "aggregation" + ]["data"] + + assert maps_agg_bold_data.ndim == 2 + assert_array_equal( + nifti_maps_masked_bold_pick_0.shape, maps_agg_bold_data.shape + ) + assert_array_equal(nifti_maps_masked_bold_pick_0, maps_agg_bold_data) + + +def test_MapsAggregation_errors() -> None: + """Test errors for MapsAggregation.""" + with pytest.raises(ValueError, match="can only be used with BOLD data"): + MapsAggregation( + maps="Smith_rsn_10", + time_method="select", + time_method_params={"pick": [0]}, + on="VBM_GM", + ) + + with pytest.raises( + ValueError, match="can only be used with `time_method`" + ): + MapsAggregation( + maps="Smith_rsn_10", + time_method_params={"pick": [0]}, + on="VBM_GM", + ) + + +def test_MapsAggregation_warning( + maps_datagrabber: PatternDataladDataGrabber, +) -> None: + """Test warning for MapsAggregation.""" + with maps_datagrabber as dg: + element = dg[("sub-01", "sub-001", "rest", "1")] + element_data = DefaultDataReader().fit_transform(element) + with pytest.warns( + RuntimeWarning, match="No time dimension to aggregate" + ): + marker = MapsAggregation( + maps="Smith_rsn_10", + time_method="select", + time_method_params={"pick": [0]}, + on="BOLD", + ) + element_data["BOLD"]["data"] = element_data["BOLD"]["data"].slicer[ + ..., 0:1 + ] + marker.fit_transform(element_data) diff --git a/junifer/typing/_typing.py b/junifer/typing/_typing.py index dbdc40226..df5fd646a 100644 --- a/junifer/typing/_typing.py +++ b/junifer/typing/_typing.py @@ -65,9 +65,7 @@ ConditionalDependencies = Sequence[ ] ExternalDependencies = Sequence[MutableMapping[str, Union[str, Sequence[str]]]] MarkerInOutMappings = MutableMapping[str, MutableMapping[str, str]] -DataGrabberPatterns = dict[ - str, Union[dict[str, str], list[dict[str, str]]] -] +DataGrabberPatterns = dict[str, Union[dict[str, str], list[dict[str, str]]]] ConfigVal = Union[bool, int, float, str] Element = Union[str, tuple[str, ...]] Elements = Sequence[Element]