wilcoxon DE, vs_all_perturbed

@ilan-gold

What kind of feature would you like to request?

Additional function parameters / changed functionality / changed defaults?

Please describe your wishes

@ilan-gold another part of my optimizations I added while doing perturb seq benchmarks was including a mode in wilcoxon DE for vs_all_perturbed. This has emerged as a technique in perturb seq: when performing differential expression on each target perturbation, instead of using all other cells as a reference, this approach uses all other cells except the controls (so only cells that had a perturbation on a coding gene).

e.g. it's introduced here: https://arxiv.org/pdf/2506.22641 in section 3.3 and 4.2 , and you can see it in the corresponding code: https://github.com/shiftbioscience/Diversity_By_Design/blob/main/data/norman19/get_data.py#L40 . This is also used in the more recent Dynamic Range Fraction paper in perturb seq, and probably some other papers derived from these works.

Practically I think this is most relevant when non-targeting (control) cell count is high relative to total cell count. Which does tend to happen, because every batch tends to include non-targeting cells to compensate for batch effects.

That said, it seems it might be a pattern that is here to stay, e.g. it's superior under certain conditions.

So my first question is -- does this deserve support in the scanpy wilcoxon DE implementation? My alternative was to filter data beforehand, but that is clunky.

If so:

This should use a one-versus-rest approach where a single ranking is performed on the 'rest' cells, rather than re-doing ranking on every perturbation. That was the hoist issue I fixed in the other PR (the 2 line change). So this would just need to hit that path.
I can spearhead getting this implemented end-to-end in scanpy and illico.
As for how, I'm thinking it might make sense if there's a parameter that allows specifying a single group to exclude from 'rest' when using 'rest' as the the reference? It should be clear from the parameter name that this only impacts the reference cells and is unrelated to 'groups' (which perturbations DE is computed for)... This seemed a better option than adding another mode (e.g. vs_all_perturbed) because there are non-perturb-seq uses of wilcoxon DE. It also didn't seem necessary to allow multiple exclusion groups, I can't think of any other cases where it'd be helpful to offer to allow excluding more than one group. And this exclusion approach should be a lightweight change.

Thoughts?

If you think this is right, I can draft a PR to show what it looks like in scanpy and ping the illico author(s) to coordinate on that end.

@remydubois

As for how, I'm thinking it might make sense if there's a parameter that allows specifying a single group to exclude from 'rest' when using 'rest' as the the reference? It should be clear from the parameter name that this only impacts the reference cells and is unrelated to 'groups' (which perturbations DE is computed for)... This seemed a better option than adding another mode (e.g. vs_all_perturbed) because there are non-perturb-seq uses of wilcoxon DE. It also didn't seem necessary to allow multiple exclusion groups, I can't think of any other cases where it'd be helpful to offer to allow excluding more than one group. And this exclusion approach should be a lightweight change.

I think a general exclude: list[str] parameter could solve this right" Because they aren't compared either as reference or the other side (perturbed in this case), they are completely excluded. @remydubois thoughts?

Definitely belongs in both scanpy but might be able to get away without illico since we construct the AnnData object ourselves

@ilan-gold

Hey, interesting use case.

I think a general exclude: list[str] parameter could solve this right

I agree it would solve it. I'm wondering if users wouldn't get confused by how this interacts with groups. Also, it would ask a bit of guardrailing on scanpy's side so as to make sure values inputted by the user for those two don't contradict each other.

As for whether illico should handle some (or all) of it, here is how i see things:

The groups argument (the functionality it currently has in scanpy) would shorten execution time of the OVO (because there would be less comparisons to perform) but not the OVR path (because the costly operation is sorting, and we would sort the same values as if groups=None), if implemented in illico.
However, the exclude argument that you mention would speed up the OVO path (same reasons as above) and the OVR path (because there would be less values to sort).

I could be open to adding the exclude functionality in illico in the short term, which would not be contradictory with what's currently in the PR @ilan-gold IIUC (where the filtering of resulting p values is handled by scanpy because the groups keyword is not supported). But maybe there's a better way (like filtering beforehand on scanpy's side).

@remydubois

@remydubois

it seems to me this would only impact the case where reference is 'rest' -- the idea is to do OVR with control cells excluded from 'rest' / some reduction in the size of 'rest'. So it wouldn't impact OVO.

@ilan-gold @remydubois

that said, maybe in scanpy it could be reference_exclusion: list[str] ? listed right after 'reference' in the docs, with a comment like so: "if reference of 'rest' is used, any groups specified here are excluded from the reference union. Note that comparisons will still be performed for all groups in 'groups'; this parameter does not impact 'groups'."

and I assume we'd have a warning if it's specified when reference != 'rest', and also probably clarify it's role in stats computation in docs as well. And a little bit of guardrails to prevent specifying non-existing groups, or entirely removing the reference (all groups specified)... Would stick to scanpy patterns for these types of things.

and yes as for implementation, at least with the current pattern, the data could be filtered when the new AnnData object is created in scanpy, or something like that. I just wasn't sure if it would make sense to stop at that layer since it seems illico and scanpy are on a path to being tightly intertwined, and not passing this thru leads to later compatibility work (it's work either way I guess haha).

it seems to me this would only impact the case where reference is 'rest' -- the idea is to do OVR with control cells excluded from 'rest' / some reduction in the size of 'rest'. So it wouldn't impact OVO.

Yes, agreed, sorry my answer wasn't clear. My point was more: we could also use this argument for the OVR path (at least on illico's side) as a way to filter groups. The motivation would be that one single argument (on illico's side) would serve the goal of two (both groups and exclude), considering illico does not implement group. But I'm aware that's kind of an api change (or at least usage change) from what you described and I dont have a strong opinion on that (nor would implement it if it didnt make sense for scanpy). My angle of reading is more "what could illico do to speed things up" considering that's its main reason to be, with a reasonable implementation budget (typically, one argument instead of two).

The only drawback I could see in rebuilding the AnnData manually is that it would maybe break the lazy loading capacity ? I'm not so sure about that though.

The motivation would be that one single argument (on illico's side) would serve the goal of two (both groups and exclude), considering illico does not implement group

Taking a step back, maybe a return_dict argument would be applicable? What about this:

adata.uns["rank_genes_groups"] = scanpy.tl.rank_genes_groups(adata[subset_without_control], **kwargs, return_dict=True)

This is basically what is done in https://github.com/shiftbioscience/Diversity_By_Design/blob/main/data/norman19/get_data.py#L40

illico and scanpy are on a path to being tightly intertwined

I put some more thoughts into this and I tend to agree (now) with your statement above. At first, I was a bit reluctant to implement those exact groups and exclude into illico because I was worried of increased codebase complexity for little added value, but I threw an attempt at it and was happy to get something functional with not too much changes, so I'm okay with having them in illico (PR is here, only targeting the Numba backend to get an idea).

Now, whether scanpy or illico should handle it, I see pros and cons for both:

One reason I see to keep both handled by scanpy is that it remains more adaptative. Typically, scanpy might see some new use-cases similar to the one you raise here coming up now and then, every time a bit different but likely handleable with some AnnData slicing tricks or results slicing.
One reason to rely on illico for exclude is that slicing at scanpy's level breaks the lazy-loading of AnnData: ad.read_h5ad(…, backed="r")[indices_of_selected_groups,].X load the entire indices_of_selected_groups slice of .X into memory. In the use case of excluding just one group, it boils down to loading almost all of it in RAM.
One reason to rely on illico for groups is that for the OVO path, it would bring a speedup. On huge datasets like TAHOE, this can be quite welcome, as even if illico is fast, it still takes ~1 minute to process a plate.

One reason to rely on illico for exclude is that slicing at scanpy's level breaks the lazy-loading of AnnData: ad.read_h5ad(…, backed="r")[indices_of_selected_groups,].X load the entire indices_of_selected_groups slice of .X into memory. In the use case of excluding just one group, it boils down to loading almost all of it in RAM.

I agree about this. My only question here would be if people are dumping log1p-ed data to disk or if we should be encouraging them to use dask. In that case, this becomes less critical because dask is lazy.

One reason to rely on illico for groups is that for the OVO path, it would bring a speedup. On huge datasets like TAHOE, this can be quite welcome, as even if illico is fast, it still takes ~1 minute to process a plate

Why?

@ilan-gold

dask is lazy

Very true, Dask-read Anndata wouldn't suffer from that.

Why?

Woops, forgot to add "it would bring a speedup if the user wants to only compare a small subset of groups to controls". But after second thoughts (for the groups argument):

For the OVO case, scanpy could filter the AnnData beforehand (instead of filtering results after), and the same (as if it was implemented in illico) speedup would (theoretically) be achieved
For the OVR case, there is no speedup (or very little) anyways, so filtering results afterwards is good enough.

All in all in don't have a strong opinion on either. I'm happy to see illico enriched with new features but I wouldn't want this to become a blocking point for a scanpy release. One way to ensure no blocking happens would be to implement those on scanpy's side - if needed, (groups already is, yet non-optimally but fully functionally) and when those two get into illico (if they ever do) we could rediscuss.

wdyt @ilan-gold ?

@remydubois

I'm happy to see illico enriched with new features but I wouldn't want this to become a blocking point for a scanpy release.

My general feeling is that if it can be done by scanpy, it should be done by scanpy to keep the burden off illico's API. But if you want to implement this @remydubois, that would of course be great. I'm quite happy with the level of integration now in any case.

I've asked @zboldyga to take a look at our dask situation with rank_genes_groups - I'll be on holiday for 1.5 weeks, but when I'm back, now that the anndata release is out, this should be my top priority.

Finally, returning to my original suggestion of

adata.uns["rank_genes_groups"] = scanpy.tl.rank_genes_groups(adata[subset_without_control], **kwargs, return_dict=True)

as long as dask is fast, does this seem reasonable? The subsetting / lop1p would be lazily applied so everything should be efficient (theoretically).

@ilan-gold

Hey,

as long as dask is fast, does this seem reasonable?

Yes, I agree it should be fine. I added both groups and exclude in illico (cf this PR) bc I needed groups for my personal use, but I agree with your rationale.

One last thing: I realized that regardless of whether those are handled by the scanpy layer or by illico, we still need minimalistic support for Dask in illico (which would be achieved by simply adding a new DataHandler). I think adding such a handler would be underoptimal (cf the discussion we had here "Is what is implemented currently dask-ifiable" @ilan-gold ) but would be Ok as a first step.

@zboldyga if you want to chip in on illico's PR or extend Dask support, I welcome all contributions ! Something I could definitely get help on is triple check that the sanity checks ran on those two args to make sure they don't interact badly with each others (illico.utils.groups.sanitize_group_args) are fine.

@remydubois

Nice plan @remydubois - it seems then that illico needs dask support (for sparse and dense subchunks?) and then we can assess here. But that can proceed independent of the PR open here ATM. I'm back from vacation so we can get moving again.

I was happily surprised how easy it was to integrate a dask backend into a joblib's parallel loop (following this simple recipe) so I gave it a try. I wanted to see how far we could get beyond the naive dask integration I described earlier.
I tested the three implementations of the 'groups' functionality (selected 500 groups out of 2058 in k562):

Input is a sliced (masked) dask array but the parallelization remains thread-based (labeled subset-dask-adata)
Input as a sliced dask array and the parallelization is handled by dask (following the recipe above) (subset-dask-adata-over-dask-backend)
Input is the (full) lazy h5 AnnData, and groups are specified in the groups arguments implemented in illico´s PR I linked above (h5-groups). Parallelization is thread-based.

The result is that option 2 is the fastest by a significant margin. In other words, slicing (masking) the anndata before running illico ´s kernels is the fastest option (with this backend configuration, which would be embarassingly simple to support in illico).

Shall we get the Dask related discussion in a dedicated issue/PR ?

Note: I assume those results hold for 'exclude' as it works similarly.
Note2: joblib's recipe only works for single-machine setup, but that's already a good start IMO.

The result is that option 2 is the fastest by a significant margin. In other words, slicing (masking) the anndata before running illico ´s kernels is the fastest option (with this backend configuration, which would be embarassingly simple to support in illico). Shall we get the Dask related discussion in a dedicated issue/PR ?

Yes please, tag me! I am not sure I 100% follow so code/pseudocode would be much easier.

Yes, I agree it should be fine. I added both groups and exclude in illico (cf remydubois/illico#20) bc I needed groups for my personal use, but I agree with your rationale.

BTW The other PR looks good. So maybe scanpy doesn't need a new kwarg then. Thanks!

Navigation Menu

Repository navigation

wilcoxon DE, vs_all_perturbed #4114

What kind of feature would you like to request?

Please describe your wishes

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

wilcoxon DE, vs_all_perturbed #4114

Description

What kind of feature would you like to request?

Please describe your wishes

Activity

ilan-gold commented on May 12, 2026

remydubois commented on May 12, 2026

zboldyga commented on May 13, 2026

remydubois commented on May 13, 2026

ilan-gold commented on May 13, 2026

remydubois commented on May 16, 2026

ilan-gold commented on May 16, 2026

remydubois commented on May 18, 2026

ilan-gold commented on May 20, 2026

remydubois commented on Jun 2, 2026

ilan-gold commented on Jun 3, 2026

remydubois commented on Jun 5, 2026

ilan-gold commented on Jun 6, 2026

ilan-gold commented on Jun 6, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions