feat: add illicofor rank_genes_groups

docs/release-notes/4037.feat.md

@@ -0,0 +1 @@
+Add `mean_in_log_space` argument to {func}`scanpy.tl.rank_genes_groups` for customizing how log-fold-change is calculated {user}`ilan-gold`

pyproject.toml

@@ -51,7 +51,7 @@ classifiers = [
 ]
 dynamic = [ "version" ]
 dependencies = [
-    "anndata>=0.10.8",
+    "anndata>=0.11",
     "certifi",
     "fast-array-utils[accel,sparse]>=1.4",
     "h5py>=3.11",
@@ -97,6 +97,7 @@ scanorama = [ "scanorama" ]
 scrublet = [ "scikit-image>=0.23.1" ]
 # highly_variable_genes method 'seurat_v3'
 skmisc = [ "scikit-misc>=0.5.1" ]
+illico = [ "illico>=0.5.1" ]
 scanpy2 = [ "igraph>=0.10.8", "scikit-misc>=0.5.1" ]
 
 [dependency-groups]
@@ -107,6 +108,7 @@ dev = [
 test = [
     "scanpy[dask-ml]",
     "scanpy[dask]",
+    "scanpy[illico]",
     "scanpy[leiden]",
     "scanpy[plotting]",
     "scanpy[scrublet]",

src/scanpy/_settings/presets.py

@@ -31,7 +31,9 @@
 ]
 
 
-type DETest = Literal["logreg", "t-test", "wilcoxon", "t-test_overestim_var"]
+type DETest = Literal[
+    "logreg", "t-test", "wilcoxon", "wilcoxon_illico", "t-test_overestim_var"
+]
 type HVGFlavor = Literal["seurat", "cell_ranger", "seurat_v3", "seurat_v3_paper"]
 type LeidenFlavor = Literal["leidenalg", "igraph"]
 
@@ -81,6 +83,7 @@ class PcaPreset(NamedTuple):
 class RankGenesGroupsPreset(NamedTuple):
     method: DETest
     mask_var: str | None
+    mean_in_log_space: bool
 
 
 class ScalePreset(NamedTuple):
@@ -185,9 +188,11 @@ def pca() -> Mapping[Preset, PcaPreset]:
     def rank_genes_groups() -> Mapping[Preset, RankGenesGroupsPreset]:
         """Correlation method for :func:`~scanpy.tl.rank_genes_groups`."""
         return {
-            Preset.ScanpyV1: RankGenesGroupsPreset(method="t-test", mask_var=None),
+            Preset.ScanpyV1: RankGenesGroupsPreset(
+                method="t-test", mask_var=None, mean_in_log_space=True
+            ),
             Preset.ScanpyV2Preview: RankGenesGroupsPreset(
-                method="wilcoxon", mask_var=None
+                method="wilcoxon", mask_var=None, mean_in_log_space=False
             ),
         }
 

src/scanpy/tools/_rank_genes_groups.py

@@ -7,6 +7,7 @@
 import numba
 import numpy as np
 import pandas as pd
+from anndata import AnnData
 from fast_array_utils.numba import njit
 from fast_array_utils.stats import mean_var
 from scipy import sparse
@@ -28,7 +29,6 @@
     from collections.abc import Generator, Iterable
     from typing import Literal
 
-    from anndata import AnnData
     from numpy.typing import NDArray
 
 
@@ -141,6 +141,7 @@ def __init__(
             self.expm1_func = lambda x: np.expm1(x * np.log(base))
         else:
             self.expm1_func = np.expm1
+        self.group_col = adata.obs[groupby].array
 
         self.groups_order, self.groups_masks_obs = _utils.select_groups(
             adata, groups, groupby
@@ -202,7 +203,7 @@ def __init__(
         self.grouping_mask = adata.obs[groupby].isin(self.groups_order)
         self.grouping = adata.obs.loc[self.grouping_mask, groupby]
 
-    def _basic_stats(self) -> None:
+    def _basic_stats(self, *, exponentiate_values: bool = False) -> None:
         """Set self.{means,vars,pts}{,_rest} depending on X."""
         n_genes = self.X.shape[1]
         n_groups = self.groups_masks_obs.shape[0]
@@ -218,6 +219,8 @@ def _basic_stats(self) -> None:
         else:
             mask_rest = self.groups_masks_obs[self.ireference]
             x_rest = self.X[mask_rest]
+            if exponentiate_values:
+                x_rest = self.expm1_func(x_rest)
             self.means[self.ireference], self.vars[self.ireference] = mean_var(
                 x_rest, axis=0, correction=1
             )
@@ -231,6 +234,8 @@ def _basic_stats(self) -> None:
 
         for group_index, mask_obs in enumerate(self.groups_masks_obs):
             x_mask = self.X[mask_obs]
+            if exponentiate_values:
+                x_mask = self.expm1_func(x_mask)
 
             if self.comp_pts:
                 self.pts[group_index] = get_nonzeros(x_mask) / x_mask.shape[0]
@@ -245,6 +250,8 @@ def _basic_stats(self) -> None:
             if self.ireference is None:
                 mask_rest = ~mask_obs
                 x_rest = self.X[mask_rest]
+                if exponentiate_values:
+                    x_rest = self.expm1_func(x_rest)
                 (
                     self.means_rest[group_index],
                     self.vars_rest[group_index],
@@ -260,8 +267,6 @@ def t_test(
     ) -> Generator[tuple[int, NDArray[np.floating], NDArray[np.floating]], None, None]:
         from scipy import stats
 
-        self._basic_stats()
-
         for group_index, (mask_obs, mean_group, var_group) in enumerate(
             zip(self.groups_masks_obs, self.means, self.vars, strict=True)
         ):
@@ -313,8 +318,6 @@ def wilcoxon(
     ) -> Generator[tuple[int, NDArray[np.floating], NDArray[np.floating]], None, None]:
         from scipy import stats
 
-        self._basic_stats()
-
         n_genes = self.X.shape[1]
         # First loop: Loop over all genes
         if self.ireference is not None:
@@ -422,27 +425,88 @@ def logreg(
             if len(self.groups_order) <= 2:
                 break
 
-    def compute_statistics(  # noqa: PLR0912
+    def compute_statistics(  # noqa: PLR0912, PLR0915
         self,
         method: DETest,
         *,
-        corr_method: _CorrMethod = "benjamini-hochberg",
-        n_genes_user: int | None = None,
-        rankby_abs: bool = False,
-        tie_correct: bool = False,
+        corr_method: _CorrMethod,
+        n_genes_user: int | None,
+        rankby_abs: bool,
+        tie_correct: bool,
+        mean_in_log_space: bool,
         **kwds,
     ) -> None:
         if method in {"t-test", "t-test_overestim_var"}:
+            self._basic_stats(exponentiate_values=False)
             generate_test_results = self.t_test(method)
-        elif method == "wilcoxon":
-            generate_test_results = self.wilcoxon(tie_correct=tie_correct)
+            if not mean_in_log_space:
+                # If we are not exponentiating after the mean aggregation, we need to recalculate the stats.
+                self._basic_stats(exponentiate_values=True)
+        elif "wilcoxon" in method:
+            if "illico" in method:
+                from illico import asymptotic_wilcoxon
+
+                illico_df = asymptotic_wilcoxon(
+                    AnnData(
+                        X=self.X,
+                        var=pd.DataFrame(index=self.var_names),
+                        obs=pd.DataFrame(
+                            index=pd.RangeIndex(self.X.shape[0]).astype("str"),
+                            # This self.group_col means illico will run tests against *all* data
+                            # instead of what's in self.groups_order as controlled by the `groups` arg.
+                            # TODO: Only run the subset once illico supports a `groups` argument
+                            data={"group": self.group_col},
+                        ),
+                    ),
+                    reference=self.groups_order[self.ireference]
+                    if self.ireference is not None
+                    else None,
+                    group_keys="group",
+                    return_as_scanpy=False,
+                    is_log1p=True,
+                    tie_correct=tie_correct,
+                    use_continuity=False,
+                    alternative="two-sided",
+                    use_rust=False,
+                )
+                # Generate a lookup of category -> result excluding the refernece if it is present.
+                generate_test_results_map = {
+                    group_cat: (
+                        group["z_score"].to_numpy(copy=True),
+                        group["p_value"].to_numpy(copy=True),
+                    )
+                    for (_, group) in illico_df.groupby(level="pert")
+                    if (
+                        group_cat := np.unique(
+                            group.index.get_level_values("pert").to_numpy(copy=True)
+                        ).item()
+                    )
+                    != (
+                        None
+                        if self.ireference is None
+                        else self.groups_order[self.ireference]
+                    )
+                }
+                # Create the iterator that is expected by the other method-branches.
+                groups_order_list = self.groups_order.tolist()
+                generate_test_results = (
+                    (
+                        groups_order_list.index(group_cat),
+                        *generate_test_results_map[group_cat],
+                    )
+                    for group_cat in self.groups_order
+                    if group_cat in generate_test_results_map
+                )
+            else:
+                generate_test_results = self.wilcoxon(tie_correct=tie_correct)
+            # If we're not exponentiating after the mean aggregation, then do it now.
+            self._basic_stats(exponentiate_values=not mean_in_log_space)
         elif method == "logreg":
             generate_test_results = self.logreg(**kwds)
 
         self.stats = None
 
         n_genes = self.X.shape[1]
-
         for group_index, scores, pvals in generate_test_results:
             group_name = str(self.groups_order[group_index])
 
@@ -482,9 +546,12 @@ def compute_statistics(  # noqa: PLR0912
                     mean_rest = self.means_rest[group_index]
                 else:
                     mean_rest = self.means[self.ireference]
-                foldchanges = (self.expm1_func(mean_group) + 1e-9) / (
-                    self.expm1_func(mean_rest) + 1e-9
-                )  # add small value to remove 0's
+                foldchanges = (
+                    (self.expm1_func(mean_group) + 1e-9)
+                    / (self.expm1_func(mean_rest) + 1e-9)
+                    if mean_in_log_space
+                    else (mean_group + 1e-9) / (mean_rest + 1e-9)
+                )  # add small value to avoid zeros
                 self.stats[group_name, "logfoldchanges"] = np.log2(
                     foldchanges[global_indices]
                 )
@@ -512,9 +579,12 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
     corr_method: _CorrMethod = "benjamini-hochberg",
     tie_correct: bool = False,
     layer: str | None = None,
+    mean_in_log_space: bool | Default = Default(
+        preset=("rank_genes_groups", "mean_in_log_space")
+    ),
     **kwds,
 ) -> AnnData | None:
-    """Rank genes for characterizing groups.
+    r"""Rank genes for characterizing groups.
 
     Expects logarithmized data.
 
@@ -575,6 +645,11 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
         The key in `adata.uns` information is saved to.
     copy
         Whether to copy `adata` or modify it inplace.
+    mean_in_log_space
+        Whether to do :math:`\log(\operatorname{mean}(e^x))` (`False`)
+        or :math:`\log(e^{\operatorname{mean}(x)})` (`True`).
+        The former is accurate, while the latter is a faster approximation
+        that underestimates this accurate result in the presence of many outliers.
     kwds
         Are passed to test methods. Currently this affects only parameters that
         are passed to :class:`sklearn.linear_model.LogisticRegression`.
@@ -597,7 +672,7 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
         Structured array to be indexed by group id storing the log2
         fold change for each gene for each group. Ordered according to
         scores. Only provided if method is 't-test' like.
-        Note: this is an approximation calculated from mean-log values.
+        Note: if `mean_in_log_space=True`, this is an approximation calculated from mean-log values.
     `adata.uns['rank_genes_groups' | key_added]['pvals']` : structured :class:`numpy.ndarray` (dtype `float`)
         p-values.
     `adata.uns['rank_genes_groups' | key_added]['pvals_adj']` : structured :class:`numpy.ndarray` (dtype `float`)
@@ -627,6 +702,8 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
 
     if isinstance(mask_var, Default):
         mask_var = settings.preset.rank_genes_groups.mask_var
+    if isinstance(mean_in_log_space, Default):
+        mean_in_log_space = settings.preset.rank_genes_groups.mean_in_log_space
     if method is None or isinstance(method, Default):
         method = settings.preset.rank_genes_groups.method
 
@@ -716,6 +793,7 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
             n_genes_user=n_genes_user,
             rankby_abs=rankby_abs,
             tie_correct=tie_correct,
+            mean_in_log_space=mean_in_log_space,
             **kwds,
         )
 

src/testing/scanpy/_pytest/marks.py

@@ -41,6 +41,7 @@ def _generate_next_value_(
     skimage = "scikit-image"
     skmisc = "scikit-misc"
     zarr = auto()
+    illico = auto()
     # external
     bbknn = auto()
     harmony = "harmonyTS"