Make Milo differential abundance testing match R Milo - #1109
Merged
Merged
Conversation
The edger solver called glmQLFit without legacy=TRUE, so under edgeR 4 it ran the new quasi-likelihood method while miloR pins the legacy one. On identical neighbourhood counts logFC and p-values now agree with miloR::testNhoods to machine precision. The pydeseq2 solver rebuilt contrasts from the formula string and fell back to the last against the first level whenever model_contrasts had no "-". COVID_severityAsymptomatic therefore tested Critical against Healthy, and a continuous covariate reported its effect over its whole range instead of per unit. It now tests the coefficient vector from _contrast_vector like the edger and mixed model paths, including their "+ 0" rule for contrasts. make_nhoods left every index cell out of its own neighbourhood because scanpy connectivities have an empty diagonal, while miloR includes it. Each neighbourhood thereby lost one cell of the index cell's sample, which biased the counts towards the null. Verified against miloR 2.6.0 and edgeR 4.8.2 on the Stephenson COVID data (113 donors, 4307 neighbourhoods). A benchmark with permuted labels and planted depletion shows the pydeseq2 solver is better calibrated than edgeR (3.2 against 14.8 false calls at SpatialFDR < 0.1 under permuted labels) but less powerful (TPR 0.31 against 0.50 with ten donors per group). Cook's outlier handling, size factors, min_mu and fit_type did not change this, so the pydeseq2 settings are left as they were.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1109 +/- ##
==========================================
+ Coverage 79.97% 79.99% +0.01%
==========================================
Files 55 55
Lines 7536 7518 -18
==========================================
- Hits 6027 6014 -13
+ Misses 1509 1504 -5
🚀 New features to boost your workflow:
|
The pydeseq2 solver used to accept StatusCovid-StatusHealthy in ~Site+Status, where the reference level Healthy has no coefficient; the contrast vector rewrite rejected it. A term naming a reference level now weighs zero, which is its value under treatment coding, so the old contrasts keep working for the pydeseq2 solver and the mixed model.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This aligns
Milo.da_nhoodsandmake_nhoodswith miloR 2.8.1, the current Bioconductor release: the edger solver uses the legacy quasi-likelihood fit that miloR pins, the pydeseq2 solver tests the same coefficient as the other solvers instead of falling back to the last level for contrasts without-and scaling continuous covariates by their range, and index cells now belong to their own neighbourhood.The table compares pertpy with miloR on identical neighbourhood counts of the Stephenson COVID data (113 donors, 4307 neighbourhoods), plus one end-to-end run for the index cell change.
The remaining pydeseq2 differences come from the more conservative Wald test of DESeq2, e.g. 128 of the 144 Asymptomatic neighbourhoods that miloR calls draw at least 80% of their Asymptomatic cells from a single donor.
~Status~Status~Site+Status~Site+Status~Site+COVID_severity,COVID_severityAsymptomatic~Site+COVID_severity,COVID_severityAsymptomatic~Site+COVID_severity_continuous~Site+COVID_severity_continuous~Status, end to end on 15k cells with k = 30