Remove lrfeatures and put into lr_results. Fix for #369.

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)}")
 

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

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:

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)

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),

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)

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")

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:

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:

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=."""