fix(tools): raise when filter_rank_genes_groups genes are missing from adata

docs/release-notes/4497.fix.md

@@ -0,0 +1 @@
+Fix {func}`~scanpy.tl.filter_rank_genes_groups` silently scoring genes missing from `adata` (e.g. removed after {func}`~scanpy.tl.rank_genes_groups`) with the values of the last gene; it now raises an informative error. Also fix re-filtering an already filtered result when `pts=True` was used {smaller}`K R Roy`

src/scanpy/tools/_rank_genes_groups.py

@@ -1089,6 +1089,45 @@ def _calc_frac(x: NDArray[np.number] | CSBase, /) -> NDArray[np.float64]:
     return n_nonzero / x.shape[0]
 
 
+def _get_var_idx(
+    adata: AnnData,
+    var_names: pd.api.extensions.ExtensionArray,
+    *,
+    key: str,
+    group: str,
+    use_raw: bool,
+) -> NDArray[np.intp]:
+    """Get positions of ranked genes, raising if any are missing (see #1487).
+
+    `get_indexer` returns -1 for missing names, which would silently select
+    the last gene. NaN names (from an already filtered result) are skipped.
+    """
+    var_idx = (adata.raw if use_raw else adata).var_names.get_indexer(var_names)
+    if not (missing := (var_idx < 0) & pd.notna(var_names)).any():
+        return var_idx
+    msg = (
+        f"Genes {list(var_names[missing][:5])!r} from `adata.uns[{key!r}]` "
+        f"(group {group!r}) are not in `adata{'.raw' if use_raw else ''}.var_names`. "
+        "Were genes removed after `rank_genes_groups`, or was it run with a different "
+        "`use_raw` or layer? Rerun `rank_genes_groups` on this object, or pass the "
+        "same `use_raw` / `use` it was run with."
+    )
+    raise KeyError(msg)
+
+
+def _get_cached_frac(
+    pts: pd.DataFrame, group: str, var_names: pd.api.extensions.ExtensionArray
+) -> NDArray[np.floating]:
+    """Look up cached `pts`/`pts_rest` fractions for ranked genes.
+
+    NaN names (from an already filtered result) get NaN, so they stay filtered out.
+    """
+    frac = np.full(len(var_names), np.nan)
+    present = pd.notna(var_names)
+    frac[present] = pts[group].loc[var_names[present]].to_numpy()
+    return frac
+
+
 @_doc_params(
     use=doc_use(
         "Which matrix to compute fractions and fold changes from. "
@@ -1257,18 +1296,20 @@ def filter_rank_genes_groups(  # noqa: PLR0912, PLR0913
         var_names = gene_names[cluster].array
 
         if not use_logfolds or not use_fraction:
-            var_idx = (adata.raw if use_raw else adata).var_names.get_indexer(var_names)
+            var_idx = _get_var_idx(
+                adata, var_names, key=key, group=cluster, use_raw=use_raw
+            )
             sub_x = x[:, var_idx]
             in_group = (adata.obs[groupby] == cluster).to_numpy()
             x_in = sub_x[in_group]
             x_out = sub_x[~in_group]
 
         if use_fraction:
-            fraction_in_cluster_matrix.loc[:, cluster] = (
-                adata.uns[key]["pts"][cluster].loc[var_names].array
+            fraction_in_cluster_matrix.loc[:, cluster] = _get_cached_frac(
+                adata.uns[key]["pts"], cluster, var_names
             )
-            fraction_out_cluster_matrix.loc[:, cluster] = (
-                adata.uns[key]["pts_rest"][cluster].loc[var_names].array
+            fraction_out_cluster_matrix.loc[:, cluster] = _get_cached_frac(
+                adata.uns[key]["pts_rest"], cluster, var_names
             )
         else:
             fraction_in_cluster_matrix.loc[:, cluster] = _calc_frac(x_in)

tests/test_filter_rank_genes_groups.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import numpy as np
+import pandas as pd
 import pytest
 
 from scanpy.tools import filter_rank_genes_groups, rank_genes_groups
@@ -134,3 +135,34 @@ def test_filter_rank_genes_groups_out() -> None:
             np.array(adata.uns[key]["names"].tolist()),
             np.array(res["names"].tolist()),
         )
+
+
+def test_filter_rank_genes_groups_missing_genes() -> None:
+    """Genes removed after ranking must raise, not be scored as the last gene (#1487)."""
+    adata = pbmc68k_reduced().raw.to_adata()
+    rank_genes_groups(adata, "bulk_labels", method="wilcoxon", n_genes=5)
+    removed = list(adata.uns["rank_genes_groups"]["names"]["CD14+ Monocyte"][:3])
+    adata = adata[:, ~adata.var_names.isin(removed)].copy()
+
+    with pytest.raises(KeyError, match=r"are not in `adata\.var_names`"):
+        filter_rank_genes_groups(adata)
+
+
+@pytest.mark.parametrize("pts", [False, True], ids=["computed", "cached_pts"])
+def test_filter_rank_genes_groups_refilter(*, pts: bool) -> None:
+    """Filtering an already filtered result (with NaN names) still works.
+
+    With ``pts=True`` the fractions come from the cached ``pts``/``pts_rest``.
+    """
+    adata = pbmc68k_reduced()
+    rank_genes_groups(adata, "bulk_labels", method="wilcoxon", n_genes=5, pts=pts)
+    filter_rank_genes_groups(adata)
+    filter_rank_genes_groups(
+        adata, key="rank_genes_groups_filtered", key_added="refiltered"
+    )
+    first = pd.DataFrame(adata.uns["rank_genes_groups_filtered"]["names"])
+    second = pd.DataFrame(adata.uns["refiltered"]["names"])
+    # the first filter must leave NaN names, or this test checks nothing
+    assert first.isna().to_numpy().any()
+    # re-applying the same thresholds keeps the same genes
+    pd.testing.assert_frame_equal(first, second)