feat: add illicofor rank_genes_groups

pyproject.toml

@@ -98,7 +98,8 @@ scanorama = [ "scanorama" ]
 scrublet = [ "scikit-image>=0.25" ]
 # highly_variable_genes method 'seurat_v3'
 skmisc = [ "scikit-misc>=0.5.1" ]
-scanpy2 = [ "igraph>=0.10.8", "scikit-misc>=0.5.1" ]
+illico = [ "illico>=0.6" ]
+scanpy2 = [ "igraph>=0.10.8", "scanpy[illico]", "scikit-misc>=0.5.1" ]
 
 [dependency-groups]
 dev = [
@@ -108,6 +109,7 @@ dev = [
 test = [
     "scanpy[dask-ml]",
     "scanpy[dask]",
+    "scanpy[illico]",
     "scanpy[leiden]",
     "scanpy[plotting]",
     "scanpy[scrublet]",

src/scanpy/_settings/presets.py

@@ -33,7 +33,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"]
 

src/scanpy/tools/_rank_genes_groups.py

@@ -7,14 +7,15 @@
 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
 
 from .. import _utils
 from .. import logging as logg
 from .._compat import CSBase
-from .._settings import Default
+from .._settings import Default, Preset
 from .._settings.presets import DETest
 from .._utils import (
     _numba_thread_limit,
@@ -28,7 +29,6 @@
     from collections.abc import Generator, Iterable
     from typing import Literal
 
-    from anndata import AnnData
     from numpy.typing import NDArray
 
 
@@ -47,6 +47,35 @@ def _select_top_n(scores: NDArray, n_top: int):
     return global_indices
 
 
+def _illico_results_to_iter(
+    illico_df: pd.DataFrame,
+    groups_order: NDArray,
+    ireference: int | None,
+    *,
+    copy_pvalues: bool,
+) -> Generator[tuple[int, NDArray[np.floating], NDArray[np.floating]]]:
+    """Yield per-group ``(index, z, p)`` tuples from illico's long-form output.
+
+    illico returns a DataFrame with a 2-level MultiIndex ``(pert, feature)``
+    and columns including ``z_score`` and ``p_value``. We stream one group
+    at a time via `pandas.Series.loc`, trusting illico_df groups are ordered
+    by ``var_name``.
+    """
+    ref_label = None if ireference is None else groups_order[ireference]
+    z_series = illico_df["z_score"]
+    p_series = illico_df["p_value"]
+    illico_groups = set(illico_df.index.unique(level="pert"))
+    return (
+        (
+            group_index,
+            z_series.loc[group_name].to_numpy(),
+            p_series.loc[group_name].to_numpy(copy=copy_pvalues),
+        )
+        for group_index, group_name in enumerate(groups_order)
+        if group_name != ref_label and group_name in illico_groups
+    )
+
+
 @njit
 def rankdata(data: NDArray[np.number]) -> NDArray[np.float64]:
     """Parallelized version of scipy.stats.rankdata."""
@@ -141,6 +170,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
@@ -424,6 +454,40 @@ def logreg(
             if len(self.groups_order) <= 2:
                 break
 
+    def illico(
+        self, *, tie_correct: bool, corr_method: _CorrMethod
+    ) -> Generator[tuple[int, NDArray[np.floating], NDArray[np.floating]], None, None]:
+        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"),
+                    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,
+            groups=self.groups_order,
+        )
+        return _illico_results_to_iter(
+            illico_df,
+            self.groups_order,
+            self.ireference,
+            # p-values are altered by this correction method
+            copy_pvalues=corr_method == "benjamini-hochberg",
+        )
+
     def compute_statistics(  # noqa: PLR0912
         self,
         method: DETest,
@@ -441,8 +505,12 @@ def compute_statistics(  # noqa: PLR0912
             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 method == "wilcoxon":
-            generate_test_results = self.wilcoxon(tie_correct=tie_correct)
+        elif "wilcoxon" in method:
+            generate_test_results = (
+                self.illico(tie_correct=tie_correct, corr_method=corr_method)
+                if "illico" in method
+                else 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":
@@ -451,7 +519,6 @@ def compute_statistics(  # noqa: PLR0912
         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])
 
@@ -651,6 +718,8 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
         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
+        if settings.preset is Preset.ScanpyV2Preview:
+            method = "wilcoxon_illico"
 
     mask_var = _check_mask(adata, mask_var, "var")
 

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"

tests/test_rank_genes_groups.py

@@ -17,13 +17,14 @@
 from scanpy._utils.random import _LegacyRng
 from scanpy.get import rank_genes_groups_df
 from scanpy.tools import rank_genes_groups
-from scanpy.tools._rank_genes_groups import _RankGenes
+from scanpy.tools._rank_genes_groups import _illico_results_to_iter, _RankGenes
 from testing.scanpy._helpers import random_mask
 from testing.scanpy._helpers.data import pbmc68k_reduced
+from testing.scanpy._pytest.marks import needs
 from testing.scanpy._pytest.params import ARRAY_TYPES, ARRAY_TYPES_MEM
 
 if TYPE_CHECKING:
-    from collections.abc import Callable
+    from collections.abc import Callable, Sequence
     from typing import Any, Literal
 
     from numpy.lib.npyio import NpzFile
@@ -78,6 +79,25 @@ def get_true_scores(method: Literal["t-test", "wilcoxon"]) -> Expected:
     return Expected(names=expected["names"].astype("T"), scores=expected["scores"])
 
 
+def get_illico_results_df(
+    n_groups: int, n_genes: int, *, seed: int = 0
+) -> pd.DataFrame:
+    """Synthetic illico-shaped output used by the _illico_results_to_iter tests.
+
+    Features are in deliberately non-ascending to test ``var_name`` ordering.
+    """
+    rng = np.random.default_rng(seed)
+    groups = [f"g{i}" for i in range(n_groups)]
+    features = [f"f{n_genes - 1 - j}" for j in range(n_genes)]  # reversed -> not sorted
+    return pd.DataFrame(
+        {
+            "z_score": rng.normal(size=n_groups * n_genes),
+            "p_value": rng.uniform(size=n_groups * n_genes),
+        },
+        index=pd.MultiIndex.from_product([groups, features], names=["pert", "feature"]),
+    )
+
+
 # TODO: Make dask compatible
 @pytest.mark.parametrize("method", ["t-test", "wilcoxon"])
 @pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM)
@@ -352,6 +372,113 @@ def test_mask_not_equal():
     assert not np.array_equal(no_mask, with_mask)
 
 
+@pytest.mark.parametrize(
+    ("groups_order", "ireference", "expected_indices"),
+    [
+        (["g0", "g1", "g2"], None, [0, 1, 2]),
+        (["g0", "g1", "g2"], 1, [0, 2]),
+    ],
+    ids=["vs_rest", "vs_reference"],
+)
+@pytest.mark.parametrize("corr_method", ["benjamini-hochberg", "bonferroni"])
+def test_illico_iter(
+    groups_order: list[str],
+    ireference: int | None,
+    expected_indices: list[int],
+    corr_method: Literal["benjamini-hochberg", "bonferroni"],
+):
+    df = get_illico_results_df(n_groups=3, n_genes=4)
+    feature_order = df.index.unique(level="feature")
+    out = list(
+        _illico_results_to_iter(
+            df,
+            np.array(groups_order),
+            ireference,
+            copy_pvalues=corr_method == "benjamini-hochberg",
+        )
+    )
+    assert sorted(t[0] for t in out) == sorted(expected_indices)
+    for gi, z, p in out:
+        sub = df.xs(groups_order[gi], level=0).reindex(feature_order)
+        np.testing.assert_array_equal(z, sub["z_score"].to_numpy())
+        np.testing.assert_array_equal(p, sub["p_value"].to_numpy())
+
+
+@pytest.mark.parametrize("corr_method", ["benjamini-hochberg", "bonferroni"])
+@pytest.mark.parametrize("test", ["ovo", "ovr"])  # pairwise or vs. rest
+@pytest.mark.parametrize(
+    "mean_in_log_space", [True, False], ids=["log_space_mean", "linear_space_mean"]
+)
+@pytest.mark.parametrize(
+    "tie_correct", [True, False], ids=["tie_correct", "no_tie_correct"]
+)
+@pytest.mark.parametrize("groups", [["CD14+ Monocyte", "Dendritic"], "all"])
+@pytest.mark.filterwarnings("ignore:invalid value encountered:RuntimeWarning")
+@needs.illico
+def test_illico(
+    test: Literal["ovo", "ovr"],
+    corr_method: Literal["benjamini-hochberg", "bonferroni"],
+    subtests: pytest.Subtests,
+    groups: Literal["all"] | Sequence[str],
+    *,
+    mean_in_log_space: bool,
+    tie_correct: bool,
+):
+
+    pbmc = pbmc68k_reduced()
+    pbmc.raw.X.sum_duplicates()
+    pbmc.raw.X.sort_indices()
+    pbmc_illico = pbmc.copy()
+
+    reference = pbmc.obs["bulk_labels"].iloc[0] if test == "ovo" else "rest"
+    sc.tl.rank_genes_groups(
+        pbmc_illico,
+        groupby="bulk_labels",
+        method="wilcoxon_illico",
+        reference=reference if test == "ovo" else "rest",
+        n_genes=pbmc.n_vars,
+        tie_correct=tie_correct,
+        corr_method=corr_method,
+        mean_in_log_space=mean_in_log_space,
+        groups=groups,
+    )
+
+    sc.tl.rank_genes_groups(
+        pbmc,
+        groupby="bulk_labels",
+        method="wilcoxon",
+        reference=reference if test == "ovo" else "rest",
+        n_genes=pbmc.n_vars,
+        tie_correct=tie_correct,
+        corr_method=corr_method,
+        mean_in_log_space=mean_in_log_space,
+        groups=groups,
+    )
+    scanpy_results = pbmc.uns["rank_genes_groups"]
+    illico_results = pbmc_illico.uns["rank_genes_groups"]
+    assert set(illico_results.keys()) == set(scanpy_results.keys()), (
+        "Output keys do not match Scanpy's output format."
+    )
+
+    for k, ref in scanpy_results.items():
+        with subtests.test(k):
+            if k in ["params", "names"]:
+                # We can skip names ordering check as if incorrect, other values will mismatch
+                continue
+            res = np.array(illico_results[k].tolist())
+            ref_arr = np.array(ref.tolist())
+            mask = np.isfinite(ref_arr) * np.isfinite(
+                res
+            )  # Mask to ignore inf values in the comparison
+            np.testing.assert_allclose(
+                ref_arr[mask],
+                res[mask],
+                rtol=0,
+                atol=1e-6,
+                err_msg=f"Mismatch in '{k}' values between asymptotic_wilcoxon and Scanpy outputs.",
+            )
+
+
 @pytest.mark.parametrize(
     ("mean_in_log_space", "expected_logfc"),
     [
@@ -363,10 +490,18 @@ def test_mask_not_equal():
         (False, -1.0),
     ],
 )
-@pytest.mark.parametrize("method", ["wilcoxon", "t-test", "t-test_overestim_var"])
+@pytest.mark.parametrize(
+    "method",
+    [
+        "wilcoxon",
+        "t-test",
+        "t-test_overestim_var",
+        pytest.param("wilcoxon_illico", marks=needs.illico),
+    ],
+)
 def test_mean_in_log_space(
     expected_logfc: float,
-    method: Literal["wilcoxon", "t-test", "t-test_overestim_var"],
+    method: Literal["wilcoxon", "wilcoxon_illico", "t-test", "t-test_overestim_var"],
     *,
     mean_in_log_space: bool,
 ):