diff --git a/CHANGELOG.md b/CHANGELOG.md index d06f530..882cc97 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 27fb84a..203c92f 100644 --- a/pyproject.toml +++ b/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" } diff --git a/src/decoupler/mt/_ora.py b/src/decoupler/mt/_ora.py index 58dc93b..5f3d2d4 100644 --- a/src/decoupler/mt/_ora.py +++ b/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 ) diff --git a/src/decoupler/mt/_query_set.py b/src/decoupler/mt/_query_set.py index 3941bb8..1f284bc 100644 --- a/src/decoupler/mt/_query_set.py +++ b/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]) diff --git a/tests/mt/test_ora.py b/tests/mt/test_ora.py index f0f41cb..6b328e4 100644 --- a/tests/mt/test_ora.py +++ b/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() diff --git a/tests/mt/test_query_set.py b/tests/mt/test_query_set.py index cb99691..26d924a 100644 --- a/tests/mt/test_query_set.py +++ b/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)