diff --git a/README.md b/README.md index ecffd7d..d77f7c8 100644 --- a/README.md +++ b/README.md @@ -249,7 +249,7 @@ is the reference. Counts are aggregated by taxon after applicable read-count bounds to the differential-abundance request, so those bounds do not exclude samples from this analysis. -All `analysis.differential` options below (`method`, `offset`, +All `analysis.differential` options below (`method`, `offset`, `dispersion`, `min_prevalence` and `replacement_delta`) are resolved at study scope. Set them in the study configuration; group and run overrides are not used by this endpoint, including when called from a group Analysis workspace. @@ -259,7 +259,15 @@ by this endpoint, including when called from a group Analysis workspace. - `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`, with a library-size offset. The effect is log2 fold change; positive values indicate greater abundance in the second condition. `offset` selects `tss` - (default) or `rle` scaling. + (default) or `rle` scaling. `dispersion` selects the fits and the test: + `glm_nb` (default) fits `MASS::glm.nb` per taxon with its own + maximum-likelihood theta and Wald-tests the group coefficient; `glmGamPoi` + fits one `glmGamPoi::glm_gp` across the taxa passing `min_prevalence`, + sharing a dispersion trend across them, and reports `glmGamPoi::test_de`'s + quasi-likelihood F-test. Under `glmGamPoi` the p-value comes from that + F-test, the standard error is `test_de`'s `lfc_se`, and theta is each + taxon's own estimate before shrinkage. It requires the pinned glmGamPoi R + package and is refused for `clr_welch`. - `clr_welch`: Welch's t-test on centred log-ratios after multiplicative zero replacement with `zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change. `replacement_delta` controls the replacement fraction diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 2181a23..debcff0 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -363,9 +363,10 @@ analysis: # the group coefficient. glmGamPoi fits one glmGamPoi::glm_gp across the # taxa passing min_prevalence, with overdispersion shrinkage toward the # cross-taxon trend, and reports glmGamPoi::test_de's quasi-likelihood - # F-test of the group contrast with its log2 fold changes. A zero shrunken - # dispersion withholds inference for its taxon. Requires the pinned - # glmGamPoi R package. Refused with clr_welch. + # F-test of the group contrast with its log2 fold changes; the reported + # theta inverts glm_gp's negative-binomial overdispersion, and a taxon + # whose overdispersion is zero is marked at the Poisson boundary, as under + # glm_nb. Requires the pinned glmGamPoi R package. Refused with clr_welch. dispersion: glm_nb # Fraction of samples, 0 to 1, in which a taxon must have reads to be # tested. Taxa below it are listed as filtered, without a p-value. For diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e6ee284..fff182b 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -439,7 +439,7 @@ export interface DifferentialResult { min_prevalence: number /** nb_glm only. */ offset?: 'tss' | 'rle' - /** nb_glm only: where each taxon's theta comes from. */ + /** nb_glm only: which fits and test produced the rows (per-taxon MASS::glm.nb with a Wald test, or one glmGamPoi::glm_gp with test_de's quasi-likelihood F-test). */ dispersion?: 'glm_nb' | 'glmGamPoi' /** clr_welch only. */ replacement_delta?: number diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index 5e46c3f..f87d7bc 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -291,27 +291,32 @@ da_result <- data.frame(status = rep("ok", da_n), note = rep("", da_n), estimate = rep(NA_real_, da_n), se = rep(NA_real_, da_n), statistic = rep(NA_real_, da_n), pvalue = rep(NA_real_, da_n), theta = rep(NA_real_, da_n), stringsAsFactors = FALSE) -# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa supplies each taxon's -# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa, with each taxon's -# overdispersion shrunk toward the cross-taxon trend, supplies glmGamPoi's own -# published inference -- glmGamPoi::test_de's quasi-likelihood F-test of the -# group coefficient. Its genes x samples convention, so -# the counts are transposed, and size_factors takes the linear size factor -# (the per-taxon fits below use its log as their offset). A fit- or test-level -# error fails every taxon. No fallback to another dispersion estimator or -# significance test is permitted. +# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa supplies glmGamPoi's +# own published inference -- glmGamPoi::test_de's quasi-likelihood F-test of +# the group coefficient. glm_gp refits the coefficients at its cross-taxon +# dispersion trend and shrinks each taxon's quasi-likelihood dispersion by +# empirical Bayes. glm_gp works genes x samples, so the counts are transposed, +# and size_factors takes the linear size factor (the per-taxon fits below use +# its log as their offset). The contrast reaches test_de as numeric weights on +# the design's columns, never as text: test_de parses and evaluates a text +# contrast as R code, so a condition label such as "run-B" or "grpB/run_2" +# would stop it, and a label that reads as R code would run. A fit- or +# test-level error fails every taxon. No fallback to another dispersion +# estimator or significance test is permitted. da_glp <- NULL da_glp_failed <- "" if (da_use_glp) { + da_design <- model.matrix(~ da_g) + da_contrast <- as.numeric(colnames(da_design) == da_term) da_gp <- tryCatch( - glmGamPoi::glm_gp(t(da_counts), design = model.matrix(~ da_g), + glmGamPoi::glm_gp(t(da_counts), design = da_design, size_factors = exp(da_offset), overdispersion_shrinkage = TRUE), error = function(e) e) if (inherits(da_gp, "error")) { da_glp_failed <- paste("glm_gp stopped with an error:", conditionMessage(da_gp)) } else { da_glp <- tryCatch( - glmGamPoi::test_de(da_gp, contrast = da_term, max_lfc = Inf, compute_lfc_se = TRUE), + glmGamPoi::test_de(da_gp, contrast = da_contrast, max_lfc = Inf, compute_lfc_se = TRUE), error = function(e) e) if (inherits(da_glp, "error")) { da_glp_failed <- paste("test_de stopped with an error:", conditionMessage(da_glp)) @@ -349,20 +354,20 @@ for (da_j in seq_len(da_n)) { # glmGamPoi's own test on its own fit. `lfc` is the log2 fold change of # the group contrast and `lfc_se` its standard error; both are converted # to the natural-log scale so every method reports one effect scale. The - # per-taxon scale the test used is the shrunken quasi-likelihood - # dispersion `ql_disp_shrunken`, whose inverse is reported as the - # dispersion theta. Zero dispersion is the Poisson boundary, where the F - # ratio would report p = 0 for the wrong reason, so inference is withheld; - # an unusable dispersion fails only its own taxon. - da_ql <- da_gp[["overdispersion_shrinkage_list"]][["ql_disp_shrunken"]][da_j] - if (!is.finite(da_ql) || da_ql < 0) { + # reported theta is the inverse of glm_gp's negative-binomial + # overdispersion for the taxon (`overdispersions`: variance = mu + mu^2 * + # overdispersion), the same parameter glm.nb's theta inverts. It is the + # taxon's own estimate before shrinkage: the coefficients and the F-test's + # reduced fit use the cross-taxon trend, and the F ratio divides by the + # shrunken quasi-likelihood dispersion; only lfc_se, from predict(), is + # weighted by this estimate. Zero overdispersion is the Poisson boundary: + # the row is marked, as under glm_nb, and keeps its statistics, since the + # test does not read it. An unusable overdispersion fails only its own + # taxon. + da_od <- da_gp[["overdispersions"]][da_j] + if (!is.finite(da_od) || da_od < 0) { da_result[da_j, "status"] <- "failed" - da_result[da_j, "note"] <- "the shrunken dispersion for this taxon is not usable" - next - } - if (da_ql == 0) { - da_result[da_j, "status"] <- "boundary" - da_result[da_j, "note"] <- "zero dispersion: Poisson boundary; quasi-likelihood inference withheld" + da_result[da_j, "note"] <- "the negative-binomial overdispersion for this taxon is not usable" next } da_est <- da_glp[["lfc"]][da_j] * log(2) @@ -376,7 +381,11 @@ for (da_j in seq_len(da_n)) { next } da_result[da_j, c("estimate", "se", "statistic", "pvalue", "theta")] <- - c(da_est, da_se, da_stat, da_p, 1 / da_ql) + c(da_est, da_se, da_stat, da_p, 1 / da_od) + if (da_od == 0) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "the taxon's negative-binomial overdispersion estimate is zero: these counts show no overdispersion beyond Poisson, so theta is unbounded; the test used the cross-taxon trend, so its statistics stand" + } } else { # "glm_nb" estimates theta per taxon inside MASS::glm.nb's alternating # fit, then Wald-tests the group coefficient at fixed scale 1, which is @@ -453,17 +462,21 @@ toward the cross-taxon trend and reports `glmGamPoi::test_de`'s quasi-likelihood F-test of the group contrast — its `lfc` (log2 fold change) and `lfc_se`, both converted to the natural-log scale so the columns carry one effect scale. In that path `statistic` is the quasi-likelihood F, `pvalue` -its p-value, and `dispersion_theta` the inverse of the shrunken -quasi-likelihood dispersion the test used for the taxon. Names are stripped -and matched without regard to case. +its p-value, and `theta` the inverse of `glm_gp`'s negative-binomial +overdispersion for the taxon (its `overdispersions` element, the taxon's own +estimate before shrinkage), the same parameter `MASS::glm.nb`'s theta inverts; +it is `nothing` where that overdispersion is zero. The contrast is passed to +`test_de` as numeric weights, so condition labels are never parsed as R code. +Names are stripped and matched without regard to case. Returns one row per taxon with `status` (`ok`, `boundary` or `failed`), `note`, `estimate`, `se`, `statistic`, `pvalue` and `theta`; a failed taxon has `nothing` for every statistic. Boundary fits retain their statistics: under `"glm_nb"` theta at or above 1e7 marks the Poisson end and theta at or below -1e-8 marks extreme variance; under `"glmGamPoi"` a zero shrunken dispersion -marks the Poisson boundary, where the F ratio would report p = 0 for the wrong -reason, so inference is withheld. Constant counts, complete separation and +1e-8 marks extreme variance; under `"glmGamPoi"` a zero negative-binomial +overdispersion marks the Poisson boundary, and the row keeps its statistics, +since the F-test divides by the shrunken quasi-likelihood dispersion, which +stays defined there. Constant counts, complete separation and unusable coefficients produce failed rows. Errors caught from the model-fitting calls also produce failed rows; a caught `glm_gp` or `test_de` error fails every taxon, with no fallback to `"glm_nb"`. @@ -500,8 +513,9 @@ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:Abst finally RCall.reval("rm(list = intersect(ls(), c('da_counts', 'da_group', 'da_levels', " * "'da_offset', 'da_theta_upper', 'da_theta_lower', 'da_use_glp', 'da_g', " * - "'da_term', 'da_n', 'da_result', 'da_glp', 'da_glp_failed', 'da_gp', " * - "'da_j', 'da_y', 'da_absent', 'da_ql', 'da_est', 'da_se', 'da_stat', " * + "'da_term', 'da_n', 'da_result', 'da_design', 'da_contrast', 'da_glp', " * + "'da_glp_failed', 'da_gp', " * + "'da_j', 'da_y', 'da_absent', 'da_od', 'da_est', 'da_se', 'da_stat', " * "'da_p', 'da_warns', 'da_fit', 'da_co', 'da_row', 'da_theta'))); invisible(gc())") end end diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 46a6e79..2c686a2 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -295,7 +295,9 @@ end # quasi-likelihood F-test, computed independently. lfc is log2 and # lfc_se its standard error, so both are scaled by log(2) to meet # the natural-log convention of every method's estimate and se. - # One matrix column per quantity, one row per taxon. + # theta inverts glm_gp's negative-binomial overdispersion, not the + # quasi-likelihood scale the F-test divides by. One matrix column + # per quantity, one row per taxon. direct = rcopy(R""" local({ g <- factor($g, levels = c("A", "B")) @@ -304,7 +306,7 @@ end td <- glmGamPoi::test_de(gp, contrast = "gB", max_lfc = Inf, compute_lfc_se = TRUE) cbind(td$pval, td$f_statistic, td$lfc * log(2), td$lfc_se * log(2), - 1 / gp$overdispersion_shrinkage_list$ql_disp_shrunken) + 1 / gp$overdispersions) }) """) for (j, t) in enumerate(effect_taxa) @@ -320,43 +322,86 @@ end end end - @testset "glmGamPoi unusable dispersions are isolated per taxon" begin + @testset "glmGamPoi unusable overdispersions are isolated per taxon" begin if !HAVE_GLMGAMPOI # Issue #37's contract: these comparisons fail, not skip, in CI. @test !_GLMGAMPOI_IN_CI @test_skip false else + # test_de refits the reduced model from the shrinkage list's + # dispersion trend and divides by ql_disp_shrunken, so altering + # `overdispersions` after the fit changes theta, the boundary check + # and lfc_se (predict() weights by it), never the F statistic or + # its p-value. MetaManifold.RRuntime.with_r_lock() do reval(raw""" test_original_gp <- get("glm_gp", asNamespace("glmGamPoi")) - test_ql <- 0 + test_od <- 0 assignInNamespace("glm_gp", function(...) { fit <- test_original_gp(...) - fit$overdispersion_shrinkage_list$ql_disp_shrunken[1] <- test_ql + fit$overdispersions[1] <- test_od fit }, ns="glmGamPoi") """) try off = log.(tss_factors(effect, effect_samples)) - for (ql, status) in ((0.0, "boundary"), (-1.0, "failed"), (NaN, "failed")) - RCall.globalEnv[:test_ql] = ql + # A zero overdispersion is the Poisson boundary: marked, with + # the statistics kept and no finite theta, as under glm_nb. + # A negative or non-finite one fails its taxon alone. + for (od, status, tested) in ((0.0, "boundary", true), (-1.0, "failed", false), + (NaN, "failed", false)) + RCall.globalEnv[:test_od] = od rows = fit_nb(effect, effect_groups, off; levels=("A", "B"), dispersion="glmGamPoi") @test rows.status[1] == status - @test isnothing(rows.pvalue[1]) - @test occursin("dispersion", rows.note[1]) + @test isnothing(rows.pvalue[1]) == !tested + @test isnothing(rows.theta[1]) + @test occursin("overdispersion", rows.note[1]) @test all(rows.status[2:end] .== "ok") end finally reval(raw""" assignInNamespace("glm_gp", test_original_gp, ns="glmGamPoi") - rm(test_original_gp, test_ql) + rm(test_original_gp, test_od) """) end end end end + @testset "Condition labels never change the fits or reach R as code" begin + if !HAVE_GLMGAMPOI + # Issue #37's contract: these comparisons fail, not skip, in CI. + @test !_GLMGAMPOI_IN_CI + @test_skip false + else + # Condition labels join group, run and sub-group names with '/', and + # those names may carry '-' and '.'. test_de parses and evaluates a + # text contrast as R code, so the contrast goes to it as numeric + # weights. These labels, and one that reads as R code, must give the + # A/B results under both dispersions, and that code must never run. + # The coded label opens with '[' so that the text the old path built, + # "da_g[...]", indexed the live factor da_g and so ran the code. + off = log.(tss_factors(effect, effect_samples)) + evaluated = "mm_contrast_label_evaluated" + labels = (("run-A", "run-B"), ("grpA/run_1", "grpB/run_2"), ("2024-01", "2024-02"), + ("A", "[assign(\"$evaluated\", TRUE, envir = globalenv())]")) + for dispersion in ("glm_nb", "glmGamPoi") + base = fit_nb(effect, effect_groups, off; levels=("A", "B"), dispersion) + @test all(base.status .== "ok") + for (ref, con) in labels + groups = [g == "A" ? ref : con for g in effect_groups] + rows = fit_nb(effect, groups, off; levels=(ref, con), dispersion) + @test rows.status == base.status + for col in (:estimate, :se, :statistic, :pvalue, :theta) + @test all(isapprox.(rows[!, col], base[!, col]; rtol=1e-10)) + end + end + end + @test !rcopy(R"exists($evaluated, envir = globalenv())") + end + end + @testset "Degenerate and filtered taxa are reported, not dropped" begin if !HAVE_MASS @test_skip false