Skip to content

Add warning for rank_genes_groups #3792

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add warning for rank_genes_groups #3792

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ac8ae0f
Add warning for rank_genes_groups
grst Sep 10, 2025
4bd458e
Add changelog
grst Sep 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions docs/release-notes/3792.fix.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Add a warning message to :func:`~scanpy.tl.rank_genes_groups` that it is not
appropriate for most use-cases.
6 changes: 6 additions & 0 deletions src/scanpy/tools/_rank_genes_groups.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -665,6 +665,12 @@ def rank_genes_groups( # noqa: PLR0912, PLR0913, PLR0915
msg = f"reference = {reference} needs to be one of groupby = {cats}."
raise ValueError(msg)

logg.warning(
"Comparing between cells leads to highly inflated p-values, "
"since cells are not independent observations (doi:10.1038/s41467-021-25960-2). "
"Consider using more appropriate methods such as pseudobulk+PyDESeq2 instead."
)
Comment on lines +668 to +672
Copy link
Contributor

@ilan-gold ilan-gold Sep 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a scenario in which you would want to use rank_genes_groups?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • if you have independent observations in your AnnData (e.g. bulk), it would not be wrong to use it (but other methods like PyDESeq2 would still be preferred)
  • It can be useful as a "first glance"/qualitative analysis to check cell-type markers or top genes between conditions.

I don't think it should be included in scanpy 2.0 in its current form.

Happy to hear other opinions, e.g. from @Zethson

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes and no. I think it is valuable for scanpy to have something built in, that allows you to check DEGs without requiring other packages. E.g., I sometimes use it together with dotplots to quickly check that my clusters capture some cell types I expect in the data as a sort of soft quality control step when just quickly exploring some data in a notebook, before running a proper analysis. Theoretically, there is also the use case of having an AnnData object where the observations are samples, not cells. In that case, it would still not be an ideal solution for single-cell RNA-sequencing, but there are other modalities that could be represented as AnnData objects for which the assumptions of PyDESeq2 and other algorithms don't hold. I am however unconvinced that anybody actually uses it like that and this benefit may not outway the risk of having a incredibly high false positivity rate on many scRNAseq datasets if this function is left in and people end up using it.

Copy link
Member

@Zethson Zethson Sep 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use case 2 that Gregor mentioned immediately came to my mind and there's also benchmarking papers that back that up (Wilcoxon does well).

I agree with Gregor that something has to be done. A warning can be pretty intrusive though because lots of downstream scverse ecosystem tools use it and some even loop over it which would lead to tons of clutter and warnings. I wonder whether an alternative would be to have a rudimentary optional and simple pydeseq2 interface while adding a fat warning in our API docs & tutorials about the usage of the simpler rank_genes_groups options.

I think that scanpy users are fine with optional dependencies. Leiden & igraph etc also have to be installed explicitly.

Copy link
Contributor

@maltekuehl maltekuehl Sep 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Having scanpy not cover DEGs (arguably the main output of any scRNAseq analysis) feels like an important limitation of the package. From reading raw data up to clustering and annotating clusters, everything can be done in pure scanpy/AnnData (at least a basic version of it). But the main desired result cannot be produced and it is quite unclear from the scanpy documentation how it can be achieved. Searching for "pseudobulk" in the scanpy docs yields just one uninformative result. The best practices guide does not cover the current state-of-the-art. Decoupler's documentation does cover it, but that is a different package. Since scanpy will be used by many beginners and biologists venturing into computational biology, this feels a bit inadequate.

Edit: Limited wifi on train, apologies for repeated points and messed up message order. I did not see anything before sending my comments.

Copy link
Member

@Zethson Zethson Sep 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would a notebook in scanpy as a stopgap make sense? Then we can move it to sc best practices? The notebook could use PyDESeq2 and/or Pertpy.

Overhauling the DE content in BP has been on our radar for a while! It will happen soonish, trust me haha. The best practice notebook also has a lot more details and references that shouldn't matter for a scanpy tutorial. I'd keep them separate for now.

Although I just merged a PR that removes scvi-tools from pertpy which makes it much leaner, I would not build a scanpy tutorial on pertpy. I agree with Malte that a native scanpy solution is useful and our users kind of expect it.

I think a future with a pure and powerful DE package, which Gregor envisions, and a simple scanpy solution, is a good future.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we'd need decoupler in this equation. Pseudobulk is now handled by scanpy and then use PyDESeq2 from there.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we all agree what a good solution looks like, the issue here is more how we get the users to adopt them.

Copy link
Contributor

@ilan-gold ilan-gold Sep 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The best practice notebook also has a lot more details and references that shouldn't matter for a scanpy tutorial.

+

I would not build a scanpy tutorial on pertpy

Sounds like we should do a scanpy tutorial then built on PyDESeq2 + scanpy?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sounds good. And add warning boxes with references to that tutorial in the 3k PBMC tutorial and the rank_genes_groups function


if key_added is None:
key_added = "rank_genes_groups"
adata.uns[key_added] = {}
Expand Down
Loading