Repository navigation
Zb/markers de api - #4400
Zb/markers de api#4400zboldyga wants to merge 5 commits into
Conversation
Split the 1200-line tools/_rank_genes_groups.py into the tools/markers
package (comparison, stats, kernels, scorers, results). Add
sc.tl.markers.{wilcoxon,ttest,logreg}, which return a long DataFrame
(one row per group and gene) and do not modify adata. rank_genes_groups
now runs on the same internals and keeps its adata.uns output.
- new markers preset (mean_in_log_space, wilcoxon_backend)
- wilcoxon vs reference allocates scores per group (yielded arrays were reused)
- logreg with two groups reports both groups with correct signs (new API only)
- ttest scores are computed on the data as given; mean_in_log_space only affects log_fc
Sums counts per sample with sc.get.aggregate and tests each condition level against a reference with PyDESeq2, optionally per cell group and adjusting for sample-level covariates (e.g. a paired donor design). Requires a sample key, so cells are never treated as replicates. Adds the optional extra scanpy[pydeseq2], installed in the stable test env only, and the needs.pydeseq2 test mark.
The seven sc.pl.rank_genes_groups* plots accept results=<table from a sc.tl.markers function> (plus groupby) as an alternative to reading adata.uns[key]. Reading results now goes through one internal adapter; the uns path is unchanged. Plotting from results reproduces the existing reference images. sc.get.rank_genes_groups_df delegates to a private _rank_genes_groups_df.
- new tutorial: pseudobulk DE with sc.tl.de.deseq2 on interferon-stimulated PBMCs (Kang et al. 2018), including a null comparison showing why cell-level tests should not be used across conditions - API docs: a question-to-tool table for marker genes vs. DE - clustering tutorial: marker section uses sc.tl.markers (only the affected cells were re-executed)
…rank_genes_groups_df Following sc.external (scverse#3645), they emit a FutureWarning from 1.13.0 that names the replacements (sc.tl.markers, sc.tl.de.deseq2, filtering the returned DataFrame) and move to the deprecated API page; removal is for 2.0. The sc.pl.rank_genes_groups* plots stay and accept results=. Internal callers (queries.enrich, legacy plots) use the private _rank_genes_groups_df so users only see warnings for their own calls. Tests that still exercise the deprecated functions ignore the warning. The results-based plot test now compares against the uns-based plot directly instead of sharing reference images (avoids a parallel race).
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
|
Not a maintainer of scanpy, so just my 2c. Overall, I think we definitely need to push people more towards using PyDESeq2. However, I am unsure whether this should be implemented again as a binding in scanpy - how much boilerplate code is actually saved by this? The scverse already maintains PyDESeq2 direclty and |
Hey there -- @grst @ilan-gold @maltekuehl @flying-sheep @Zethson @Intron7
I'm tagging you all because you've discussed or worked on scanpy's rank_gene_groups. I would like to clean up this code and this draft PR is a proposal of how this might look.
There are a few overlapping problems I think are worth solving:
Even though this has been discussed at length, e.g. docs: add warning for rank_genes_groups #3792 , a significant part of the single cell field keeps using rank_genes_groups (and t-tests or wilcoxon) for differential expression. The 2025 ARC VCC used a MWU-based package (pdex). And the paper that introduces the dynamic range fraction (and it's supporting code) are based on wilcoxon: https://www.nature.com/articles/s41587-026-03307-w ... Newcomers need more guidance and scanpy is exactly the place to do this, in my opinion.
There was mention of breaking down the rank_genes_groups API: Split up rank_genes_groups #4305
The existing code is messy: shared state, repeated boilerplate, patterns not matching pythonic practices / rest of scanpy code. I worked on rank_genes_groups code in other PRs earlier this year and the state of the code made reviews lengthy.
All said, I've drafted a PR here to suggest a direction that addresses each of these. This is designed as a scanpy 2.0 release, with deprecation built in where relevant. My intent is to solicit feedback -- any strong feelings about why this isn't right, or something you'd prefer to adjust? (No hard feelings, this is only a day's work, so please share if you disagree with something).
sc.tl.markers.wilcoxon
sc.tl.markers.ttest
sc.tl.markers.logreg
sc.tl.de.deseq2
(and sc.tl.rank_genes_groups and sc.tl.filter_rank_genes_groups deprecated, to be dropped in 2.0)
could also use more explicit namespaces like marker_genes and differential_expression... But the idea is to organize by use case.
A light facade supports pydeseq2 (optional dependency). DE is such a core part of work in this field, it seems logical it should be built into scanpy. This is lightweight code. And this pattern can later be expanded to include other robust DE tools.
I will add a solid tutorial specifically for differential expression, which had been alluded to in the prior discussion but not added yet. The initial draft I included here is AI generated and likely nonsense, but the final version will be correct and straightforward.
Some underlying code organization and minor refactoring to support this structure. Probably reducing some parameters. Returning outputs directly instead of inserting them to 'uns'. Once we generally agree on the direction I can refine this more to minimize the scope while cleaning up what is needed for a stable 2.0 API.
Some minor docs and tutorial adjustments to capture these changes. And usability details like how to filter (1-2 lines of pandas / python instead of the prior builtin filter_rank_genes_groups function).
I will likely also implement further refactors and minor bug fixes after this work. But broken down across multiple PRs to limit scope.
Look forward to hearing thoughts! And feel free to tag anyone else you'd like to include.