feat: add illicofor rank_genes_groups (#4038)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Philipp A. <flying-sheep@web.de>
Co-authored-by: Copilot <copilot@github.com>
Co-authored-by: Ilay Kavitzky <ilay.kavitzky@gmail.com>
Co-authored-by: Cuiwei Gao <48gaocuiwei@gmail.com>
Co-authored-by: Jhonatan Felix <108437587+JhonatanFelix@users.noreply.github.com>
Co-authored-by: Zach Boldyga <zboldyga@users.noreply.github.com>

Author: 8 people

docs/release-notes/3653.feat.md

@@ -0,0 +1 @@
+Add {attr}`scanpy.settings.preset` setting with the new preset {attr}`~scanpy.Preset.ScanpyV2Preview` {smaller}`P Angerer`

docs/release-notes/3653.feature.md

@@ -0,0 +1,3 @@
+Introduce blazing-fast [`illico`][] as a new `wilcoxon` `method` in {func}`~scanpy.tl.rank_genes_groups`. Use it via either `method="wilcoxon_illico"` or via {attr}`~scanpy.Preset.ScanpyV2Preview` (preferred) {smaller}`I Gold`
+
+[`illico`]: https://github.com/remydubois/illico

docs/release-notes/4038.feat.md

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

pyproject.toml

@@ -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"]
 
@@ -233,7 +235,11 @@ def draw_graph() -> Mapping[Preset, DrawGraphPreset]:
 
     @preset_property
     def rank_genes_groups() -> Mapping[Preset, RankGenesGroupsPreset]:
-        """Correlation method for :func:`~scanpy.tl.rank_genes_groups`."""
+        """
+        Correlation method for :func:`~scanpy.tl.rank_genes_groups`.
+
+        V2 `method` preset is for `illico` implementation of `wilcoxon` test.
+        """
         return {
             Preset.ScanpyV1: RankGenesGroupsPreset(
                 method="t-test", mask_var=None, mean_in_log_space=True

src/scanpy/_settings/presets.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 .._compat import CSBase, warn
+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])
 
@@ -649,8 +716,18 @@ def rank_genes_groups(  # noqa: PLR0912, PLR0913, PLR0915
         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 scanpy presets are used for v2, use illico - prevents the presets from showing the `wilcoxon_illico` method and allows us to silently replace `wilcoxon`'s implementation.
     if method is None or isinstance(method, Default):
         method = settings.preset.rank_genes_groups.method
+        if settings.preset is Preset.ScanpyV2Preview:
+            method = "wilcoxon_illico"
+    # Otherwise, nudge people to use the presets.
+    elif "illico" in method:
+        msg = (
+            "`wilcoxon_illico` flavor will be removed in scanpy 2.0 and be simply the new `wilcoxon` implementation."
+            "To remove theis warning, you can locally do `with sc.settings.override(preset=sc.Preset.ScanpyV2Preview)`."
+        )
+        warn(msg, DeprecationWarning)
 
     mask_var = _check_mask(adata, mask_var, "var")
 

src/scanpy/tools/_rank_genes_groups.py

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