feat: add illicofor rank_genes_groups

[ilan-gold]

TODOs:

See: https://github.com/scverse/scanpy/actions/runs/24088645078/job/70268566419?pr=4038

• Investigate seemingly mismatched ovo "column" in scores (cc @remydubois)
• Maybe get a compat layer for the 1e-9 helper in illico (see the filterwarnings on the new test)
• Investigate why ovr results are so far off
• Once we can match our current implementation (along with the new exp_post_agg feature feat: allow exponentiation post agg for log-fold-change in rank_genes_groups #4037) to illico, we will add the illico backend (or all changes have been documented)

---

• Closes Integration with illico #4012
• Tests included or not required because:

• Release notes not necessary because:

[codecov]

⚠️ JUnit XML file not found

The CLI was unable to find any JUnit XML files to upload.
For more help, visit our troubleshooting guide.

[remydubois]

Hey, good catch for the OVO column ordering. I have no idea where the bug comes from, it seems to be related to group labels for PBMC that numpy sorts in a weird way. I will try to get a fix in the coming days.

[remydubois]

So it seems to be due to np.argsort and np.sort not sorting the weird characters (" ", "+", "/", "#", etc) of the bulk_labels the same way. I will ship a fix sometime today or this weekend.

NB: I spent a bit of time looking into the PBMC dataset and it seems rather unusual: a lot of values are negative leading to non-defined logfoldchanges, on top of which there seems to be very little value diversity in the dataset. As a result, a lot of genes end up with identical ranksum, but because illico and scanpy do not compute z-score the exact same way (mathematically equivalent but programmatically different), gene ordering is impacted. The gene ordering impacts the position of NaN values in the results, but non-NaN values do match because they are very similar. I am not sure as if this dataset is a good test case.

NB2: details of the programmatic differences:
illico does:

U = ranksum - n_tgt * (n_tgt+ 1)/2
z = U - n_ref * n_tgt / 2

scanpy does:

z = 
ranksum - n_tgt * (n_tgt + 1 + n_ref) / 2

Those are mathematically equivalent but result in differences of the order of 1.e-9 approximately, changing gene orders when ranksums are equal.

[ilan-gold]

I am not sure as if this dataset is a good test case.

Interesting, good observation but glad we are aware of this now. This could be another target for scanpy 2.0.

I will try a different dataset but glad to have this documented.

[remydubois]

Forget my previous message which was actually kind of off topic.

This PBMC dataset was actually a very good test case.

1. It unveiled a silent bug in illico as it seems np.argsort and np.sort are not sorting the weird characters (" ", "+", "/", "#", etc) of the bulk_labels the same way. I will ship a fix in the next few days.
2. Because it has very limited diversity in the data (the sheer values in .X are not very diverse), a lot of genes end up with identical ranksums (see below example), hence, identical z-scores.
3. The sorting method used by scanpy is more elaborate as it allows user to select only top n genes. As a result, it does not sort identical values in the same order as illico (even if all genes are returned). That explains the genes ordering difference that I was still facing after fixing the name sorting issue described in 1. I believe I can add the functionality to return n_genes only like scanpy does, and implement the same sorting routine as scanpy, which appears to solve that issue.
4. Not related to PBMC, but z-scores also mismatched because scanpy casts them to float32, and illico keeps thems as float64. I will fix that in the next patch as well.
5. For the 1.e-9 adjustment, I don't know what's the best way to go. Currently, illico handles this with a np.where(mu_ref == 0, np.inf, mu_tgt / mu_ref) but I'm quite open to even change what's currently in illico if this e-9 offset is the de-facto standard in other softwares like Seurat or so.

Example: let's take CD14+ Monocyte as control group, CD34+ as perturbed group, and look at the genes CDK6 and TCAEL8. Although they don't exactly have the same values, both CDK6 and TCEAL8 have equal ranksums, hence, equal z-scores. Due to the sorting methodology difference, they end up not sorted in the same order between illico and scanpy.
From the next release, the test suite in illico will not only test closeness of the values but also explicitly test matching ordering of the genes.

NB: what is the situation with the Dask issue ? Current version is not dask-compatible.

[ilan-gold]

NB: what is the situation with the Dask issue ? Current version is not dask-compatible.

That should come next, I want to keep this scoped for now to the in-memory stuff.

Another thing if you're doing a batch of fixes this weekend - supporting cs{r,c}_array would be awesome. Should only be one line of code.

Also for your numba types, only if you want instead of named tuples, FAU has a plug-in for handling putting cs{c,r}_{array,matrix} into numba kernels: https://github.com/scverse/fast-array-utils/blob/main/src/fast_array_utils/_plugins/numba_sparse.py. But it might not be worth it to bring on the dependency just for this purpose

[ilan-gold]

For the 1.e-9 adjustment, I don't know what's the best way to go. Currently, illico handles this with a np.where(mu_ref == 0, np.inf, mu_tgt / mu_ref) but I'm quite open to even change what's currently in illico if this e-9 offset is the de-facto standard in other softwares like Seurat or so.

I'm seeing some instances of https://github.com/satijalab/seurat/blob/main/R/differential_expression.R#L1089 i.e., they calculate the mean expression with this "+1" offset so don't need the correction.

So we have three different things here, I guess. I think just giving the option to match scanpy (since that seems a project goal) is good, and then we can maybe revisit what seurat does as a new parameter for scanpy 2.0.

Not related to PBMC, but z-scores also mismatched because scanpy casts them to float32, and illico keeps thems as float64. I will fix that in the next patch as well.

For this, I would be open to also making this part of a scanpy 2.0 set of default i.e., using float64, as there are other places we might benefit from this. So if you were to keep float64, I don't think that would be a blocker. Again, an option could be good so we can smooth out the transition, although this is so small as to maybe not warrant it.

Overall, my ideal scenario would be that illico and scanpy with wilcoxon match, and then we use illico by default in scanpy 2.0. I'm pushing to minimize result changes to make the transition smoother (since results changing are always no fun) and so far the changes don't seem particularly destructiv.

Short of exact matches, like I said we can just document changes. It sounds like the biggest things that would close the gap are:

1. argsort bug
2. 1e-9 option
3. The option to return n_genes
4. cs{r,c}_array support

I think we can patch the float{32,64} stuff in a general "numerical accuracy" scanpy 2.0 preset (which we have on main right now).

Overall, this is really cool and I'm super happy that things seem to be pretty close!

[remydubois]

Hey,

I just released version 0.5.0rc1 which should fix the issues identified on this test case:

1. Perturbation (or group) names are no longer re-sorted when outputs are formatted for scanpy ensuring they are ordered the same everywhere.
2. Fold change is now computed by adding the same 1.e-9 factor, on top of being accumulated into a f64 placeholder (regardless of the original data dtype) ensuring a better matching.
3. New argument n_genes allowing to return only top n DE genes per perturbation (which, even if not specified, results in genes sorting methodology being identical. This solves the issue raised by genes having identical scores).
4. Support for cs[cr]_array. I did not explicitely add test cases for those as the test suite is already quite heavy and illico does nothing more than accessing .data, .indices, and .indptr of those objects. I did test it manually quickly and it ran with no issue.

Regarding the numerical precision question: so far I force the recarrays to have the same dtype as what's currently implemented in scanpy, see FC for instance. I believe we could change that when the need comes as you say.

Testing locally, it seems everything now runs smoothly for PMBC. Let me know how it goes with the rest of the CI.

[ilan-gold]

@remydubois I think the issue is that you are sorting / cutting off by n_top genes globally but scanpy does this per group, but that is just a guess.

In any case, I think I got a little lost-in-the-sauce. I think matching p-values and scores should be enough because scanpy can internally handle LFC when it wraps illico. I just pushed a commit with this change. Of course, if you want to match LFC, no argument from me, but I realized it's probably not essential for internal consistency.

I'm pushing the results of scores and pvals_adj which all match except the two cases highlighted in the tests.

With this in mind, I am going to push ahead with integrating illico, sorry for the churn on the log fold changes, but I think we can just let scanpy do it as along as I can get scores/pvals out of illico.

BTW: https://remydubois.github.io/illico/api.html looks both incomplete and like it's leaking internals.

[ilan-gold]

Ok everything passing locally with just getting P-value and z-score from illico. This is looking good now to me at least conceptually. Last thing would be to make sure mean-var calculation is fast since we won't rely on illico (at least for now) to do logfc, but I have something for that in #4041.

I assume the LFC calculation (i.e., actually calculating the means) is the less computationally intensive part?

[remydubois]

With this in mind, I am going to push ahead with integrating illico

Ok ! That's good news to me.

Because this was bugging me out a little bit: I checked out on your ig/illico branch and indeed reproduced test failure, until I added a pbmc = ad.AnnData(pbmc.X.copy(), obs=pbmc.obs.copy(), var=pbmc.var.copy()) (in order to make sure pmbc gets stripped out of all of its metadata or attributes):

def test_illico(test, corr_method, exp_post_agg):
    from illico.asymptotic_wilcoxon import asymptotic_wilcoxon

    pbmc = pbmc68k_reduced()
    pbmc = ad.AnnData(pbmc.X.copy(), obs=pbmc.obs.copy(), var=pbmc.var.copy())
    # ... Rest of the test. They all pass on my machine

Then, all tests passed (LFC, pvals, scores, pvals_adj, up to atol=1.e-9) which I still consider good news, even with in mind that LFC will remain scanpy-computed (non-matching LFC might mean non-matching names, which would have been a real issue). The reason why my "tests" passed locally before shipping 0.5.0rc1 for PBMC was exactly because I was stripping it entirely before running and comparing outputs, just to be sure.

By the way, wouldn't it be safer for scanpy's test suite to add an explicit check on the ordering of the genes, like I did here (which did fail, before stripping pbmc) ?

This is looking good now to me at least conceptually

Very happy to read it's moving forward 🚀

I assume the LFC calculation (i.e., actually calculating the means) is the less computationally intensive part?

Yes it's neglectable, I could not measure its impact on the overall runtime.

[zboldyga]

@ilan-gold curious why this is hardcoded?

[ilan-gold]

From what @remydubois said on zulip: https://scverse.zulipchat.com/#narrow/channel/315570-random/topic/Article.20on.20speeding.20up.20computation/near/581587388 it is basically not worth the headache. Furthermore, if we were to ever adopt the codebase, it would be the numba version, not the rust version. His observations about Rust largely match ours (the performance is almost always identical to numba)

[ilan-gold]

@remydubois Everything seems to pass with groups - I think exclude_from_ovr + dask should be follow-up PRs because both apply to all tests, not just wilcoxon.

[remydubois]

@ilan-gold good news !

I think exclude_from_ovr + dask should be follow-up PRs

That makes sense. Does it mean any specifics for illico releases (i.e: should I release/not release some features in the v0.6.0) ?

[ilan-gold]

That makes sense. Does it mean any specifics for illico releases (i.e: should I release/not release some features in the v0.6.0) ?

I think you should be pretty good releasing what you want - we can always catch up here e.g., when we do a proper assessment of dask.