Fix mt.ora selecting the bottom instead of the top n_up features

CHANGELOG.md

@@ -8,6 +8,14 @@ and this project adheres to [Semantic Versioning][].
 [keep a changelog]: https://keepachangelog.com/en/1.0.0/
 [semantic versioning]: https://semver.org/spec/v2.0.0.html
 
+## 2.2.1
+
+### Bugfixes
+- Fixed `mt.ora` selecting the bottom `nvar - n_up` features instead of the top `n_up` ones. Ranks are ascending, so the threshold is now `nvar - n_up`, restoring the behaviour of `decoupler<2`. With the default `n_up` (top 5%) this selected 95% of the features as observed, which also made `mt.ora` fail with `invalid contingency table` whenever `nvar - n_up` exceeded `n_bg` (#346)
+- Fixed `n_bm` in `mt.ora` selecting the bottom `n_bm - 1` features instead of `n_bm`
+- `mt.ora` now validates that `n_up` and `n_bm` do not overlap and that `n_bg` is large enough for the number of selected features, instead of failing inside the Fisher exact test
+- `mt.query_set` now raises an informative error when `n_bg` is smaller than the number of features to test, instead of failing inside `scipy.stats.fisher_exact`
+
 ## 2.2.0
 
 ### Added

pyproject.toml

@@ -4,7 +4,7 @@ requires = [ "hatchling" ]
 
 [project]
 name = "decoupler"
-version = "2.2.0"
+version = "2.2.1"
 description = "Python package to perform enrichment analysis from omics data."
 readme = "README.md"
 license = { file = "LICENSE" }

src/decoupler/mt/_ora.py

@@ -231,7 +231,7 @@ def _func_ora(
 
     n_up
         Number of top-ranked features, based on their magnitude, to select as observed features.
-        If ``None``, the top 5% of positive features are selected.
+        If ``None``, the top 5% of features are selected.
     n_bm
         Number of bottom-ranked features, based on their magnitude, to select as observed features.
     %(n_bg)s
@@ -261,6 +261,13 @@ def _func_ora(
     assert isinstance(n_up, int | float) and n_up > 0, "n_up must be numeric and > 0"
     assert isinstance(n_bm, int | float) and n_bm >= 0, "n_bm must be numeric and positive"
     assert isinstance(n_bg, int | float) and n_bg >= 0, "n_bg must be numeric and positive"
+    n_up, n_bm = int(np.ceil(n_up)), int(np.ceil(n_bm))
+    assert n_up + n_bm <= nvar, f"For nvar={nvar}, n_up={n_up} and n_bm={n_bm} overlap, decrease any of them"
+    assert n_bg == 0 or n_bg >= n_up + n_bm, (
+        f"n_bg={n_bg} must be larger or equal than the number of selected features n_up + n_bm={n_up + n_bm}, "
+        "otherwise the contingency table is invalid. Increase n_bg, decrease n_up or n_bm, "
+        "or set n_bg=None to use a feature specific background"
+    )
     m = f"ora - calculating {nsrc} scores across {nobs} observations with n_up={n_up}, n_bm={n_bm}, n_bg={n_bg}"
     _log(m, level="info", verbose=verbose)
     es = np.zeros((nobs, nsrc))
@@ -273,7 +280,8 @@ def _func_ora(
             row = mat[i]
         # Find ranks
         row = sts.rankdata(row, method="ordinal")
-        row = ranks[(row > n_up) | (row < n_bm)]
+        # Ranks are ascending, the top n_up features are the ones with the largest ranks
+        row = ranks[(row > (nvar - n_up)) | (row <= n_bm)]
         es[i], pv[i] = _runora(
             row=set(row), ranks=set(ranks), cnct=cnct, starts=starts, offsets=offsets, n_bg=n_bg, ha_corr=ha_corr
         )

src/decoupler/mt/_query_set.py

@@ -77,6 +77,11 @@ def query_set(
             d = len(set_d)
         else:
             d = int(n_bg - a - b - c)
+            assert d >= 0, (
+                f"n_bg={n_bg} must be larger or equal than the number of features to test for source={source} "
+                f"({a + b + c}), otherwise the contingency table is invalid. Increase n_bg or set n_bg=None "
+                "to use a feature specific background"
+            )
         od = _oddsr(a=a, b=b, c=c, d=d, ha_corr=ha_corr, log=True)
         _, pv = sts.fisher_exact([[a, b], [c, d]], alternative=alternative)
         df.append([source, od, pv])

tests/mt/test_ora.py

@@ -1,6 +1,7 @@
 import math
 
 import numpy as np
+import pandas as pd
 import pytest
 import scipy.sparse as sps
 import scipy.stats as sts
@@ -40,7 +41,7 @@ def test_runora(
     cnct, starts, offsets = idxmat
     row = sts.rankdata(X[0], method="ordinal")
     ranks = np.arange(row.size, dtype=np.int_)
-    row = ranks[(row > 2) | (row < 0)]
+    row = ranks[row > (row.size - 2)]
     es, pv = dc.mt._ora._runora.py_func(
         row=set(row),
         ranks=set(ranks),
@@ -77,7 +78,7 @@ def test_func_ora(
     rnk = set(ranks)
     for i in range(st_es.shape[0]):
         row = sts.rankdata(X[i], method="ordinal")
-        row = set(ranks[row > n_up])
+        row = set(ranks[row > (X.shape[1] - n_up)])
         for j in range(st_es.shape[1]):
             fset = dc.pp.net._getset(cnct=cnct, starts=starts, offsets=offsets, j=j)
             fset = set(fset)
@@ -100,3 +101,70 @@ def test_func_ora(
             st_es[i, j], _ = np.log(es)
     assert np.isclose(dc_es, st_es).all()
     assert np.isclose(dc_pv, st_pv).all()
+
+
+@pytest.mark.parametrize(
+    "n_up,n_bm,obs_idxs",
+    [
+        [5, 0, [15, 16, 17, 18, 19]],
+        [1, 0, [19]],
+        [20, 0, list(range(20))],
+        [3, 2, [0, 1, 17, 18, 19]],
+        [0.5, 0, [19]],
+    ],
+)
+def test_func_ora_selection(
+    n_up,
+    n_bm,
+    obs_idxs,
+):
+    # Row values are ascending, so the top n_up features are the last ones
+    nvar = 20
+    X = np.arange(nvar, dtype=float).reshape(1, nvar)
+    var = np.array([f"G{i:02d}" for i in range(nvar)])
+    net = pd.DataFrame(
+        {
+            "source": ["S1"] * 5,
+            "target": var[[15, 16, 17, 18, 19]],
+            "weight": [1.0] * 5,
+        }
+    )
+    sources, cnct, starts, offsets = dc.pp.idxmat(features=var, net=net, verbose=False)
+    dc_es, dc_pv = dc.mt._ora._func_ora(
+        mat=X, cnct=cnct, starts=starts, offsets=offsets, n_up=n_up, n_bm=n_bm, n_bg=None
+    )
+    # Build the expected contingency table from the features that should be selected
+    row = set(obs_idxs)
+    fset = {15, 16, 17, 18, 19}
+    a = len(row & fset)
+    b = len(fset - row)
+    c = len(row - fset)
+    d = nvar - len(row | fset)
+    st_es = np.log(((a + 0.5) * (d + 0.5)) / ((b + 0.5) * (c + 0.5)))
+    st_pv = sts.fisher_exact([[a, b], [c, d]])[1]
+    assert np.isclose(dc_es[0, 0], st_es)
+    assert np.isclose(dc_pv[0, 0], st_pv)
+
+
+def test_func_ora_validate(
+    mat,
+    idxmat,
+):
+    X, obs, var = mat
+    cnct, starts, offsets = idxmat
+    nvar = X.shape[1]
+    kwargs = {"mat": X, "cnct": cnct, "starts": starts, "offsets": offsets}
+    with pytest.raises(AssertionError, match="overlap"):
+        dc.mt._ora._func_ora(**kwargs, n_up=nvar, n_bm=1, n_bg=None)
+    with pytest.raises(AssertionError, match="contingency table is invalid"):
+        dc.mt._ora._func_ora(**kwargs, n_up=5, n_bm=0, n_bg=4)
+
+
+def test_ora_wide():
+    # Selecting the top 5% of a wide matrix must not exceed n_bg
+    adata, net = dc.ds.toy(nobs=5, nvar=1_000, seed=42, verbose=False)
+    dc.mt.ora(adata, net, tmin=3, n_bg=100)
+    es = adata.obsm["score_ora"].values
+    pv = adata.obsm["padj_ora"].values
+    assert np.isfinite(es).all()
+    assert ((pv >= 0) & (pv <= 1)).all()

tests/mt/test_query_set.py

@@ -1,4 +1,5 @@
 import pandas as pd
+import pytest
 
 import decoupler as dc
 
@@ -12,3 +13,11 @@ def test_query_set(
     cols = {"source", "stat", "pval", "padj"}
     assert cols.issubset(df.columns)
     df = dc.mt.query_set(features=ft, net=net, n_bg=None, tmin=0)
+
+
+def test_query_set_n_bg(
+    net,
+):
+    ft = set(net["target"])
+    with pytest.raises(AssertionError, match="contingency table is invalid"):
+        dc.mt.query_set(features=ft, net=net, n_bg=2, tmin=0)