diff --git a/docs/changes/latest.inc b/docs/changes/latest.inc index 56a8cfe82..007443848 100644 --- a/docs/changes/latest.inc +++ b/docs/changes/latest.inc @@ -59,6 +59,7 @@ Enhancements - Modify the behaviour of the ``collect`` parameter in HTCondor ``queue`` function to run a collect job even if some of the previous jobs fail. This is useful to collect the results of a pipeline even if some of the jobs fail (:gh:`190` by `Fede Raimondo`_). +- Allow for empty parcels in :class:`junifer.markers.ParcelAggregation`, that will result in NaNs (:gh:`194` by `Fede Raimondo`_). Bugs ~~~~ @@ -76,6 +77,8 @@ Bugs - Fix an issue with datalad cache and locks in which the overriden settings in Junifer were not propagated to subprocesses, resulting in using the default settings (:gh:`199` by `Fede Raimondo`_). +- Fix a bug in which :class:`junifer.markers.ParcelAggregation` could yield duplicated column names if two or more parcels were used and label names were not unique (:gh:`194` by `Fede Raimondo`_). + API changes ~~~~~~~~~~~ diff --git a/junifer/markers/parcel_aggregation.py b/junifer/markers/parcel_aggregation.py index f7c9c3774..075464125 100644 --- a/junifer/markers/parcel_aggregation.py +++ b/junifer/markers/parcel_aggregation.py @@ -13,7 +13,7 @@ from nilearn.maskers import NiftiMasker from ..api.decorators import register_marker from ..data import get_mask, load_parcellation from ..stats import get_aggfunc_by_name -from ..utils import logger +from ..utils import logger, warn_with_log from .base import BaseMarker @@ -159,6 +159,22 @@ class ParcelAggregation(BaseMarker): labels = all_labels[0] else: # Merge the parcellations + + # Check for duplicated labels + all_labels_flat = [ + item for sublist in all_labels for item in sublist + ] + if len(all_labels_flat) != len(set(all_labels_flat)): + warn_with_log( + "The parcellations have duplicated labels. " + "Each label will be prefixed with the parcellation name." + ) + for i_parcellation, t_labels in enumerate(all_labels): + all_labels[i_parcellation] = [ + f"{self.parcellation[i_parcellation]}_{t_label}" + for t_label in t_labels + ] + overlapping_voxels = False parc_data = all_parcelations[0].get_fdata() labels = all_labels[0] for t_parc, t_labels in zip(all_parcelations[1:], all_labels[1:]): @@ -171,9 +187,19 @@ class ParcelAggregation(BaseMarker): # This makes sure that the voxels that are in multiple # parcellations are assigned to the parcellation that was # first in the list. + if np.any(parc_data[t_parc_data != 0] != 0): + overlapping_voxels = True + parc_data[parc_data == 0] += t_parc_data[parc_data == 0] labels.extend(t_labels) + if overlapping_voxels: + warn_with_log( + "The parcellations have overlapping voxels. " + "The overlapping voxels will be assigned to the " + "parcellation that was first in the list." + ) + parcellation_img_res = new_img_like( all_parcelations[0], parc_data, @@ -208,17 +234,14 @@ class ParcelAggregation(BaseMarker): # Get the values for each parcel and apply agg function logger.debug("Computing ROI means") - parcellation_roi_vals = sorted(np.unique(parcellation_values)) - out_labels = [] out_values = [] # Iterate over the parcels (existing) - for t_v in parcellation_roi_vals: + for t_v in range(1, len(labels) + 1): t_values = agg_func(data[:, parcellation_values == t_v], axis=-1) out_values.append(t_values) # Update the labels just in case a parcel has no voxels # in it - out_labels.append(labels[t_v - 1]) out_values = np.array(out_values).T - out = {"data": out_values, "col_names": out_labels} + out = {"data": out_values, "col_names": labels} return out diff --git a/junifer/markers/tests/test_parcel_aggregation.py b/junifer/markers/tests/test_parcel_aggregation.py index 93f878714..e54bba8ef 100644 --- a/junifer/markers/tests/test_parcel_aggregation.py +++ b/junifer/markers/tests/test_parcel_aggregation.py @@ -4,6 +4,7 @@ # Synchon Mandal # License: AGPL +import warnings from pathlib import Path import nibabel as nib @@ -328,8 +329,12 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: nib.save(parcellation1_img, parcellation1_path) nib.save(parcellation2_img, parcellation2_path) - register_parcellation("Schaefer100x7_low", parcellation1_path, labels1) - register_parcellation("Schaefer100x7_high", parcellation2_path, labels2) + register_parcellation( + "Schaefer100x7_low", parcellation1_path, labels1, overwrite=True + ) + register_parcellation( + "Schaefer100x7_high", parcellation2_path, labels2, overwrite=True + ) # Use the ParcelAggregation object on the original parcellation marker_original = ParcelAggregation( @@ -355,7 +360,11 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: on="VBM_GM", ) # Test passing "on" as a keyword argument input = {"VBM_GM": {"data": img, "meta": {}}} - split_mean = marker_split.fit_transform(input)["VBM_GM"] + + # No warnings should be raised + with warnings.catch_warnings(): + warnings.simplefilter("error", category=UserWarning) + split_mean = marker_split.fit_transform(input)["VBM_GM"] split_mean_data = split_mean["data"] assert split_mean_data.ndim == 2 @@ -408,8 +417,12 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: nib.save(parcellation1_img, parcellation1_path) nib.save(parcellation2_img, parcellation2_path) - register_parcellation("Schaefer100x7_low2", parcellation1_path, labels1) - register_parcellation("Schaefer100x7_high2", parcellation2_path, labels2) + register_parcellation( + "Schaefer100x7_low2", parcellation1_path, labels1, overwrite=True + ) + register_parcellation( + "Schaefer100x7_high2", parcellation2_path, labels2, overwrite=True + ) # Use the ParcelAggregation object on the original parcellation marker_original = ParcelAggregation( @@ -435,7 +448,101 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: on="VBM_GM", ) # Test passing "on" as a keyword argument input = {"VBM_GM": {"data": img, "meta": {}}} - split_mean = marker_split.fit_transform(input)["VBM_GM"] + with pytest.warns(RuntimeWarning, match="overlapping voxels"): + split_mean = marker_split.fit_transform(input)["VBM_GM"] + split_mean_data = split_mean["data"] + + assert split_mean_data.ndim == 2 + assert split_mean_data.shape[0] == 1 + assert split_mean_data.shape[1] == 105 + + # Overlapping voxels should be NaN + assert np.isnan(split_mean_data[:, 50:55]).all() + + non_nan = split_mean_data[~np.isnan(split_mean_data)] + # Data should be the same + assert_array_equal(orig_mean_data, non_nan[None, :]) + + # Labels should be "low" for the first 50 and "high" for the second 50 + assert all(x.startswith("low") for x in split_mean["col_names"][:50]) + assert all(x.startswith("high") for x in split_mean["col_names"][50:]) + + +def test_ParcelAggregation_3D_multiple_duplicated_labels( + tmp_path: Path, +) -> None: + """Test ParcelAggregation with two parcellations with duplicated labels. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + + # Get the testing parcellation + parcellation, labels, _ = load_parcellation("Schaefer100x7") + + assert parcellation is not None + + # Get the oasis VBM data + oasis_dataset = datasets.fetch_oasis_vbm(n_subjects=1) + vbm = oasis_dataset.gray_matter_maps[0] + img = nib.load(vbm) + + # Create two parcellations from it + parcellation_data = parcellation.get_fdata() + parcellation1_data = parcellation_data.copy() + parcellation1_data[parcellation1_data > 50] = 0 + parcellation2_data = parcellation_data.copy() + parcellation2_data[parcellation2_data <= 50] = 0 + parcellation2_data[parcellation2_data > 0] -= 50 + labels1 = labels[:50] + labels2 = labels[49:-1] # One label is duplicated + + parcellation1_img = new_img_like(parcellation, parcellation1_data) + parcellation2_img = new_img_like(parcellation, parcellation2_data) + + parcellation1_path = tmp_path / "parcellation1.nii.gz" + parcellation2_path = tmp_path / "parcellation2.nii.gz" + + nib.save(parcellation1_img, parcellation1_path) + nib.save(parcellation2_img, parcellation2_path) + + register_parcellation( + "Schaefer100x7_low", parcellation1_path, labels1, overwrite=True + ) + register_parcellation( + "Schaefer100x7_high", parcellation2_path, labels2, overwrite=True + ) + + # Use the ParcelAggregation object on the original parcellation + marker_original = ParcelAggregation( + parcellation="Schaefer100x7", + method="mean", + name="gmd_schaefer100x7_mean", + on="VBM_GM", + ) # Test passing "on" as a keyword argument + input = {"VBM_GM": {"data": img, "meta": {}}} + orig_mean = marker_original.fit_transform(input)["VBM_GM"] + + orig_mean_data = orig_mean["data"] + assert orig_mean_data.ndim == 2 + assert orig_mean_data.shape[0] == 1 + assert orig_mean_data.shape[1] == 100 + # assert_array_almost_equal(auto, jun_values3d_mean) + + # Use the ParcelAggregation object on the two parcellations + marker_split = ParcelAggregation( + parcellation=["Schaefer100x7_low", "Schaefer100x7_high"], + method="mean", + name="gmd_schaefer100x7_mean", + on="VBM_GM", + ) # Test passing "on" as a keyword argument + input = {"VBM_GM": {"data": img, "meta": {}}} + + with pytest.warns(RuntimeWarning, match="duplicated labels."): + split_mean = marker_split.fit_transform(input)["VBM_GM"] split_mean_data = split_mean["data"] assert split_mean_data.ndim == 2 @@ -445,6 +552,7 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: # Data should be the same assert_array_equal(orig_mean_data, split_mean_data) - # Labels should be "low" for the first 50 and "high" for the second 50 - assert all(x.startswith("low") for x in split_mean["col_names"][:50]) - assert all(x.startswith("high") for x in split_mean["col_names"][50:]) + # Labels should be prefixed with the parcellation name + col_names = [f"Schaefer100x7_low_{x}" for x in labels1] + col_names += [f"Schaefer100x7_high_{x}" for x in labels2] + assert col_names == split_mean["col_names"]