Add warning for rank_genes_groups

[grst]

And once again I saw I scientific presentation with "1000s of DE genes" and a "surprising correlation between #cells and p-values". I really thought we were over this, but I think we should make more noise about the fact that tl.rank_genes_groups is a bad idea for most use-cases.

Therefore, I suggest to add a warning message when using it.

CC @maltekuehl

• [ ] Closes #

• Tests included or not required because: no change to program logic

• Release notes not necessary because:

[codecov]

Codecov Report

❌ Patch coverage is 0% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 20.49%. Comparing base (deb46bf) to head (4bd458e).

Files with missing lines	Patch %	Lines
src/scanpy/tools/_rank_genes_groups.py	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3792      +/-   ##
==========================================
- Coverage   20.49%   20.49%   -0.01%     
==========================================
  Files         116      116              
  Lines       12436    12437       +1     
==========================================
  Hits         2549     2549              
- Misses       9887     9888       +1     

Files with missing lines	Coverage Δ
src/scanpy/tools/_rank_genes_groups.py	8.13% <0.00%> (-0.03%)	⬇️

[ilan-gold]

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

[grst]

• 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

[maltekuehl]

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.

[Zethson]

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.

[maltekuehl]

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.

[Zethson]

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.

[grst]

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

[grst]

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

[ilan-gold]

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?

[grst]

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