[flying-sheep]

Note

since people are still providing input, this PR is pretty vibe-coded at the moment. I’ll go over it in detail once everything is hashed out.

Current design:

• deprecate inplace everywhere, copy almost everywhere: pp.sample, pp.subsample, pp.filter_cells and pp.filter_genes subset the object (they change its shape).
• functions get use (where to read from) and out (where to write to). out=False (default) writes to the default place (for matrix transforms: back to use), out=True returns the result instead of writing it, and an accessor writes there.
• a str in use/out is always A.resolved (unless the v1 preset is set and the parameter took strings with another meaning in scanpy 1.12.x). No sc.Prefix, and no guessing the meaning from the string’s content.
• functions with several outputs: a single accessor names the primary output k, and the others follow it: k in kinds of slots no other output goes to, k_{name} where several share a kind. Outputs that aren’t written by default stay unwritten. A mapping (out={"connectivities": A.obsp["c"]}) writes only what it names, out=True returns a NamedTuple (PcaOut, NeighborsOut, …). Naming the default places is identical to out=False. Functions without a primary output (scrublet, score_genes_cell_cycle, dpt) only take a mapping.
• use takes what was passed to the producer’s out. Consumers of neighbors take the graph (use=A.obsp["nn"]), replacing neighbors_key/obsp/adjacency; distances and params are found where pp.neighbors stored them next to it. umap only needs data to lay out more than 2 * n_components connected components, and then asks for it explicitly (use=A.obsp["connectivities"], rep=A.obsm["pca"]) instead of silently guessing.
• groupings take references: batch, covariates, and groupby/groups accept A.obs["x"].
• functions writing only to .uns:
  • rank_genes_groups, filter_rank_genes_groups, dendrogram, paga) get out: bool and keep key_added; their readers keep key. These switch to A.uns[k] once anndata has .uns accessors.
  • Likewise experimental.pp.normalize_pearson_residuals_pca and recipe_pearson_residuals (out=True replaces inplace=False; PCA destination via kwargs_pca)
  • normalize_pearson_residuals_pca and recipe_pearson_residuals return an AnnData object and otherwise write to .uns
• any accessor is only stringified for display purposes, not to create data from the string

Examples

# in-place on X
sc.pp.log1p(adata)
# read counts, write a layer
sc.pp.log1p(adata, use=A.layers["counts"], out=A.layers["log1p"])
# return the matrix instead of writing it
logged = sc.pp.log1p(adata, out=True)

# several outputs, chained: writes obsm["p"], varm["p"], uns["p"]
sc.pp.pca(adata, use="layers.log1p", out="obsm.p")
# writes obsp["nn"], obsp["nn_distances"], uns["nn"]
sc.pp.neighbors(adata, use="obsm.p", out="obsp.nn")
sc.tl.umap(adata, use="obsp.nn")
sc.tl.leiden(adata, use="obsp.nn", out="obs.clusters")

Example for how this improves code:

-logged = sc.pp.log1p(adata_obs, copy=True)
-sc.pp.highly_variable_genes(logged)
-adata_obs = adata_obs[:, logged.var["highly_variable"]]
+sc.pp.log1p(adata_obs, out="layers.log1p")
+sc.pp.highly_variable_genes(adata_obs, use="layers.log1p")
+adata_obs = adata_obs[:, adata_obs.var["highly_variable"]]

Blocked on anndata (follow-ups)

• .uns accessors (A.uns accessors anndata#2657): out=A.uns[k] for .uns-only functions
• column refs into obsm/varm DataFrames: writing several outputs as one frame (out=A.obsm["cc"]), and chaining their columns
• storable accessors (Allow storing heterogeneous lists in .uns (as JSON or so) anndata#1979): store use in .uns params instead of use_rep/layer strings
• min anndata with anndata.acc: accessor-only internals, and decorators replacing the inline legacy-parameter handling

Follow up after this:

• Default destinations

Follow-ups once anndata ships …

… .uns accessors (scverse/anndata#2657)

• out=A.uns[k] for rank_genes_groups, filter_rank_genes_groups, dendrogram, paga; their readers take use=A.uns[k]. Deprecate key_added/key there. (See “Accessors stored in .uns” below for pp.scrublet, tl.score_genes_cell_cycle and tl.dpt)
• ingest(use=A.uns[k]) for the neighbors model; deprecate neighbors_key.
• normalize_pearson_residuals_pca and recipe_pearson_residuals have out: bool and return AnnData. Choosing A.uns[...] should write an AnnData there?
• Make params a real output field in the …Out tuples, and decide whether out=True returns it.
• Let groupby/groups of rank_genes_groups, filter_rank_genes_groups, dendrogram, embedding_density and paga take any obs-aligned vector (e.g. A.obsm["x"][:, 0]), not only .obs columns: store the accessor in params so readers (pl.rank_genes_groups_*, pl.dendrogram, pl.paga) can fetch it. Delete _obs_col.

… references into DataFrames in obsm/varm (and scverse/anndata#2722)

• Let score_genes_cell_cycle, scrublet, dpt and hvg (both) write a single frame to obsm[k]/varm[k]; this also brings back hvg's out.
• groupby, covariates, mask and color accept frame columns.
• Chaining tests: hvg → pca(mask=…), cell cycle → regress_out, qc → color=A.obsm["qc"]["total_counts"].

… Accessors stored in .uns (scverse/anndata#2657)

• Store use accessors in params; store accessors in NeighborsDict. Delete _use_to_uns/_use_from_uns/_ref_to_uns, keeping a reader for the old format.
• pp.scrublet, tl.score_genes_cell_cycle and tl.dpt will allow out=A.uns[k] to store accessors in there (see “.uns accessors” above)

Follow-up once minimum anndata ≥ 0.13.3

• Make the internals accessor-only: delete _dest, _SetRep, _Rep and the **choices branches.
• Add @legacy_use/@legacy_out decorators to replace the inline handling of old parameters; this is the PR's "decorators" TODO.
• Delete the fallback classes in get/_types.py and most of the own_deprecations() call sites.

Separate PRs

• Defaults for the v2 preset, e.g. log1p writing to layers["log1p"] as grst suggested.
• split hvg by flavor** (like Split up rank_genes_groups and highly_variable_genes #4305 for rank_genes_groups; not blocked on anndata)
  • Create sc.pp.hvg.{seurat,cell_ranger,seurat_v3,pearson_residuals} with flavor-specific parameters. Keep highly_variable_genes(flavor=…) as a dispatcher.
  • Give each its own …Out tuple, with a shared helper for the batch fields.
  • Give each out with an accessor (out=A.var["hv"] writes hv_* columns and uns["hv"]), a mapping, and out=True returning the tuple.
  • Give pl.highly_variable_genes its use back.

Notes

• functions that haven’t been released in a stable version get no backwards compat code

---

• Closes anndata.acc replacement for functions taking obsm, layer, … #4334, closes Remove inplace/copy #2583, closes don’t return dicts #17
• Tests included or not required because:

• Release notes not necessary because:

[codecov]

Codecov Report

❌ Patch coverage is 95.31388% with 53 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.05%. Comparing base (c1f5350) to head (191fa32).
✅ All tests successful. No failed tests found.

Files with missing lines	Patch %	Lines
src/scanpy/get/_use.py	90.56%	20 Missing ⚠️
src/scanpy/get/_out.py	90.79%	15 Missing ⚠️
src/scanpy/neighbors/_common.py	93.61%	3 Missing ⚠️
src/scanpy/plotting/_v2/_tl.py	72.72%	3 Missing ⚠️
src/scanpy/preprocessing/_simple.py	90.00%	3 Missing ⚠️
src/scanpy/experimental/pp/_normalization.py	87.50%	2 Missing ⚠️
src/scanpy/tools/_rank_genes_groups.py	83.33%	2 Missing ⚠️
src/scanpy/_keys.py	0.00%	1 Missing ⚠️
src/scanpy/metrics/_common.py	94.73%	1 Missing ⚠️
src/scanpy/plotting/_v2/_core.py	0.00%	1 Missing ⚠️
... and 2 more

Additional details and impacted files

@@            Coverage Diff            @@
##           main    #4359       +/-   ##
=========================================
+ Coverage      0   85.05%   +85.05%     
=========================================
  Files         0      138      +138     
  Lines         0    12263    +12263     
=========================================
+ Hits          0    10430    +10430     
- Misses        0     1833     +1833     

Flag	Coverage Δ
hatch-test.low-vers	78.57% <74.97%> (?)
hatch-test.pre	84.84% <94.87%> (?)

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
src/scanpy/_docs.py	100.00% <100.00%> (ø)
src/scanpy/_settings/presets.py	87.24% <100.00%> (ø)
src/scanpy/_utils/__init__.py	77.19% <100.00%> (ø)
src/scanpy/experimental/_docs.py	100.00% <100.00%> (ø)
...c/scanpy/experimental/pp/_highly_variable_genes.py	100.00% <100.00%> (ø)
src/scanpy/experimental/pp/_recipes.py	100.00% <100.00%> (ø)
src/scanpy/external/pp/__init__.py	100.00% <ø> (ø)
src/scanpy/get/__init__.py	100.00% <100.00%> (ø)
src/scanpy/get/_aggregated.py	96.37% <100.00%> (ø)
src/scanpy/get/_types.py	100.00% <100.00%> (ø)
... and 53 more

... and 75 files with indirect coverage changes

[grst]

functions operating on a main matrix get use and out – default for out is to operate in place, out=None means returning the matrix instead of writing it into AnnData.

If we want to discourage use of inplace we could also make the out param explicit:

• not specified -> fails
• out=False -> inplace
• out=True -> return

[flying-sheep]

Yeah, I think bool instead of None|Default is better. I thought that users probably don’t need to pass a dynamic value (since it completely changes how the function works), but out=False for inplace and out=True for return is just more intuitive.

[grst]

Or sentinels, e.g. A.INPLACE | A.RETURN? But I guess the boolean is nicer to type.

[flying-sheep]

Yeah, an enum would of course be clearest, but I think the bool is intuitive enough, don’t you think?

I currently left the default to be False, i.e. “overwrite what use points to”, but it doesn’t need to stay that, we could:

• force people to specify out=... as you suggest, or
• use function-dependent conventions like e.g. the defaults for log1p could be use="layers.counts" and out="layers.logged" or something.

[grst]

use function-dependent conventions like e.g. the defaults for log1p could be use="layers.counts" and out="layers.logged" or something.

I like that, this would establish some conventions

[flying-sheep]

Yeah! I’ll meet with Ilan today about this, and we’ll probably hash something out. I intentionally left it out of this PR so far since the design is still changing, so I don’t want to muddle the discussion with bikeshedding about which conventional names to settle on.

[flying-sheep]

Plan:

• remove key_added
• for multi-out, add global lookup and validation table for each function’s slots, sth. like sc.OUTS["pp.neighbors"] == dict(distances=A.obsm["distances"], connectivities=A.obsm["connectivities"])
• either:
  • add prefix parameter to all use/out functions, as mutually exclusive replacement to out
  • allow to use a class like sc.Prefix("foo") for out, e.g. sc.pp.neighbors(..., use=sc.Prefix("foo"))

[flying-sheep]

OK, multi-out is blocked on not having .uns accessors then: scverse/anndata#2657

Storing accessors for e.g. neighbors in anndata is scverse/anndata#1979 (and the follow-up of storing accessors)

[ilan-gold]

Got through _simple.py - my main two things are use_rep's deprecation (could punt on this) and whether or not the plan is to store accessors in uns right now. The current state of things has that, it seems, but I thought that that was not a goal.

[flying-sheep]

@selmanozleyen any insights from your return value efforts in squidpy?

[selmanozleyen]

I will share share some things that are no go to me, some edge cases and in the end my opinion on this.

No go's

I am really against this kind of behaviour:

• out="obsp.connectivities" writes to obsp["connectivities"], obsp["distances"], uns["neighbors"]
• but out="obsp.connectivities2" writes to obsp["connectivities2"], obsp["connectivities2_distances"], uns["connectivities2"]

For example what would this mean then? ( also I don't think this is an edge case, it might happen to someone)

sc.pp.neighbors(ad)                                 # obsp connectivities/distances, uns["neighbors"]
sc.pp.neighbors(ad, n_neighbors=5, out="obsp.distances")

Before I share my opinions here are some edge cases we might need to think about:

other edge cases

1- same ref's to different keys. Do we allow?

sc.pp.neighbors(ad, out={"connectivities": A.obsp["x"], "distances": A.obsp["x"]}) 

2- a stupid case

sc.pp.pca(ad); sc.pp.neighbors(ad, out="obsp.pca") 
# primary outputs not overwritten but the uns are
# this particular example is unlikely but it implies a rule we'd need to follow
# whenever there is sc.pp.neighbors(ad, out="obsp.x") we'd need to ensure there is no function x since
# it would write to uns["x"] as sc.pp.neighbours would do. I.e., by default there is an implicit write here that might cause complications

my opinion

I would treat out more as an advanced tooling, I'd let the defaults be the convinience for the "users". We can have default presets and things like that for them to customize etc.

But out should be more clear and explicit imo. If one key is given other's shouldn't be written explicitly with it. I don't think it's worse UX if it prevents bugs that the old UX tricked the users into not noticing those bugs.

another point

why don't we use typed dicts for out mappings? It would give nice autocomplete it's extra code but it's a plus on the UX and typing I'd say.

[flying-sheep]

I am really against this kind of behaviour:

Agreed, there should be no magic strings. neighbors really is the only instance where things are muddled. I ignored it since the distance matrix can be input or output, but there’s still no reason to do magic strings.

same ref's to different keys. Do we allow?

no, makes no sense. good catch for user-friendliness.

we'd need to ensure there is no function x since it would write to uns["x"] as sc.pp.neighbours would do. I.e., by default there is an implicit write here that might cause complications

that was the motivation behind the initial design with sc.Prefix – always adding it adds some kind of namespaces and prevents clashes, but people didn’t like that it needs a new class.

I really don’t know if there’s an ideal fix, both options here have their trade-offs.

If someone has a better idea today I’ll do it, but if nobody comes forth, I can’t wait forever.

why don't we use typed dicts for out mappings?

no reason. could just live next to the FuncNameOut NamedTuples.

[flying-sheep]

all done.

• renamed FuncOut NamedTuples to FuncResult and added FuncOut TypedDicts to type out
• no more magic strings, out=False is defaults, all accessors are handled the same
• duplicates in out now raise an error

[grst]

Are there any constraints on out?

For instance, let's take pp.log1p(adata):

• typically, one would write to out=A.layers["log1p"]
• could I write to A.obsm["log1p"]? It would fit, but layers would be preferable, since it is aligned to both axes
• for the sake of the argument, let's assume adata.shape[0] == adata.shape[1]. Could I write to A.varm["log1p"]?

[flying-sheep]

It’s currently a bit overconstrained (mainly around obsp and dim), my plan was to relax these as people come up with use cases.

could I write to A.obsm["log1p"]?

yup, that works now as well.

let's assume adata.shape[0] == adata.shape[1]. Could I write to A.varm["log1p"]?

Pretty sure that not, I’m checking that the dimensions match in many cases (but might have missed a spot).

What will be enabled in the future though is doing e.g. sc.pp.neighbors(adata, use=A.varm["foo"]) or sc.pp.neighbors(adata, out="varp.xyz") to measure distances and connectivities between genes. Basically when a dimension can be unambiguously inferred from the input, it should work if it can work.

[Intron7]

I can agree with @grst that you want to discourage people from using inplace operations. However copying count matrices around can also be very expensive. If there was a good way to encode linear operators in layers we could get rid of this until than we need out=False