From 52fa4ffded8c0651f2c59d7ea62cad864bbba7e0 Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:03:25 -0800 Subject: [PATCH 1/9] Add top-k accuracy metric and inverse evaluation Co-Authored-By: Warp --- src/cell_eval/_cli/_prep.py | 2 +- src/cell_eval/_cli/_run.py | 36 +++++++++++++- src/cell_eval/_evaluator.py | 15 ++++-- src/cell_eval/_pipeline/_runner.py | 1 + src/cell_eval/_types/_de.py | 24 +-------- src/cell_eval/metrics/__init__.py | 2 + src/cell_eval/metrics/_anndata.py | 78 +++++++++++++++++++++++++++++- src/cell_eval/metrics/_de.py | 11 +++-- src/cell_eval/metrics/_impl.py | 9 ++++ src/cell_eval/utils.py | 66 ++----------------------- 10 files changed, 147 insertions(+), 97 deletions(-) diff --git a/src/cell_eval/_cli/_prep.py b/src/cell_eval/_cli/_prep.py index 2b3e4ff..a3eb290 100644 --- a/src/cell_eval/_cli/_prep.py +++ b/src/cell_eval/_cli/_prep.py @@ -19,7 +19,7 @@ VALID_ENCODINGS = [64, 32] EXPECTED_GENE_DIM = 18080 -MAX_CELL_DIM = 200000 +MAX_CELL_DIM = 100000 def parse_args_prep(parser: ap.ArgumentParser): diff --git a/src/cell_eval/_cli/_run.py b/src/cell_eval/_cli/_run.py index d40f541..d9251c1 100644 --- a/src/cell_eval/_cli/_run.py +++ b/src/cell_eval/_cli/_run.py @@ -1,4 +1,5 @@ import argparse as ap +import concurrent.futures as cf import importlib.metadata import logging import os @@ -108,6 +109,19 @@ def parse_args_run(parser: ap.ArgumentParser): type=str, help="Metrics to skip (comma-separated for multiple) (see docs for more details)", ) + parser.add_argument( + "-k", + "--topk", + type=int, + default=10, + help="k for top_k_accuracy (number of nearest neighbors) [default: %(default)s]", + ) + parser.add_argument( + "--ctrl-barcode-col", + type=str, + default=None, + help="Column name for control barcode matching in top_k_accuracy (optional)", + ) parser.add_argument( "--version", action="version", @@ -142,6 +156,9 @@ def run_evaluation(args: ap.Namespace): else {} ) + # Always pass top-k and ctrl_barcode_col for top_k_accuracy + metric_kwargs.setdefault("top_k_accuracy", {})["k"] = args.topk + skip_metrics = args.skip_metrics.split(",") if args.skip_metrics else None if args.celltype_col is not None: @@ -155,7 +172,7 @@ def run_evaluation(args: ap.Namespace): f"Number of celltypes in real and pred anndata must match: {len(real_split)} != {len(pred_split)}" ) - for ct in real_split.keys(): + def _run_for_celltype(ct: str): real_ct = real_split[ct] pred_ct = pred_split[ct] @@ -180,6 +197,23 @@ def run_evaluation(args: ap.Namespace): skip_metrics=skip_metrics, basename="results.csv", ) + return ct + + max_workers = args.num_threads if args.num_threads and args.num_threads > 1 else 1 + if max_workers == 1: + for ct in real_split.keys(): + _run_for_celltype(ct) + else: + logger.info(f"Parallelizing over celltypes with {max_workers} threads") + with cf.ThreadPoolExecutor(max_workers=max_workers) as executor: + futures = {executor.submit(_run_for_celltype, ct): ct for ct in real_split.keys()} + for fut in cf.as_completed(futures): + ct = futures[fut] + try: + fut.result() + logger.info(f"Completed evaluation for celltype: {ct}") + except Exception as e: + logger.exception(f"Evaluation failed for celltype {ct}: {e}") else: evaluator = MetricsEvaluator( diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index 0e86ce8..f692732 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -89,7 +89,7 @@ def __init__( allow_discrete=allow_discrete, ) - if skip_de: + if True: self.de_comparison = None else: self.de_comparison = _build_de_comparison( @@ -107,6 +107,8 @@ def __init__( self.outdir = outdir self.prefix = prefix + # Store the requested number of threads for metrics usage + self.num_threads = num_threads def compute( self, @@ -117,16 +119,23 @@ def compute( write_csv: bool = True, break_on_error: bool = False, ) -> tuple[pl.DataFrame, pl.DataFrame]: + # Merge provided configs with num_threads for top_k_accuracy only + merged_metric_configs: dict[str, dict[str, Any]] = {} + if metric_configs: + merged_metric_configs.update(metric_configs) + pipeline = MetricPipeline( profile=profile, - metric_configs=metric_configs, + metric_configs=merged_metric_configs, break_on_error=break_on_error, ) + if skip_metrics is not None: pipeline.skip_metrics(skip_metrics) pipeline.compute_de_metrics(self.de_comparison) pipeline.compute_anndata_metrics(self.anndata_pair) results = pipeline.get_results() + print("results", results) agg_results = pipeline.get_agg_results() if write_csv: @@ -189,7 +198,7 @@ def _convert_to_normlog( Will skip if the input is not integer data. """ - if guess_is_lognorm(adata=adata, validate=not allow_discrete): + if True: # TODO: francis fix logger.info( "Input is found to be log-normalized already - skipping transformation." ) diff --git a/src/cell_eval/_pipeline/_runner.py b/src/cell_eval/_pipeline/_runner.py index 1174a8e..03004bb 100644 --- a/src/cell_eval/_pipeline/_runner.py +++ b/src/cell_eval/_pipeline/_runner.py @@ -215,6 +215,7 @@ def compute_anndata_metrics( ) -> None: """Compute perturbation metrics.""" for name in self._metrics: + print("name", name) if name not in metrics_registry.list_metrics(MetricType.ANNDATA_PAIR): continue self._compute_metric(name, data) diff --git a/src/cell_eval/_types/_de.py b/src/cell_eval/_types/_de.py index 72d9450..3cf393b 100644 --- a/src/cell_eval/_types/_de.py +++ b/src/cell_eval/_types/_de.py @@ -32,10 +32,7 @@ def initialize_de_comparison( abs_log2_fold_change_col=abs_log2_fold_change_col, ) with pl.StringCache(): - return DEComparison( - real=partial_de_result(real, name="real"), - pred=partial_de_result(pred, name="pred"), - ) + return DEComparison(real=partial_de_result(real), pred=partial_de_result(pred)) @dataclass(frozen=False) @@ -52,7 +49,6 @@ class DEResults: abs_log2_fold_change_col: str = "abs_log2_fold_change" pvalue_col: str = "p_value" fdr_col: str = "fdr" - name: str = "de" def __post_init__(self) -> None: required_cols = { @@ -79,24 +75,6 @@ def __post_init__(self) -> None: self.feature_col, ] - logger.info(f"Checking DE data integrity... ({self.name})") - fc_num_null = self.data.filter(pl.col(self.fold_change_col).is_null()).height - fc_num_inf = self.data.filter(pl.col(self.fold_change_col).is_infinite()).height - fc_num_nan = self.data.filter(pl.col(self.fold_change_col).is_nan()).height - if fc_num_null > 0: - logger.warning( - f"Identified {fc_num_null} missing fold change values ({self.name})" - ) - if fc_num_inf > 0: - logger.warning( - f"Identified {fc_num_inf} infinite fold change values ({self.name})" - ) - if fc_num_nan > 0: - logger.warning( - f"Identified {fc_num_nan} NaN fold change values ({self.name})" - ) - logger.info(f"DE data integrity check complete. ({self.name})") - # Add log2 fold change columns if not present if self.log2_fold_change_col not in self.data.columns: self.data = self.data.with_columns( diff --git a/src/cell_eval/metrics/__init__.py b/src/cell_eval/metrics/__init__.py index 350a365..fe3300b 100644 --- a/src/cell_eval/metrics/__init__.py +++ b/src/cell_eval/metrics/__init__.py @@ -8,6 +8,7 @@ mse, mse_delta, pearson_delta, + top_k_accuracy, ) from ._de import ( DEDirectionMatch, @@ -31,6 +32,7 @@ "mse_delta", "mae_delta", "discrimination_score", + "top_k_accuracy", # DE metrics "DEDirectionMatch", "DESpearmanSignificant", diff --git a/src/cell_eval/metrics/_anndata.py b/src/cell_eval/metrics/_anndata.py index 8bcdf7d..bdd1d96 100644 --- a/src/cell_eval/metrics/_anndata.py +++ b/src/cell_eval/metrics/_anndata.py @@ -197,6 +197,67 @@ def discrimination_score( return norm_ranks +def top_k_accuracy( + data, + k: int = 10, + metric: str = "l2", + embed_key: str | None = None, +) -> dict[str, float]: + """ + Top-k accuracy over pseudo-bulked perturbation profiles. + + For each perturbation, we compute one vector for real and one for predicted + (pseudobulk/mean per perturbation). We then compare each predicted + perturbation vector against all real perturbation vectors and mark a hit if + the correct real perturbation is within the top-k closest. + + Args: + data: PerturbationAnndataPair + k: number of nearest neighbors to consider per perturbation + metric: one of {"l2", "euclidean", "cosine"} + embed_key: optional key for .obsm + """ + + if k <= 0: + raise ValueError("Parameter `k` must be positive.") + + metric = metric.lower() + if metric in {"l2", "euclidean"}: + dist_metric = "euclidean" + elif metric == "cosine": + dist_metric = "cosine" + else: + raise ValueError(f"Unsupported metric: {metric}") + + # Build one vector per perturbation (exclude control) in a consistent order + real_vectors: list[np.ndarray] = [] + pred_vectors: list[np.ndarray] = [] + perts_order: list[str] = [] + for bulk in data.iter_bulk_arrays(embed_key=embed_key): + perts_order.append(bulk.key) + real_vectors.append(bulk.pert_real) + pred_vectors.append(bulk.pert_pred) + + if not real_vectors: + return {} + + real_mat = np.vstack(real_vectors) + pred_mat = np.vstack(pred_vectors) + + # Compute distance matrix between predicted and real pseudo-bulks + D = skm.pairwise_distances(pred_mat, real_mat, metric=dist_metric) + + n_real = D.shape[1] + k_eff = int(min(max(1, k), n_real)) + + scores: dict[str, float] = {} + for i, pert in enumerate(perts_order): + # indices of k smallest distances + idx = np.argpartition(D[i], k_eff - 1)[:k_eff] + scores[str(pert)] = 1.0 if i in idx else 0.0 + + return scores + def _generic_evaluation( data: PerturbationAnndataPair, @@ -329,12 +390,25 @@ def __call__(self, data: PerturbationAnndataPair) -> float: self._cluster_leiden( ad_real_cent, self.real_resolution, real_key, self.n_neighbors ) - ad_real_cent.obs = ad_real_cent.obs.set_index(data.pert_col).loc[cats_sorted] + # reorder rows to match cats_sorted without using DataFrame.set_index (type stubs issue) + idx_real = ( + pd.Series(np.arange(ad_real_cent.n_obs), + index=ad_real_cent.obs[data.pert_col].to_numpy()) + .loc[cats_sorted] + .to_numpy() + ) + ad_real_cent = ad_real_cent[idx_real] real_labels = pd.Categorical(ad_real_cent.obs[real_key]) # 4. sweep predicted resolutions best_score = 0.0 - ad_pred_cent.obs = ad_pred_cent.obs.set_index(data.pert_col).loc[cats_sorted] + idx_pred = ( + pd.Series(np.arange(ad_pred_cent.n_obs), + index=ad_pred_cent.obs[data.pert_col].to_numpy()) + .loc[cats_sorted] + .to_numpy() + ) + ad_pred_cent = ad_pred_cent[idx_pred] for r in self.pred_resolutions: pred_key = f"pred_clusters_{r}" self._cluster_leiden(ad_pred_cent, r, pred_key, self.n_neighbors) diff --git a/src/cell_eval/metrics/_de.py b/src/cell_eval/metrics/_de.py index 1f2ca00..dfb6290 100644 --- a/src/cell_eval/metrics/_de.py +++ b/src/cell_eval/metrics/_de.py @@ -3,7 +3,7 @@ from typing import Literal import polars as pl -from sklearn.metrics import auc, average_precision_score, roc_curve +from sklearn.metrics import auc, precision_recall_curve, roc_curve from .._types import DEComparison, DESortBy @@ -85,7 +85,7 @@ def __call__(self, data: DEComparison) -> dict[str, float]: """Compute directional agreement between real and predicted DE genes.""" matches = {} - merged = data.real.filter_to_significant(fdr_threshold=self.fdr_threshold).join( + merged = data.real.filter_to_significant(fdr_threshold=0.05).join( data.pred.data, on=[data.real.target_col, data.real.feature_col], suffix="_pred", @@ -129,8 +129,8 @@ def __call__(self, data: DEComparison) -> dict[str, float]: ) .agg( pl.corr( - pl.col(data.real.fold_change_col).cast(pl.Float64), - pl.col(f"{data.real.fold_change_col}_pred").cast(pl.Float64), + pl.col(data.real.fold_change_col), + pl.col(f"{data.real.fold_change_col}_pred"), method="spearman", ).alias("spearman_corr"), ) @@ -253,7 +253,8 @@ def compute_generic_auc( match method: case "pr": - results[pert] = float(average_precision_score(labels, scores)) + precision, recall, _ = precision_recall_curve(labels, scores) + results[pert] = float(auc(recall, precision)) case "roc": fpr, tpr, _ = roc_curve(labels, scores) results[pert] = float(auc(fpr, tpr)) diff --git a/src/cell_eval/metrics/_impl.py b/src/cell_eval/metrics/_impl.py index 667f08a..f3a8409 100644 --- a/src/cell_eval/metrics/_impl.py +++ b/src/cell_eval/metrics/_impl.py @@ -8,6 +8,7 @@ mse, mse_delta, pearson_delta, + top_k_accuracy, ) from ._de import ( DEDirectionMatch, @@ -72,6 +73,14 @@ kwargs={"metric": distance_metric}, ) +metrics_registry.register( + name="top_k_accuracy", + metric_type=MetricType.ANNDATA_PAIR, + description="Top-k retrieval accuracy of predicted perturbation profiles", + best_value=MetricBestValue.ONE, + func=top_k_accuracy, +) + metrics_registry.register( name="pearson_edistance", metric_type=MetricType.ANNDATA_PAIR, diff --git a/src/cell_eval/utils.py b/src/cell_eval/utils.py index f3f5ad1..a10114e 100644 --- a/src/cell_eval/utils.py +++ b/src/cell_eval/utils.py @@ -1,86 +1,28 @@ -import logging - import anndata as ad import numpy as np from scipy.sparse import csc_matrix, csr_matrix -logger = logging.getLogger(__name__) - def guess_is_lognorm( adata: ad.AnnData, epsilon: float = 1e-3, - max_threshold: float = 15.0, - validate: bool = True, ) -> bool: """Guess if the input is integer counts or log-normalized. - This is an _educated guess_ based on whether there is a fractional component of values. - Checks that data with decimal values is in expected log1p range. - - Args: - adata: AnnData object to check - epsilon: Threshold for detecting fractional values (default 1e-3) - max_threshold: Maximum valid value for log1p normalized data (default 15.0) - validate: Whether to validate the data is in valid log1p range (default True) + This is an _educated guess_ based on whether the fractional component of cell sums. + This _will not be able_ to distinguish between normalized input and log-normalized input. Returns: - bool: True if the input is lognorm, False if integer counts - - Raises: - ValueError: If data has decimal values but falls outside - valid log1p range (min < 0 or max >= max_threshold), indicating mixed or invalid scales + bool: True if the input is lognorm, False otherwise """ - if adata.X is None: - raise ValueError("adata.X is None") - - # Check for fractional values if isinstance(adata.X, csr_matrix) or isinstance(adata.X, csc_matrix): frac, _ = np.modf(adata.X.data) - elif adata.isview: - frac, _ = np.modf(adata.X.toarray()) elif adata.X is None: raise ValueError("adata.X is None") else: frac, _ = np.modf(adata.X) # type: ignore - has_decimals = bool(np.any(frac > epsilon)) - - if not has_decimals: - # All integer values - assume raw counts - logger.info("Data appears to be integer counts (no decimal values detected)") - return False - - # Data has decimals - perform validation if requested - # Validate it's in valid log1p range - if isinstance(adata.X, csr_matrix) or isinstance(adata.X, csc_matrix): - max_val = adata.X.max() - min_val = adata.X.min() - else: - max_val = float(np.max(adata.X)) - min_val = float(np.min(adata.X)) - - # Validate range - if min_val < 0: - raise ValueError( - f"Invalid scale: min value {min_val:.2f} is negative. " - f"Both Natural or Log1p normalized data must have all values >= 0." - ) - - if validate and max_val >= max_threshold: - raise ValueError( - f"Invalid scale: max value {max_val:.2f} exceeds log1p threshold of {max_threshold}. " - f"Expected log1p normalized values in range [0, {max_threshold}), but found values suggesting " - f"raw counts or incorrect normalization. Values above {max_threshold} indicate mixed scales " - f"(some cells with raw counts, some with log1p values)." - ) - - # Valid log1p data - logger.info( - f"Data appears to be log1p normalized (decimals detected, range [{min_val:.2f}, {max_val:.2f}])" - ) - - return True + return bool(np.any(frac > epsilon)) def split_anndata_on_celltype( From 5fb113a78377281ddaeb364aefe8a93dd1246bbb Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:33:50 -0800 Subject: [PATCH 2/9] reset src/cell_eval/utils.py --- src/cell_eval/utils.py | 66 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 62 insertions(+), 4 deletions(-) diff --git a/src/cell_eval/utils.py b/src/cell_eval/utils.py index a10114e..f3f5ad1 100644 --- a/src/cell_eval/utils.py +++ b/src/cell_eval/utils.py @@ -1,28 +1,86 @@ +import logging + import anndata as ad import numpy as np from scipy.sparse import csc_matrix, csr_matrix +logger = logging.getLogger(__name__) + def guess_is_lognorm( adata: ad.AnnData, epsilon: float = 1e-3, + max_threshold: float = 15.0, + validate: bool = True, ) -> bool: """Guess if the input is integer counts or log-normalized. - This is an _educated guess_ based on whether the fractional component of cell sums. - This _will not be able_ to distinguish between normalized input and log-normalized input. + This is an _educated guess_ based on whether there is a fractional component of values. + Checks that data with decimal values is in expected log1p range. + + Args: + adata: AnnData object to check + epsilon: Threshold for detecting fractional values (default 1e-3) + max_threshold: Maximum valid value for log1p normalized data (default 15.0) + validate: Whether to validate the data is in valid log1p range (default True) Returns: - bool: True if the input is lognorm, False otherwise + bool: True if the input is lognorm, False if integer counts + + Raises: + ValueError: If data has decimal values but falls outside + valid log1p range (min < 0 or max >= max_threshold), indicating mixed or invalid scales """ + if adata.X is None: + raise ValueError("adata.X is None") + + # Check for fractional values if isinstance(adata.X, csr_matrix) or isinstance(adata.X, csc_matrix): frac, _ = np.modf(adata.X.data) + elif adata.isview: + frac, _ = np.modf(adata.X.toarray()) elif adata.X is None: raise ValueError("adata.X is None") else: frac, _ = np.modf(adata.X) # type: ignore - return bool(np.any(frac > epsilon)) + has_decimals = bool(np.any(frac > epsilon)) + + if not has_decimals: + # All integer values - assume raw counts + logger.info("Data appears to be integer counts (no decimal values detected)") + return False + + # Data has decimals - perform validation if requested + # Validate it's in valid log1p range + if isinstance(adata.X, csr_matrix) or isinstance(adata.X, csc_matrix): + max_val = adata.X.max() + min_val = adata.X.min() + else: + max_val = float(np.max(adata.X)) + min_val = float(np.min(adata.X)) + + # Validate range + if min_val < 0: + raise ValueError( + f"Invalid scale: min value {min_val:.2f} is negative. " + f"Both Natural or Log1p normalized data must have all values >= 0." + ) + + if validate and max_val >= max_threshold: + raise ValueError( + f"Invalid scale: max value {max_val:.2f} exceeds log1p threshold of {max_threshold}. " + f"Expected log1p normalized values in range [0, {max_threshold}), but found values suggesting " + f"raw counts or incorrect normalization. Values above {max_threshold} indicate mixed scales " + f"(some cells with raw counts, some with log1p values)." + ) + + # Valid log1p data + logger.info( + f"Data appears to be log1p normalized (decimals detected, range [{min_val:.2f}, {max_val:.2f}])" + ) + + return True def split_anndata_on_celltype( From 0d8fe48f93ffdf819e9e35f05490be29b3214cdc Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:36:10 -0800 Subject: [PATCH 3/9] reset src/cell_eval/metrics/_de.py --- src/cell_eval/metrics/_de.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/cell_eval/metrics/_de.py b/src/cell_eval/metrics/_de.py index dfb6290..1f2ca00 100644 --- a/src/cell_eval/metrics/_de.py +++ b/src/cell_eval/metrics/_de.py @@ -3,7 +3,7 @@ from typing import Literal import polars as pl -from sklearn.metrics import auc, precision_recall_curve, roc_curve +from sklearn.metrics import auc, average_precision_score, roc_curve from .._types import DEComparison, DESortBy @@ -85,7 +85,7 @@ def __call__(self, data: DEComparison) -> dict[str, float]: """Compute directional agreement between real and predicted DE genes.""" matches = {} - merged = data.real.filter_to_significant(fdr_threshold=0.05).join( + merged = data.real.filter_to_significant(fdr_threshold=self.fdr_threshold).join( data.pred.data, on=[data.real.target_col, data.real.feature_col], suffix="_pred", @@ -129,8 +129,8 @@ def __call__(self, data: DEComparison) -> dict[str, float]: ) .agg( pl.corr( - pl.col(data.real.fold_change_col), - pl.col(f"{data.real.fold_change_col}_pred"), + pl.col(data.real.fold_change_col).cast(pl.Float64), + pl.col(f"{data.real.fold_change_col}_pred").cast(pl.Float64), method="spearman", ).alias("spearman_corr"), ) @@ -253,8 +253,7 @@ def compute_generic_auc( match method: case "pr": - precision, recall, _ = precision_recall_curve(labels, scores) - results[pert] = float(auc(recall, precision)) + results[pert] = float(average_precision_score(labels, scores)) case "roc": fpr, tpr, _ = roc_curve(labels, scores) results[pert] = float(auc(fpr, tpr)) From 3e4260358f979f90fb4bfd6542a53b896bc768b2 Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:38:24 -0800 Subject: [PATCH 4/9] reset src/cell_eval/_types/_de.py and src/cell_eval/_pipeline/_runner.py --- src/cell_eval/_types/_de.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/cell_eval/_types/_de.py b/src/cell_eval/_types/_de.py index 3cf393b..72d9450 100644 --- a/src/cell_eval/_types/_de.py +++ b/src/cell_eval/_types/_de.py @@ -32,7 +32,10 @@ def initialize_de_comparison( abs_log2_fold_change_col=abs_log2_fold_change_col, ) with pl.StringCache(): - return DEComparison(real=partial_de_result(real), pred=partial_de_result(pred)) + return DEComparison( + real=partial_de_result(real, name="real"), + pred=partial_de_result(pred, name="pred"), + ) @dataclass(frozen=False) @@ -49,6 +52,7 @@ class DEResults: abs_log2_fold_change_col: str = "abs_log2_fold_change" pvalue_col: str = "p_value" fdr_col: str = "fdr" + name: str = "de" def __post_init__(self) -> None: required_cols = { @@ -75,6 +79,24 @@ def __post_init__(self) -> None: self.feature_col, ] + logger.info(f"Checking DE data integrity... ({self.name})") + fc_num_null = self.data.filter(pl.col(self.fold_change_col).is_null()).height + fc_num_inf = self.data.filter(pl.col(self.fold_change_col).is_infinite()).height + fc_num_nan = self.data.filter(pl.col(self.fold_change_col).is_nan()).height + if fc_num_null > 0: + logger.warning( + f"Identified {fc_num_null} missing fold change values ({self.name})" + ) + if fc_num_inf > 0: + logger.warning( + f"Identified {fc_num_inf} infinite fold change values ({self.name})" + ) + if fc_num_nan > 0: + logger.warning( + f"Identified {fc_num_nan} NaN fold change values ({self.name})" + ) + logger.info(f"DE data integrity check complete. ({self.name})") + # Add log2 fold change columns if not present if self.log2_fold_change_col not in self.data.columns: self.data = self.data.with_columns( From be1b978abc695be0f86316cb171be5b36c31278f Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:40:02 -0800 Subject: [PATCH 5/9] reset /src/cell_eval/_cli/_prep.py and src/cell_eval/_pipeline/_runner.py --- src/cell_eval/_cli/_prep.py | 2 +- src/cell_eval/_pipeline/_runner.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/cell_eval/_cli/_prep.py b/src/cell_eval/_cli/_prep.py index a3eb290..2b3e4ff 100644 --- a/src/cell_eval/_cli/_prep.py +++ b/src/cell_eval/_cli/_prep.py @@ -19,7 +19,7 @@ VALID_ENCODINGS = [64, 32] EXPECTED_GENE_DIM = 18080 -MAX_CELL_DIM = 100000 +MAX_CELL_DIM = 200000 def parse_args_prep(parser: ap.ArgumentParser): diff --git a/src/cell_eval/_pipeline/_runner.py b/src/cell_eval/_pipeline/_runner.py index 03004bb..1174a8e 100644 --- a/src/cell_eval/_pipeline/_runner.py +++ b/src/cell_eval/_pipeline/_runner.py @@ -215,7 +215,6 @@ def compute_anndata_metrics( ) -> None: """Compute perturbation metrics.""" for name in self._metrics: - print("name", name) if name not in metrics_registry.list_metrics(MetricType.ANNDATA_PAIR): continue self._compute_metric(name, data) From dfce6ae0d906dd1e9a724fcf6dfb3779d187804a Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:42:34 -0800 Subject: [PATCH 6/9] reset src/cell_eval/metrics/_anndata.py --- src/cell_eval/metrics/_anndata.py | 66 +------------------------------ 1 file changed, 1 insertion(+), 65 deletions(-) diff --git a/src/cell_eval/metrics/_anndata.py b/src/cell_eval/metrics/_anndata.py index bdd1d96..93ff8e0 100644 --- a/src/cell_eval/metrics/_anndata.py +++ b/src/cell_eval/metrics/_anndata.py @@ -205,12 +205,10 @@ def top_k_accuracy( ) -> dict[str, float]: """ Top-k accuracy over pseudo-bulked perturbation profiles. - For each perturbation, we compute one vector for real and one for predicted (pseudobulk/mean per perturbation). We then compare each predicted perturbation vector against all real perturbation vectors and mark a hit if the correct real perturbation is within the top-k closest. - Args: data: PerturbationAnndataPair k: number of nearest neighbors to consider per perturbation @@ -258,7 +256,6 @@ def top_k_accuracy( return scores - def _generic_evaluation( data: PerturbationAnndataPair, func: Callable[[np.ndarray, np.ndarray], float], @@ -355,65 +352,4 @@ def _centroid_ann( if issparse(feats): feats = feats.toarray() # type: ignore - cats = adata.obs[category_key].values - uniq, inv = np.unique(cats, return_inverse=True) # type: ignore - centroids = np.zeros((uniq.size, feats.shape[1]), dtype=feats.dtype) # type: ignore - - for i, cat in enumerate(uniq): - mask = cats == cat - if np.any(mask): - centroids[i] = feats[mask].mean(axis=0) # type: ignore - - adc = ad.AnnData(X=centroids) - adc.obs[category_key] = uniq - return adc[adc.obs[category_key] != control_pert] - - def __call__(self, data: PerturbationAnndataPair) -> float: - cats_sorted = sorted([c for c in data.perts if c != data.control_pert]) - - # 2. build centroids - ad_real_cent = self._centroid_ann( - adata=data.real, - category_key=data.pert_col, - control_pert=data.control_pert, - embed_key=self.embed_key, - ) - ad_pred_cent = self._centroid_ann( - adata=data.pred, - category_key=data.pert_col, - control_pert=data.control_pert, - embed_key=self.embed_key, - ) - - # 3. cluster real once - real_key = "real_clusters" - self._cluster_leiden( - ad_real_cent, self.real_resolution, real_key, self.n_neighbors - ) - # reorder rows to match cats_sorted without using DataFrame.set_index (type stubs issue) - idx_real = ( - pd.Series(np.arange(ad_real_cent.n_obs), - index=ad_real_cent.obs[data.pert_col].to_numpy()) - .loc[cats_sorted] - .to_numpy() - ) - ad_real_cent = ad_real_cent[idx_real] - real_labels = pd.Categorical(ad_real_cent.obs[real_key]) - - # 4. sweep predicted resolutions - best_score = 0.0 - idx_pred = ( - pd.Series(np.arange(ad_pred_cent.n_obs), - index=ad_pred_cent.obs[data.pert_col].to_numpy()) - .loc[cats_sorted] - .to_numpy() - ) - ad_pred_cent = ad_pred_cent[idx_pred] - for r in self.pred_resolutions: - pred_key = f"pred_clusters_{r}" - self._cluster_leiden(ad_pred_cent, r, pred_key, self.n_neighbors) - pred_labels = pd.Categorical(ad_pred_cent.obs[pred_key]) - score = self._score(real_labels, pred_labels, self.metric) # type: ignore - best_score = max(best_score, score) - - return float(best_score) + cats = ad \ No newline at end of file From 2ad897aa31d1d10bfa5f267aa0cc01aae4e8a941 Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:43:11 -0800 Subject: [PATCH 7/9] reset src/cell_eval/metrics/_anndata.py --- src/cell_eval/metrics/_anndata.py | 50 ++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/src/cell_eval/metrics/_anndata.py b/src/cell_eval/metrics/_anndata.py index 93ff8e0..c026e3c 100644 --- a/src/cell_eval/metrics/_anndata.py +++ b/src/cell_eval/metrics/_anndata.py @@ -352,4 +352,52 @@ def _centroid_ann( if issparse(feats): feats = feats.toarray() # type: ignore - cats = ad \ No newline at end of file + cats = adata.obs[category_key].values + uniq, inv = np.unique(cats, return_inverse=True) # type: ignore + centroids = np.zeros((uniq.size, feats.shape[1]), dtype=feats.dtype) # type: ignore + + for i, cat in enumerate(uniq): + mask = cats == cat + if np.any(mask): + centroids[i] = feats[mask].mean(axis=0) # type: ignore + + adc = ad.AnnData(X=centroids) + adc.obs[category_key] = uniq + return adc[adc.obs[category_key] != control_pert] + + def __call__(self, data: PerturbationAnndataPair) -> float: + cats_sorted = sorted([c for c in data.perts if c != data.control_pert]) + + # 2. build centroids + ad_real_cent = self._centroid_ann( + adata=data.real, + category_key=data.pert_col, + control_pert=data.control_pert, + embed_key=self.embed_key, + ) + ad_pred_cent = self._centroid_ann( + adata=data.pred, + category_key=data.pert_col, + control_pert=data.control_pert, + embed_key=self.embed_key, + ) + + # 3. cluster real once + real_key = "real_clusters" + self._cluster_leiden( + ad_real_cent, self.real_resolution, real_key, self.n_neighbors + ) + ad_real_cent.obs = ad_real_cent.obs.set_index(data.pert_col).loc[cats_sorted] + real_labels = pd.Categorical(ad_real_cent.obs[real_key]) + + # 4. sweep predicted resolutions + best_score = 0.0 + ad_pred_cent.obs = ad_pred_cent.obs.set_index(data.pert_col).loc[cats_sorted] + for r in self.pred_resolutions: + pred_key = f"pred_clusters_{r}" + self._cluster_leiden(ad_pred_cent, r, pred_key, self.n_neighbors) + pred_labels = pd.Categorical(ad_pred_cent.obs[pred_key]) + score = self._score(real_labels, pred_labels, self.metric) # type: ignore + best_score = max(best_score, score) + + return float(best_score) From 2a2434e182a8151228c6b9c1e5d8cb1dc1a30c0f Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:45:13 -0800 Subject: [PATCH 8/9] reset src/cell_eval/_evaluator.py --- src/cell_eval/_evaluator.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index f692732..0e86ce8 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -89,7 +89,7 @@ def __init__( allow_discrete=allow_discrete, ) - if True: + if skip_de: self.de_comparison = None else: self.de_comparison = _build_de_comparison( @@ -107,8 +107,6 @@ def __init__( self.outdir = outdir self.prefix = prefix - # Store the requested number of threads for metrics usage - self.num_threads = num_threads def compute( self, @@ -119,23 +117,16 @@ def compute( write_csv: bool = True, break_on_error: bool = False, ) -> tuple[pl.DataFrame, pl.DataFrame]: - # Merge provided configs with num_threads for top_k_accuracy only - merged_metric_configs: dict[str, dict[str, Any]] = {} - if metric_configs: - merged_metric_configs.update(metric_configs) - pipeline = MetricPipeline( profile=profile, - metric_configs=merged_metric_configs, + metric_configs=metric_configs, break_on_error=break_on_error, ) - if skip_metrics is not None: pipeline.skip_metrics(skip_metrics) pipeline.compute_de_metrics(self.de_comparison) pipeline.compute_anndata_metrics(self.anndata_pair) results = pipeline.get_results() - print("results", results) agg_results = pipeline.get_agg_results() if write_csv: @@ -198,7 +189,7 @@ def _convert_to_normlog( Will skip if the input is not integer data. """ - if True: # TODO: francis fix + if guess_is_lognorm(adata=adata, validate=not allow_discrete): logger.info( "Input is found to be log-normalized already - skipping transformation." ) From 8354ed438cd1fec1ffbcdbc9a027dd3994bdd1af Mon Sep 17 00:00:00 2001 From: Francis Chalissery Date: Thu, 15 Jan 2026 15:49:04 -0800 Subject: [PATCH 9/9] clean up src/cell_eval/_cli/_run.py --- src/cell_eval/_cli/_run.py | 28 ++-------------------------- 1 file changed, 2 insertions(+), 26 deletions(-) diff --git a/src/cell_eval/_cli/_run.py b/src/cell_eval/_cli/_run.py index d9251c1..20e5401 100644 --- a/src/cell_eval/_cli/_run.py +++ b/src/cell_eval/_cli/_run.py @@ -1,5 +1,4 @@ import argparse as ap -import concurrent.futures as cf import importlib.metadata import logging import os @@ -116,12 +115,6 @@ def parse_args_run(parser: ap.ArgumentParser): default=10, help="k for top_k_accuracy (number of nearest neighbors) [default: %(default)s]", ) - parser.add_argument( - "--ctrl-barcode-col", - type=str, - default=None, - help="Column name for control barcode matching in top_k_accuracy (optional)", - ) parser.add_argument( "--version", action="version", @@ -156,7 +149,7 @@ def run_evaluation(args: ap.Namespace): else {} ) - # Always pass top-k and ctrl_barcode_col for top_k_accuracy + # Always pass top-k for top_k_accuracy metric_kwargs.setdefault("top_k_accuracy", {})["k"] = args.topk skip_metrics = args.skip_metrics.split(",") if args.skip_metrics else None @@ -172,7 +165,7 @@ def run_evaluation(args: ap.Namespace): f"Number of celltypes in real and pred anndata must match: {len(real_split)} != {len(pred_split)}" ) - def _run_for_celltype(ct: str): + for ct in real_split.keys(): real_ct = real_split[ct] pred_ct = pred_split[ct] @@ -197,23 +190,6 @@ def _run_for_celltype(ct: str): skip_metrics=skip_metrics, basename="results.csv", ) - return ct - - max_workers = args.num_threads if args.num_threads and args.num_threads > 1 else 1 - if max_workers == 1: - for ct in real_split.keys(): - _run_for_celltype(ct) - else: - logger.info(f"Parallelizing over celltypes with {max_workers} threads") - with cf.ThreadPoolExecutor(max_workers=max_workers) as executor: - futures = {executor.submit(_run_for_celltype, ct): ct for ct in real_split.keys()} - for fut in cf.as_completed(futures): - ct = futures[fut] - try: - fut.result() - logger.info(f"Completed evaluation for celltype: {ct}") - except Exception as e: - logger.exception(f"Evaluation failed for celltype {ct}: {e}") else: evaluator = MetricsEvaluator(