Repository navigation
Add function for calculating niches - #831
Merged
Merged
Conversation
LLehner
marked this pull request as draft
May 27, 2024 16:24
…into niche_definitions
for more information, see https://pre-commit.ci
…into niche_definitions
for more information, see https://pre-commit.ci
for more information, see https://pre-commit.ci
This was referenced Jun 12, 2024
for more information, see https://pre-commit.ci
for more information, see https://pre-commit.ci
LLehner
commented
Mar 13, 2025
| inplace, | ||
| ) | ||
|
|
||
| if resolutions is None: |
Member
Author
There was a problem hiding this comment.
Why not set this as default right away?
| if resolutions is None: | ||
| resolutions = [0.5] | ||
|
|
||
| if distance is None: |
| elif libraries is not None and len(libraries) > 0: | ||
| return [f"{base_column}_{lib}" for lib in libraries] | ||
|
|
||
| # For neighborhood and utag, we need to handle resolutions |
Member
Author
There was a problem hiding this comment.
Would do this with if-else block maybe
| groups: str | None, | ||
| n_neighbors: int | None, | ||
| resolutions: float | list[float], | ||
| subset_groups: list[str] | None, |
Member
Author
There was a problem hiding this comment.
why remove this completely? It serves a different purpose than masking. This was also described in the original tutorial. E.g. a mask allows to restrict niche calculation to specific cell types but they can have neighbors outside the masked region. Subsetting means only observations within that subset are considered as neighbors.
|
|
||
| def _get_utag_niches( | ||
| adata: AnnData, | ||
| subset_groups: list[str] | None, |
Member
Author
There was a problem hiding this comment.
Here subsetting serves the same purpose as before.
| if library_key not in data.tables[table_key].obs.columns: | ||
| raise ValueError(f"'library_key' must be a column in 'adata.obs', got {library_key}") | ||
|
|
||
| if n_neighbors is not None and not isinstance(n_neighbors, int): |
Member
Author
There was a problem hiding this comment.
Why check these separately?
| for param_name in param_specs["unused"]: | ||
| param_value = param_dict.get(param_name) | ||
|
|
||
| # Special handling for boolean parameters with default values |
selmanozleyen
added a commit
to shashkat/squidpy
that referenced
this pull request
Sep 15, 2026
UTAG, CellCharter and SpatialLeiden get `:cite:` entries; the bibliographic details come from Crossref via DOI content negotiation rather than from recall. `calculate_niche_neighborhood` gets no citation because it has none: PR scverse#831 introduced it as "based on neighborhood profiles similar to" monkeybread's `_neighborhood_profile.py`, and monkeybread cites nobody for the method, so the source link is the provenance. Two comments restore provenance the rewrite dropped, phrased as the relationship rather than as "adapted from", since neither is a copy any more: - the hop rings say what CellCharter does (`adj_hop @ adj`, then `adj_hop > adj_visited`, in scipy) and why this does not: a numba BFS gives the same disjoint rings without a comparison that is degree-dependent on a weighted graph. - the neighborhood profile says monkeybread counts per cell with a `Counter` where this is a one-hot product, and that the scaling follows theirs. CellCharter's `tl/_gmm.py` is deliberately not restored. It is torchgmm with `init_strategy="kmeans"`; this uses scikit-learn's `GaussianMixture` with `init_params="random_from_data"`, so the old "adapted from" was no longer true. The docstring says so instead, along with `use_rep` being the recommended input and PCA the fallback. Checked: pybtex parses the file, a minimal sphinx build with -W resolves all three keys and renders the accented names, and all four docstrings parse as RST after docrep substitution.
selmanozleyen
added a commit
that referenced
this pull request
Sep 25, 2026
* Moved embedding step out of library_key loop in calculate_niche_custom and store the niche embedding in adata.obsm now. Not sure how much computational benefit this offers, as now, instead of many small embeddings (for each library-subsetted adata), one big embedding is being computed (for the whole adata). But irrespective of computational advantage, I think it can be useful to store the embedding in adata.obsm to inspect for end user. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * incorporated suggestions in utag flavor: i) now can supply a different layer of adata and ii) no longer creating temporary adata object in _UtagEmbedder * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fixed docs * improved implementation of calculate_nhood_profile in NhoodEmbedder * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove import * fixed _hop function to account for node degree when ensuring visited cells are not re-visited * removed unnecessary .tolil() and .tocsr() in _setdiag function * dont convert output to df in _calculate_neighborhood_profile() and other cleaning of _NhoodProfileEmbedder * improved implementation of _CellcharterEmbedder by formalizing computation of hop adj matrices (_compute_hop_adjacency_matrices function) and removing unnecessary adata creation * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove import * added tests for _compute_hop_adjacency_matrices * avoid duplication of matrix when calling _compute_hop_adjacency_matrices from _CellcharterEmbedder * use standardized hopping function _compute_hop_adjacency_matrices() in _NhoodProfileEmbedder * keep the warning style to consistently use logg instead of warnings.warn * throw error when n_hop_weights < distance in NhoodProfileEmbedder * i) made only data as positional arg in user-facing funcs, ii) convert the result columns back to categoricals in _postprocess_niche_results, iii) explained resolutions arg better in calculate_niche_neighborhood * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * updated documentation of user-facing functions to be more descriptive * move use_rep from CellcharterEmbedder to calculate_niche_cellcharter. Also renamed the CellcharterEmbedder to a more generic name * remove n_components arg from _NHopPCAEmbedder * made it not necessary to have sparse expression matrix to use calculate_niche_cellcharter * replaced sc.tl.pca with sc.pp.pca as recommended API * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix hop semantics, build niche features from one primitive Fixes - Shell exclusion compared path counts, not sets, so it depended on node degree. CellCharter's `_aggregate_neighbors` starts with `adj = adj.astype(bool)`, which the port dropped. - The search walks `indices` and never reads `data`, so an explicitly stored zero was traversed as a real edge while hop 1 correctly treated it as a non-edge. - Short `n_hop_weights` raised `UnboundLocalError`; `distance=0` silently returned an embedding with no spatial content. Both now fail under the names the caller passed. - The `calculate_niche` deprecation is a `FutureWarning` rather than an unfilterable log line, and its stacklevel points at user code. - Niche labels are categorical for every flavor; `_LeidenClusterer` assigned `list(...)`, which pandas flattens to object, where cellcharter's `pd.Categorical` kept the dtype. - `calculate_niche_cellcharter(aggregation="sum")` was accepted once `_aggregate_over` became shared with the neighborhood flavor. Restricted back to mean/variance. Revert - `neighborhood` goes back to matrix powers from disjoint shells. Rings suit concatenation, where blocks must not restate each other; powers suit the weighted sum, where multiplicity is the distance weighting. Bit-identical to main again. Restructure - A numba breadth-first search reaches each cell once at its shortest distance, so the rings are disjoint by construction and there is no visited set to subtract. One shared `_expand` advances a level over a FIFO, and `GridBuilder`'s ring expansion uses the same primitive: measured 8.7-20.9x faster at n_rings >= 3, bit-identical across 48 configurations. - `nhood_aggregate` is one primitive behind all three flavors. `_nhood_blocks` is gone: what was shared was verification, so that is `_assert_hop_request`, and each caller builds its own hops. - Embedders are functions, not classes, and the ABC is gone. - Sparse in, sparse out: the concatenation is no longer densified before PCA. - `Array | CSBase` from fast_array_utils instead of `NDArrayA`. Dtypes, no longer hardcoded - Index scratch follows the graph's own index dtype, which fixes a 16->32 thread regression and takes `_shell_adjacencies` from 439 to 358 ms on 200k cells. - `n_jobs` caps the search and its per-thread scratch, resolved with `get_n_numba_threads`. - `niche_embedding` is float32 for every flavor, where `neighborhood` was float64. This changes Leiden partitions on near-degenerate graphs (12 of 24 probed configurations) and needs a release note. `test_niche.py` goes 23 -> 21 functions covering more cases, and imports no private squidpy name. * Give clusterers a protocol The niche flavors took a _NicheClusterer subclass whose cluster() both fitted and wrote to adata.obs, and which looped over resolutions internally. They now take a mapping of obs column name to a structural Clusterer, so one clusterer labels one column and a resolution sweep is a mapping rather than a loop. - squidpy/types.py is new, with the runtime-checkable Clusterer protocol. Structural on purpose: sklearn's GaussianMixture is a DensityMixin rather than a ClusterMixin, and cuML does not import sklearn at all, so a base-class check would reject both. - gr/_clusterers.py holds LeidenClusterer as a sklearn estimator, so clone/get_params/set_params come for free and each fit goes to a fresh clone. Public name in a private module, absent from __all__ and api.md, as nhood_aggregate and compute_hop_adjacency_matrices already are. - _GMMClusterer is gone. GaussianMixture already satisfies the protocol, so cellcharter passes it directly, which is the point of making the protocol structural. - _NicheClusterer and the ABC are gone, matching what this branch already did to the embedders. This mirrors scverse/feat/cluster-auto-k-seed's design without its dependencies. It deviates in one place: there LeidenClusterer.fit takes a prebuilt graph and an embedder builds it, so a resolution sweep builds one graph; here fit takes the embedding and builds its own, so N resolutions build N graphs. That is the embedder split, and it is the follow-up. Behaviour: 79 of 82 captured records identical across 10 configurations, including every label for neighborhood, utag and spatialleiden, every embedding and every dtype. The three that differ are cellcharter's labels, because _GMMClusterer drew its seed directly from the root generator while the pipeline now spawns one generator per clusterer. Keeping the spawn rule uniform matches Leiden exactly at any number of resolutions and moves the GMM seed; the reverse choice would have moved every multi-resolution Leiden seed instead. * Drop the BFS kernel docstring Removed, along with a stale claim: scratch is now two buffers per thread, not the three the text described. * Trim the hop rings docstring 28 lines of numpydoc sections for a 27-statement private helper, down to 3. * Trim the clusterer docstring 18 lines with a Parameters block for five self-describing arguments, down to 3. * Keep library niche labels categorical The per-library loop seeds a column with the string "not_a_niche" and writes each library's labels in as strings, so the aggregate stayed object where the non-library path gave a category. Cast once every library has been seen; casting inside the loop would reject the labels the next library brings. The aggregation tail was duplicated verbatim between calculate_niche_spatialleiden and _calculate_niche_custom, so it is now _merge_library_columns and the fix lands once. Verified: cc_lib and nbhd_lib label dtypes go str -> category, utag and spatialleiden likewise, and no label or embedding value changes. * Act on the review findings Independent reviewers went over f564b49 in two rounds. This is everything that was mine to fix. Correctness - The per-library loop assigned `added_columns` only under `if itr == 0:`, so a `library_key` whose first value is NaN skipped library 0, left the list empty, and silently discarded every later library's labels. Both loops now take the names from what was written. - Re-running a `library_key` call in place kept the previous run's labels: the first run leaves a categorical column and `.loc` cannot add unseen categories to one. - With `library_key` the embedding was fitted once on the pooled object and sliced, where upstream/main fitted it per library. Pooling leaves each section's own offset inside the rows that are then clustered, which is what stratifying exists to remove. Now per library; no embedding is written for a stratified run, since the blocks are in different spaces and need not share a width or a column basis. - `use_rep` stopped truncating to the first `n_components` columns and stopped rejecting a narrower embedding. Both behaviours are in v1.8.3 and in main, so that was a released-contract break; restored, message verbatim. - The BFS sized its scratch from `indptr` but indexed it with `indices`, so a rectangular input wrote out of bounds and killed the process with SIGBUS. Now a ValueError, via `np.shape` so an array-like still works. - `nhood_aggregate` accepted any `aggregation` when every hop was 0, and returned an all-NaN matrix for weights summing to zero. Both checked up front, along with the weight count, which previously ran after the whole aggregation. - `embedding_key_added=None` wrote `obsm[None]`, which only failed later at `write_h5ad`. - Ring 0 kept whatever `.tocsr()` returned while rings 1+ were `csr_matrix`, so a `csr_array` input came back as a mixed list where `*` means different things per element. Performance - The counting pass handed `_bfs_shells` int64 dummies while the filling pass handed it an `indices.dtype` `out`, so the `parallel=True` kernel compiled twice. Both dummies now match: one specialization, cold start 1077 -> 660 ms, reproduced independently at 1076 -> 653. Interfaces - `Clusterer` validates what it is handed: a deterministic estimator satisfies the protocol and then rejects the seed the pipeline sets, so `runtime_checkable` now has a caller. - Dropped `__all__` from `_clusterers.py`, no other private module declares one. - The weights warning pointed inside squidpy; `n_jobs` was undocumented; `generators` is `rngs`, matching `library_rngs` in the same file. Tests - A mutation run showed the suite caught none of the guards. Twenty added, one per guard, plus the per-library and `use_rep` contracts. - `test_niche_copy_semantics` compared two unseeded runs for equality: 1 disagreement in 250 runs with the real fixture, 0 with `rng=0`. Needs a release note: per-library embedding, no `.obsm` when stratified, and the `use_rep` truncation are all user-visible changes. * Close four test gaps 3.2 The BFS stripes sources across threads with range(thread, n, n_threads), and every existing case ran on fewer nodes than NUMBA_NUM_THREADS, where that is a no-op. Mutating the stride to n_threads + 1 previously left 83 tests green while producing 385 wrong entries at n=100. Two tests now cover it: rings compared against scipy.sparse.csgraph.dijkstra at n = 50/137/200, directed and undirected, and thread-count invariance across n_jobs 1/2/3/8. The same mutation now fails 9 of them. 3.3 The cellcharter path had no numeric assertion: flipping the sign in the variance formula and dropping the hop-0 raw-feature block both passed. aggregation='variance' is now checked against E[x^2] - E[x]^2 computed per neighborhood, and the concatenation width handed to PCA is asserted as (distance + 1) * n_vars. Those mutations now fail 1 and 2 tests. 3.8 n_jobs uses the shared %(n_jobs_threads)s docrep template instead of hand-written text, so it documents what get_n_numba_threads enforces: None/-1 mean NUMBA_NUM_THREADS, over-asking warns and clamps, 0 and below -1 raise. 3.10 Restores the coverage the deleted seeding test took with it: cellcharter on a sparse X, and cellcharter through the library_key path, which had none at all. * Validate resolutions once 3.6 `_validate_niche_args` has one caller: the deprecated `calculate_niche`. Everything it checks is therefore unreachable from the four functions that replace it in v1.9.0, so the new API let bad input through to scanpy: resolutions='high' calculate_niche(flavor='utag') TypeError: 'resolutions' must be a float, a tuple of floats, ... calculate_niche_utag TypeError: must be real number, not str 3.7 The only container test was `isinstance(resolutions, list)`, so anything else counted as a single resolution: tuple (0.5, 1.0) -> TypeError: must be real number, not tuple, from inside scanpy ndarray [0.5, 1.0] -> TypeError: only 0-dimensional arrays can be converted list [0.5, 0.5] -> one column, not two, because the names collide in the dict A tuple is legitimate, but only for spatialleiden, where it is the (latent, spatial) pair; the shared validator accepted it for every flavor. Duplicates were the worst of the three: fewer columns than asked for, no warning. `_resolution_values` now normalises and checks in one place, called from `_leiden_clusterers` (neighborhood, utag) and once at the top of `calculate_niche_spatialleiden` (pairs_ok=True). Both are before any embedding work, so a typo fails immediately rather than after the PCA. ndarrays and ints now work; the dead block in `_validate_niche_args` is gone, which also drops its rejection of an int resolution that the flavors accept. Six tests. Restoring the `isinstance(..., list)` line fails all six. * Document the utag PCA The documented steps said the aggregated matrix is the niche embedding, but `_utag_embedding` returns `sc.pp.pca(aggregated)`, so `.obsm[embedding_key_added]` holds PCA scores. Main is equally undocumented here; it only carried an inline comment, which this branch's rewrite dropped. The docstring is the right place for it. * Drop four explanatory comments * Consolidate the niche test fixtures The four `nhood_aggregate` rejections and the four bad `resolutions` become two tables; the positive case each one carried is now its own named test rather than a tail assertion. The two `library_key` embedding tests built the same two-section object by hand, so that is `_two_sections` now, and `_with_embedding` collapses into a keyword on `_tiny`. Same coverage: restoring `isinstance(resolutions, list)` still fails 7, and the BFS thread stride mutation still fails 9. * Stop recursing for library niches `calculate_niche_spatialleiden` was the only flavor with a public `prefix`, and passing it together with `library_key` silently did nothing: the stratified path hardcodes `prefix=f"lib={lib_id}_"` in the recursive call, so a caller's value is the thing being overwritten. The parameter was never a feature. Spatialleiden stratifies by calling itself with `library_key=None`, so the argument the library loop sets had to sit on the public signature to be reachable; the other three flavors pass the same argument to the private `_run_niche_pipeline` and never exposed it. It is absent from v1.8.3, where prefixing was hardcoded. The `else:` branch was the whole single-library algorithm and the recursion existed only to reach it, so it is `_spatialleiden_once` now, with the pass-through arguments bound once and both branches calling it. `prefix` comes off the public signature. The re-entry was also redoing `extract_adata_if_sdata`, the rng normalisation, `_resolution_values` and the SpatialData `sanitize_table` check once per library; that happens once now. `copy=False` and `library_key=None` go with it, since they only existed to tame the recursion. The `try: import spatialleiden` guard moves into the helper, because the outer binding goes unused once the loop moves. Same message, now raised after the `library_key` assertion rather than before it. Labels and dtypes are identical on 18 of 18 captured records across seven configurations: single, multi-resolution, a (latent, spatial) pair, `min_niche_size`, `mask`, and each of those again with `library_key`. * Use pytest.param for hop cases The ids were re-derived from the cases by index, `[case[1:] ...]` against `ids=[case[0] ...]`, so a new case had to be sliced correctly in two places to keep its label. `pytest.param(..., id=...)` keeps the id with the case it names, and `HOP_RING_CASES` had no other reader, so it goes inline. Same eight ids, verified by collection. * Cite the niche publications UTAG, CellCharter and SpatialLeiden get `:cite:` entries; the bibliographic details come from Crossref via DOI content negotiation rather than from recall. `calculate_niche_neighborhood` gets no citation because it has none: PR #831 introduced it as "based on neighborhood profiles similar to" monkeybread's `_neighborhood_profile.py`, and monkeybread cites nobody for the method, so the source link is the provenance. Two comments restore provenance the rewrite dropped, phrased as the relationship rather than as "adapted from", since neither is a copy any more: - the hop rings say what CellCharter does (`adj_hop @ adj`, then `adj_hop > adj_visited`, in scipy) and why this does not: a numba BFS gives the same disjoint rings without a comparison that is degree-dependent on a weighted graph. - the neighborhood profile says monkeybread counts per cell with a `Counter` where this is a one-hot product, and that the scaling follows theirs. CellCharter's `tl/_gmm.py` is deliberately not restored. It is torchgmm with `init_strategy="kmeans"`; this uses scikit-learn's `GaussianMixture` with `init_params="random_from_data"`, so the old "adapted from" was no longer true. The docstring says so instead, along with `use_rep` being the recommended input and PCA the fallback. Checked: pybtex parses the file, a minimal sphinx build with -W resolves all three keys and renders the accented names, and all four docstrings parse as RST after docrep substitution. * Drop the underscore from the niche pipeline `calculate_niche_custom` is the name #1245's description and #1285 both use for the high-level entry point, and it has never existed under that spelling. A public name in a private module, absent from `__all__` and api.md, as `nhood_aggregate` and `compute_hop_adjacency_matrices` already are here. Nothing is exposed and nothing renders in the docs; exporting it later becomes one line instead of a rename in a release. * Separate n_clusters from n_components `n_components` was doing two jobs on the cellcharter flavor: it set the number of mixture components, and it also truncated a `use_rep` embedding to that many columns. The second reads it as a width. CellCharter's own clusterer names it `n_clusters` and documents that "the dimensionality of each component is automatically inferred from the data"; it has no width parameter and never truncates. The trap is sklearn's own naming, where `PCA(n_components=k)` is a width and `GaussianMixture(n_components=k)` is a cluster count. So `n_clusters` counts clusters, and `n_components` sizes the PCA, which is where a component count belongs and where nothing used it before: the PCA path ran `sc.pp.pca(aggregated)` at scanpy's default, so asking for 4 gave a 4-wide embedding through `use_rep` and a 50-wide one without it. The truncation and its `ValueError` are gone; a k-cluster GMM is well posed in any dimensionality. `use_rep` now feeds the aggregation instead of replacing it, which is the other half of the same misreading. It was bypassing the hop rings entirely, so the aggregated matrix was computed and discarded and the flavor had no spatial component at all on the path its own warning recommends. CellCharter selects features the same way, `X = adata.X if use_rep is None else adata.obsm[use_rep]`, and their `sample_key` handles per-sample neighborhoods while the representation stays shared. use_rep='X_scvi' identical after scrambling the graph: 100% -> 82% width to GMM 4 -> 24 no use_rep identical after scrambling the graph: 81% width to GMM 35 `n_components` with `use_rep` now raises rather than being silently meaningless. The deprecated `calculate_niche` keeps its released contract: its `n_components` still counts clusters, forwarded as `n_clusters`. Verified identical on v1.8.3, main and here — `n_components=7` gives 7 mixture components and 7 niches on all three. 111 tests. Four replace the two that asserted the truncation. * Name the PCA width explicitly `n_pca_components` rather than `n_components`, so neither the signature nor the internals can be read as the cluster count again. `GaussianMixture` spells its own cluster count `n_components`, which is what caused the confusion in the first place; leaving the same word on a width keeps the trap open for the next reader. The deprecated `calculate_niche` keeps `n_components` with its released meaning and forwards it as `n_clusters`. The only remaining uses of the word are that function, its validator, sklearn's own keyword, and the docstring line that explains the collision. Existing tests passed `n_components` meaning the cluster count; they now pass `n_clusters`. 111 pass. * Initialise the CellCharter GMM with k-means random_from_data lets EM settle in a local optimum that merges well-separated niches. k-means is scikit-learn's default and what CellCharter's torchgmm uses. * Reduce before aggregating in the CellCharter flavor CellCharter aggregates an already reduced representation over the hop rings and clusters the concatenation directly. Run PCA on X first (10 components by default, scVI's latent width) instead of aggregating every gene and reducing afterwards, so the PCA fallback and use_rep take the same path and no cells x genes x hops intermediate is built. The use_rep docstring still described it as replacing the aggregation; it is aggregated like the PCA it replaces. Add a regression test where two cell types are mixed in one region and in pure blocks in another: niches must follow the neighborhood, which clustering the representation alone cannot. * Harden the CellCharter PCA fallback Mask by var['highly_variable'] as sc.pp.pca(adata) does, so a user who ran highly_variable_genes first does not silently get an all-gene PCA; masking here rather than through the AnnData overload keeps the function free of side effects. Clamp the default width to what X can give (a one-marker panel is its own reduction), reject an explicit n_pca_components against the same ceiling with a message that names the parameter, and fill a preallocated block instead of np.hstack, which keeps the features' dtype and halves peak memory. * Restore mask to the niche fit, and follow-ups ## mask `mask` came from monkeybread's `cellular_niches`, where rows are dropped before the clustering graph is built. `6638075b` ported that; `b8dff5de` replaced it with a post-hoc relabel, so everything is clustered now and masked rows are renamed afterwards. Labels on the kept cells, against an unmasked run: v1.8.3 35% masked cells absent from the fit main 100% the mask has no effect on the model at all here 20% 100% is the tell. The docstring has promised the exclusion since 2024 with no flavor qualifier, and on v1.8.3 it held for neighborhood (0/20 masked cells got a niche) and was false for the other three (20/20) -- silently, because `mask` is in no flavor's `unused` list. Renamed `cluster_mask`, because the name should say which stage it restricts: the profile is built from every observation, so a masked one still reaches its neighbors and is absent only from the clustering. Applied by the three flavors that build an embedding. spatialleiden raises rather than accepting one it cannot honour: it clusters the graphs, so an observation either takes part or loses its edges, and `is_membership_fixed` is never passed to `optimise_partition_multiplex`. v1.8.3 did nothing there and did not warn, so no working call breaks. `_fitted_on` aligns the mask in one place and fixes two things: a mask that omits observations now keeps them, which is how the documented three-entry example reads and which used to raise `IndexingError` on main and `IndexError` on v1.8.3, so the documented usage has never worked; and a mask sharing no index value, or excluding everything, raises. ## Follow-ups One library loop instead of two. `_stratify` is the loop both pipelines had their own copy of, called with a `run_one` callback; 104 records byte-identical across the refactor. spatialleiden keeps its own `run_one` rather than becoming a Clusterer, because it takes two graphs and no feature matrix. `key_added` on all four flavors, defaulting to the existing column names. This retires the private `prefix`, which existed only to plumb the per-library label through the recursion `_stratify` replaced. A warning when a graph has edges across libraries under `library_key`. The loop slices `obsp`, which only preserves a block-diagonal graph, and nothing checked. spatialleiden checks both its graphs, and its `latent_connectivities_key` defaults to whatever `scanpy.pp.neighbors` left in `connectivities`, which knows nothing of libraries. `min_niche_size` was also listed as unused by every flavor while all four apply it. `use_rep` on utag, which could emit an embedding but not take one back, so its PCA sat inside the per-library loop and each library got its own basis -- widths need not even match (50 columns in one library, 7 in another). `'X'` now spells `adata.X` on both `use_rep` arguments. `_aggregate_over` scales the rows after the sum rather than normalizing the adjacency first, which rounded 1/k to float32 and left a representable mean inexact. A non-floating product is promoted first, since the CellCharter rings are bool and an integer X summed to an integer that cannot hold the quotient. Dividing by the signed row sum also fixes the weighted mean for negative edge weights. 153 tests. * Document what library_key does to niche labels Stratifying gives every library its own embedding and its own clustering, so labels carry a lib=<id>_ prefix and a niche in one library is unrelated to the same-numbered niche in another. Say so where users read it, and point at the alternative: no library_key, a batch-corrected use_rep and a graph built with spatial_neighbors(library_key=...). The Returns sections promised the embedding in .obsm even when stratifying, where nothing is written because the per-library embeddings share no axes. * Drop library_key from every flavor but spatialleiden It cannot build a graph. The flavors take spatial_connectivities_key, never spatial_key, and call no builder, so all library_key can do is subset a graph somebody else built. On a graph already built with spatial_neighbors (library_key=...) that subsetting is a no-op, measured to 0.00e+00; on one built pooled it drops the crossing edges without replacing them, leaving 198/200 cells short at least one neighbour where a genuine per-library run gives every cell its full degree. So it either does nothing or degrades, and never does the thing its name suggests. What is left is fitting a separate model per library, which makes the labels of one library unrelated to another's: three identical sections at n_clusters=4 give 12 labels with none shared, against 4 shared without it. The point of running across samples is a shared vocabulary, so that is now the only behaviour for neighborhood, utag and cellcharter. spatialleiden keeps it. It has no embedding step and both its graphs are inputs, so none of the above applies, and one multiplex run over pooled sections under-clusters through the global null model. calculate_niche is released with the parameter, so it raises for the three flavors that no longer take one rather than silently pooling, and still forwards it for spatialleiden. The loop leaves the pipeline with it. _stratify was unwrapping the SpatialData, looping over libraries and putting the table back; the first and last are what every path needs, so they are _on_table now and _stratify is only the loop on top. calculate_niche_custom runs once over everything and takes neither library_key nor graph_keys. spatialleiden calls the loop when library_key is given and the plain pipeline when it is not, so the stratified shape is opt-in rather than the frame every flavor is built inside. Two aggregation bugs surfaced on the way, both older than this change: Integer features lost their ring means. The block the rings are written into took the features' dtype, so an integer X through use_rep truncated every mean back to a count -- 0 of 240 ring entries kept a fractional part. Bool features aggregated to 1/k in every column, max error 0.75 against a float64 reference, because `bool @ bool` saturates to True instead of summing and the widening happened to the product rather than the operand. * Drop library_key, address review No flavor takes library_key now. spatialleiden builds neither of its graphs, so per-library graphs were already the caller's job, and what library_key added on top was one Leiden run per library, a separate label space for each: what the other three flavors already dropped. The deprecated calculate_niche keeps the parameter and raises for every flavor, naming the alternative. _stratify, the cross-library warning and the lib=<id>_ prefix go with it. The neighborhood profile counts each reached cell once again. Past hop 1 it read walk counts off the matrix power, so a cell reached by two paths counted twice and a cell took a flat 1/d share of its own 2-hop profile. Hop 1 keeps the edge weights; past it only reachability is defined. Both flavors warn on a weighted graph. Pins that a lower-degree cell gets its own composition, which v1.8.3's padding got wrong. From the review: - fast-array-utils>=1.5, for types.HasArrayNamespace - LeidenClusterer.labels_ is .cat.codes; the niche labels are unchanged - _bfs_shells reads get_num_threads() instead of taking n_threads - unlabelled observations come from groups.notna(), not the one-hot sums Also: - the neighbor count takes that mask instead of every hop matrix being copied with it, including under aggregation="sum" where it changed nothing; bit-identical, and the 3-hop profile is ~25% faster - the BFS no longer writes ring 0 only for it to be replaced by the input - min_niche_size applies before the categorical is built; spatialleiden reads its labels back inside its own loop, and they are strings now as on every other flavor, instead of integers unless min_niche_size was set - niche categories sort numerically, 0, 1, 2 ... 10, as sc.tl.leiden has them, instead of as text * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Deprecate library_key on calculate_niche calculate_niche is released with library_key, and raising broke every v1.8.x script that passed it. It now warns and keeps v1.8.3's behaviour until the function itself goes in v1.9.0: each library fitted on its own subset, labels prefixed lib=<id>_, observations with no library left as not_a_niche. The loop is its own function, _niches_per_library, marked for deletion with calculate_niche, which is its only caller. The flavor functions do not take library_key and do not go through it; the umbrella runs the flavor on each library's subset and merges the labels back. The columns it merges come from each flavor's key_added and _resolution_values, so a rerun overwrites them rather than keeping the previous run's labels. * Show the per-library slicing The library_key warning and docstring said to call the flavor on each library's subset without showing how, and did not say which suggestion answers which need. They now separate the two: graphs built per library either way, library_key unset for niches shared across libraries, and a slice per library for niches per library. The docstring example copies the labels back, which a loop over slices otherwise discards; run on two libraries it gives the same 300/300 labels as the deprecated path. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Selman Özleyen <32667648+selmanozleyen@users.noreply.github.com> Co-authored-by: selmanozleyen <syozleyen@gmail.com> Co-authored-by: anon <tim.treis@helmholtz-munich.de>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
Adds a function that calculates niches using different strategies. The initial function calculates niches based on neighborhood profiles similar to here.
This PR will get updated with methods discussed in #789.