From 8969dc930d931cd3c42bc9ebeae51e0aac8aef0d Mon Sep 17 00:00:00 2001 From: Andrew Newman Date: Sun, 8 Nov 2026 20:46:25 +1000 Subject: [PATCH] Remove lrfeatures and put into lr_results. Fix for #369. --- stlearn/pl/cci_plot_helpers.py | 6 +-- stlearn/tl/cci/__init__.py | 4 +- stlearn/tl/cci/_results.py | 38 +++++++++----- stlearn/tl/cci/analysis.py | 23 ++++++--- stlearn/tl/cci/perm_utils.py | 12 ++--- stlearn/tl/cci/permutation.py | 2 + tests/acceptance/run_cci.py | 9 ++-- tests/io/test_cci_round_trip.py | 2 +- tests/tl/test_cci.py | 1 + tests/tl/test_lr_table.py | 87 +++++++++++++++++++++++++++++---- 10 files changed, 141 insertions(+), 43 deletions(-) diff --git a/stlearn/pl/cci_plot_helpers.py b/stlearn/pl/cci_plot_helpers.py index 83b3b9f6..09946cfa 100644 --- a/stlearn/pl/cci_plot_helpers.py +++ b/stlearn/pl/cci_plot_helpers.py @@ -37,7 +37,8 @@ def lr_scatter( figsize: tuple | None = None, show_all: bool = False, ): - """`feature` of the LR pairs against their rank; `data` is the LR table.""" + """`feature` of the LR pairs against their rank; `data` is the LR table, whose + var holds the summary and the expression features side by side.""" lr_df = data.var if max_text > len(lr_df): @@ -55,8 +56,7 @@ def lr_scatter( highlight_lrs = highlight_lrs[0:max_text] lrs = lr_df.index.to_numpy().astype(str)[0:n_top] - lr_features = data.uns["lrfeatures"] - lr_df = pd.concat([lr_df, lr_features], axis=1).loc[lrs, :] + lr_df = lr_df.loc[lrs, :] if feature not in lr_df.columns: raise ValueError(f"Inputted {feature}; must be one of {list(lr_df.columns)}") diff --git a/stlearn/tl/cci/__init__.py b/stlearn/tl/cci/__init__.py index f74cb1e9..4d4270a4 100644 --- a/stlearn/tl/cci/__init__.py +++ b/stlearn/tl/cci/__init__.py @@ -1,10 +1,12 @@ -from ._results import LR_RESULTS, lr_result +from ._results import LR_FEATURES, LR_QUANTILES, LR_RESULTS, lr_result from .analysis import adj_pvals, run, run_cci, run_lr_go from .edges import get_edges from .gridding import create_grid, grid from .ligand_receptors import LRPairs, load_lrs, lr_pairs, read_lr_pairs __all__ = [ + "LR_FEATURES", + "LR_QUANTILES", "LR_RESULTS", "LRPairs", "adj_pvals", diff --git a/stlearn/tl/cci/_results.py b/stlearn/tl/cci/_results.py index b38c2b3b..3578f45e 100644 --- a/stlearn/tl/cci/_results.py +++ b/stlearn/tl/cci/_results.py @@ -1,11 +1,15 @@ """st.tl.cci.run returns an LR table: the obs/spots/cells as rows and LR pairs as -columns with the most significant first. Its var is the LR summary, each of LR_RESULTS -is a layer, and its uns holds lrfeatures. +columns with the most significant first. Its var is the LR summary, with each +pair's expression features as the LR_FEATURES columns, and each of LR_RESULTS is a +layer. The layers are CSC matrices with one structure: an entry for each spot tested for an LR, those with a score above 0 for it. Elsewhere a result has its default, a score of 0 or a p-value of 1, which lr_result fills in. +varm[LR_QUANTILES] holds the expression quantiles of each pair's ligand and +receptor. The LR pair on the var axis stays aligned with the layers' columns. + st.tl.cci.run_cci returns cell type interactions and st.tl.cci.run_lr_go's GO terms. The table it ran on keeps only what is per spot: the neighbours and cell type counts. """ @@ -19,7 +23,19 @@ from statsmodels.stats.multitest import multipletests LR_RESULTS = ("lr_scores", "p_vals", "p_adjs", "neg_log10_p_adjs", "lr_sig_scores") -DEFAULTS = {"p_vals": 1.0, "p_adjs": 1.0} # a result's value where an LR was not tested +# A result's value where an LR was not tested +DEFAULTS = {"p_vals": 1.0, "p_adjs": 1.0} +# Statistics per LR pair: the average of the L and R non-zero medians and +# the proportion of spots with a count of zero. +# Rank is an order of values n-1 to 0 (highest to lowest). +# Median rank is ordered (ascending) using the pair's non-zero median. +# Proportion rank is ordered (descending) using the pair's proportions. +# Mean rank is the average of the median and proportion rank. +LR_FEATURES = ("nonzero-median", "zero-prop", "median_rank", "prop_rank", "mean_rank") +# lr_quantiles is stored in varm - the LR pairs by quantile, of each ligand's ('L_q') +# and receptor's ('R_q') expression quantiles. +LR_QUANTILES = "lr_quantiles" + type Matrix = sparse.spmatrix | np.ndarray @@ -171,14 +187,12 @@ def store_lr_results( padjs: Matrix, pval_adj_cutoff: float, ) -> None: - """Writes the LR summary to the LR table's var and the LR_RESULTS layers, CSC - on the tested entries, with LRs ordered by their number of significant spots. - - The inputs' columns follow `lr_results.var_names` on the way in. Previous - results are replaced; the rest of the table (X, varm, other layers) moves with - the summary, as scanpy does. 'source', where each pair came from, is kept; any - other summary columns, such as run_cci's, are dropped: they describe the - previous results. + """Writes the LR summary to the LR table's var and the LR_RESULTS layers, and + orders the LRs by their number of significant spots. + + The inputs' columns are in `lr_results.var_names` order. Reordering var + reorders the rest of the table with it (X, varm, other layers), as scanpy does. + LR_FEATURES and 'source' var columns are kept - other var columns are dropped. """ scores = scored_entries(lr_scores) adj, ps = on_structure(padjs, scores), on_structure(pvals, scores) @@ -197,6 +211,8 @@ def store_lr_results( ) if "source" in lr_results.var: summary.insert(0, "source", lr_results.var["source"]) + features = [col for col in LR_FEATURES if col in lr_results.var] + summary = summary.join(lr_results.var[features]) order: npt.NDArray[np.intp] = np.argsort(-summary["n_spots_sig"].to_numpy()) for key in LR_RESULTS: diff --git a/stlearn/tl/cci/analysis.py b/stlearn/tl/cci/analysis.py index 9a37440b..34677fb2 100644 --- a/stlearn/tl/cci/analysis.py +++ b/stlearn/tl/cci/analysis.py @@ -13,6 +13,8 @@ from ...utils import _set_n_threads from ._results import ( + LR_FEATURES, + LR_QUANTILES, check_lr_results, correct_pvals, lr_result, @@ -149,6 +151,11 @@ def run( Summary of significant spots detected per LR: 'n_spots', 'n_spots_sig' and 'n_spots_sig_pval'; and with LRPairs, 'source', where each pair came from. + varm + 'lr_quantiles': the expression quantiles of each pair's ligand + ('L_0.5' ... 'L_1.0') and receptor ('R_0.5' ... 'R_1.0'), from + which the test picked its expression-matched background genes and + storing statistics in LR_FEATURES. layers 'lr_scores', 'lr_sig_scores', 'p_vals', 'p_adjs' and 'neg_log10_p_adjs', spots by LR pairs. 'lr_scores' is the raw @@ -160,7 +167,7 @@ def run( matrices leave out; read a layer with st.tl.cci.lr_result, which fills them in, one LR at a time for a large sample. uns - Returns 'lr_species' and 'lrfeatures' + 'lr_species', the species of lrs. With save_bg 'lrs_to_bg' and 'lr_spot_indices'. Relevant information stored in adata: adata.obsp['spot_neighbours'] @@ -274,14 +281,14 @@ def run( neg_binom=neg_binom, random_state=random_state, ) - lr_results = anndata.AnnData( - obs=pandas.DataFrame(index=adata.obs_names), - var=pandas.DataFrame(index=tested.names()), - ) - - lr_results.var["source"] = tested.to_frame()["source"].to_numpy(dtype=str) + features = tests.lr_features.loc[tested.names()] + var = pandas.DataFrame( + {"source": tested.to_frame()["source"].to_numpy(dtype=str)}, + index=tested.names(), + ).join(features[list(LR_FEATURES)]) + lr_results = anndata.AnnData(obs=pandas.DataFrame(index=adata.obs_names), var=var) + lr_results.varm[LR_QUANTILES] = features.drop(columns=list(LR_FEATURES)) lr_results.uns["lr_species"] = lrs.species - lr_results.uns["lrfeatures"] = tests.lr_features padjs = correct_pvals(lr_scores, tests.pvals, "spot", adj_method) store_lr_results(lr_results, lr_scores, tests.pvals, padjs, pval_adj_cutoff) diff --git a/stlearn/tl/cci/perm_utils.py b/stlearn/tl/cci/perm_utils.py index c40a3239..50a149f8 100644 --- a/stlearn/tl/cci/perm_utils.py +++ b/stlearn/tl/cci/perm_utils.py @@ -8,6 +8,7 @@ from numba import njit, prange from scipy import sparse +from ._results import LR_FEATURES from .ligand_receptors import LRPairs @@ -218,11 +219,9 @@ def _gen_rand_pairs(genes1, genes2, n_pairs, seed): def get_lr_features(lr_expr: pd.DataFrame | AnnData, lrs: LRPairs, quantiles): - """Gets expression features of LR pairs; nonzero-median, zero-prop, quantiles. - A row per pair, indexed by its name. - - lr_expr is the LR genes' expression, spots by genes: a DataFrame, or an - AnnData of those genes, whose sparse X is used as it is, never densified. + """Get expression features for a given LR pair (per row). The LR_FEATURES columns + then the 'L_q' and 'R_q' quantiles. Rows are ordered by mean_rank; lowest + expressed first, not as in lrs. """ quantiles = np.array(quantiles) if isinstance(lr_expr, AnnData): @@ -268,8 +267,7 @@ def get_lr_features(lr_expr: pd.DataFrame | AnnData, lrs: LRPairs, quantiles): prop_ranks = [np.where(prop_order == i)[0][0] for i in range(len(lrs))] mean_ranks = np.array([median_ranks, prop_ranks]).mean(axis=0) - # Saving the lrfeatures... - cols = ["nonzero-median", "zero-prop", "median_rank", "prop_rank", "mean_rank"] + cols = LR_FEATURES lr_features_data = { cols[0]: np.array(lr_median_means, dtype=np.float64), cols[1]: np.array(lr_prop_means, dtype=np.float64), diff --git a/stlearn/tl/cci/permutation.py b/stlearn/tl/cci/permutation.py index a06cef23..0123e6b0 100644 --- a/stlearn/tl/cci/permutation.py +++ b/stlearn/tl/cci/permutation.py @@ -66,6 +66,8 @@ def perform_spot_testing( get_lrs_scores returns with sparse_result) or dense. lr_scores has a column for each of `lrs`, in the order of lrs.names(), as get_lrs_scores gives them. The p-values come back sparse, an entry for each. + + Returns lr_features in expression-rank order. """ quantiles = np.array(quantiles) scores = scored_entries(lr_scores) diff --git a/tests/acceptance/run_cci.py b/tests/acceptance/run_cci.py index a89ba28b..33b7d061 100644 --- a/tests/acceptance/run_cci.py +++ b/tests/acceptance/run_cci.py @@ -195,9 +195,12 @@ def main() -> None: result.uns["neighbours_indices"] = graph.indices if args.use_label: result.obs["cci_het"] = numpy.asarray(adata.obsm["cci_het"]).ravel() - # The layout of baselines written before the LR table. - result.uns["lr_summary"] = lr_results.var - result.uns["lrfeatures"] = lr_results.uns["lrfeatures"] + # For backward compatibility - layout of baselines before the LR table. + features = list(stlearn.tl.cci.LR_FEATURES) + result.uns["lr_summary"] = lr_results.var.drop(columns=features) + result.uns["lrfeatures"] = lr_results.var[features].join( + lr_results.varm[stlearn.tl.cci.LR_QUANTILES] + ) result.uns["metrics"] = json.dumps(metrics) pathlib.Path(args.result).parent.mkdir(parents=True, exist_ok=True) result.write_h5ad(args.result, compression="gzip") diff --git a/tests/io/test_cci_round_trip.py b/tests/io/test_cci_round_trip.py index 16f4f316..4409acf6 100644 --- a/tests/io/test_cci_round_trip.py +++ b/tests/io/test_cci_round_trip.py @@ -109,7 +109,7 @@ def test_every_key(self): pd.testing.assert_frame_equal( getattr(read[name], slot), getattr(table, slot) ) - for slot in ("obsm", "obsp", "layers", "uns"): + for slot in ("obsm", "varm", "obsp", "layers", "uns"): written, back = getattr(table, slot), getattr(read[name], slot) self.assertEqual(sorted(back), sorted(written), msg=f"{name}.{slot}") for key in written: diff --git a/tests/tl/test_cci.py b/tests/tl/test_cci.py index 6e12855c..ae21b64d 100644 --- a/tests/tl/test_cci.py +++ b/tests/tl/test_cci.py @@ -22,6 +22,7 @@ class TestCCI(unittest.TestCase): """Tests for `stlearn` CCI capability.""" + all_edges: list[Any] def setUp(self) -> None: diff --git a/tests/tl/test_lr_table.py b/tests/tl/test_lr_table.py index 8c73ef52..15c466ea 100644 --- a/tests/tl/test_lr_table.py +++ b/tests/tl/test_lr_table.py @@ -1,12 +1,6 @@ #!/usr/bin/env python -"""st.tl.cci.run's results in a table of their own. - -st.tl.cci.run returns an AnnData with the same spots as rows and the LR pairs as -columns: var is the LR summary, each per-spot result is a layer. The expression -table keeps only what is per spot, so a store can open it without the results. -Everything that reads or adds to the results takes the table as lr_results=. -""" +"""st.tl.cci.run's results in a table of their own.""" import unittest from unittest import mock @@ -24,7 +18,11 @@ from stlearn.tl.cci import LR_RESULTS, lr_result from tests.utils import CCI_LRS, cci_table -SUMMARY = ["source", "n_spots", "n_spots_sig", "n_spots_sig_pval"] +FEATURES = ["nonzero-median", "zero-prop", "median_rank", "prop_rank", "mean_rank"] +SUMMARY = ["source", "n_spots", "n_spots_sig", "n_spots_sig_pval", *FEATURES] +# varm key of the ligands' and receptors' expression quantiles, run's defaults. +QUANTILES = "lr_quantiles" +RUN_QUANTILES = [0.5, 0.75, 0.85, 0.9, 0.95, 0.97, 0.98, 0.99, 0.995, 0.9975, 0.999, 1] RUN_CCI_UNS = [ "lr_cci_cell_type", "lr_cci_raw_cell_type", @@ -81,7 +79,62 @@ def test_run_returns_the_lr_table(self): np.testing.assert_array_equal( lr.var["n_spots_sig"], (lr_result(lr, "p_adjs") < 0.05).sum(axis=0) ) - self.assertIn("lrfeatures", lr.uns) + self.assertEqual(list(lr.uns), ["lr_species"]) + self.assertEqual(list(lr.varm), [QUANTILES]) + + def _genes(self, name): + """The expression of pair `name`'s ligand and receptor, spots by 2.""" + ligand, receptor = name.split("_") + return self.adata[:, [ligand, receptor]].X.toarray() + + def test_expression_features_are_var_columns(self): + """Each pair's expression features sit beside its summary in var, so they + are read by name like any other column and can never be out of step with + the layers' columns (what uns['lrfeatures'] was).""" + lr = self.lr + features = lr.var[FEATURES] + self.assertFalse(features.isna().any().any()) + for name in lr.var_names: + with self.subTest(lr=name): + expr = self._genes(name) + nonzero_medians = [ + np.quantile(gene[gene > 0], 0.5, method="nearest") + for gene in expr.T + ] + self.assertEqual( + features.loc[name, "nonzero-median"], np.mean(nonzero_medians) + ) + self.assertEqual(features.loc[name, "zero-prop"], (expr == 0).mean()) + # The ranks run over the tested pairs, 0 the lowest expressed; so the most + # significant pair, LigA_RecA, with its co-expressed corner, ranks top. + for rank in ("median_rank", "prop_rank"): + self.assertCountEqual(features[rank], range(lr.n_vars)) + np.testing.assert_array_equal( + features["mean_rank"], features[["median_rank", "prop_rank"]].mean(axis=1) + ) + self.assertEqual(features["mean_rank"].idxmax(), "LigA_RecA") + + def test_expression_quantiles_are_in_varm(self): + """The quantiles of each pair's ligand and receptor expression, which the + test matched its background genes on, are a DataFrame in varm: a row per + pair in the table's order, 'L_q' columns then 'R_q'.""" + lr = self.lr + quantiles = lr.varm[QUANTILES] + self.assertIsInstance(quantiles, pd.DataFrame) + self.assertTrue(quantiles.index.equals(lr.var_names)) + expected_columns = [f"L_{float(q)}" for q in RUN_QUANTILES] + expected_columns += [f"R_{float(q)}" for q in RUN_QUANTILES] + self.assertEqual(list(quantiles.columns), expected_columns) + self.assertFalse(quantiles.isna().any().any()) + for name in lr.var_names: + with self.subTest(lr=name): + expr = self._genes(name) + expected = np.quantile(expr, RUN_QUANTILES, axis=0, method="nearest") + np.testing.assert_array_equal( + quantiles.loc[name].to_numpy(), expected.T.ravel() + ) + # Not all zero: the fixture's co-expressed corner expresses LigA at 6. + self.assertEqual(quantiles.loc["LigA_RecA", "L_1.0"], 6) def test_layers_are_sparse_with_the_tested_spots(self): """Each layer is a CSC matrix with an entry for each spot scored for an @@ -143,6 +196,22 @@ def test_adj_pvals_reorders_the_table(self): lr_result(lr, "lr_scores", name), lr_result(self.lr, "lr_scores", name), ) + self.assertEqual(list(lr.var.columns), SUMMARY) + pd.testing.assert_frame_equal( + lr.var[FEATURES], self.lr.var.loc[lr.var_names, FEATURES] + ) + pd.testing.assert_frame_equal( + lr.varm[QUANTILES], self.lr.varm[QUANTILES].loc[lr.var_names] + ) + + def test_adj_pvals_drops_run_cci_columns_but_keeps_features(self): + lr = self.with_cci.copy() + self.assertTrue(set(lr.var.columns) - set(SUMMARY), "run_cci added columns") + st.tl.cci.adj_pvals(lr) + self.assertEqual(list(lr.var.columns), SUMMARY) + pd.testing.assert_frame_equal( + lr.var[FEATURES], self.lr.var.loc[lr.var_names, FEATURES] + ) def test_plots(self): """Every plot of LR results draws from the table given as lr_results=."""