diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index 01b0ce5f2..8b9c0985f 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -62,7 +62,7 @@ jobs: cache-to: type=gha,mode=max tags: wonky:conn load: true - # always use light mode for testing container since the full version was tested during build + # Always use light mode for testing container since the full version was tested during build - run: | datalad get data/halfpipe data/atlases docker run --rm \ diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 606017587..8ec1b24ab 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -65,7 +65,7 @@ jobs: environments: test frozen: true - run: pixi run unittest - - run: pixi run smoketestlight + - run: pixi run smoketestlight --cov-append - uses: codecov/codecov-action@v6 if: ${{ always() }} with: diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 3a7f90b2f..93eead49e 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -20,7 +20,13 @@ repos: rev: v1.19.1 hooks: - id: mypy - additional_dependencies: [pandas-stubs, types-tqdm, types-setuptools, types-Jinja2] + additional_dependencies: + - pandas-stubs + - pytest + - types-tqdm + - types-setuptools + - types-Jinja2 + - textual args: [--config-file=pyproject.toml] - repo: https://github.com/codespell-project/codespell rev: v2.4.2 diff --git a/codecov.yml b/codecov.yml index 7b9ef93df..e4af71446 100644 --- a/codecov.yml +++ b/codecov.yml @@ -14,13 +14,16 @@ # # noisy red coverage status on github PRs. # target: auto # threshold: 1% -comment: # this is a top-level key +comment: # this is a top-level key layout: reach, diff, flags, files behavior: default - require_changes: false # if true: only post the comment if coverage changes - require_base: no # [yes :: must have a base report to post] - require_head: yes # [yes :: must have a head report to post] + require_changes: false # if true: only post the comment if coverage changes + require_base: false # if true: must have a base report to post + require_head: false # if true: must have a head report to post ignore: -- '*/tests/' # ignore folders related to testing -- '*/data/' + - "**/tests/**" + - "**/test_*.py" + - "**/conftest.py" + - "**/data/**" + - "**/setup.py" diff --git a/data/halfpipe/participants.tsv b/data/halfpipe/participants.tsv index 9d371793b..9f8d33065 120000 --- a/data/halfpipe/participants.tsv +++ b/data/halfpipe/participants.tsv @@ -1 +1 @@ -../../.git/annex/objects/gG/5m/MD5E-s16387--b0d7908c25384a3ecd7c317238748596.tsv/MD5E-s16387--b0d7908c25384a3ecd7c317238748596.tsv \ No newline at end of file +../../.git/annex/objects/Fj/K7/MD5E-s16372--3c23f937d8662e3bbd2f588f3d992e8b.tsv/MD5E-s16372--3c23f937d8662e3bbd2f588f3d992e8b.tsv \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index f162fb4c4..fc1dd056d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -164,13 +164,18 @@ warn_unused_ignores = true ignore_missing_imports = true module = [ "bids.*", + "brainspace.*", "datalad.*", + "joblib", "matplotlib.*", + "nilearn", + "nilearn.*", "numba.*", "patsy.*", "rich.*", "scipy.*", "seaborn.*", + "sklearn.*", "statsmodels.*", "templateflow.*", "nibabel.*", diff --git a/wonkyconn/atlas.py b/wonkyconn/atlas.py index 3c28807f3..6131d4a10 100644 --- a/wonkyconn/atlas.py +++ b/wonkyconn/atlas.py @@ -7,8 +7,8 @@ import numpy as np import pandas as pd import scipy -from nilearn.image import iter_img, load_img, math_img, resample_to_img # type: ignore[import-not-found] -from nilearn.maskers import NiftiLabelsMasker, NiftiMasker # type: ignore[import-not-found] +from nilearn.image import iter_img, load_img, math_img, resample_to_img +from nilearn.maskers import NiftiLabelsMasker, NiftiMasker from numpy import typing as npt from .logger import logger @@ -29,7 +29,7 @@ class Atlas(ABC): """ seg: str - image: nib.nifti1.Nifti1Image + image: nib.nifti1.Nifti1Image # pyright: ignore[reportAttributeAccessIssue] structure: npt.NDArray[np.bool_] = field(default_factory=lambda: np.ones((3, 3, 3), dtype=bool)) @@ -51,7 +51,7 @@ def get_centroids(self) -> npt.NDArray[np.float64]: npt.NDArray[np.float64]: An array of centroid coordinates. """ centroid_points = self.get_centroid_points() - centroid_coordinates = nib.affines.apply_affine(self.image.affine, centroid_points) + centroid_coordinates = nib.affines.apply_affine(self.image.affine, centroid_points) # pyright: ignore[reportAttributeAccessIssue] return centroid_coordinates def get_distance_matrix(self) -> npt.NDArray[np.float64]: @@ -65,7 +65,7 @@ def get_distance_matrix(self) -> npt.NDArray[np.float64]: centroids = self.get_centroids() return scipy.spatial.distance.squareform(scipy.spatial.distance.pdist(centroids)) - def load_yeo7_network(self) -> nib.nifti1.Nifti1Image: + def load_yeo7_network(self) -> nib.nifti1.Nifti1Image: # pyright: ignore[reportAttributeAccessIssue] """ Load and resample the yeo 7 networks to the atlas's space. @@ -102,10 +102,10 @@ def create(seg: str, path: Path) -> "Atlas": None """ - image = nib.nifti1.load(path) + image = nib.nifti1.load(path) # pyright: ignore[reportAttributeAccessIssue] if image.ndim <= 3 or image.shape[3] == 1: - return DsegAtlas(seg, nib.funcs.squeeze_image(image)) + return DsegAtlas(seg, nib.funcs.squeeze_image(image)) # pyright: ignore[reportAttributeAccessIssue] else: return ProbsegAtlas(seg, image) @@ -142,7 +142,8 @@ def get_yeo7_membership(self) -> pd.DataFrame: for n in network_labels: cur_region = math_img(f"img=={n}", img=yeo7_nii) masker = NiftiMasker(cur_region) - atlas_parcel_in_network = masker.fit_transform(self.image) + masker.fit() + atlas_parcel_in_network = masker.transform(self.image) atlas_parcel_in_network = np.unique(atlas_parcel_in_network)[1:] region_membership.loc[atlas_parcel_in_network, f"yeo7-{int(n)}"] = 1 return region_membership @@ -161,7 +162,7 @@ def _get_centroid_point(self, i: int, array: npt.NDArray[np.float64]) -> tuple[f def get_centroid_points(self) -> npt.NDArray[np.float64]: return np.asarray( - [self._get_centroid_point(i, image.get_fdata()) for i, image in enumerate(nib.funcs.four_to_three(self.image))] + [self._get_centroid_point(i, image.get_fdata()) for i, image in enumerate(nib.funcs.four_to_three(self.image))] # pyright: ignore[reportAttributeAccessIssue] ) def get_yeo7_membership(self) -> pd.DataFrame: diff --git a/wonkyconn/config.py b/wonkyconn/config.py index b3bd4b535..54ee2d198 100644 --- a/wonkyconn/config.py +++ b/wonkyconn/config.py @@ -3,7 +3,7 @@ import argparse from dataclasses import dataclass, field from pathlib import Path -from typing import Iterable, Sequence +from typing import Literal, TypeAlias def _coerce_path(value: str | Path | None) -> Path | None: @@ -13,8 +13,14 @@ def _coerce_path(value: str | Path | None) -> Path | None: return Path(value).expanduser().resolve() +Metric: TypeAlias = Literal["motion", "analytic-insights", "gradients", "prediction"] + +light_mode_metrics: set[Metric] = {"motion", "analytic-insights"} +all_metrics: set[Metric] = {"motion", "analytic-insights", "gradients", "prediction"} + + @dataclass -class WonkyConnConfig: +class WonkyconnConfig: """Shared configuration for CLI and GUI.""" bids_dir: Path | None = None @@ -22,57 +28,38 @@ class WonkyConnConfig: analysis_level: str = "group" phenotypes: Path | None = None atlas: list[tuple[str, Path]] = field(default_factory=list) - verbosity: int = 2 + log_level: str | None = None debug: bool = False - light_mode: bool = False + metrics: set[Metric] | None = None theme: str | None = None # GUI-only suppress_warnings: bool = False + site_correction: bool = False + + @property + def light_mode(self) -> bool: + return self.metrics == light_mode_metrics @classmethod - def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyConnConfig": + def from_cli_args(cls, args: argparse.Namespace | None) -> "WonkyconnConfig": """Create a config from argparse args (may be partial when GUI is requested).""" if args is None: return cls() - verbosity = args.verbosity - if isinstance(verbosity, Sequence) and not isinstance(verbosity, (str, bytes)): - verbosity = verbosity[0] - atlas_entries: list[tuple[str, Path]] = list() - for label, atlas_path in getattr(args, "atlas", []) or []: + for label, atlas_path in args.atlas or []: atlas_entries.append((label, Path(atlas_path).expanduser().resolve())) + metrics: set[Metric] = light_mode_metrics if args.light_mode else set(args.metrics) + return cls( - bids_dir=_coerce_path(getattr(args, "bids_dir", None)), - output_dir=_coerce_path(getattr(args, "output_dir", None)), - analysis_level=getattr(args, "analysis_level", "group"), - phenotypes=_coerce_path(getattr(args, "phenotypes", None)), + bids_dir=_coerce_path(args.bids_dir), + output_dir=_coerce_path(args.output_dir), + analysis_level=args.analysis_level, + phenotypes=_coerce_path(args.phenotypes), atlas=atlas_entries, - verbosity=int(verbosity) if verbosity is not None else 2, - debug=bool(getattr(args, "debug", False)), - light_mode=bool(getattr(args, "light_mode", False)), - suppress_warnings=bool(getattr(args, "suppress_warnings", False)), - ) - - def to_namespace(self) -> argparse.Namespace: - """Convert to argparse.Namespace expected by workflow.""" - if self.bids_dir is None: - raise ValueError("bids_dir is required") - if self.output_dir is None: - raise ValueError("output_dir is required") - if self.phenotypes is None: - raise ValueError("phenotypes is required") - if not self.atlas: - raise ValueError("At least one atlas entry is required") - - atlas_as_str: Iterable[tuple[str, str]] = ((label, str(path)) for label, path in self.atlas) - return argparse.Namespace( - bids_dir=self.bids_dir, - output_dir=self.output_dir, - analysis_level=self.analysis_level, - phenotypes=str(self.phenotypes), - atlas=list(atlas_as_str), - verbosity=self.verbosity, - debug=self.debug, - light_mode=self.light_mode, + log_level=args.log_level, + debug=bool(args.debug), + metrics=metrics, + suppress_warnings=bool(args.suppress_warnings), + site_correction=bool(args.site_correction), ) diff --git a/wonkyconn/features/age_sex_prediction.py b/wonkyconn/features/age_sex_prediction.py index d2d82420c..7447f0a92 100644 --- a/wonkyconn/features/age_sex_prediction.py +++ b/wonkyconn/features/age_sex_prediction.py @@ -1,23 +1,84 @@ from __future__ import annotations +from dataclasses import dataclass, field from typing import TYPE_CHECKING, Dict, List import numpy as np import pandas as pd -from joblib import parallel_backend # type: ignore[import-not-found] -from nilearn.connectome import sym_matrix_to_vec # type: ignore[import-not-found] +from joblib import parallel_backend +from nilearn.connectome import sym_matrix_to_vec from numpy.typing import NDArray -from sklearn.decomposition import PCA # type: ignore[import-not-found] -from sklearn.impute import SimpleImputer # type: ignore[import-not-found] -from sklearn.linear_model import LogisticRegression, Ridge # type: ignore[import-not-found] -from sklearn.model_selection import StratifiedShuffleSplit, cross_validate # type: ignore[import-not-found] -from sklearn.pipeline import Pipeline # type: ignore[import-not-found] -from sklearn.preprocessing import LabelEncoder, StandardScaler # type: ignore[import-not-found] +from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.decomposition import PCA +from sklearn.impute import SimpleImputer +from sklearn.linear_model import LinearRegression, LogisticRegression, Ridge +from sklearn.model_selection import StratifiedShuffleSplit, cross_validate +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import LabelEncoder, StandardScaler + +from ..logger import logger if TYPE_CHECKING: from ..base import ConnectivityMatrix +@dataclass +class SiteRegressor(BaseEstimator, TransformerMixin): + sites: NDArray[np.str_] + + model: LinearRegression = field(default_factory=LinearRegression) + + def _get_dummies(self: SiteRegressor, X: pd.DataFrame) -> pd.DataFrame: # noqa: N803 + """ + Convert site labels to dummy variables. + + Args: + X: A DataFrame containing the site labels. + + Returns: + A DataFrame with dummy variables for each site. + """ + return pd.get_dummies(self.sites[X.index], drop_first=False, dtype=np.float32) + + def fit(self: SiteRegressor, X: pd.DataFrame, y: pd.DataFrame | None = None) -> SiteRegressor: # noqa: N803 + """ + Fit the site regressor to the connectivity data. + + Args: + X: Connectivity data with one row per subject. Its index is used to + look up the corresponding site labels in ``self.sites``. + y: Ignored. This parameter exists for compatibility with the scikit-learn API. + + Returns: + self: Returns the instance itself. + + Raises: + ValueError: If the training data contains fewer than two sites. + """ + fit_sites = np.asarray(self.sites[X.index]) + if np.unique(fit_sites).size < 2: + raise ValueError("SiteRegressor requires at least two sites in the training data.") + + y = X + + # Estimate coefficients on the training fold only + self.model.fit(self._get_dummies(X), y) + return self + + def transform(self: SiteRegressor, X: pd.DataFrame) -> pd.DataFrame: # noqa: N803 + """ + Transform the connectivity data by regressing out site effects. + + Args: + X: Connectivity data to correct, with one row per subject. Its index + is used to look up the corresponding site labels in ``self.sites``. + + Returns: + pd.DataFrame: Connectivity data with the fitted site effects removed. + """ + return X - self.model.predict(self._get_dummies(X)) + + def training_pipeline( connectivity_data: NDArray[np.float32], target_labels: NDArray[np.float64] | NDArray[np.str_], @@ -26,6 +87,7 @@ def training_pipeline( n_pca: int, n_jobs: int = 4, random_state: int = 1, + sites: NDArray[np.str_] | None = None, ) -> pd.DataFrame: """Runs a cross-validation pipeline for age or sex prediction. @@ -37,45 +99,83 @@ def training_pipeline( n_pca (int): Number of principal components to extract. n_jobs (int): Number of cores for parallel calculation. random_state (int): Seed for reproducibility. + sites (NDArray[np.str_] | None): Site labels for the data. Returns: pd.DataFrame: Statistics (mean, 95% CI) of the scores obtained. """ - connectivity_data = np.asarray(connectivity_data, dtype=np.float32, order="C") + connectivity_data_frame = pd.DataFrame(connectivity_data, dtype=np.float32) if task_type == "classification": - y_train = LabelEncoder().fit_transform(target_labels) - estimator = LogisticRegression(max_iter=5000, solver="saga", penalty="l2", n_jobs=n_jobs, random_state=random_state) - cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) + y_train = pd.Series(LabelEncoder().fit_transform(target_labels)) # pyright: ignore[reportArgumentType, reportCallIssue] + estimator = LogisticRegression(max_iter=5000, solver="lbfgs", random_state=random_state) + + bins = pd.Series(target_labels) + scoring_metrics = {"accuracy": "accuracy", "roc_auc": "roc_auc"} else: - y_train = np.asarray(target_labels) + y_train = pd.Series(target_labels) estimator = Ridge(alpha=1.0) - bins = pd.qcut(y_train, q=5, labels=False, duplicates="drop") - cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) - splits = list(cv_strategy.split(np.zeros_like(bins), bins)) - cv_strategy = splits + bins, edges = pd.qcut(y_train, q=5, labels=False, retbins=True, duplicates="drop") + labels = pd.Series([f"age-group-{int(round(edges[i]))}-{int(round(edges[i + 1]))}" for i in range(len(edges) - 1)]) + bins = bins.map(labels) scoring_metrics = {"mae": "neg_mean_absolute_error", "r2": "r2"} - pipe = Pipeline( + if sites is not None: + data_frame = pd.DataFrame({"site": sites, "bins": bins}).astype(str) + # Combine bins and sites + bins = data_frame.agg("_".join, axis=1) + + counts = bins.value_counts() + singletons = counts.index[counts == 1] + mask = bins.isin(singletons).to_numpy() + if mask.any(): + count = mask.sum().item() + logger.warning(f"Excluding {count} {'subject'} that are the only subjects in their subgroups: {singletons.tolist()}") + keep = np.flatnonzero(np.logical_not(mask)) + + cv_strategy = StratifiedShuffleSplit(n_splits=n_splits, test_size=0.2, random_state=random_state) + splits = [ + (keep[train_indices], keep[test_indices]) + for train_indices, test_indices in cv_strategy.split(np.zeros(keep.size), bins.to_numpy()[keep]) + ] + + steps: list[tuple[str, BaseEstimator]] = [ + # keep_empty_features=True avoids sklearn's "Skipping features without any + # observed values" warning. + ("imputer", SimpleImputer(strategy="median", keep_empty_features=True).set_output(transform="pandas")), + ] + + if sites is not None: + steps.append(("site_regression", SiteRegressor(sites))) + + steps.extend( [ - ("imputer", SimpleImputer(strategy="median")), ("scaler", StandardScaler()), - ("pca", PCA(n_components=n_pca, svd_solver="randomized", random_state=random_state)), + ( + "pca", + PCA( + n_components=n_pca, + svd_solver="randomized", + random_state=random_state, + ), + ), ("estimator", estimator), ] ) + pipe = Pipeline(steps) with parallel_backend("threading", n_jobs=n_jobs): cv_results = cross_validate( pipe, - connectivity_data, + connectivity_data_frame, y_train, - cv=cv_strategy, + cv=splits, scoring=scoring_metrics, n_jobs=n_jobs, + error_score="raise", ) scores_df = pd.DataFrame({k.replace("test_", ""): v for k, v in cv_results.items() if k.startswith("test_")}) @@ -89,6 +189,7 @@ def age_sex_scores( connectivity_matrices: List[ConnectivityMatrix], ages: NDArray[np.float64], genders: NDArray[np.str_], + sites: NDArray[np.str_] | None, n_splits: int, n_pca: int, n_jobs: int = 4, @@ -100,6 +201,7 @@ def age_sex_scores( connectivity_matrices (List[ConnectivityMatrix]): List of matrix objects. ages: Vector of subject ages. genders: Vector of subject genders. + sites: Vector of subject sites. n_splits (int): Number of splits for cross-validation. n_pca (int): Number of PCA components. n_jobs (int): Number of joblib threads. @@ -119,6 +221,7 @@ def age_sex_scores( n_pca=n_pca, n_jobs=n_jobs, random_state=random_state, + sites=sites, ) age_summary = training_pipeline( @@ -129,6 +232,7 @@ def age_sex_scores( n_pca=n_pca, n_jobs=n_jobs, random_state=random_state, + sites=sites, ) return { diff --git a/wonkyconn/features/calculate_degrees_of_freedom.py b/wonkyconn/features/calculate_degrees_of_freedom.py index cfb9772fd..a023f0aa6 100644 --- a/wonkyconn/features/calculate_degrees_of_freedom.py +++ b/wonkyconn/features/calculate_degrees_of_freedom.py @@ -98,5 +98,5 @@ def calculate_for_key( else: raise ValueError(f"Unexpected value for `{keys}`: {value}") - percentages = pd.Series(proportions) * 100 + percentages: pd.Series = pd.Series(proportions) * 100 return percentages.mean() diff --git a/wonkyconn/features/calculate_gradients_correlation.py b/wonkyconn/features/calculate_gradients_correlation.py index f8bb087ea..f92140501 100644 --- a/wonkyconn/features/calculate_gradients_correlation.py +++ b/wonkyconn/features/calculate_gradients_correlation.py @@ -1,13 +1,15 @@ import glob +import warnings from pathlib import Path from typing import Iterable, List, Tuple import nibabel as nib import numpy as np -from brainspace.gradient import GradientMaps # type: ignore[import-not-found] -from nilearn import image # type: ignore[import-not-found] -from nilearn.connectome import sym_matrix_to_vec, vec_to_sym_matrix # type: ignore[import-not-found] -from nilearn.maskers import NiftiLabelsMasker # type: ignore[import-not-found] +from brainspace.gradient import GradientMaps +from nibabel.nifti1 import Nifti1Image +from nilearn import image +from nilearn.connectome import sym_matrix_to_vec, vec_to_sym_matrix +from nilearn.maskers import NiftiLabelsMasker from scipy import stats from ..base import ConnectivityMatrix @@ -38,7 +40,7 @@ def remove_nan_from_matrix(matrix: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: return conn_clean, kept_idx -def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> nib.Nifti1Image: +def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> Nifti1Image: """Remove ROIs from an atlas that are not present in a connectivity matrix. Args: @@ -46,8 +48,8 @@ def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> nib.Ni kept_idx (np.ndarray): Indices of rows/columns kept after NaN removal. Returns: - tuple[nib.Nifti1Image, list[int]]: Atlas image with only kept ROIs and - the list of kept label values. + nib.Nifti1Image: Atlas image containing only the kept ROIs, with removed + regions set to zero. """ # Load atlas @@ -62,10 +64,10 @@ def remove_nan_roi_atlas(atlas: nib.Nifti1Image, kept_idx: np.ndarray) -> nib.Ni keep_mask = np.isin(atlas_data, kept_labels) kept_atlas_data = atlas_data.copy() kept_atlas_data[~keep_mask] = 0 - return nib.Nifti1Image(kept_atlas_data, atlas.affine, atlas.header), kept_labels + return Nifti1Image(kept_atlas_data, atlas.affine, atlas.header) -def overlapping_atlas_with_mask(subject_atlas: nib.Nifti1Image, group_mask: nib.Nifti1Image) -> nib.Nifti1Image: +def overlapping_atlas_with_mask(subject_atlas: Nifti1Image, group_mask: Nifti1Image) -> Nifti1Image: """Create a new atlas containing only regions that overlap with the group gradient mask. Args: @@ -76,9 +78,9 @@ def overlapping_atlas_with_mask(subject_atlas: nib.Nifti1Image, group_mask: nib. nib.Nifti1Image: Atlas with non-overlapping regions zeroed out. """ - mask_gradient_resampled = image.resample_to_img( + mask_gradient_resampled: Nifti1Image = image.resample_to_img( group_mask, subject_atlas, interpolation="nearest", copy_header=True, force_resample=True - ) + ) # pyright: ignore[reportAssignmentType] # Get arrays atlas_data = subject_atlas.get_fdata() @@ -127,7 +129,10 @@ def group_mean_connectivity( matrices = [np.asarray(cm.load(), dtype=np.float64) for cm in connectivity_matrices] matrices_vec = [sym_matrix_to_vec(mat, discard_diagonal=False) for mat in matrices] - mean_vec = np.nanmean(matrices_vec, axis=0) + with warnings.catch_warnings(): + # Edges that are NaN across all subjects yield all-NaN slices; the NaN is intended. + warnings.simplefilter("ignore", RuntimeWarning) + mean_vec = np.nanmean(matrices_vec, axis=0) mean_matrix = vec_to_sym_matrix(mean_vec, diagonal=None) return mean_matrix @@ -153,7 +158,7 @@ def process_single_matrix( """ matrix = np.asarray(connectivity_matrix, dtype=np.float64) conn_clean, kept_idx = remove_nan_from_matrix(matrix) - atlas_mask_without_nan, _ = remove_nan_roi_atlas(atlas, kept_idx) + atlas_mask_without_nan = remove_nan_roi_atlas(atlas, kept_idx) # filter out labels > 400 atlas_data = atlas_mask_without_nan.get_fdata() @@ -180,7 +185,7 @@ def process_single_matrix( gm = GradientMaps(approach="pca", n_components=5, alignment="procrustes", kernel="normalized_angle") ind_gradient = gm.fit(masked_matrix, reference=group_gradients_np) - return ind_gradient.aligned_, group_gradients_np + return ind_gradient.aligned_, group_gradients_np # pyright: ignore[reportReturnType] def extract_gradients( @@ -201,11 +206,11 @@ def extract_gradients( repo_root = Path(__file__).resolve().parent.parent path_gradients = repo_root / "data" / "gradients" - gradient_mask = nib.load(path_gradients / "gradientmask_cortical.nii.gz") + gradient_mask = nib.nifti1.load(path_gradients / "gradientmask_cortical.nii.gz") # pyright: ignore[reportAttributeAccessIssue] # Load all group gradient templates gradient_files = sorted(glob.glob(str(path_gradients / "templates" / "gradient*_cortical_only.nii.gz"))) - gradient_imgs = [nib.load(fname) for fname in gradient_files] + gradient_imgs = [nib.nifti1.load(fname) for fname in gradient_files] # pyright: ignore[reportAttributeAccessIssue] mean_connectome = group_mean_connectivity(connectivity_matrices) gradient_aligned, template_gradient = process_single_matrix(mean_connectome, atlas, gradient_mask, gradient_imgs) diff --git a/wonkyconn/features/distance_dependence.py b/wonkyconn/features/distance_dependence.py index 03279a996..a16982743 100644 --- a/wonkyconn/features/distance_dependence.py +++ b/wonkyconn/features/distance_dependence.py @@ -19,5 +19,5 @@ def calculate_distance_dependence(qcfc: pd.DataFrame, distance_matrix: npt.NDArr """ i, j = map(np.asarray, zip(*qcfc.index, strict=False)) distance_vector = distance_matrix[i, j] - r, _ = spearmanr(distance_vector, qcfc.correlation, nan_policy="omit") - return np.abs(r) + statistic, _ = spearmanr(distance_vector, qcfc.correlation, nan_policy="omit") + return np.abs(statistic).item() # pyright: ignore[reportArgumentType, reportCallIssue] diff --git a/wonkyconn/features/network.py b/wonkyconn/features/network.py index 5b5b56160..ffa76e1ad 100644 --- a/wonkyconn/features/network.py +++ b/wonkyconn/features/network.py @@ -4,6 +4,7 @@ annotations, ) +import warnings from typing import Tuple import numpy as np @@ -45,8 +46,8 @@ def single_subject_within_network_connectivity( roi_index -= 1 # calculate similarity of the thresholded individual level seed based connectivity with the binary mask - subj_average_connectivity_within_network = [] - subj_variance_connectivity_within_network = [] + subj_average_connectivity_within_network_list = [] + subj_variance_connectivity_within_network_list = [] subj_corr_with_network = [] for idx_roi in roi_index: seed_based_map = connectivity_matrix.load()[int(idx_roi), :] @@ -57,8 +58,10 @@ def single_subject_within_network_connectivity( ] networks = np.asarray(isolate_parcel) networks[networks == 0] = np.nan - mean_within_network_connection = np.nanmean(networks, axis=1) - std_within_network_connection = np.nanstd(networks, axis=1) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + mean_within_network_connection = np.nanmean(networks, axis=1) + std_within_network_connection = np.nanstd(networks, axis=1) # Remove rows with NaN values mask = np.array( @@ -68,21 +71,23 @@ def single_subject_within_network_connectivity( ] ) - correlation_with_given_network = np.corrcoef( - seed_based_map[mask], region_membership[f"yeo7-{yeo_network_index}"].to_numpy()[mask] - ) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + correlation_with_given_network = np.corrcoef( + seed_based_map[mask], region_membership[f"yeo7-{yeo_network_index}"].to_numpy()[mask] + ) - subj_average_connectivity_within_network.append(mean_within_network_connection) - subj_variance_connectivity_within_network.append(std_within_network_connection) + subj_average_connectivity_within_network_list.append(mean_within_network_connection) + subj_variance_connectivity_within_network_list.append(std_within_network_connection) subj_corr_with_network.append(correlation_with_given_network[1, 0]) # summarise of the given subject - subj_average_connectivity_within_network = np.asarray(subj_average_connectivity_within_network).mean(axis=0) - subj_variance_connectivity_within_network = np.asarray(subj_variance_connectivity_within_network).mean(axis=0) + subj_average_connectivity_within_network = np.asarray(subj_average_connectivity_within_network_list).mean(axis=0) + subj_variance_connectivity_within_network = np.asarray(subj_variance_connectivity_within_network_list).mean(axis=0) return ( subj_average_connectivity_within_network, subj_variance_connectivity_within_network, - np.nanmean(subj_corr_with_network), + np.nanmean(subj_corr_with_network), # pyright: ignore[reportReturnType] ) diff --git a/wonkyconn/features/quality_control_connectivity.py b/wonkyconn/features/quality_control_connectivity.py index dcf9fb02e..f79903dc5 100644 --- a/wonkyconn/features/quality_control_connectivity.py +++ b/wonkyconn/features/quality_control_connectivity.py @@ -20,6 +20,7 @@ def calculate_qcfc( data_frame: pd.DataFrame, connectivity_matrices: Iterable[ConnectivityMatrix], metric_key: str = "MeanFramewiseDisplacement", + site_correction: bool = False, ) -> pd.DataFrame: """ metric calculation: quality control / functional connectivity @@ -27,13 +28,16 @@ def calculate_qcfc( For each edge, we then computed the correlation between the weight of that edge and the mean relative RMS motion. QC-FC relationships were calculated as partial correlations that - accounted for participant age and sex + accounted for participant age and sex (and site when ``site_correction`` is enabled). Parameters: data_frame (pd.DataFrame): The data frame containing the covariates "age" and "gender". + "site" is also required if site correction is applied. It needs to have one row for each connectivity matrix. connectivity_matrices (Iterable[ConnectivityMatrix]): The connectivity matrices to calculate QCFC for. metric_key (str, optional): The key of the metric to use for QCFC calculation. Defaults to "MeanFramewiseDisplacement". + site_correction (bool, optional): If True, include site as an additional covariate + (requires a "site" column in ``data_frame``). Defaults to False. Returns: pd.DataFrame: The QCFC values between connectivity matrices and the metric. @@ -44,7 +48,10 @@ def calculate_qcfc( ) if np.isnan(metrics).all(): raise ValueError(f"None of the connectivity matrices have a metric with key '{metric_key}'") - covariates = np.asarray(dmatrix("age + gender", data_frame)) + if site_correction: + covariates = np.asarray(dmatrix("age + C(gender) + C(site)", data_frame)) + else: + covariates = np.asarray(dmatrix("age + C(gender)", data_frame)) connectivity_arrays = [ connectivity_matrix.load() @@ -65,7 +72,8 @@ def calculate_qcfc( axis=1, ) - correlation, count = partial_correlation(connectivity_array, metrics, covariates) + with np.errstate(divide="ignore", invalid="ignore"): + correlation, count = partial_correlation(connectivity_array, metrics, covariates) # pyright: ignore[reportCallIssue] p_value = correlation_p_value(correlation, count) qcfc = pd.DataFrame(dict(i=i, j=j, correlation=correlation, p_value=p_value)) @@ -95,7 +103,7 @@ def significant_level(x: "pd.Series[float]", alpha: float = 0.05, correction: st res, _, _, _ = multipletests(x, alpha=alpha, method=correction) else: res = x < alpha - return res + return np.asarray(res) def calculate_qcfc_percentage(qcfc: pd.DataFrame) -> float: diff --git a/wonkyconn/features/tests/test_network.py b/wonkyconn/features/tests/test_network.py index a6e314a81..1efa365b8 100644 --- a/wonkyconn/features/tests/test_network.py +++ b/wonkyconn/features/tests/test_network.py @@ -22,9 +22,9 @@ def test_single_subject_within_network_connectivity(data_path: Path) -> None: # / "halfpipe/derivatives/halfpipe/" / "sub-10171/func/task-rest/sub-10171_task-rest_feature-cCompCor_atlas-Schaefer2018Combined_timeseries.json" ) - dl.get(str(dseg_path)) - dl.get(str(relmat_path)) - dl.get(str(metadata_path)) + dl.get(str(dseg_path)) # pyright: ignore[reportAttributeAccessIssue] + dl.get(str(relmat_path)) # pyright: ignore[reportAttributeAccessIssue] + dl.get(str(metadata_path)) # pyright: ignore[reportAttributeAccessIssue] atlas = Atlas.create("Schaefer2018Combined", dseg_path) with metadata_path.open("r") as file: metadata = json.load(file) diff --git a/wonkyconn/features/tests/test_prediction.py b/wonkyconn/features/tests/test_prediction.py new file mode 100644 index 000000000..5a0e464fc --- /dev/null +++ b/wonkyconn/features/tests/test_prediction.py @@ -0,0 +1,226 @@ +from dataclasses import dataclass +from pathlib import Path +from typing import Callable, Literal + +import numpy as np +import pandas as pd +import pytest +from numpy.typing import NDArray + +from wonkyconn.base import ConnectivityMatrix +from wonkyconn.features.age_sex_prediction import SiteRegressor, age_sex_scores, training_pipeline + +subject_count = 192 +feature_count = 4_096 +random_state = 1 + +site_labels = np.array(["site-a", "site-b"] * (subject_count // 2)) + + +@dataclass(frozen=True) +class Config: + task: Literal["classification", "regression"] + # How the target relates to the connectivity features: + # "random" - labels are independent noise (nothing to learn) + # "predictable" - a latent factor drives both features and labels + # "site-confounded" - the site drives both features and labels + signal: Literal["random", "predictable", "site-confounded"] + + +@dataclass(frozen=True) +class Dataset: + config: Config + connectivity_data: NDArray[np.float32] + labels: NDArray[np.float64] | NDArray[np.str_] + + +def sex(score: NDArray[np.float64]) -> NDArray[np.str_]: + """Convert a continuous score to a binary sex label based on the mean.""" + return np.where(score > np.mean(score), "male", "female") + + +def age(score: NDArray[np.float64]) -> NDArray[np.float64]: + """Convert a continuous score to an age label based on a linear transformation.""" + return 18.0 + (score - score.min()) / (score.max() - score.min()) * 62.0 + + +@pytest.fixture( + params=[ + pytest.param(Config(task="classification", signal="random"), id="classification-random"), + pytest.param(Config(task="classification", signal="predictable"), id="classification-predictable"), + pytest.param(Config(task="classification", signal="site-confounded"), id="classification-site-confounded"), + pytest.param(Config(task="regression", signal="random"), id="regression-random"), + pytest.param(Config(task="regression", signal="predictable"), id="regression-predictable"), + pytest.param(Config(task="regression", signal="site-confounded"), id="regression-site-confounded"), + ], +) +def dataset(request: pytest.FixtureRequest) -> Dataset: + config: Config = request.param + rng = np.random.default_rng(random_state) + + loadings = rng.normal(size=feature_count) + + connectivity_data: NDArray[np.float32] + labels: NDArray[np.float64] | NDArray[np.str_] + + if config.signal == "site-confounded": + latent = np.where(site_labels == "site-a", 1.0, -1.0) + else: + latent = rng.normal(size=subject_count) + noise = rng.normal(scale=0.5, size=subject_count) + + func: Callable[[NDArray[np.float64]], NDArray[np.float64] | NDArray[np.str_]] + match config.task: + case "classification": + func = sex + case "regression": + func = age + + match config.signal: + case "random": + labels = func(rng.normal(size=subject_count)) + case "predictable": + labels = func(latent + noise) + case "site-confounded": + labels = func(latent) + + noise = rng.normal(scale=0.3, size=(subject_count, feature_count)) + connectivity_data = (np.outer(latent, loadings) + noise).astype(np.float32) + return Dataset( + config=config, + connectivity_data=connectivity_data, + labels=labels, + ) + + +@pytest.mark.parametrize( + "sites", + [ + pytest.param(None, id="without-site-correction"), + pytest.param(site_labels, id="with-site-correction"), + ], +) +def test_training_pipeline( + dataset: Dataset, + sites: NDArray[np.str_] | None, +) -> None: + if dataset.config.task == "classification": + expected_metrics = frozenset({"accuracy", "roc_auc"}) + primary_metric = "roc_auc" + learnable_threshold = 0.8 + chance_ceiling = 0.65 + else: + expected_metrics = frozenset({"mae", "r2"}) + primary_metric = "r2" + learnable_threshold = 0.3 + chance_ceiling = 0.2 + + summary = training_pipeline( + dataset.connectivity_data, + dataset.labels, + task_type=dataset.config.task, + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=random_state, + sites=sites, + ) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == expected_metrics + assert list(summary.columns) == ["mean", "ci_lower", "ci_upper"] + assert np.isfinite(summary.to_numpy()).all() + assert (summary["ci_lower"] <= summary["mean"]).all() + assert (summary["mean"] <= summary["ci_upper"]).all() + + score = float(summary.loc[primary_metric, "mean"]) # type: ignore[arg-type] + + if dataset.config.signal == "site-confounded": + if sites is None: + # Without correction the site confound is fully exploitable. + assert score > learnable_threshold + else: + # Site-effect correction must strip the confound, dropping the + # otherwise-perfect prediction back to chance level. + assert score < chance_ceiling + elif dataset.config.signal == "predictable": + # Labels that carry signal should be recovered above the chance baseline. + assert score > learnable_threshold + else: + # Random labels should not be recoverable above chance level. + assert score < chance_ceiling + + +def test_age_sex_scores(tmp_path: Path) -> None: + rng = np.random.default_rng(random_state) + + region_count = np.sqrt(feature_count).astype(int).item() + + matrices = [] + for i in range(subject_count): + square = rng.normal(size=(region_count, region_count)).astype(np.float32) + path = tmp_path / f"sub-{i}.tsv" + np.savetxt(path, square + square.T, delimiter="\t") + matrices.append(ConnectivityMatrix(path, metadata=dict())) + + ages = rng.uniform(18.0, 80.0, subject_count) + genders = np.array(["m", "f"] * (subject_count // 2)) + sites = np.array(["site-a", "site-b"] * (subject_count // 2)) + + scores = age_sex_scores(matrices, ages, genders, sites, n_splits=3, n_pca=5, n_jobs=1) + assert len(scores) == 8 + assert all(np.isfinite(value) for value in scores.values()) + + +def test_site_regressor_requires_two_sites() -> None: + regressor = SiteRegressor(np.array(["site-a", "site-a"])) + with pytest.raises(ValueError, match="at least two sites"): + regressor.fit(pd.DataFrame(np.zeros((2, 3), dtype=np.float32))) + + +def test_training_pipeline_exclude_singleton_classes(caplog: pytest.LogCaptureFixture) -> None: + labels = np.array(["male"] * 10 + ["female"] * 10) + sites = np.array(["site-a"] * 5 + ["site-b"] * 5 + ["site-a"] * 9 + ["site-b"] * 1) + + rng = np.random.default_rng(random_state) + connectivity_data = rng.normal(size=(len(labels), feature_count)).astype(np.float32) + + with caplog.at_level("WARNING"): + summary = training_pipeline( + connectivity_data, + labels, + task_type="classification", + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=random_state, + sites=sites, + ) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == {"accuracy", "roc_auc"} + assert np.isfinite(summary.to_numpy()).all() + assert "site-b_female" in caplog.text + + +def test_training_pipeline_numeric_sites() -> None: + labels = np.array(["male", "female"] * 10) + sites = np.array([1.0, 2.0] * 10) + + rng = np.random.default_rng(random_state) + connectivity_data = rng.normal(size=(len(labels), feature_count)).astype(np.float32) + + summary = training_pipeline( + connectivity_data, + labels, + task_type="classification", + n_splits=3, + n_pca=5, + n_jobs=1, + random_state=random_state, + sites=sites, + ) + + assert isinstance(summary, pd.DataFrame) + assert set(summary.index) == {"accuracy", "roc_auc"} + assert np.isfinite(summary.to_numpy()).all() diff --git a/wonkyconn/file_index/base.py b/wonkyconn/file_index/base.py index adf98c14b..9c52c62a0 100644 --- a/wonkyconn/file_index/base.py +++ b/wonkyconn/file_index/base.py @@ -5,33 +5,49 @@ from __future__ import annotations from collections import defaultdict -from hashlib import sha1 +from dataclasses import dataclass, field +from functools import cached_property from pathlib import Path -from typing import Mapping +from typing import Any, Mapping -def create_defaultdict_of_set() -> defaultdict[str, set[Path]]: - """Factory for a ``defaultdict`` whose values are sets of ``Path``.""" +def _defaultdict_of_sets() -> dict[str, set[Any]]: return defaultdict(set) -class FileIndex: - def __init__(self) -> None: - self.paths_by_tags: dict[str, dict[str, set[Path]]] = defaultdict(create_defaultdict_of_set) - self.tags_by_paths: dict[Path, dict[str, str]] = defaultdict(dict) +def _defaultdict_of_defaultdict_of_sets() -> dict[str, dict[str, set[Path]]]: + return defaultdict(_defaultdict_of_sets) + + +def _defaultdict_of_dict() -> dict[Path, dict[str, str]]: + return defaultdict(dict) - @property - def hexdigest(self) -> str: - """ - A forty character hash code of the paths in the index, obtained using the `sha1` algorithm. - """ - hash_algorithm = sha1() - for path in sorted(self.tags_by_paths.keys(), key=str): - path_bytes = str(path).encode() - hash_algorithm.update(path_bytes) +@dataclass +class FileIndex: + paths_by_tags: dict[str, dict[str, set[Path]]] = field(default_factory=_defaultdict_of_defaultdict_of_sets) + tags_by_paths: dict[Path, dict[str, str]] = field(default_factory=_defaultdict_of_dict) - return hash_algorithm.hexdigest() + @cached_property + def paths(self) -> set[Path]: + return set(self.tags_by_paths.keys()) + + def _get_phenotypes_without_key(self, key: str) -> set[Path]: + if key not in self.paths_by_tags: + return self.paths.copy() + else: + cache = self.__dict__.setdefault("_phenotypes_without_key", dict()) + phenotypes = cache.get(key) + if phenotypes is None: + phenotypes = self.paths.difference(*self.paths_by_tags[key].values()) + cache[key] = phenotypes + return phenotypes + + def _invalidate_caches(self) -> None: + if "phenotypes" in self.__dict__: + del self.__dict__["phenotypes"] + if "_phenotypes_without_key" in self.__dict__: + del self.__dict__["_phenotypes_without_key"] def get(self, **tags: str | None) -> set[Path]: """ @@ -47,21 +63,20 @@ def get(self, **tags: str | None) -> set[Path]: """ matches: set[Path] | None = None - for key, value in tags.items(): + for key, query in tags.items(): if key not in self.paths_by_tags: return set() values = self.paths_by_tags[key] - if value is not None: - if value not in values: - return set() - paths: set[Path] = values[value] + if query is None: + paths = self._get_phenotypes_without_key(key) + elif query in values: + paths = values[query] else: - paths_in_index = set(self.tags_by_paths.keys()) - paths = paths_in_index.difference(*values.values()) + return set() if matches is not None: - matches &= paths + matches.intersection_update(paths) else: matches = paths.copy() diff --git a/wonkyconn/file_index/bids.py b/wonkyconn/file_index/bids.py index 96e01a99b..802a87a9b 100644 --- a/wonkyconn/file_index/bids.py +++ b/wonkyconn/file_index/bids.py @@ -102,6 +102,7 @@ def parse(path: Path) -> dict[str, str] | None: class BIDSIndex(FileIndex): def put(self, root: Path) -> None: + self._invalidate_caches() for path in root.glob("**/*"): tags = parse(path) diff --git a/wonkyconn/logger.py b/wonkyconn/logger.py index 5bcee33d3..b54fd6789 100644 --- a/wonkyconn/logger.py +++ b/wonkyconn/logger.py @@ -20,17 +20,3 @@ def _setup_logger(log_level: str = "INFO") -> logging.Logger: logger = _setup_logger() - - -def set_verbosity(verbosity: int | list[int]) -> None: - """Set the logger verbosity level (0=ERROR, 1=WARNING, 2=INFO, 3=DEBUG).""" - if isinstance(verbosity, list): - verbosity = verbosity[0] - if verbosity == 0: - logger.setLevel("ERROR") - elif verbosity == 1: - logger.setLevel("WARNING") - elif verbosity == 2: - logger.setLevel("INFO") - elif verbosity == 3: - logger.setLevel("DEBUG") diff --git a/wonkyconn/run.py b/wonkyconn/run.py index c89820016..4aeb1e948 100644 --- a/wonkyconn/run.py +++ b/wonkyconn/run.py @@ -3,14 +3,25 @@ import argparse import sys from pathlib import Path -from typing import Sequence +from typing import Any, Sequence from . import __version__ -from .config import WonkyConnConfig +from .config import WonkyconnConfig from .logger import logger from .workflow import workflow +class VerbosityAction(argparse.Action): + def __call__( + self, + parser: argparse.ArgumentParser, + namespace: argparse.Namespace, + values: Any, + option_string: str | None = None, + ) -> None: + namespace.log_level = {0: "ERROR", 1: "WARNING", 2: "INFO", 3: "DEBUG"}[values] + + def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: """Build the CLI argument parser for wonkyconn.""" parser = argparse.ArgumentParser( @@ -57,25 +68,46 @@ def global_parser(exit_on_error: bool = True) -> argparse.ArgumentParser: help="Specify the atlas label and the path to the atlas file (for example --atlas Schaefer2018 /path/to/atlas.nii.gz)", ) - parser.add_argument("-v", "--version", action="version", version=__version__) - parser.add_argument("--debug", action="store_true", default=False) - parser.add_argument( + metrics_group = parser.add_mutually_exclusive_group(required=False) + metrics_group.add_argument( "--light-mode", required=False, action="store_true", default=False, help="Disable sex and age prediction to reduce runtime.", ) - parser.add_argument( - "--verbosity", - help=""" - Verbosity level. - """, + metrics = ["motion", "analytic-insights", "gradients", "prediction"] + metrics_group.add_argument( + "--metrics", required=False, + nargs="+", + choices=metrics, + default=metrics, + ) + + logging_group = parser.add_mutually_exclusive_group(required=False) + logging_group.add_argument( + "--verbosity", choices=[0, 1, 2, 3], - default=2, type=int, - nargs=1, + action=VerbosityAction, + ) + logging_group.add_argument( + "--log-level", + required=False, + choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], + default="INFO", + type=str, + ) + + parser.add_argument("-v", "--version", action="version", version=__version__) + parser.add_argument("--debug", action="store_true", default=False) + parser.add_argument( + "--site-correction", + required=False, + action="store_true", + default=False, + help="Apply site correction to the data.", ) parser.add_argument( "--textual", @@ -101,11 +133,11 @@ def _loosen_parser_for_gui(parser: argparse.ArgumentParser) -> None: action.required = False -def _build_initial_config(args: argparse.Namespace | None) -> WonkyConnConfig: - return WonkyConnConfig.from_cli_args(args) +def _build_initial_config(args: argparse.Namespace | None) -> WonkyconnConfig: + return WonkyconnConfig.from_cli_args(args) -def _run_textual_ui(config: WonkyConnConfig) -> WonkyConnConfig | None: +def _run_textual_ui(config: WonkyconnConfig) -> WonkyconnConfig | None: if not sys.stdout.isatty(): print("The Textual UI requires an interactive terminal; run without --textual instead.", file=sys.stderr) sys.exit(2) @@ -140,30 +172,28 @@ def main(argv: None | Sequence[str] = None) -> None: if use_gui: _loosen_parser_for_gui(parser) - parsed_args: argparse.Namespace | None + args: argparse.Namespace | None try: - parsed_args = parser.parse_args(raw_args) + args = parser.parse_args(raw_args) except SystemExit as exc: if use_gui and exc.code != 0: - parsed_args = None + args = None else: raise except argparse.ArgumentError: if not use_gui: raise - parsed_args = None + args = None if use_gui: - initial_config = _build_initial_config(parsed_args) - result_config = _run_textual_ui(initial_config) - if result_config is None: + initial_config = _build_initial_config(args) + config = _run_textual_ui(initial_config) + if config is None: return - args_for_workflow = result_config.to_namespace() - debug_enabled = result_config.debug - suppress_warnings = result_config.suppress_warnings + debug_enabled = config.debug + suppress_warnings = config.suppress_warnings else: - config = _build_initial_config(parsed_args) - args_for_workflow = config.to_namespace() + config = _build_initial_config(args) debug_enabled = config.debug suppress_warnings = config.suppress_warnings @@ -173,7 +203,7 @@ def main(argv: None | Sequence[str] = None) -> None: warnings.filterwarnings("ignore", category=RuntimeWarning) try: - workflow(args_for_workflow) + workflow(config) except Exception as e: logger.exception("Exception: %s", e, exc_info=True) if debug_enabled: diff --git a/wonkyconn/tests/test_atlas.py b/wonkyconn/tests/test_atlas.py index cf72ece16..e0c04e143 100644 --- a/wonkyconn/tests/test_atlas.py +++ b/wonkyconn/tests/test_atlas.py @@ -9,10 +9,7 @@ from wonkyconn.atlas import Atlas -def test_dseg_atlas(data_path: Path) -> None: - # atlas_path = data_path / YEO_NETWORK_MAP - # dl.get(str(atlas_path)) - +def test_dseg_atlas() -> None: url = ( "https://raw.githubusercontent.com/ThomasYeoLab/CBIG/master/" "stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/" @@ -29,7 +26,7 @@ def test_dseg_atlas(data_path: Path) -> None: resolution=2, suffix="dseg", extension=".nii.gz", - ) + ) # pyright: ignore[reportCallIssue] assert isinstance(path, Path) atlas = Atlas.create("Schaefer2018400Parcels7Networks", path) @@ -60,10 +57,8 @@ def _get_centroids(path: Path): return centroids -def test_probseg_atlas(data_path: Path) -> None: - # "TODO: @haoting wants to revisit this test, to check if the assertion values make sense" - # atlas_path = data_path / YEO_NETWORK_MAP - # dl.get(str(atlas_path)) +def test_probseg_atlas() -> None: + # TODO: @haoting wants to revisit this test, to check if the assertion values make sense path = get_template( template="MNI152NLin2009cAsym", @@ -72,7 +67,7 @@ def test_probseg_atlas(data_path: Path) -> None: suffix="probseg", resolution=3, # matches “res-03” extension=".nii.gz", - ) + ) # pyright: ignore[reportCallIssue] assert isinstance(path, Path) _centroids = _get_centroids(path) diff --git a/wonkyconn/tests/test_cli.py b/wonkyconn/tests/test_cli.py index 25343b29b..e6edf7a2e 100644 --- a/wonkyconn/tests/test_cli.py +++ b/wonkyconn/tests/test_cli.py @@ -16,13 +16,13 @@ from tqdm.auto import tqdm from wonkyconn import __version__ -from wonkyconn.config import WonkyConnConfig +from wonkyconn.config import WonkyconnConfig from wonkyconn.file_index.bids import BIDSIndex from wonkyconn.run import global_parser, main from wonkyconn.workflow import workflow -def test_version(capsys): +def test_version(capsys: pytest.CaptureFixture[str]) -> None: try: main(["-v"]) except SystemExit: @@ -31,7 +31,7 @@ def test_version(capsys): assert __version__ == captured.out.split()[0] -def test_help(capsys): +def test_help(capsys: pytest.CaptureFixture[str]) -> None: try: main(["-h"]) except SystemExit: @@ -40,68 +40,6 @@ def test_help(capsys): assert "Evaluating the residual motion in fMRI connectome and visualize reports" in captured.out -def test_cli_and_textual_namespace_consistency(tmp_path: Path): - """Ensure CLI (run.py) and Textual UI produce namespaces with the same attributes. - - Both interfaces should produce namespaces that workflow() can consume, - meaning they must have identical attribute names. - """ - # Create minimal valid paths for testing - bids_dir = tmp_path / "bids" - bids_dir.mkdir() - output_dir = tmp_path / "output" - output_dir.mkdir() - phenotypes = tmp_path / "participants.tsv" - phenotypes.touch() - atlas_path = tmp_path / "atlas.nii.gz" - atlas_path.touch() - - # Get namespace from CLI parser - parser = global_parser() - cli_args = parser.parse_args( - [ - str(bids_dir), - str(output_dir), - "group", - "--phenotypes", - str(phenotypes), - "--atlas", - "TestAtlas", - str(atlas_path), - ] - ) - - # Get namespace from WonkyConnConfig (used by Textual UI) - config = WonkyConnConfig( - bids_dir=bids_dir, - output_dir=output_dir, - analysis_level="group", - phenotypes=phenotypes, - atlas=[("TestAtlas", atlas_path)], - verbosity=2, - debug=False, - ) - config_namespace = config.to_namespace() - - # Get the attribute names from both namespaces - cli_attrs = set(vars(cli_args).keys()) - config_attrs = set(vars(config_namespace).keys()) - - # Attributes that are interface-specific and not passed to workflow() - # These are handled separately before calling workflow() - interface_specific_attrs = {"textual", "version", "suppress_warnings"} - - # The workflow-relevant CLI attributes (excluding interface-specific ones) - workflow_cli_attrs = cli_attrs - interface_specific_attrs - - # Both should have the same attributes for workflow consumption - assert workflow_cli_attrs == config_attrs, ( - f"Namespace mismatch!\n" - f"CLI-only attrs (not in config): {workflow_cli_attrs - config_attrs}\n" - f"Config-only attrs (not in CLI): {config_attrs - workflow_cli_attrs}" - ) - - def _copy_file(path: Path, new_path: Path, sub: str) -> None: new_path = Path(re.sub(r"sub-\d+", sub, str(new_path))) new_path.parent.mkdir(parents=True, exist_ok=True) @@ -128,9 +66,9 @@ def _copy_file(path: Path, new_path: Path, sub: str) -> None: @pytest.mark.heavy_smoke -def test_giga_connectome(data_path: Path, tmp_path: Path): +def test_giga_connectome(data_path: Path, tmp_path: Path) -> None: data_path = data_path / "giga_connectome" / "connectome_Schaefer20187Networks_dev" - dl.get(str(data_path)) + dl.get(str(data_path)) # pyright: ignore[reportAttributeAccessIssue] bids_dir = tmp_path / "bids" bids_dir.mkdir() @@ -178,16 +116,18 @@ def test_giga_connectome(data_path: Path, tmp_path: Path): ] args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) assert (output_dir / "metrics.tsv").is_file() assert (output_dir / "metrics.png").is_file() @pytest.mark.smoke -def test_halfpipe(data_path: Path, tmp_path: Path): +@pytest.mark.parametrize("site_correction", [False, True], ids=["no-site-correction", "site-correction"]) +def test_halfpipe(data_path: Path, tmp_path: Path, site_correction: bool) -> None: bids_dir = data_path / "halfpipe" - dl.get(str(bids_dir)) + dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] index = BIDSIndex() index.put(bids_dir) @@ -198,7 +138,7 @@ def test_halfpipe(data_path: Path, tmp_path: Path): phenotypes_path = bids_dir / "participants.tsv" atlas_path = data_path / "atlases" - dl.get(str(atlas_path)) + dl.get(str(atlas_path)) # pyright: ignore[reportAttributeAccessIssue] atlas_args: list[str] = list() atlas_args.append("--atlas") @@ -216,9 +156,12 @@ def test_halfpipe(data_path: Path, tmp_path: Path): str(output_dir), "group", ] + if site_correction: + argv.insert(0, "--site-correction") args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) # Add persistent storage to extract figure as artifact persistent_dir = Path("figures_artifacts") @@ -235,9 +178,10 @@ def test_halfpipe(data_path: Path, tmp_path: Path): @pytest.mark.heavy_smoke -def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): +@pytest.mark.parametrize("site_correction", [False, True], ids=["no-site-correction", "site-correction"]) +def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path, site_correction: bool) -> None: bids_dir = data_path / "halfpipe" - dl.get(str(bids_dir)) + dl.get(str(bids_dir)) # pyright: ignore[reportAttributeAccessIssue] index = BIDSIndex() index.put(bids_dir) @@ -248,7 +192,7 @@ def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): phenotypes_path = bids_dir / "participants.tsv" atlas_path = data_path / "atlases" - dl.get(str(atlas_path)) + dl.get(str(atlas_path)) # pyright: ignore[reportAttributeAccessIssue] atlas_args: list[str] = list() atlas_args.append("--atlas") @@ -265,9 +209,12 @@ def test_halfpipe_with_full_metrics(data_path: Path, tmp_path: Path): str(output_dir), "group", ] + if site_correction: + argv.insert(0, "--site-correction") args = parser.parse_args(argv) - workflow(args) + config = WonkyconnConfig.from_cli_args(args) + workflow(config) # Add persistent storage to extract figure as artifact persistent_dir = Path("figures_artifacts") diff --git a/wonkyconn/tests/test_correlation.py b/wonkyconn/tests/test_correlation.py index 0eb2cec23..0ea7320bf 100644 --- a/wonkyconn/tests/test_correlation.py +++ b/wonkyconn/tests/test_correlation.py @@ -15,7 +15,7 @@ def test_correlation() -> None: y = np.random.normal(size=(m,)) cov = np.random.normal(size=(m, 2)) - correlation, count = partial_correlation(x, y, cov) + correlation, count = partial_correlation(x, y, cov) # pyright: ignore[reportCallIssue] p_value = correlation_p_value(correlation, count) assert np.all(np.isfinite(correlation)) diff --git a/wonkyconn/tests/test_gradients_correlation.py b/wonkyconn/tests/test_gradients_correlation.py index 7226d83ec..b4e840fb2 100644 --- a/wonkyconn/tests/test_gradients_correlation.py +++ b/wonkyconn/tests/test_gradients_correlation.py @@ -36,7 +36,7 @@ def create_fake_connectivity( return connectivity_matrices -def test_gradients(): +def test_gradients() -> None: repo_root = Path(__file__).resolve().parent.parent path = repo_root / "data" / "gradients" @@ -45,7 +45,8 @@ def test_gradients(): connectivity_matrices = create_fake_connectivity(n_regions=434, n_subjects=3) random_gradient, group_gradients = calculate_gradients_correlation.extract_gradients( - connectivity_matrices, atlas=nib.load(atlas_file) + connectivity_matrices, + atlas=nib.nifti1.load(atlas_file), # pyright: ignore[reportAttributeAccessIssue] ) # Calculate similarity for random gradient diff --git a/wonkyconn/tests/test_textual_app.py b/wonkyconn/tests/test_textual_app.py new file mode 100644 index 000000000..17ad2e242 --- /dev/null +++ b/wonkyconn/tests/test_textual_app.py @@ -0,0 +1,7 @@ +import textual.app + + +def test_import_textual_app() -> None: + from wonkyconn import textual_app + + assert issubclass(textual_app.WonkyConnApp, textual.app.App) diff --git a/wonkyconn/tests/test_workflow.py b/wonkyconn/tests/test_workflow.py new file mode 100644 index 000000000..b60787c1e --- /dev/null +++ b/wonkyconn/tests/test_workflow.py @@ -0,0 +1,27 @@ +from pathlib import Path + +import pandas as pd +import pytest + +from wonkyconn.config import WonkyconnConfig +from wonkyconn.workflow import load_data_frame + + +@pytest.mark.parametrize("column", ["participant_id", "age", "gender", "site"]) +def test_load_data_frame_missing_column(tmp_path: Path, column: str) -> None: + """Enabling site correction requires a 'site' column in the phenotypes file.""" + + phenotypes = dict( + participant_id=["sub-1", "sub-2"], + age=[30.0, 40.0], + gender=["m", "f"], + site=["site-a", "site-b"], + ) + del phenotypes[column] + + phenotypes_path = tmp_path / "participants.tsv" + pd.DataFrame(phenotypes).to_csv(phenotypes_path, sep="\t", index=False) + + config = WonkyconnConfig(phenotypes=phenotypes_path, site_correction=True) + with pytest.raises(ValueError, match=column): + load_data_frame(config) diff --git a/wonkyconn/textual_app.py b/wonkyconn/textual_app.py index 0e149e50c..d61369cce 100644 --- a/wonkyconn/textual_app.py +++ b/wonkyconn/textual_app.py @@ -3,10 +3,10 @@ from pathlib import Path from typing import Iterable -from textual.app import App, ComposeResult # type: ignore[import-not-found] -from textual.containers import Center, Container, Horizontal, Vertical # type: ignore[import-not-found] -from textual.events import DescendantFocus # type: ignore[import-not-found] -from textual.widgets import ( # type: ignore[import-not-found] +from textual.app import App, ComposeResult +from textual.containers import Center, Container, Horizontal, Vertical +from textual.events import DescendantFocus +from textual.widgets import ( Button, Checkbox, DirectoryTree, @@ -17,14 +17,14 @@ Select, Static, ) -from textual.widgets._directory_tree import DirEntry # type: ignore[import-not-found] -from textual.widgets._select import NoSelection # type: ignore[import-not-found] -from textual.widgets._tree import Tree # type: ignore[import-not-found] +from textual.widgets._directory_tree import DirEntry +from textual.widgets._select import NoSelection +from textual.widgets._tree import Tree -from .config import WonkyConnConfig +from .config import WonkyconnConfig, all_metrics, light_mode_metrics -class WonkyConnApp(App[WonkyConnConfig | None]): +class WonkyConnApp(App[WonkyconnConfig | None]): """Textual UI for configuring wonkyconn.""" CSS = """ @@ -107,7 +107,7 @@ class WonkyConnApp(App[WonkyConnConfig | None]): BINDINGS = [("escape", "cancel", "Cancel"), ("ctrl+s", "run", "Run")] - def __init__(self, initial_config: WonkyConnConfig): + def __init__(self, initial_config: WonkyconnConfig): super().__init__() self.initial_config = initial_config self.selected_path: Path | None = None @@ -195,12 +195,12 @@ def compose(self) -> ComposeResult: with Horizontal(): yield Select( options=[ - ("Errors only (0)", "0"), - ("Warnings (1)", "1"), - ("Info (2)", "2"), - ("Debug (3)", "3"), + ("Errors only", "ERROR"), + ("Warnings and errors", "WARNING"), + ("Info", "INFO"), + ("Debug", "DEBUG"), ], - id="verbosity", + id="log_level", ) yield Checkbox("Debug", id="debug") yield Checkbox("Skip age/sex prediction and gradient", id="light_mode") @@ -289,7 +289,7 @@ def _load_initial_values(self) -> None: self._path_values["atlas_path"] = str(path) self.query_one("#atlas_path_display", Button).label = f"Atlas Path: {path}" - self.query_one("#verbosity", Select).value = str(self.initial_config.verbosity) + self.query_one("#log_level", Select).value = str(self.initial_config.log_level) self.query_one("#debug", Checkbox).value = bool(self.initial_config.debug) self.query_one("#light_mode", Checkbox).value = bool(self.initial_config.light_mode) self.query_one("#suppress_warnings", Checkbox).value = bool(self.initial_config.suppress_warnings) @@ -301,7 +301,7 @@ def _set_status(self, message: str, error: bool = False) -> None: if error: status.add_class("status-error") - def _validate_and_build_config(self) -> WonkyConnConfig | None: + def _validate_and_build_config(self) -> WonkyconnConfig | None: errors: list[str] = list() bids_str = self._path_values.get("bids_dir", "").strip() @@ -337,8 +337,8 @@ def _validate_and_build_config(self) -> WonkyConnConfig | None: elif not (atlas_path.exists() or atlas_path.is_symlink()) or (atlas_path.exists() and not atlas_path.is_file()): errors.append(f"Atlas file must exist: {atlas_path}") - raw_verbosity = self.query_one("#verbosity", Select[str]).value - verbosity = int("2" if isinstance(raw_verbosity, NoSelection) else raw_verbosity) + raw_log_level = self.query_one("#log_level", Select[str]).value + log_level = "INFO" if isinstance(raw_log_level, NoSelection) else raw_log_level debug = self.query_one("#debug", Checkbox).value light_mode = self.query_one("#light_mode", Checkbox).value suppress_warnings = self.query_one("#suppress_warnings", Checkbox).value @@ -348,15 +348,15 @@ def _validate_and_build_config(self) -> WonkyConnConfig | None: return None assert atlas_path is not None # validated above - config = WonkyConnConfig( + config = WonkyconnConfig( bids_dir=bids_dir, output_dir=output_dir, analysis_level="group", phenotypes=phenotypes, atlas=[(atlas_label, atlas_path)], - verbosity=verbosity, + log_level=log_level, debug=debug, - light_mode=light_mode, + metrics=light_mode_metrics if light_mode else all_metrics, theme="dark" if self.dark else "light", suppress_warnings=suppress_warnings, ) diff --git a/wonkyconn/visualization/plot.py b/wonkyconn/visualization/plot.py index 6be664699..5d565d3a7 100644 --- a/wonkyconn/visualization/plot.py +++ b/wonkyconn/visualization/plot.py @@ -39,23 +39,6 @@ def plot(records: list[dict[str, Any]], group_by: list[str], output_dir: Path) - group_by (list[str]): The list of columns that the results are grouped by. output_dir (Path): The directory to save the plot image into as "metrics.png". """ - # separate dmn similarity from the rest of the metrics - dmn_sim_array = [] - for record in records: - df_dmn_similarity = record.pop("dmn_similarity") - for g in group_by: - df_dmn_similarity[g] = record[g] - dmn_sim_array.append(df_dmn_similarity[["corr_with_dmn"] + group_by]) - - df_dmn_sim_array = pd.concat(dmn_sim_array, ignore_index=True) - if len(group_by) == 2: - df_dmn_sim_array["group_labels"] = df_dmn_sim_array[group_by].apply(lambda x: "-".join(x.astype(str)), axis=1) - else: - df_dmn_sim_array["group_labels"] = df_dmn_sim_array[group_by[0]] - - # summarize the info - for record, dmn_sim in zip(records, dmn_sim_array, strict=True): - record["dmn_similarity"] = dmn_sim["corr_with_dmn"].mean() result_frame = pd.DataFrame.from_records(records, index=group_by) data_frame = result_frame.reset_index() if len(group_by) == 2: # halfpipe @@ -113,7 +96,7 @@ def plot(records: list[dict[str, Any]], group_by: list[str], output_dir: Path) - gcor_axes.set_title("Global correlation (GCOR)") gcor_axes.set_xlabel("Mean correlation") - sns.barplot(data=df_dmn_sim_array, y="group_labels", x="corr_with_dmn", color=palette[4], ax=dmn_mean_axes, errorbar="sd") + sns.barplot(data=data_frame, y="group_labels", x="dmn_similarity_mean", color=palette[4], ax=dmn_mean_axes, errorbar="sd") dmn_mean_axes.set_title("Similarity with DMN") dmn_mean_axes.set_xlabel("Mean correlation") @@ -181,7 +164,7 @@ def plot_degrees_of_freedom_loss( result_frame: pd.DataFrame, degrees_of_freedom_loss_axes: Axes, legend_axes: Axes, - colors: list[str], + colors: list[Any], ) -> None: """Plot stacked bars showing degrees-of-freedom loss by source.""" sns.barplot( diff --git a/wonkyconn/workflow.py b/wonkyconn/workflow.py index ef556225e..d2a78c285 100644 --- a/wonkyconn/workflow.py +++ b/wonkyconn/workflow.py @@ -2,7 +2,6 @@ Process fMRIPrep outputs to timeseries based on denoising strategy. """ -import argparse import sys from collections import defaultdict, namedtuple from pathlib import Path @@ -13,6 +12,8 @@ from numpy import typing as npt from tqdm.auto import tqdm +from wonkyconn.config import Metric, WonkyconnConfig + from .atlas import Atlas from .base import ConnectivityMatrix from .features.age_sex_prediction import age_sex_scores @@ -29,7 +30,7 @@ calculate_qcfc_percentage, ) from .file_index.bids import BIDSIndex -from .logger import logger, set_verbosity +from .logger import logger from .visualization.plot import plot _DEFAULT_N_SPLITS = 20 @@ -50,17 +51,16 @@ def is_halfpipe(index: BIDSIndex) -> bool: return False -def workflow(args: argparse.Namespace) -> None: +def workflow(config: WonkyconnConfig) -> None: """Run the group-level connectivity quality-control pipeline.""" - if "pytest" not in sys.modules: - set_verbosity(args.verbosity) - logger.debug(vars(args)) - - # check if light mode is enabled - if so, it will not run the age and sex prediction and gradient similarity - disable_prediction_gradient = getattr(args, "light_mode", False) + if "pytest" not in sys.modules and config.log_level: + logger.setLevel(config.log_level) + logger.debug(vars(config)) # Check BIDS path - bids_dir = args.bids_dir + bids_dir = config.bids_dir + if bids_dir is None: + raise ValueError("BIDS directory is not specified in the configuration") index = BIDSIndex() index.put(bids_dir) @@ -78,14 +78,16 @@ def workflow(args: argparse.Namespace) -> None: has_header = False # Check output path - output_dir = args.output_dir + output_dir = config.output_dir + if output_dir is None: + raise ValueError("Output directory is not specified in the configuration") output_dir.mkdir(parents=True, exist_ok=True) # Load data frame (participants: age, gender, etc.) - data_frame = load_data_frame(args) + data_frame = load_data_frame(config) # Load atlases - atlases: dict[str, Atlas] = {name: Atlas.create(name, Path(atlas_path_str)) for name, atlas_path_str in args.atlas} + atlases: dict[str, Atlas] = {name: Atlas.create(name, Path(atlas_path_str)) for name, atlas_path_str in config.atlas} logger.debug(f"Atlas dictionary contains: {list(atlases.keys())}") Group = namedtuple("Group", group_by) # type: ignore[misc] @@ -93,7 +95,7 @@ def workflow(args: argparse.Namespace) -> None: grouped_connectivity_matrix: defaultdict[tuple[str, ...], list[ConnectivityMatrix]] = defaultdict(list) segs: set[str] = set() - for timeseries_path in index.get(suffix="timeseries", extension=".tsv"): + for timeseries_path in tqdm(index.get(suffix="timeseries", extension=".tsv"), unit="connectivity matrices"): query = dict(**index.get_tags(timeseries_path)) del query["suffix"] @@ -128,6 +130,9 @@ def workflow(args: argparse.Namespace) -> None: records: list[dict[str, Any]] = list() for key, connectivity_matrices in tqdm(grouped_connectivity_matrix.items(), unit="groups"): + suffix = "_".join(f"{key}-{value}" for key, value in zip(group_by, key, strict=False)) + dmn_similarity_path = output_dir / "dmn-similarity" / f"{suffix}.tsv" + record = make_record( index, data_frame, @@ -137,26 +142,16 @@ def workflow(args: argparse.Namespace) -> None: metric_key, seg_key, atlases, - disable_prediction_gradient, + dmn_similarity_path, + site_correction=config.site_correction, + metrics=config.metrics, ) record.update(dict(zip(group_by, key, strict=False))) - if len(group_by) == 2: - record["dmn_similarity"].to_csv(output_dir / f"dmn_similarity_{'-'.join(group_by)}.tsv", sep="\t") - else: - record["dmn_similarity"].to_csv(output_dir / f"dmn_similarity_{group_by[0]}.tsv", sep="\t") - - dmn_similarity_std = record["dmn_similarity"].loc[:, "corr_with_dmn"].std() - dmn_similarity_avg = record["dmn_similarity"].loc[:, "corr_with_dmn"].mean() - record["dmn_similarity_std"] = dmn_similarity_std - record["dmn_similarity_mean"] = dmn_similarity_avg records.append(record) plot(records, group_by, output_dir) - for record in records: - record.pop("dmn_similarity") - result_frame = pd.DataFrame.from_records(records, index=group_by) result_frame.to_csv(output_dir / "metrics.tsv", sep="\t") @@ -170,10 +165,14 @@ def make_record( metric_key: str, seg_key: str, atlases: dict[str, Atlas], - disable_prediction_gradient: bool, + dmn_similarity_path: Path, + metrics: set[Metric] | None = None, + site_correction: bool = False, ) -> dict[str, Any]: """Compute all QC metrics for a single group of connectivity matrices.""" - seg_subjects: list[str] = list() + subjects: list[str] = list() + missing_subjects: list[str] = list() + filtered: list[ConnectivityMatrix] = list() for c in connectivity_matrices: @@ -187,107 +186,110 @@ def make_record( found = next((s for s in candidates if s in data_frame.index), None) if found: - seg_subjects.append(found) + subjects.append(found) filtered.append(c) else: - logger.info(f"Skipping subject {sub}: not found in phenotype file.") + missing_subjects.append(sub) - # Renaming for consistency - connectivity_matrices[:] = filtered + if missing_subjects: + logger.info( + f"Found {len(subjects)} subjects, skipped {len(missing_subjects)} subjects not in " + f"phenotype file: {', '.join(missing_subjects)}" + ) + else: + logger.info(f"Found {len(subjects)} subjects") - # Slice phenotypes (age, gender, etc.) for just this group - seg_data_frame = data_frame.loc[seg_subjects] - qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key) + # Renaming for consistency + connectivity_matrices[:] = filtered + seg_data_frame = data_frame.loc[subjects] (seg,) = index.get_tag_values(seg_key, {c.path for c in connectivity_matrices}) distance_matrix = distance_matrices[seg] - - gcor = calculate_gcor(connectivity_matrices) - - dmn_similarity_summary, t_stats_dmn_vis_fpn = network_similarity(connectivity_matrices, region_memberships[seg]) atlas = atlases[seg].image - record = dict( - median_absolute_qcfc=calculate_median_absolute(qcfc.correlation), - percentage_significant_qcfc=calculate_qcfc_percentage(qcfc), - distance_dependence=calculate_distance_dependence(qcfc, distance_matrix), - gcor=gcor, - dmn_similarity=dmn_similarity_summary, - dmn_vis_distance_vs_dmn_fpn=t_stats_dmn_vis_fpn, + record: dict[str, Any] = dict( + # Motion + median_absolute_qcfc=np.nan, + percentage_significant_qcfc=np.nan, + distance_dependence=np.nan, + gcor=np.nan, + # Analytic insights + dmn_similarity_std=np.nan, + dmn_similarity_mean=np.nan, + dmn_vis_distance_vs_dmn_fpn=np.nan, + # Gradients + gradients_similarity=np.nan, + # Prediction + sex_auc=np.nan, + sex_auc_ci_lower=np.nan, + sex_auc_ci_upper=np.nan, + sex_accuracy=np.nan, + age_mae=np.nan, + age_mae_ci_lower=np.nan, + age_mae_ci_upper=np.nan, + age_r2=np.nan, + # Degrees of freedom loss **calculate_degrees_of_freedom_loss(connectivity_matrices)._asdict(), ) - if disable_prediction_gradient: - logger.info("Light mode enabled - skipping age and sex prediction, gradient similarity.") + if metrics is None: + raise ValueError("Metrics set is not specified") + + if "motion" in metrics: + qcfc = calculate_qcfc(seg_data_frame, connectivity_matrices, metric_key, site_correction) record.update( - dict( - sex_auc=np.nan, - sex_auc_ci_lower=np.nan, - sex_auc_ci_upper=np.nan, - sex_accuracy=np.nan, - age_mae=np.nan, - age_mae_ci_lower=np.nan, - age_mae_ci_upper=np.nan, - age_r2=np.nan, - gradients_similarity=np.nan, - ) - ) # place holders - return record - # Gradient similarity - gradients, gradients_group = extract_gradients(connectivity_matrices, atlas) - record["gradients_similarity"] = calculate_gradients_similarity(gradients, gradients_group) - - # age / sex predictability metrics - try: - ages = seg_data_frame["age"].to_numpy() - genders = seg_data_frame["gender"].to_numpy() - - scores = age_sex_scores( - connectivity_matrices, - ages=ages, - genders=genders, - n_splits=_DEFAULT_N_SPLITS, - random_state=42, - n_pca=_DEFAULT_N_PCA, - n_jobs=_DEFAULT_N_JOBS, + median_absolute_qcfc=calculate_median_absolute(qcfc.correlation), + percentage_significant_qcfc=calculate_qcfc_percentage(qcfc), + distance_dependence=calculate_distance_dependence(qcfc, distance_matrix), ) - # scores is: - # { - # "sex_auc": float, - # "sex_auc_ci_lower": float, - # "sex_auc_ci_upper": float, - # "sex_accuracy": float, - # "age_mae": float, - # "age_mae_ci_lower": float, - # "age_mae_ci_upper": float, - # "age_r2": float, - # } - record.update(scores) - - except (ValueError, np.linalg.LinAlgError) as exc: - logger.warning(f"[age_sex_prediction] Skipping age/sex prediction for this group due to error: {exc!r}") - # If it fails, we still want consistent columns in the output. + if "analytic-insights" in metrics: + gcor = calculate_gcor(connectivity_matrices) + dmn_similarity_summary, t_stats_dmn_vis_fpn = network_similarity(connectivity_matrices, region_memberships[seg]) + + dmn_similarity_path.parent.mkdir(parents=True, exist_ok=True) + dmn_similarity_summary.to_csv(dmn_similarity_path, sep="\t") + record.update( - dict( - sex_auc=np.nan, - sex_auc_ci_lower=np.nan, - sex_auc_ci_upper=np.nan, - sex_accuracy=np.nan, - age_mae=np.nan, - age_mae_ci_lower=np.nan, - age_mae_ci_upper=np.nan, - age_r2=np.nan, - ) + gcor=gcor, + dmn_similarity_std=dmn_similarity_summary.loc[:, "corr_with_dmn"].std(), + dmn_similarity_mean=dmn_similarity_summary.loc[:, "corr_with_dmn"].mean(), + dmn_vis_distance_vs_dmn_fpn=t_stats_dmn_vis_fpn, ) + if "gradients" in metrics: + gradients, gradients_group = extract_gradients(connectivity_matrices, atlas) + record["gradients_similarity"] = calculate_gradients_similarity(gradients, gradients_group) + + if "prediction" in metrics: + try: + record.update( + age_sex_scores( + connectivity_matrices, + ages=seg_data_frame["age"].to_numpy(), + genders=seg_data_frame["gender"].to_numpy(), + sites=seg_data_frame["site"].to_numpy() if site_correction else None, + n_splits=_DEFAULT_N_SPLITS, + random_state=42, + n_pca=_DEFAULT_N_PCA, + n_jobs=_DEFAULT_N_JOBS, + ) + ) + except (ValueError, np.linalg.LinAlgError) as exc: + logger.warning(f"[age_sex_prediction] Skipping age/sex prediction for this group due to error: {exc!r}") + return record -def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: - """Load a phenotype TSV with ``participant_id``, ``gender``, and ``age`` columns.""" +def load_data_frame(config: WonkyconnConfig) -> pd.DataFrame: + """Load a phenotype TSV with ``participant_id``, ``gender``, and ``age`` columns. + If site correction is enabled, the ``site`` column is also required. + """ + path = config.phenotypes + if path is None: + raise ValueError("Phenotypes file path is not specified in the configuration") data_frame = pd.read_csv( - args.phenotypes, + path, sep="\t", index_col="participant_id", dtype={"participant_id": str}, @@ -296,4 +298,8 @@ def load_data_frame(args: argparse.Namespace) -> pd.DataFrame: raise ValueError('Phenotypes file is missing the "gender" column') if "age" not in data_frame.columns: raise ValueError('Phenotypes file is missing the "age" column') + if config.site_correction: + if "site" not in data_frame.columns: + raise ValueError('Phenotypes file is missing the "site" column required for site correction') + logger.info("Site correction is enabled") return data_frame