Conversation
…m 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.
for more information, see https://pre-commit.ci
…t layer of adata and ii) no longer creating temporary adata object in _UtagEmbedder
for more information, see https://pre-commit.ci
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #1284 +/- ##
==========================================
+ Coverage 78.44% 78.63% +0.19%
==========================================
Files 63 64 +1
Lines 9532 9642 +110
Branches 1594 1626 +32
==========================================
+ Hits 7477 7582 +105
- Misses 1489 1498 +9
+ Partials 566 562 -4
🚀 New features to boost your workflow:
|
for more information, see https://pre-commit.ci
…cells are not re-visited
…her cleaning of _NhoodProfileEmbedder
…ation of hop adj matrices (_compute_hop_adjacency_matrices function) and removing unnecessary adata creation
for more information, see https://pre-commit.ci
…ces from _CellcharterEmbedder
…n _NhoodProfileEmbedder
… the result columns back to categoricals in _postprocess_niche_results, iii) explained resolutions arg better in calculate_niche_neighborhood
for more information, see https://pre-commit.ci
… Also renamed the CellcharterEmbedder to a more generic name
…te_niche_cellcharter
for more information, see https://pre-commit.ci
shashkat
marked this pull request as ready for review
September 10, 2026 04:36
Contributor
Author
|
Hi, @grst and @selmanozleyen, this is ready for review. I think the suggestions were all quite useful and some of them I also thought about during the refactor, but didn't implement at that time. Here are some questions/points-of-discussion from me about the same:
|
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.
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.
This was referenced Sep 15, 2026
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.
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.
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.
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.
`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`.
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.
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.
`calculate_niche_custom` is the name scverse#1245's description and scverse#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.
`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.
`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.
Member
|
Hi @shashkat , I'd like your feedback on my changes. So here is what I changed on top of your reviews. Note: some of these items are from AI generated output.
besides the changes I made on your commits here are the changes vs main and the latest released version: #1293 (still in progress because I found some more bugs in 1.8.3) |
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.
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.
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.
Member
Member
|
I think it'd be nice to add a "share_niches: bool" parameter to the entire function that, if library_key is preset, decides whether things are jointly modeled. But that'd go beyond the scope of this PR. That's iirc also what official Cellcharter does |
## 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.
selmanozleyen
force-pushed
the
niche_refactor_suggestions
branch
from
September 18, 2026 11:32
8e59c35 to
9f01a12
Compare
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.
Contributor
|
IMO the neighborhood graph should respect niches, but the embedding+clustering should be across samples. That's the whole point, to classify cells into shared labels that can be used to compare between groups of samples. |
Contributor
Author
|
I think with the Just to reiterate, even though |
Member
Where does cellcharter provide something like this?
Is this in general? If so I don't see a reason to require library_key then (for other flavors except spatialleiden)? because all the embedding + clustering is already across samples by default |
Member
It doesn't provide the option, it just jointly models it by default.
Yeah, I agree. I also wouldn't require it, the baseline assumption should be the thing that makes the most sense (joint niches here) |
Member
|
Ok any reason not to remove it then for utag, cellcharter, neighbours? Because slice and run is all we do there and we can just document that if someone really wants it. Because |
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
4 participants
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.
Fixes: #1277
Description
Big thanks to @grst for invaluable feedback (#1277) on the niche refactor in #1245. Incorporating those suggestions in this PR.