Merge pull request #234 from LeonHafner/leonhafner/fix_log2fc

fix: consume pdex's log2_fold_change directly (stop double-logging)

Author: LeonHafner

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "cell-eval"
-version = "0.7.0"
+version = "0.7.1"
 description = "Evaluation metrics for single-cell perturbation predictions"
 readme = "README.md"
 authors = [
@@ -11,7 +11,7 @@ authors = [
 requires-python = ">=3.11"
 dependencies = [
     "igraph>=0.11.8",
-    "pdex>=0.2.0",
+    "pdex>=0.2.2",
     "polars>=1.30.0",
     "pyyaml>=6.0.2",
     "scanpy>=1.10.3",

src/cell_eval/_baseline.py

@@ -77,7 +77,7 @@ def build_base_mean_adata(
 
     if output_path is not None:
         logger.info(f"Saving baseline data to {output_path}")
-        baseline_adata.write_h5ad(output_path)  # type: ignore[invalid-argument-type]
+        baseline_adata.write_h5ad(output_path)  # ty: ignore[invalid-argument-type]
 
     if output_de_path is not None:
         logger.info("Calculating differential expression")

src/cell_eval/_cli/_prep.py

@@ -226,7 +226,7 @@ def strip_anndata(
 
         # Write the h5ad file
         logger.info(f"Writing h5ad output to {tmp_h5ad}")
-        minimal.write_h5ad(tmp_h5ad)  # type: ignore[invalid-argument-type]
+        minimal.write_h5ad(tmp_h5ad)  # ty: ignore[invalid-argument-type]
 
         # Zstd compress the h5ad file (will create pred.h5ad.zst)
         logger.info(f"Zstd compressing {tmp_h5ad}")

src/cell_eval/_types/_anndata.py

@@ -125,9 +125,9 @@ def _bulk_anndata(
 
         # Create a polars dataframe with the groupby key
         frame = pl.DataFrame(
-            matrix,
+            matrix,  # ty: ignore[invalid-argument-type]
         ).with_columns(
-            groupby_key=adata.obs[groupby_key].to_numpy(str),  # type: ignore
+            groupby_key=adata.obs[groupby_key].to_numpy(str),
         )
 
         # Pseudobulk (mean) the dataframe by the groupby key

src/cell_eval/_types/_de.py

@@ -17,7 +17,6 @@ def initialize_de_comparison(
     pred: pl.DataFrame,
     target_col: str = "target",
     feature_col: str = "feature",
-    fold_change_col: str = "fold_change",
     log2_fold_change_col: str = "log2_fold_change",
     abs_log2_fold_change_col: str = "abs_log2_fold_change",
     pvalue_col: str = "p_value",
@@ -27,7 +26,6 @@ def initialize_de_comparison(
         DEResults,
         target_col=target_col,
         feature_col=feature_col,
-        fold_change_col=fold_change_col,
         log2_fold_change_col=log2_fold_change_col,
         abs_log2_fold_change_col=abs_log2_fold_change_col,
     )
@@ -47,7 +45,6 @@ class DEResults:
     # Column names configuration
     target_col: str = "target"
     feature_col: str = "feature"
-    fold_change_col: str = "fold_change"
     log2_fold_change_col: str = "log2_fold_change"
     abs_log2_fold_change_col: str = "abs_log2_fold_change"
     pvalue_col: str = "p_value"
@@ -58,7 +55,7 @@ def __post_init__(self) -> None:
         required_cols = {
             self.target_col,
             self.feature_col,
-            self.fold_change_col,
+            self.log2_fold_change_col,
             self.pvalue_col,
             self.fdr_col,
         }
@@ -67,7 +64,6 @@ def __post_init__(self) -> None:
             raise ValueError(f"Missing required columns: {missing}")
 
         numeric_cols = [
-            self.fold_change_col,
             self.pvalue_col,
             self.fdr_col,
             self.log2_fold_change_col,
@@ -80,31 +76,32 @@ def __post_init__(self) -> None:
         ]
 
         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:
+        lfc_num_null = self.data.filter(
+            pl.col(self.log2_fold_change_col).is_null()
+        ).height
+        lfc_num_inf = self.data.filter(
+            pl.col(self.log2_fold_change_col).is_infinite()
+        ).height
+        lfc_num_nan = self.data.filter(
+            pl.col(self.log2_fold_change_col).is_nan()
+        ).height
+        if lfc_num_null > 0:
             logger.warning(
-                f"Identified {fc_num_null} missing fold change values ({self.name})"
+                f"Identified {lfc_num_null} missing log2 fold change values ({self.name})"
             )
-        if fc_num_inf > 0:
+        if lfc_num_inf > 0:
             logger.warning(
-                f"Identified {fc_num_inf} infinite fold change values ({self.name})"
+                f"Identified {lfc_num_inf} infinite log2 fold change values ({self.name})"
             )
-        if fc_num_nan > 0:
+        if lfc_num_nan > 0:
             logger.warning(
-                f"Identified {fc_num_nan} NaN fold change values ({self.name})"
+                f"Identified {lfc_num_nan} NaN log2 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:
+        # Derive abs(log2_fold_change) if not already provided.
+        if self.abs_log2_fold_change_col not in self.data.columns:
             self.data = self.data.with_columns(
-                pl.col(self.fold_change_col)
-                .log(base=2)
-                .alias(self.log2_fold_change_col)
-                .fill_nan(0.0)
-            ).with_columns(
                 pl.col(self.log2_fold_change_col)
                 .abs()
                 .alias(self.abs_log2_fold_change_col)
@@ -153,7 +150,10 @@ def get_top_genes(
         # Set FDR threshold if not provided
         fdr_threshold = fdr_threshold if fdr_threshold is not None else 0.05
 
-        descending = sort_by in {DESortBy.FOLD_CHANGE, DESortBy.ABS_FOLD_CHANGE}
+        descending = sort_by in {
+            DESortBy.LOG2_FOLD_CHANGE,
+            DESortBy.ABS_LOG2_FOLD_CHANGE,
+        }
 
         # Create a rank matrix where each row is the ordinal rank of a gene and each column is a perturbation.
         # The rank is sensitive to the sort-by column and is computed post-filtering for FDR.
@@ -219,7 +219,7 @@ def compute_overlap(
         k: int | None,
         metric: Literal["overlap", "precision"] = "overlap",
         fdr_threshold: float | None = None,
-        sort_by: DESortBy = DESortBy.ABS_FOLD_CHANGE,
+        sort_by: DESortBy = DESortBy.ABS_LOG2_FOLD_CHANGE,
     ) -> dict[str, float]:
         """
         Compute overlap metrics across perturbations.

src/cell_eval/_types/_enums.py

@@ -4,8 +4,8 @@
 class DESortBy(enum.Enum):
     """Sorting options for differential expression results."""
 
-    FOLD_CHANGE = "log2_fold_change"
-    ABS_FOLD_CHANGE = "abs_log2_fold_change"
+    LOG2_FOLD_CHANGE = "log2_fold_change"
+    ABS_LOG2_FOLD_CHANGE = "abs_log2_fold_change"
     PVALUE = "p_value"
     FDR = "fdr"
 

src/cell_eval/metrics/_de.py

@@ -13,7 +13,7 @@ def de_overlap_metric(
     k: int | None,
     metric: Literal["precision", "overlap"] = "overlap",
     fdr_threshold: float = 0.05,
-    sort_by: DESortBy = DESortBy.ABS_FOLD_CHANGE,
+    sort_by: DESortBy = DESortBy.ABS_LOG2_FOLD_CHANGE,
 ) -> dict[str, float]:
     """Compute overlap between real and predicted DE genes.
 
@@ -124,7 +124,9 @@ def __call__(self, data: DEComparison) -> dict[str, float]:
                 suffix="_pred",
                 how="left",
             )
-            .with_columns(pl.col(f"{data.real.fold_change_col}_pred").fill_null(0.0))
+            .with_columns(
+                pl.col(f"{data.real.log2_fold_change_col}_pred").fill_null(0.0)
+            )
         )
 
         for row in (
@@ -133,8 +135,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.log2_fold_change_col).cast(pl.Float64),
+                    pl.col(f"{data.real.log2_fold_change_col}_pred").cast(pl.Float64),
                     method="spearman",
                 ).alias("spearman_corr"),
             )

src/cell_eval/utils.py

@@ -41,7 +41,7 @@ def guess_is_lognorm(
     if isinstance(adata.X, csr_matrix) or isinstance(adata.X, csc_matrix):
         frac, _ = np.modf(adata.X.data)
     elif adata.is_view:
-        frac, _ = np.modf(adata.X.toarray())  # type: ignore[unresolved-attribute]
+        frac, _ = np.modf(adata.X.toarray())  # ty: ignore[unresolved-attribute]
     elif adata.X is None:
         raise ValueError("adata.X is None")
     else:
@@ -60,8 +60,8 @@ def guess_is_lognorm(
         max_val = adata.X.max()
         min_val = adata.X.min()
     else:
-        max_val = float(np.max(adata.X))  # type: ignore[no-matching-overload]
-        min_val = float(np.min(adata.X))  # type: ignore[no-matching-overload]
+        max_val = float(np.max(adata.X))  # ty: ignore[no-matching-overload]
+        min_val = float(np.min(adata.X))  # ty: ignore[no-matching-overload]
 
     # Validate range
     if min_val < 0:

tests/test_de_float_types.py

@@ -10,13 +10,13 @@ def test_de_spearman_lfc_mixed_float_types() -> None:
         {
             "target": ["pert1", "pert1", "pert2", "pert2"],
             "feature": ["gene1", "gene2", "gene1", "gene2"],
-            "fold_change": [1.5, 2.0, 0.5, 1.2],
+            "log2_fold_change": [1.5, 2.0, 0.5, 1.2],
             "p_value": [0.01, 0.02, 0.03, 0.04],
             "fdr": [0.01, 0.02, 0.03, 0.04],
         }
     )
 
-    pred_df = real_df.with_columns(pl.col("fold_change").cast(pl.Float64))
+    pred_df = real_df.with_columns(pl.col("log2_fold_change").cast(pl.Float64))
 
     comparison = DEComparison(
         real=DEResults(real_df, name="real"),
@@ -27,3 +27,25 @@ def test_de_spearman_lfc_mixed_float_types() -> None:
 
     assert isinstance(result, dict)
     assert all(isinstance(value, (int, float)) for value in result.values())
+
+
+def test_de_results_preserves_negative_log2_fold_change() -> None:
+    """`abs_log2_fold_change` must be |log2_fold_change|; negatives must not be zeroed.
+
+    Pre-fix, `__post_init__` applied `log(base=2).fill_nan(0.0)` to the values, which
+    silently zeroed every downregulated gene.
+    """
+    df = pl.DataFrame(
+        {
+            "target": ["p", "p", "p", "p"],
+            "feature": ["g1", "g2", "g3", "g4"],
+            "log2_fold_change": [-1.0, 0.0, 1.0, 2.0],
+            "p_value": [0.01, 0.01, 0.01, 0.01],
+            "fdr": [0.01, 0.01, 0.01, 0.01],
+        }
+    )
+    de = DEResults(df, name="real")
+    lfc = de.data["log2_fold_change"].cast(pl.Float64).to_list()
+    abs_lfc = de.data["abs_log2_fold_change"].cast(pl.Float64).to_list()
+    assert lfc == [-1.0, 0.0, 1.0, 2.0]
+    assert abs_lfc == [1.0, 0.0, 1.0, 2.0]

tests/test_eval.py

@@ -196,7 +196,7 @@ def test_unknown_alternative_de_metric():
             control_pert=CONTROL_VAR,
             pert_col=PERT_COL,
             outdir=OUTDIR,
-            de_method="unknown",  # type: ignore[unknown-argument]
+            de_method="unknown",  # ty: ignore[unknown-argument]
         ).compute()