From ad7f5f2c472b5c1447533ca178c7ce484108efad Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:13:46 +0000 Subject: [PATCH 1/2] bench(differential): null-calibration bench + small-sample caveat (refs #193) Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- README.md | 9 +- bench/null_calibration/.gitignore | 4 + bench/null_calibration/README.md | 124 ++++++ bench/null_calibration/REPORT.md | 112 +++++ bench/null_calibration/UPSTREAM_DRAFT.md | 51 +++ bench/null_calibration/calibration.R | 395 ++++++++++++++++++ .../results/calibration_summary.csv | 33 ++ .../results/power_summary.csv | 5 + .../results/provenance_summary.csv | 17 + bench/null_calibration/summarise.R | 131 ++++++ bench/null_calibration/webr_run.mjs | 104 +++++ config/defaults/pipeline.yml | 8 +- frontend/src/api/types.ts | 12 + frontend/src/components/DifferentialPanel.tsx | 9 + frontend/src/components/PipelineStages.tsx | 2 +- .../src/components/differentialTable.test.ts | 31 ++ frontend/src/components/differentialTable.ts | 8 +- src/analysis/differential.jl | 71 +++- test/unit/test_differential.jl | 32 ++ 19 files changed, 1153 insertions(+), 5 deletions(-) create mode 100644 bench/null_calibration/.gitignore create mode 100644 bench/null_calibration/README.md create mode 100644 bench/null_calibration/REPORT.md create mode 100644 bench/null_calibration/UPSTREAM_DRAFT.md create mode 100644 bench/null_calibration/calibration.R create mode 100644 bench/null_calibration/results/calibration_summary.csv create mode 100644 bench/null_calibration/results/power_summary.csv create mode 100644 bench/null_calibration/results/provenance_summary.csv create mode 100644 bench/null_calibration/summarise.R create mode 100644 bench/null_calibration/webr_run.mjs diff --git a/README.md b/README.md index 998fcc15..6bf9fc2b 100644 --- a/README.md +++ b/README.md @@ -260,7 +260,14 @@ not exclude samples from this analysis. - `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 per-taxon maximum-likelihood + theta (`ml`, default) or a fixed theta from `glmGamPoi::glm_gp`'s cross-taxon + dispersion trend (`glmGamPoi`). Both options Wald-test the group coefficient + with the dispersion treated as known, so every `ml` result — and a + `glmGamPoi` result whose smaller group has at most 12 samples — carries a + calibration caveat in its diagnostics, shown with the result and written + into its CSV export. The caveat changes no p-value; the null-simulation + evidence and full replication protocol live in `bench/null_calibration`. - `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/bench/null_calibration/.gitignore b/bench/null_calibration/.gitignore new file mode 100644 index 00000000..088565ad --- /dev/null +++ b/bench/null_calibration/.gitignore @@ -0,0 +1,4 @@ +# Raw per-taxon p-values are large, append-only working state. The seeds in +# calibration.R regenerate them dataset-for-dataset; committed evidence is +# the summarised CSVs and REPORT.md under results/. +results/raw/ diff --git a/bench/null_calibration/README.md b/bench/null_calibration/README.md new file mode 100644 index 00000000..f9598d2a --- /dev/null +++ b/bench/null_calibration/README.md @@ -0,0 +1,124 @@ +# Null-calibration bench for `nb_glm` (issue #193) + +Issue #193 measured that the default `nb_glm` path — `MASS::glm.nb` per taxon +with a fixed-scale Wald test, `analysis.differential.dispersion: ml` — rejects +true nulls too often in small studies (nominal 0.05 rejecting at 0.132 with 4 +samples per group and 0.069 with 12). The mechanism is textbook: the Wald test +treats the per-taxon maximum-likelihood theta as known, ignoring its sampling +uncertainty, which is large when groups are small. + +This bench replicates that finding at the scale #193 asked for, measures the +candidate remedies on identical datasets, and records what was actually run. +It is a bench, not a unit test: slow, stochastic, never run by CI. + +## What is compared + +Every simulated dataset is analysed by every enabled test (paired design — no +test gets an easier dataset): + +| test | what it is | status | +|---|---|---| +| `wald_ml` | the production default, through the real `_FIT_R` block extracted from `src/analysis/differential.jl` (`da_use_gp = FALSE`) | production | +| `wald_trend` | the production `dispersion: glmGamPoi` option on this branch, through the same block (`da_use_gp = TRUE`): dispersion trend plugged into a fixed-theta Wald fit. **Not** a quasi-likelihood test | production option | +| `lrt` | `MASS::glm.nb` likelihood-ratio test, full vs intercept-only, each fit estimating its own theta; `2 × ΔlogLik` against χ²(1) | bench-only candidate | +| `glmgampoi_ql` | glmGamPoi's own quasi-likelihood F-test (`glm_gp` + `test_de`, numeric contrast weights) | reference (the calibrated path in #193) | +| `edger_ql` | edgeR's quasi-likelihood F-test, fed the same offsets | outside reference only, not a dependency | + +A test whose package is unavailable is skipped loudly (stdout + provenance +file), never silently replaced. + +## Protocol (mapped to #193's checklist) + +- **Group sizes** `n ∈ {3, 4, 6, 8, 12, 20}` per group, equal groups. +- **Summaries**: rejection rates at α = 0.01 / 0.05 / 0.10 with Wilson 95% + intervals, Kolmogorov–Smirnov D and p of the null p-values against + Uniform(0, 1), and per-test fit-failure shares (`summarise.R`). +- **Data shapes** (`--shapes`): + - `original` — the #193 generator: taxon means exp(U(log 2, log 500)), + dispersion `0.05 + 2/√μ`, library sizes lognormal(0, 0.3); + - `constant` — same means, dispersion 0.5 for every taxon; + - `sparse` — amplicon-like: means exp(N(log 3, 1.6)), so most taxa are rare + and a few dominate; many zeros; + - `unequal_lib` — `original` with group B library sizes ×2.5. +- **Offsets** (`--offsets`): `tss` and `rle`, implemented to mirror + `tss_factors` / `rle_factors` in `differential.jl`. An RLE refusal (no taxon + present in every sample) is counted, never replaced by TSS. +- **Power** (`--truth=power --log2fc=L`): 25% of taxa truly change by ±L (log2), + the rest are null. Reports raw and BH power on the changed taxa, the BH + false discovery rate at 0.05 (averaged per dataset, BH applied within each + dataset over its tested family, as production does), and null-taxon + calibration within the same cells. +- **Taxa filter**: the #193 script's filter is kept for fidelity — a taxon + needs reads in ≥2 samples and non-constant counts. Taxon parameters are + redrawn for every dataset, as in #193. Datasets with <5 testable taxa are + skipped and counted. +- **Seeds**: `seed0 + cell_index × 1000003 + sim`, one per dataset, so cells + can run in any order or in parallel processes and reproduce exactly. + Default `seed0 = 20261020` (the #193 seed). + +### Reading a cell + +A test is **flagged anti-conservative** in a cell when the lower bound of the +Wilson 95% interval for its rejection rate at α = 0.05 exceeds 0.05, and +**flagged conservative** when the upper bound is below 0.05. The KS statistic +is supporting evidence, not the criterion — with ~10⁴ p-values it detects +departures too small to matter at any single α. + +## Running it + +Native R with the pinned library (the supported route; `Rscript -e +'renv::restore()'` first): + +```sh +# smoke: 10 sims, one cell, a few seconds +Rscript bench/null_calibration/calibration.R --mode=quick + +# one direct-replication cell +Rscript bench/null_calibration/calibration.R \ + --shapes=original --offsets=tss --n-per=4 --sims=200 --taxa=40 \ + --tests=wald_ml,lrt + +# the full #193 grid — 1,000 sims per cell; expect hours, use several +# processes by splitting --n-per / --shapes across terminals +Rscript bench/null_calibration/calibration.R --mode=full + +# aggregate whatever raw cells exist +Rscript bench/null_calibration/summarise.R +``` + +Outputs append to `results/raw/pvalues_.csv` (gitignored), one row per +taxon × test, so an interrupted run keeps finished simulations; each cell also +writes `provenance_.txt` with arguments, package versions and +`sessionInfo()`. `summarise.R` writes `results/calibration_summary.csv` (and +`power_summary.csv` when power cells exist). + +### No native R? The WebR fallback + +`webr_run.mjs` runs the same `calibration.R` / `summarise.R` under WebR (R in +WebAssembly) for environments with Node but no R. `glm.nb` is pure R (MASS +7.3-65's `negbin.R` makes no `.C`/`.Call`/`.Fortran` call), so sourcing that +one file provides genuine MASS code; `calibration.R` detects this mode, +strips the block's `MASS::` prefixes as a single mechanical substitution, and +records the fact in every provenance file. glmGamPoi and edgeR are compiled +Bioconductor packages and cannot run this way — their tests skip loudly. +Setup and usage are in the header of `webr_run.mjs`. Prefer native Rscript +whenever it exists. + +## What has actually been run + +`REPORT.md` records each committed run: cells, sizes, environment and results, +with the corresponding summary CSVs in `results/`. Anything not in `REPORT.md` +has not been run — in particular, a committed bench is not a result. + +## Limitations (read before quoting rates) + +- Taxa are simulated independently. Real amplicon tables are compositional and + correlated; independence makes BH behaviour here cleaner than in practice. +- Groups are equal-sized. Unequal groups change both theta estimation and + offset behaviour; #193's grid did not cross that factor either. +- The ≥2-sample presence filter mirrors #193, not production's default + `min_prevalence = 0` (under which untestable taxa instead surface as failed + rows outside the BH family). Failure shares are reported so the tested + family is never a hidden denominator. +- `original`-shape library sizes are mild (lognormal σ = 0.3). Harsher depth + variation is covered only through `unequal_lib`'s group shift. diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md new file mode 100644 index 00000000..68c1f34a --- /dev/null +++ b/bench/null_calibration/REPORT.md @@ -0,0 +1,112 @@ +# Null-calibration bench report (issue #193) + +This report records the one reduced run that has actually been executed. +Anything not listed here has not been run. The committed summary CSVs in +`results/` (`calibration_summary.csv`, `power_summary.csv`) are the source of +truth for every number below; per-cell run records are summarised in +`results/provenance_summary.csv` (the raw p-value files under `results/raw/` +are gitignored and reproducible from the recorded seeds and arguments). + +## Run identity + +- **Date:** 2026-10-22 (UTC, per-cell provenance records). +- **Executor:** WebR 0.6.0 (R 4.6.0, wasm32) through `webr_run.mjs`, because the + sandbox used for this run had Node but no native R. `MASS::glm.nb` was + provided by sourcing MASS 7.3-65's pure-R `negbin.R` (the exact MASS version + pinned in `renv.lock`); in that mode `calibration.R` strips the `MASS::` + prefixes from the extracted production block as a single mechanical + substitution and records the fact in every provenance file. Native + `Rscript` runs use the extracted block byte-identically. The sandbox R + version (4.6.0) differs from the lockfile's R 4.5.0. +- **Production path under test:** the real `_FIT_R` block extracted from + `src/analysis/differential.jl` (`wald_ml`, `dispersion: ml`), paired on + identical datasets with a bench-only likelihood-ratio test (`lrt`: full vs + intercept-only `glm.nb`, each fit estimating its own theta, `2 × ΔlogLik` + against χ²(1)). +- **Not run in this sandbox:** `wald_trend`, `glmgampoi_ql`, `edger_ql` + (glmGamPoi/edgeR are compiled Bioconductor packages and were absent — the + provenance records say so). Their code paths in `calibration.R` are + unexecuted here and need a native R/Bioconductor run. +- **Cells:** direct replication (original/TSS null) at 3/4/6/8/12/20 samples + per group, 200 sims × 40 taxa; shape/offset variations (constant, sparse, + unequal library sizes, original/RLE) at 4 and 12 per group, 150 × 40; power + cells (original/TSS, 25% of taxa at |log2FC| = 1) at 4 and 12 per group, + 150 × 40. This is a *reduced* run: the full 1,000-sim grid in `README.md` + has not been run. + +## Results + +A cell is flagged anti-conservative when the lower bound of the Wilson 95% +interval for its α = 0.05 rejection rate exceeds 0.05. **All 32 summarised +α = 0.05 test-cells were flagged**, including 20 samples per group; the +closest cell was LRT at 20/group on original/TSS (lower bound 0.0534). + +### Direct replication (original shape, TSS offset, null) + +Rejection rates at α = 0.01 / 0.05 / 0.10. Issue #193 reported Wald at +4/group: 0.060 / 0.132 / 0.198 and 12/group: 0.019 / 0.069 / 0.124 — the +replication lands essentially on those values. + +| n/group | sims | `wald_ml` (production) | `lrt` (bench-only) | Wald fit failures | LRT fit failures | +|---:|---:|---|---|---:|---:| +| 3 | 200 | 0.092 / 0.179 / 0.251 | 0.046 / 0.131 / 0.203 | 6.41% | 8.43% | +| 4 | 200 | 0.063 / 0.138 / 0.204 | 0.030 / 0.103 / 0.169 | 2.68% | 3.61% | +| 6 | 200 | 0.045 / 0.111 / 0.171 | 0.023 / 0.087 / 0.149 | 0.79% | 1.20% | +| 8 | 200 | 0.029 / 0.089 / 0.145 | 0.017 / 0.070 / 0.129 | 0.20% | 0.33% | +| 12 | 200 | 0.025 / 0.076 / 0.130 | 0.017 / 0.065 / 0.119 | 0.01% | 0.04% | +| 20 | 200 | 0.019 / 0.065 / 0.116 | 0.014 / 0.058 / 0.110 | 0.00% | 0.00% | + +### Other shapes and the RLE offset (α = 0.05, Wilson 95% interval) + +| shape / offset | n/group | sims | `wald_ml` | `lrt` | Wald failures | LRT failures | +|---|---:|---:|---|---|---:|---:| +| `constant` / `tss` | 4 | 150 | 0.135 (0.126–0.144) | 0.104 (0.096–0.112) | 4.03% | 5.35% | +| `constant` / `tss` | 12 | 150 | 0.073 (0.066–0.080) | 0.065 (0.059–0.072) | 0.18% | 0.23% | +| `sparse` / `tss` | 4 | 150 | 0.107 (0.098–0.116) | 0.082 (0.074–0.090) | 15.12% | 18.36% | +| `sparse` / `tss` | 12 | 150 | 0.068 (0.062–0.075) | 0.060 (0.054–0.067) | 4.52% | 5.07% | +| `unequal_lib` / `tss` | 4 | 150 | 0.149 (0.140–0.158) | 0.114 (0.106–0.122) | 1.68% | 3.17% | +| `unequal_lib` / `tss` | 12 | 150 | 0.074 (0.068–0.081) | 0.062 (0.056–0.068) | 0.00% | 0.08% | +| `original` / `rle` | 4 | 150 | 0.138 (0.130–0.147) | 0.100 (0.093–0.108) | 2.57% | 3.37% | +| `original` / `rle` | 12 | 150 | 0.080 (0.074–0.087) | 0.068 (0.061–0.074) | 0.02% | 0.02% | + +The inflation is not specific to the original mean-varying dispersion shape, +and the RLE offset does not rescue calibration. + +### Power and BH behaviour (25% of taxa at |log2FC| = 1) + +BH is applied within each simulated dataset over its tested family, as +production does. “Null rate” is the α = 0.05 rejection rate among the truly +null taxa in the same cells. + +| n/group | test | sims | null rate @0.05 | raw power | BH power | BH FDR | mean discoveries | +|---:|---|---:|---:|---:|---:|---:|---:| +| 4 | `lrt` | 150 | 0.114 | 0.414 | 0.147 | 0.184 | 1.89 | +| 4 | `wald_ml` | 150 | 0.150 | 0.480 | 0.290 | 0.379 | 4.70 | +| 12 | `lrt` | 150 | 0.083 | 0.668 | 0.422 | 0.101 | 4.82 | +| 12 | `wald_ml` | 150 | 0.096 | 0.688 | 0.509 | 0.129 | 5.99 | + +## Interpretation + +- #193 replicates. The production fixed-scale Wald test on per-taxon ML theta + rejects true nulls well above nominal at every measured group size, + including 20/group (0.065 at α = 0.05). +- The LRT is consistently closer to nominal but remains anti-conservative in + every cell, and fails more often than Wald at tiny n (it needs two fits). + It is a mitigation, not a calibration fix. +- Wald's extra BH “power” is bought with uncontrolled FDR (0.38 at 4/group + against a nominal 0.05). Nothing measured in this sandbox “stays + calibrated”; the calibrated comparator from #193 (glmGamPoi’s + quasi-likelihood F-test) still needs the native run. +- Because even 20/group remained flagged, the production calibration caveat + (see `_calibration_caveat` in `src/analysis/differential.jl`) uses “caution” + severity through the largest measured size and only drops to “info” beyond + the measured range. That caveat is fork-side only; any change to upstream + defaults remains an upstream decision (see `UPSTREAM_DRAFT.md`). + +## Checklist mapping (#193) + +Done in this reduced run: direct replication at larger sim counts and more +group sizes; additional shapes; RLE offset at 4/12 per group; an additional +test (LRT) on identical datasets; power and BH-FDR at |log2FC| = 1. +Not done: the full 1,000-sim grid; `wald_trend`; glmGamPoi QL F; edgeR QL; +native Rscript/Bioconductor execution; power at other effect/group sizes. diff --git a/bench/null_calibration/UPSTREAM_DRAFT.md b/bench/null_calibration/UPSTREAM_DRAFT.md new file mode 100644 index 00000000..5de5b181 --- /dev/null +++ b/bench/null_calibration/UPSTREAM_DRAFT.md @@ -0,0 +1,51 @@ +# Draft upstream report: null calibration of the `glm.nb` Wald test + +*Status: draft prepared downstream for review before anything is filed +upstream. All numbers are from `REPORT.md` in this directory, which defines +exactly what ran, where, and what did not run.* + +## Summary + +In null simulations that drive the production negative-binomial fitting code +(`MASS::glm.nb` per taxon, fixed-scale Wald test, per-taxon ML theta treated +as known), true nulls are rejected more often than nominal at every measured +group size. At α = 0.05 the rejection rate was 0.138 with 4 samples per group +and 0.076 with 12 (original/TSS shape, 200 sims × 40 taxa), and still 0.065 +at 20 per group. A likelihood-ratio test on the same fits/datasets was closer +to nominal (0.103 / 0.065 / 0.058) but still anti-conservative in every cell +(Wilson 95% lower bound > 0.05 in all 32 summarised test-cells). Under a +power design (25% of taxa at |log2FC| = 1), Benjamini–Hochberg at 0.05 +delivered FDR 0.38 (Wald, 4/group) and 0.13 (12/group). + +The mechanism is the standard one: the Wald test plugs in the per-taxon ML +theta and ignores its sampling uncertainty. The bench that produces these +numbers extracts the production fitting block directly from +`src/analysis/differential.jl`, so the `wald_ml` column is the shipped code +path, not a reimplementation. + +## What was and was not measured + +Measured (reduced WebR run, genuine MASS 7.3-65 `glm.nb`; see `REPORT.md` for +provenance limits): Wald vs LRT, paired on identical datasets; group sizes +3–20; constant/sparse/unequal-library shapes; RLE offsets; one power design. + +Not yet measured: the glmGamPoi quasi-likelihood F-test as an in-bench +comparator on this branch (the calibrated path reported in the downstream +issue), the dispersion-trend Wald option, edgeR, and the full 1,000-sim grid. +The bench (`bench/null_calibration/`, native `Rscript` route) is ready to run +those where R + Bioconductor are available. + +## Options for the upstream maintainer's consideration + +1. **Result-notes caveat for small samples** (adopted downstream): record a + calibration caveat with `ml`-dispersion results, scaled by group size. + Changes no statistics or defaults. +2. **LRT instead of Wald for the `glm.nb` path:** measurably closer to + nominal in this bench, but still anti-conservative, slower, and fails + more often at tiny n. The bench quantifies the trade-off; it is not a + calibration fix. +3. **Prefer a quasi-likelihood path below a group-size threshold:** the + downstream measurements that motivated this bench found the QL F-test at + nominal rates where Wald was inflated. Any such default change is + upstream's call; this draft only supplies the replication bench and + numbers to decide with. diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R new file mode 100644 index 00000000..a424dc4c --- /dev/null +++ b/bench/null_calibration/calibration.R @@ -0,0 +1,395 @@ +#!/usr/bin/env Rscript +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Null-calibration and power bench for the nb_glm tests in +# src/analysis/differential.jl (fork issue #193). This is a bench, not a unit +# test: it is slow and stochastic, and it never runs in CI by default. +# +# What it measures, on the SAME simulated datasets (paired comparison): +# wald_ml - the production default: the real `_FIT_R` block extracted +# from differential.jl, da_use_gp = FALSE (MASS::glm.nb per +# taxon, Wald test at fixed scale 1). +# wald_trend - the production glmGamPoi option ON THIS BRANCH: the same +# `_FIT_R` block with da_use_gp = TRUE (glm_gp dispersion +# trend plugged into a fixed-theta Wald fit). This is NOT +# glmGamPoi's quasi-likelihood test. Needs glmGamPoi. +# lrt - MASS::glm.nb likelihood-ratio test, full (~ group) against +# intercept-only, each fit estimating its own theta; +# 2 * (logLik(full) - logLik(null)) against chisq(1). +# Bench-only candidate fix; production does not use it. +# glmgampoi_ql - glmGamPoi's own quasi-likelihood F-test (glm_gp + test_de), +# as on the upstream PR 62 branch. Needs glmGamPoi. +# edger_ql - edgeR's quasi-likelihood F-test as an outside reference +# only. Needs edgeR. Not a dependency of the app. +# +# Tests whose package is unavailable are skipped LOUDLY (a line on stdout and +# in the provenance file), never silently replaced by another test. +# +# Usage (from the repo root, with the pinned R library restored): +# Rscript bench/null_calibration/calibration.R --mode=quick +# Rscript bench/null_calibration/calibration.R --mode=full # hours; see README +# Rscript bench/null_calibration/calibration.R --truth=power --log2fc=1 \ +# --n-per=4,12 --sims=200 +# Any --key=value below overrides the mode defaults. Under WebR (no Rscript +# argv), settings may instead come from the NB_CALIB_ARGS environment +# variable as one space-separated string; see webr_run.mjs. +# +# Output: one appended CSV per cell under --out (default +# bench/null_calibration/results/raw), so an interrupted long run keeps the +# simulations it finished. summarise.R aggregates them. + +## ---------------------------------------------------------------- settings + +parse_args <- function(argv) { + out <- list() + for (a in argv) { + if (!grepl("^--[a-z0-9-]+=", a)) stop("unrecognised argument '", a, "' (want --key=value)") + kv <- strsplit(sub("^--", "", a), "=", fixed = TRUE)[[1]] + out[[kv[1]]] <- paste(kv[-1], collapse = "=") + } + out +} + +argv <- commandArgs(trailingOnly = TRUE) +if (length(argv) == 0 && nzchar(Sys.getenv("NB_CALIB_ARGS"))) { + argv <- strsplit(trimws(Sys.getenv("NB_CALIB_ARGS")), "\\s+")[[1]] +} +args <- parse_args(argv) + +mode <- if (!is.null(args$mode)) args$mode else "quick" +defaults <- switch(mode, + quick = list(`n-per` = "4", sims = "10", taxa = "20", shapes = "original", + offsets = "tss", tests = "wald_ml,lrt", truth = "null", + log2fc = "1", seed = "20261020"), + full = list(`n-per` = "3,4,6,8,12,20", sims = "1000", taxa = "40", + shapes = "original,constant,sparse,unequal_lib", + offsets = "tss,rle", + tests = "wald_ml,wald_trend,lrt,glmgampoi_ql,edger_ql", + truth = "null", log2fc = "1", seed = "20261020"), + stop("unknown --mode '", mode, "' (want quick or full, plus overrides)") +) +for (k in names(defaults)) if (is.null(args[[k]])) args[[k]] <- defaults[[k]] +if (is.null(args$out)) args$out <- file.path("bench", "null_calibration", "results", "raw") +if (is.null(args[["differential-jl"]])) { + args[["differential-jl"]] <- file.path("src", "analysis", "differential.jl") +} + +csv_ints <- function(s) as.integer(strsplit(s, ",", fixed = TRUE)[[1]]) +csv_strs <- function(s) trimws(strsplit(s, ",", fixed = TRUE)[[1]]) + +n_pers <- csv_ints(args[["n-per"]]) +n_sims <- as.integer(args$sims) +n_taxa <- as.integer(args$taxa) +shapes <- csv_strs(args$shapes) +offsets <- csv_strs(args$offsets) +want <- csv_strs(args$tests) +truth <- args$truth +lfc <- as.numeric(args$log2fc) +seed0 <- as.integer(args$seed) +out_dir <- args$out +stopifnot(truth %in% c("null", "power")) +dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + +## ------------------------------------------------------------ MASS access +# Normally MASS is installed and the extracted production block runs +# byte-identical. In a minimal R with no MASS but with MASS's pure-R +# negbin.R already sourced (glm.nb is pure R; MASS 7.3-65 negbin.R makes no +# .C/.Call/.Fortran call), the block's `MASS::` prefixes are stripped as a +# single mechanical substitution and that fact is recorded in the +# provenance. The functions executed are still MASS's own source. + +have_mass_ns <- requireNamespace("MASS", quietly = TRUE) +mass_mode <- if (have_mass_ns) { + paste0("installed MASS ", as.character(packageVersion("MASS"))) +} else if (exists("glm.nb", envir = globalenv())) { + "sourced negbin.R (MASS:: stripped from the extracted block; see header)" +} else { + stop("MASS is not available and glm.nb is not sourced: restore renv.lock first") +} +glm_nb <- if (have_mass_ns) MASS::glm.nb else get("glm.nb", envir = globalenv()) + +have_gp <- requireNamespace("glmGamPoi", quietly = TRUE) +have_edger <- requireNamespace("edgeR", quietly = TRUE) + +## ------------------------------------------------------- production block + +extract_fit_r <- function(path) { + lines <- readLines(path, warn = FALSE) + start <- grep('const _FIT_R = raw"""', lines, fixed = TRUE) + if (length(start) != 1) stop("could not find exactly one _FIT_R block in ", path) + rest <- lines[(start + 1):length(lines)] + end <- which(rest == '"""')[1] + if (is.na(end)) stop("_FIT_R block in ", path, " is not terminated") + block <- paste(rest[seq_len(end - 1)], collapse = "\n") + if (!grepl("da_result", block, fixed = TRUE) || !grepl("glm.nb", block, fixed = TRUE)) { + stop("extracted _FIT_R does not look like the production fit loop") + } + block +} + +fit_code <- extract_fit_r(args[["differential-jl"]]) +if (!have_mass_ns) fit_code <- gsub("MASS::", "", fit_code, fixed = TRUE) + +# Every da_* global fit_nb sets or the block creates on the current branch +# (mirrors the cleanup list in differential.jl's fit_nb). +DA_VARS <- c("da_counts", "da_group", "da_levels", "da_offset", "da_theta_upper", + "da_theta_lower", "da_use_gp", "da_g", "da_term", "da_n", "da_result", + "da_phi", "da_gp_failed", "da_gp", "da_j", "da_y", "da_warns", + "da_fit", "da_co", "da_row", "da_theta", "da_absent") + +run_block <- function(counts, offset, use_gp, n_per) { + g <- globalenv() + assign("da_counts", counts, g) + assign("da_group", rep(c("A", "B"), each = n_per), g) + assign("da_levels", c("A", "B"), g) + assign("da_offset", offset, g) + assign("da_theta_upper", 1e7, g) + assign("da_theta_lower", 1e-8, g) + assign("da_use_gp", use_gp, g) + on.exit(rm(list = intersect(ls(g), DA_VARS), envir = g)) + suppressWarnings(eval(parse(text = fit_code), envir = g)) + res <- get("da_result", envir = g) + res +} + +## ---------------------------------------------------------------- offsets +# Mirror src/analysis/differential.jl exactly: factors centred to geometric +# mean 1, offset = log(factor). RLE refuses (returns NULL) when no taxon is +# present in every sample, as rle_factors does; the refusal is recorded, +# never replaced by TSS. + +geomean <- function(x) exp(mean(log(x))) + +make_offset <- function(counts, method) { + if (method == "tss") { + lib <- rowSums(counts) + if (any(lib <= 0)) return(NULL) + return(log(lib / geomean(lib))) + } + if (method == "rle") { + shared <- which(colSums(counts > 0) == nrow(counts)) + if (length(shared) == 0) return(NULL) + sub <- counts[, shared, drop = FALSE] + gm <- apply(sub, 2, geomean) + raw <- apply(sub, 1, function(row) median(row / gm)) + return(log(raw / geomean(raw))) + } + stop("unknown offset method '", method, "'") +} + +## ------------------------------------------------------------- simulation +# One dataset: samples x taxa integer counts plus the true log2 fold change +# per taxon (all zero under the null). Taxon parameters are redrawn per +# dataset, as in the #193 script. The data-generating library sizes are +# realisations, not observed library sizes: offsets are computed from the +# simulated counts afterwards, as production computes them from data. + +simulate <- function(shape, n_per, n_taxa, truth, lfc) { + n <- 2 * n_per + mu <- switch(shape, + original = exp(runif(n_taxa, log(2), log(500))), + constant = exp(runif(n_taxa, log(2), log(500))), + sparse = exp(rnorm(n_taxa, mean = log(3), sd = 1.6)), + unequal_lib = exp(runif(n_taxa, log(2), log(500))), + stop("unknown shape '", shape, "'") + ) + disp <- if (shape == "constant") rep(0.5, n_taxa) else 0.05 + 2 / sqrt(mu) + lib <- exp(rnorm(n, 0, 0.3)) + if (shape == "unequal_lib") lib[(n_per + 1):n] <- lib[(n_per + 1):n] * 2.5 + true_lfc <- rep(0, n_taxa) + if (truth == "power") { + n_de <- max(1L, round(0.25 * n_taxa)) + de <- seq_len(n_de) + true_lfc[de] <- rep(c(lfc, -lfc), length.out = n_de) + } + y <- vapply(seq_len(n_taxa), function(j) { + m <- mu[j] * lib + if (true_lfc[j] != 0) m[(n_per + 1):n] <- m[(n_per + 1):n] * 2^true_lfc[j] + as.integer(rnbinom(n, mu = m, size = 1 / disp[j])) + }, integer(n)) + storage.mode(y) <- "integer" + # The #193 script's filter, kept so the replication is faithful: a taxon + # needs reads in at least 2 samples and non-constant counts to be tested. + keep <- colSums(y > 0) >= 2 & apply(y, 2, function(v) length(unique(v)) > 1) + list(counts = y[, keep, drop = FALSE], true_lfc = true_lfc[keep], + zero_frac = mean(y == 0)) +} + +## ------------------------------------------------------------ bench tests + +# Production's per-taxon screens, applied to the bench-only tests so every +# test faces the same family: constant counts, and complete separation (no +# reads in any sample of one group). Returns "" when fittable. +screen <- function(yv, g) { + if (length(unique(yv)) < 2L) return("constant counts") + if (all(yv[g == "A"] == 0) || all(yv[g == "B"] == 0)) return("complete separation") + "" +} + +fit_lrt <- function(counts, offset, n_per) { + g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) + out <- data.frame(status = "failed", pvalue = NA_real_) + out <- out[rep(1, ncol(counts)), , drop = FALSE] + rownames(out) <- NULL + one <- function(formula, yv) { + warns <- character(0) + fit <- tryCatch( + withCallingHandlers(glm_nb(formula, control = glm.control(maxit = 100)), + warning = function(w) { warns <<- c(warns, conditionMessage(w)); invokeRestart("muffleWarning") }), + error = function(e) e) + list(fit = fit, warns = warns) + } + bad_fit <- function(r) { + inherits(r$fit, "error") || !isTRUE(r$fit[["converged"]]) || + !is.null(r$fit[["th.warn"]]) || + any(grepl("iteration limit|alternation limit|did not converge|NaNs produced", r$warns)) + } + off <- offset + for (j in seq_len(ncol(counts))) { + yv <- counts[, j] + if (nzchar(screen(yv, g))) next + # glm.nb looks its variables up in the formula's environment; yv, g and + # off all live in this frame (one()'s enclosing environment), which is + # where these formulas are created. + full <- one(yv ~ g + offset(off), yv) + if (bad_fit(full)) next + null <- one(yv ~ 1 + offset(off), yv) + if (bad_fit(null)) next + lr <- 2 * (as.numeric(logLik(full$fit)) - as.numeric(logLik(null$fit))) + if (!is.finite(lr) || lr < 0) next + p <- pchisq(lr, df = 1, lower.tail = FALSE) + if (!is.finite(p)) next + out$status[j] <- "ok" + out$pvalue[j] <- p + } + out +} + +fit_glmgampoi_ql <- function(counts, offset, n_per) { + g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) + out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_) + screens <- vapply(seq_len(ncol(counts)), + function(j) nzchar(screen(counts[, j], g)), logical(1)) + res <- tryCatch(suppressWarnings({ + gp <- glmGamPoi::glm_gp(t(counts), design = model.matrix(~ g), + size_factors = exp(offset), + overdispersion_shrinkage = TRUE) + # Numeric contrast weights on the design columns ((Intercept), gB): + # test_de parses a TEXT contrast as R code, so bench and production + # alike must never hand it a label. + glmGamPoi::test_de(gp, contrast = c(0, 1), sort_by = "none") + }), error = function(e) e) + if (inherits(res, "error")) return(out) + pcol <- intersect(c("pval", "pvalue", "p_value"), names(res))[1] + if (is.na(pcol) || length(res[[pcol]]) != ncol(counts)) return(out) + p <- as.numeric(res[[pcol]]) + ok <- !screens & is.finite(p) & p >= 0 & p <= 1 + out$status[ok] <- "ok" + out$pvalue[ok] <- p[ok] + out +} + +fit_edger_ql <- function(counts, offset, n_per) { + g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) + out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_) + screens <- vapply(seq_len(ncol(counts)), + function(j) nzchar(screen(counts[, j], g)), logical(1)) + res <- tryCatch(suppressWarnings({ + y <- edgeR::DGEList(counts = t(counts), group = g) + # Feed edgeR the bench offset rather than its own normalisation, so all + # tests see the same library-size handling: effective library size + # proportional to exp(offset). + y$samples$lib.size <- exp(offset) * geomean(rowSums(counts)) + y$samples$norm.factors <- 1 + design <- model.matrix(~ g) + y <- edgeR::estimateDisp(y, design) + fit <- edgeR::glmQLFit(y, design) + qlf <- edgeR::glmQLFTest(fit, coef = 2) + edgeR::topTags(qlf, n = Inf, sort.by = "none")$table + }), error = function(e) e) + if (inherits(res, "error")) return(out) + pcol <- intersect(c("pvalue", "p.value", "pvalue"), names(res))[1] + if (is.na(pcol) || length(res[[pcol]]) != ncol(counts)) return(out) + p <- as.numeric(res[[pcol]]) + ok <- !screens & is.finite(p) & p >= 0 & p <= 1 + out$status[ok] <- "ok" + out$pvalue[ok] <- p[ok] + out +} + +## ------------------------------------------------------------------- main + +enabled <- want +for (t in want) { + if (t %in% c("wald_trend", "glmgampoi_ql") && !have_gp) { + cat("SKIPPED test", t, ": glmGamPoi is not available in this R library\n") + enabled <- setdiff(enabled, t) + } + if (t == "edger_ql" && !have_edger) { + cat("SKIPPED test edger_ql: edgeR is not available in this R library\n") + enabled <- setdiff(enabled, t) + } +} +if (length(enabled) == 0) stop("no requested test can run in this R library") + +cell_index <- 0L +for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { + cell_index <- cell_index + 1L + # The taxa count is part of the cell identity: reruns append to a cell's + # CSV, so two settings must never share a file by accident. + cell <- sprintf("%s_%s_%s_n%d_t%d%s", truth, shape, off_m, n_per, n_taxa, + if (truth == "power") sprintf("_lfc%s", lfc) else "") + f <- file.path(out_dir, paste0("pvalues_", cell, ".csv")) + prov <- file.path(out_dir, paste0("provenance_", cell, ".txt")) + writeLines(c( + paste("cell:", cell), paste("date:", format(Sys.time(), tz = "UTC", usetz = TRUE)), + paste("args:", paste(argv, collapse = " ")), + paste("MASS:", mass_mode), + paste("glmGamPoi:", if (have_gp) as.character(packageVersion("glmGamPoi")) else "absent"), + paste("edgeR:", if (have_edger) as.character(packageVersion("edgeR")) else "absent"), + paste("tests run:", paste(enabled, collapse = ",")), + capture.output(sessionInfo()) + ), prov) + cat(sprintf("cell %s: %d sims x %d taxa, tests %s\n", + cell, n_sims, n_taxa, paste(enabled, collapse = "+"))) + header_done <- file.exists(f) + skipped_sims <- 0L + refused_offsets <- 0L + for (s in seq_len(n_sims)) { + # One seed per (cell, sim): cells can run in any order, in parallel + # processes, and reproduce dataset-for-dataset. + set.seed(seed0 + cell_index * 1000003L + s) + sim <- simulate(shape, n_per, n_taxa, truth, lfc) + counts <- sim$counts + if (ncol(counts) < 5) { skipped_sims <- skipped_sims + 1L; next } + offset <- make_offset(counts, off_m) + if (is.null(offset)) { refused_offsets <- refused_offsets + 1L; next } + fits <- list() + if ("wald_ml" %in% enabled) { + r <- run_block(counts, offset, FALSE, n_per) + fits$wald_ml <- data.frame(status = r$status, pvalue = r$pvalue) + } + if ("wald_trend" %in% enabled) { + r <- run_block(counts, offset, TRUE, n_per) + fits$wald_trend <- data.frame(status = r$status, pvalue = r$pvalue) + } + if ("lrt" %in% enabled) fits$lrt <- fit_lrt(counts, offset, n_per) + if ("glmgampoi_ql" %in% enabled) fits$glmgampoi_ql <- fit_glmgampoi_ql(counts, offset, n_per) + if ("edger_ql" %in% enabled) fits$edger_ql <- fit_edger_ql(counts, offset, n_per) + rows <- do.call(rbind, lapply(names(fits), function(t) { + data.frame(truth = truth, shape = shape, offset = off_m, n_per = n_per, + taxa = n_taxa, sim = s, taxon = seq_len(ncol(counts)), + true_log2fc = sim$true_lfc, test = t, + status = fits[[t]]$status, pvalue = fits[[t]]$pvalue, + stringsAsFactors = FALSE) + })) + write.table(rows, f, sep = ",", row.names = FALSE, + col.names = !header_done, append = header_done, quote = FALSE) + header_done <- TRUE + if (s %% 25 == 0) { cat(sprintf(" %s sim %d/%d\n", cell, s, n_sims)); flush.console() } + } + cat(sprintf("cell %s done: %d sims skipped (<5 testable taxa), %d offsets refused\n", + cell, skipped_sims, refused_offsets)) +} +cat("calibration run complete\n") diff --git a/bench/null_calibration/results/calibration_summary.csv b/bench/null_calibration/results/calibration_summary.csv new file mode 100644 index 00000000..a73478c8 --- /dev/null +++ b/bench/null_calibration/results/calibration_summary.csv @@ -0,0 +1,33 @@ +"truth","shape","offset","n_per_group","n_taxa_setting","test","n_sims","n_attempted","n_pvalues","n_failed","n_boundary","fail_rate","rate_01","rate_01_lo","rate_01_hi","rate_05","rate_05_lo","rate_05_hi","rate_10","rate_10_lo","rate_10_hi","ks_D","ks_p","mean_p_null" +"null","constant","tss",4,40,"lrt",150,6000,5679,321,0,0.0535,0.0343370311674591,0.0299070144118022,0.0393966009845745,0.103539355520338,0.0958817540765541,0.111732951941575,0.167635147032928,0.158145307258128,0.177574327478367,0.0843830805075691,1.50521819708956e-35,0.444054181891207 +"null","constant","tss",4,40,"wald_ml",150,6000,5758,242,147,0.0403333333333333,0.067210837096214,0.0610277913962478,0.0739709697428515,0.134769017019799,0.126192001748904,0.14383303624288,0.195206668982286,0.185173587708648,0.205646165822299,0.104078202508359,1.3341509762296e-54,0.433539982517614 +"null","constant","tss",12,40,"lrt",150,6000,5986,14,0,0.00233333333333333,0.0143668559973271,0.0116487013673261,0.0177079118301315,0.0649849649181423,0.0590152557624449,0.0715126495578783,0.115603073838958,0.107748388407623,0.123950809035566,0.0376376568717833,8.62249043869314e-08,0.485492376920367 +"null","constant","tss",12,40,"wald_ml",150,6000,5989,11,6,0.00183333333333333,0.0213725162798464,0.0180049117092347,0.0253537288632171,0.0728001335782267,0.066490416983351,0.0796575271661756,0.127233260978461,0.119031959024314,0.135912455794454,0.0445719596399606,9.25695055432286e-11,0.481410606985188 +"null","original","rle",4,40,"lrt",150,5997,5795,202,0,0.0336835084208771,0.0291630716134599,0.0251329702386332,0.0338169874810617,0.100431406384815,0.0929553445016146,0.108436858943485,0.167385677308024,0.15799492693621,0.177217110265444,0.0999852967297406,9.5741484666021e-51,0.437673327317577 +"null","original","rle",4,40,"wald_ml",150,5997,5843,154,79,0.0256795064198766,0.0641793599178504,0.0581774162795938,0.0707539843971647,0.138113982543214,0.129504908831245,0.147198584878162,0.199897313024131,0.189841646759709,0.210347322840407,0.121629858990801,1.65933724634618e-75,0.426544757195022 +"null","original","rle",12,40,"lrt",150,6000,5999,1,0,0.000166666666666667,0.0170028338056343,0.0140269075883063,0.0205969385096592,0.0675112518753126,0.0614348261387838,0.0741412113767197,0.126187697949658,0.118023356072251,0.134830474773197,0.0374388440525218,9.94059353574812e-08,0.476768173372017 +"null","original","rle",12,40,"wald_ml",150,6000,5999,1,0,0.000166666666666667,0.024837472912152,0.021192855359116,0.029090241552791,0.0801800300050008,0.0735734789459446,0.0873239003112516,0.138523087181197,0.130012534390866,0.14749628709781,0.0425443113010688,7.40678934295594e-10,0.472785456890374 +"null","original","tss",3,40,"lrt",200,7988,7315,673,0,0.0842513770655984,0.045933014354067,0.0413694275154052,0.0509732557662729,0.131237183868763,0.12369250293294,0.139168971747875,0.203417634996582,0.19434974613276,0.21279685973194,0.135182796664952,1.55071172931089e-116,0.415587740716814 +"null","original","tss",3,40,"wald_ml",200,7988,7476,512,273,0.0640961442163245,0.0918940609951846,0.0855537288027265,0.0986535791203389,0.179240235420011,0.170711249407119,0.188098689760206,0.250668806848582,0.240974270777598,0.260619443404176,0.161656443920632,4.03457010651865e-170,0.401366722446981 +"null","original","tss",4,40,"lrt",200,7997,7708,289,0,0.0361385519569839,0.0303580695381422,0.026755637175536,0.0344283824052878,0.102880124545926,0.0962945813035816,0.109861298236999,0.169434353918007,0.161224866817917,0.177973166853889,0.0933013909161196,1.0455020069296e-58,0.443086204545911 +"null","original","tss",4,40,"wald_ml",200,7997,7783,214,106,0.0267600350131299,0.0628292432224078,0.0576509835270075,0.0684388391582219,0.13825003212129,0.130760030253546,0.146096956012513,0.203520493383014,0.194723061730016,0.212610447703097,0.113643078400003,9.87146482984532e-88,0.431436851581758 +"null","original","tss",6,40,"lrt",200,8000,7904,96,0,0.012,0.0232793522267206,0.0201794197046358,0.0268424459659169,0.086918016194332,0.0809063372182911,0.09333102780675,0.148658906882591,0.14098681509377,0.156672346521709,0.0645951152356356,4.52120174669666e-29,0.461303882402388 +"null","original","tss",6,40,"wald_ml",200,8000,7937,63,28,0.007875,0.0447272269119315,0.0403957707152817,0.0494991683188881,0.111251102431649,0.104520582486221,0.118357744451632,0.170719415396245,0.162601431411705,0.179155984712773,0.0787526584971869,3.50452954348045e-43,0.453506610955192 +"null","original","tss",8,40,"lrt",200,8000,7974,26,0,0.00325,0.0171808377225984,0.0145524202088425,0.0202742256019111,0.0701028342111863,0.0647034005936096,0.0759162726113939,0.128668171557562,0.121497414397098,0.136196533199172,0.0453296175635426,1.17327512113714e-14,0.473696316973658 +"null","original","tss",8,40,"wald_ml",200,8000,7984,16,6,0.002,0.0293086172344689,0.0258291673085302,0.0332407906085308,0.0886773547094188,0.0826378982356264,0.0951124322059991,0.145415831663327,0.137853814201029,0.153318897577363,0.0552258902999433,1.41431816022389e-21,0.467708032138111 +"null","original","tss",12,40,"lrt",200,8000,7997,3,0,0.000375,0.0171314242841065,0.0145105151972706,0.0202160145372216,0.0645241965737151,0.059345808818977,0.0701207559357476,0.119294735525822,0.112372755059871,0.126582293438552,0.0329734074224069,5.6092770626337e-08,0.480625097706312 +"null","original","tss",12,40,"wald_ml",200,8000,7999,1,1,0.000125,0.0251281410176272,0.0219194190982992,0.028792751186341,0.0757594699337417,0.0701620857850788,0.0817641350549703,0.130016252031504,0.122823186861121,0.137564510376576,0.0381036842306565,1.63501153735572e-10,0.476689665536424 +"null","original","tss",20,40,"lrt",200,8000,8000,0,0,0,0.014125,0.0117624417845792,0.0169539514627402,0.058375,0.0534462764601491,0.0637276410450682,0.11025,0.103572963173447,0.117301159323374,0.020738749813529,0.00205314046936761,0.487703822259244 +"null","original","tss",20,40,"wald_ml",200,8000,8000,0,0,0,0.018625,0.0158851798715362,0.0218268938087898,0.0645,0.0593234953117845,0.0700945427828028,0.116,0.109166456166853,0.123202146883254,0.026145410693358,3.55641366754023e-05,0.48530796517549 +"null","sparse","tss",4,40,"lrt",150,5365,4380,985,0,0.183597390493942,0.0237442922374429,0.0196354002458669,0.0286878479396827,0.0817351598173516,0.0739836169914801,0.0902197342433841,0.144063926940639,0.133976279706903,0.154775371262067,0.0691867976694867,1.2302251751396e-18,0.455338729993458 +"null","sparse","tss",4,40,"wald_ml",150,5365,4554,811,236,0.15116495806151,0.0408432147562582,0.0354710856191992,0.0469893207412142,0.106719367588933,0.0980810744949493,0.116020593506725,0.170838823012736,0.160186222147154,0.182046273852014,0.0823244156861178,3.1118285502842e-27,0.44737008911558 +"null","sparse","tss",12,40,"lrt",150,5859,5562,297,0,0.0506912442396313,0.015282272563826,0.0123767061844007,0.0188569284972326,0.0604099244875944,0.054446958150994,0.0669796875450717,0.110931319669184,0.102945034249156,0.119454663618446,0.0326604118895657,1.40505153549206e-05,0.480162876414684 +"null","sparse","tss",12,40,"wald_ml",150,5859,5594,265,33,0.0452295613585936,0.0196639256346085,0.0163415019144068,0.0236456003405133,0.0682874508401859,0.0619693876083744,0.075198030309613,0.119413657490168,0.111176063604179,0.128173597988946,0.0345742276076426,3.11048126202616e-06,0.476462685961144 +"null","unequal_lib","tss",4,40,"lrt",150,5999,5809,190,0,0.0316719453242207,0.0395937338612498,0.0348758076860468,0.0449201857375891,0.113616801514891,0.105710109456138,0.122034181902751,0.17868824238251,0.169050117338897,0.188751049879355,0.0933208396676188,2.28941310680147e-44,0.440256995667013 +"null","unequal_lib","tss",4,40,"wald_ml",150,5999,5898,101,56,0.0168361393565594,0.0764665988470668,0.0699568766474143,0.0835276696829271,0.14852492370295,0.139678040812683,0.157829350917821,0.211766700576467,0.20152919358122,0.222379424818132,0.115545521434757,8.05155216532249e-69,0.429486140222383 +"null","unequal_lib","tss",12,40,"lrt",150,6000,5995,5,0,0.000833333333333333,0.0128440366972477,0.0102895388904349,0.0160224515004104,0.0618849040867389,0.0560617457461293,0.0682691711384984,0.116263552960801,0.108394130866014,0.124624439205597,0.023602226833183,0.00251353912720786,0.48526062352396 +"null","unequal_lib","tss",12,40,"wald_ml",150,6000,6000,0,1,0,0.0196666666666667,0.0164483267012244,0.0234996733360393,0.0741666666666667,0.0678051993175382,0.0810730588695571,0.13,0.121726673024468,0.138746803756492,0.0305808238880169,2.67468380464776e-05,0.481384434000172 +"power","original","tss",4,40,"lrt",150,5997,5768,229,0,0.0381857595464399,0.0325109522711552,0.0276325349035269,0.038216785218464,0.114134194143417,0.105010338234011,0.123940999484436,0.183306433018215,0.172073169738517,0.195100216790487,0.105618632297901,1.89778614361325e-42,0.434784684139481 +"power","original","tss",4,40,"wald_ml",150,5997,5841,156,88,0.0260130065032516,0.0721602549510585,0.0648766537321604,0.0801914513597557,0.150466651491008,0.140199738589643,0.161344329177655,0.218757113589802,0.206780906784307,0.23122475633943,0.127157614018003,4.02317867444891e-62,0.422783825089475 +"power","original","tss",12,40,"lrt",150,6000,5997,3,0,5e-04,0.0237777777777778,0.0197159062402261,0.0286520171896841,0.0831111111111111,0.0753968197212776,0.0915365560857739,0.150888888888889,0.140728778500017,0.161644533543095,0.0677456389208428,2.30351108889158e-18,0.456624861778725 +"power","original","tss",12,40,"wald_ml",150,6000,5998,2,1,0.000333333333333333,0.0328888888888889,0.0280635063862089,0.0385110969967455,0.096,0.0877341616927958,0.104955005267343,0.16,0.1495793559677,0.171000636028456,0.0737015111182361,1.17376844197725e-21,0.452433541509862 diff --git a/bench/null_calibration/results/power_summary.csv b/bench/null_calibration/results/power_summary.csv new file mode 100644 index 00000000..b9685f56 --- /dev/null +++ b/bench/null_calibration/results/power_summary.csv @@ -0,0 +1,5 @@ +"shape","offset","n_per_group","test","n_sims","power_raw_05","power_bh_05","bh_fdr_05","mean_discoveries" +"original","tss",4,"lrt",150,0.413696715583508,0.147449336128581,0.183793650793651,1.89333333333333 +"original","tss",4,"wald_ml",150,0.479972375690608,0.290055248618785,0.379349206349206,4.7 +"original","tss",12,"lrt",150,0.668002672010688,0.421509686038744,0.100953102453102,4.82 +"original","tss",12,"wald_ml",150,0.688251001335114,0.509345794392523,0.128748011248011,5.99333333333333 diff --git a/bench/null_calibration/results/provenance_summary.csv b/bench/null_calibration/results/provenance_summary.csv new file mode 100644 index 00000000..5cc78144 --- /dev/null +++ b/bench/null_calibration/results/provenance_summary.csv @@ -0,0 +1,17 @@ +cell,date_utc,args,mass_mode,glmGamPoi,edgeR,tests_run,r_version,platform,webr +null_constant_tss_n12_t40,2026-10-22 22:42:23 UTC,"--shapes=constant --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_constant_tss_n4_t40,2026-10-22 22:43:30 UTC,"--shapes=constant --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_rle_n12_t40,2026-10-22 22:54:06 UTC,"--shapes=original --offsets=rle --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_rle_n4_t40,2026-10-22 22:55:33 UTC,"--shapes=original --offsets=rle --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n12_t40,2026-10-22 22:32:15 UTC,"--shapes=original --offsets=tss --n-per=12 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n20_t40,2026-10-22 22:37:21 UTC,"--shapes=original --offsets=tss --n-per=20 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n3_t40,2026-10-22 22:27:17 UTC,"--shapes=original --offsets=tss --n-per=3 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n4_t40,2026-10-22 22:32:41 UTC,"--shapes=original --offsets=tss --n-per=4 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n6_t40,2026-10-22 22:38:09 UTC,"--shapes=original --offsets=tss --n-per=6 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_original_tss_n8_t40,2026-10-22 22:27:17 UTC,"--shapes=original --offsets=tss --n-per=8 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_sparse_tss_n12_t40,2026-10-22 22:46:13 UTC,"--shapes=sparse --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_sparse_tss_n4_t40,2026-10-22 22:47:35 UTC,"--shapes=sparse --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_unequal_lib_tss_n12_t40,2026-10-22 22:50:11 UTC,"--shapes=unequal_lib --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +null_unequal_lib_tss_n4_t40,2026-10-22 22:51:15 UTC,"--shapes=unequal_lib --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +power_original_tss_n12_t40_lfc1,2026-10-22 22:57:57 UTC,"--shapes=original --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=power --log2fc=1 --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 +power_original_tss_n4_t40_lfc1,2026-10-22 22:59:41 UTC,"--shapes=original --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=power --log2fc=1 --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0 diff --git a/bench/null_calibration/summarise.R b/bench/null_calibration/summarise.R new file mode 100644 index 00000000..21231182 --- /dev/null +++ b/bench/null_calibration/summarise.R @@ -0,0 +1,131 @@ +#!/usr/bin/env Rscript +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Aggregate calibration.R's per-cell p-value CSVs into the #193 summary +# tables. Base R only, so it runs anywhere calibration.R's output lands. +# +# Rscript bench/null_calibration/summarise.R \ +# --results=bench/null_calibration/results/raw \ +# --out=bench/null_calibration/results +# +# Writes calibration_summary.csv (one row per cell x test: rejection rates +# at alpha = 0.01/0.05/0.10 with Wilson 95% intervals, Kolmogorov-Smirnov D +# and p against Uniform(0,1), fit-failure share) and, when power cells are +# present, power_summary.csv (raw and BH power on truly-changed taxa, and +# the Benjamini-Hochberg false discovery rate at 0.05 averaged per dataset). +# +# A "test" here is one taxon in one simulated dataset that produced a +# p-value (production's BH family: status ok or boundary). Failures are +# counted separately and never enter a rate's denominator silently: every +# row reports both n_pvalues and n_failed. + +argv <- commandArgs(trailingOnly = TRUE) +if (length(argv) == 0 && nzchar(Sys.getenv("NB_CALIB_ARGS"))) { + argv <- strsplit(trimws(Sys.getenv("NB_CALIB_ARGS")), "\\s+")[[1]] +} +args <- list() +for (a in argv) { + if (!grepl("^--[a-z0-9-]+=", a)) stop("unrecognised argument '", a, "' (want --key=value)") + kv <- strsplit(sub("^--", "", a), "=", fixed = TRUE)[[1]] + args[[kv[1]]] <- paste(kv[-1], collapse = "=") +} +results_dir <- if (!is.null(args$results)) args$results else + file.path("bench", "null_calibration", "results", "raw") +out_dir <- if (!is.null(args$out)) args$out else + file.path("bench", "null_calibration", "results") +dir.create(out_dir, recursive = TRUE, showWarnings = FALSE) + +files <- list.files(results_dir, pattern = "^pvalues_.*[.]csv$", full.names = TRUE) +if (length(files) == 0) stop("no pvalues_*.csv under ", results_dir) +dat <- do.call(rbind, lapply(files, function(f) { + d <- read.csv(f, stringsAsFactors = FALSE) + d$file <- basename(f) + d +})) +if (!"taxa" %in% names(dat)) dat$taxa <- NA_integer_ +dat$has_p <- !is.na(dat$pvalue) & dat$status %in% c("ok", "boundary") + +# Wilson score interval for k/n at 95% (the same interval prop.test uses +# without continuity correction, written out so the bench needs nothing +# beyond base R distributions). +wilson <- function(k, n) { + if (n == 0) return(c(NA_real_, NA_real_, NA_real_)) + z <- qnorm(0.975) + p <- k / n + den <- 1 + z^2 / n + mid <- (p + z^2 / (2 * n)) / den + half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / den + c(p, max(0, mid - half), min(1, mid + half)) +} + +key_of <- function(d) paste(d$truth, d$shape, d$offset, d$n_per, d$taxa, d$test, sep = "|") +groups <- split(dat, key_of(dat)) + +summ_rows <- lapply(groups, function(d) { + p <- d$pvalue[d$has_p] + # Under a power cell only the truly-null taxa calibrate the test. + p_null <- d$pvalue[d$has_p & d$true_log2fc == 0] + ks <- if (length(p_null) >= 10) { + k <- suppressWarnings(ks.test(p_null, "punif")) + c(unname(k$statistic), k$p.value) + } else c(NA_real_, NA_real_) + w <- function(alpha) wilson(sum(p_null < alpha), length(p_null)) + w01 <- w(0.01); w05 <- w(0.05); w10 <- w(0.10) + data.frame( + truth = d$truth[1], shape = d$shape[1], offset = d$offset[1], + n_per_group = d$n_per[1], n_taxa_setting = d$taxa[1], test = d$test[1], + n_sims = length(unique(d$sim)), n_attempted = nrow(d), + n_pvalues = length(p), n_failed = sum(d$status == "failed"), + n_boundary = sum(d$status == "boundary"), + fail_rate = sum(d$status == "failed") / nrow(d), + rate_01 = w01[1], rate_01_lo = w01[2], rate_01_hi = w01[3], + rate_05 = w05[1], rate_05_lo = w05[2], rate_05_hi = w05[3], + rate_10 = w10[1], rate_10_lo = w10[2], rate_10_hi = w10[3], + ks_D = ks[1], ks_p = ks[2], mean_p_null = if (length(p_null)) mean(p_null) else NA_real_, + stringsAsFactors = FALSE + ) +}) +summary_df <- do.call(rbind, summ_rows) +summary_df <- summary_df[order(summary_df$truth, summary_df$shape, summary_df$offset, + summary_df$n_per_group, summary_df$test), ] +write.csv(summary_df, file.path(out_dir, "calibration_summary.csv"), row.names = FALSE) +cat(sprintf("wrote %s (%d rows)\n", file.path(out_dir, "calibration_summary.csv"), nrow(summary_df))) + +pow <- dat[dat$truth == "power", , drop = FALSE] +if (nrow(pow) > 0) { + pow$padj <- NA_real_ + # BH within each simulated dataset and test, over that dataset's tested + # family - the same family production adjusts over. + fams <- split(seq_len(nrow(pow)), + paste(pow$file, pow$sim, pow$test, sep = "|")) + for (idx in fams) { + fam <- idx[pow$has_p[idx]] + if (length(fam) > 0) pow$padj[fam] <- p.adjust(pow$pvalue[fam], method = "BH") + } + pg <- split(pow, key_of(pow)) + prows <- lapply(pg, function(d) { + de <- d[d$true_log2fc != 0 & d$has_p, , drop = FALSE] + by_sim <- split(d, d$sim) + fdp <- vapply(by_sim, function(s) { + disc <- s$has_p & !is.na(s$padj) & s$padj < 0.05 + if (!any(disc)) 0 else mean(s$true_log2fc[disc] == 0) + }, numeric(1)) + data.frame( + shape = d$shape[1], offset = d$offset[1], n_per_group = d$n_per[1], + test = d$test[1], n_sims = length(by_sim), + power_raw_05 = if (nrow(de)) mean(de$pvalue < 0.05) else NA_real_, + power_bh_05 = if (nrow(de)) mean(!is.na(de$padj) & de$padj < 0.05) else NA_real_, + bh_fdr_05 = mean(fdp), + mean_discoveries = mean(vapply(by_sim, function(s) + sum(s$has_p & !is.na(s$padj) & s$padj < 0.05), numeric(1))), + stringsAsFactors = FALSE + ) + }) + power_df <- do.call(rbind, prows) + power_df <- power_df[order(power_df$shape, power_df$offset, + power_df$n_per_group, power_df$test), ] + write.csv(power_df, file.path(out_dir, "power_summary.csv"), row.names = FALSE) + cat(sprintf("wrote %s (%d rows)\n", file.path(out_dir, "power_summary.csv"), nrow(power_df))) +} +cat("summarise complete\n") diff --git a/bench/null_calibration/webr_run.mjs b/bench/null_calibration/webr_run.mjs new file mode 100644 index 00000000..50eedd29 --- /dev/null +++ b/bench/null_calibration/webr_run.mjs @@ -0,0 +1,104 @@ +#!/usr/bin/env node +// SPDX-License-Identifier: AGPL-3.0-only +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +// +// Sandbox-only driver: run calibration.R / summarise.R under WebR (R in +// WebAssembly) when no native R is available. The R code executed is the +// bench's own, unmodified; MASS comes from sourcing MASS's pure-R negbin.R +// (glm.nb makes no compiled calls), which calibration.R detects and records +// in its provenance output. Prefer native Rscript with the renv.lock pin +// whenever it exists - it is faster and carries glmGamPoi/edgeR. +// +// One-time setup (outside the repo, so no repo lockfile changes): +// mkdir ~/webr-probe && cd ~/webr-probe && npm init -y && npm i webr@0.6.0 +// git clone https://github.com/cran/MASS.git mass-src +// git -C mass-src checkout b6a2274 # MASS 7.3-65, the renv.lock pin +// +// Run, from the repo root: +// WEBR_PATH=~/webr-probe/node_modules/webr \ +// node bench/null_calibration/webr_run.mjs \ +// --mass-r ~/webr-probe/mass-src/R/negbin.R \ +// --script calibration -- --n-per=4 --sims=200 --tests=wald_ml,lrt +// Files land in bench/null_calibration/results/raw (override with +// --host-out). For summarise, existing raw CSVs in --host-out are mounted +// in first, and the summary CSVs are copied back to --host-summary-out +// (default bench/null_calibration/results). + +import { existsSync, mkdirSync, readFileSync, readdirSync, writeFileSync } from 'node:fs'; +import { dirname, join, resolve } from 'node:path'; +import { fileURLToPath, pathToFileURL } from 'node:url'; + +const here = dirname(fileURLToPath(import.meta.url)); +const argv = process.argv.slice(2); +const sep = argv.indexOf('--'); +const driverArgs = sep === -1 ? argv : argv.slice(0, sep); +const rArgs = sep === -1 ? [] : argv.slice(sep + 1); + +function opt(name, fallback) { + const hit = driverArgs.find((a) => a.startsWith(`${name}=`)); + if (hit) return hit.slice(name.length + 1); + const i = driverArgs.indexOf(name); + return i === -1 ? fallback : driverArgs[i + 1]; +} + +const script = opt('--script', 'calibration'); +const massR = opt('--mass-r', ''); +const hostOut = resolve(opt('--host-out', join(here, 'results', 'raw'))); +const hostSummaryOut = resolve(opt('--host-summary-out', join(here, 'results'))); +const webrPath = process.env.WEBR_PATH || 'webr'; +// WEBR_PATH may name the package directory or its ESM entry directly. +const webrEntry = webrPath === 'webr' ? 'webr' + : pathToFileURL(resolve(webrPath.endsWith('.mjs') ? webrPath : join(webrPath, 'dist', 'webr.mjs'))).href; +const { WebR } = await import(webrEntry); + +// No stdout/stderr callbacks here: under Node the WebR worker cannot clone +// function options, and the defaults already stream R output to this +// process's stdout/stderr line by line. +const webR = new WebR({ interactive: false }); +await webR.init(); + +const put = async (vfsPath, hostPath) => { + const dir = vfsPath.slice(0, vfsPath.lastIndexOf('/')); + if (dir) await webR.evalRVoid(`dir.create("${dir}", recursive = TRUE, showWarnings = FALSE)`); + await webR.FS.writeFile(vfsPath, readFileSync(hostPath)); +}; + +await put('/work/calibration.R', join(here, 'calibration.R')); +await put('/work/summarise.R', join(here, 'summarise.R')); +await put('/work/differential.jl', join(here, '..', '..', 'src', 'analysis', 'differential.jl')); +if (massR) await put('/work/negbin.R', resolve(massR.replace(/^~/, process.env.HOME || ''))); + +// Mount any existing raw results so summarise (and resumed calibration +// cells appending to the same CSV) see them. +if (existsSync(hostOut)) { + for (const f of readdirSync(hostOut)) { + if (/\.(csv|txt)$/.test(f)) await put(`/work/results/raw/${f}`, join(hostOut, f)); + } +} + +const benchArgs = [...rArgs.filter((a) => !a.startsWith('--differential-jl=') && !a.startsWith('--out=') && !a.startsWith('--results=')), + '--differential-jl=/work/differential.jl', + // calibration.R writes raw cells to --out; summarise.R writes summaries to + // --out. They are different directories. + script === 'summarise' ? '--out=/work/results' : '--out=/work/results/raw', + '--results=/work/results/raw']; +const quote = (s) => `'${s.replace(/'/g, "'\\''")}'`; +await webR.evalRVoid(`Sys.setenv(NB_CALIB_ARGS = ${quote(benchArgs.join(' '))})`); +if (massR) await webR.evalRVoid('source("/work/negbin.R")'); +await webR.evalRVoid('setwd("/work")'); +await webR.evalRVoid(`source("/work/${script}.R")`); + +mkdirSync(hostOut, { recursive: true }); +mkdirSync(hostSummaryOut, { recursive: true }); +for (const [vfsDir, hostDir] of [['/work/results/raw', hostOut], ['/work/results', hostSummaryOut]]) { + const listing = await (await webR.evalR( + `list.files("${vfsDir}", pattern = "\\\\.(csv|txt)$")`)).toJs(); + for (const name of listing.values) { + if (vfsDir.endsWith('/raw') === false && name.startsWith('pvalues_')) continue; + const bytes = await webR.FS.readFile(`${vfsDir}/${name}`); + writeFileSync(join(hostDir, name), bytes); + console.log(`[webr_run] ${name} -> ${join(hostDir, name)}`); + } +} +await webR.close(); +process.exit(0); diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 44b2c882..4bc1d023 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -364,7 +364,13 @@ analysis: # across all taxa in one call -- sharing information across taxa, which # steadies low-count taxa -- and fits each taxon through # MASS::negative.binomial(theta = 1/dispersion), so the Wald test is - # unchanged. Requires the glmGamPoi R package. + # unchanged. Requires the glmGamPoi R package. Both fixed-theta Wald + # tests treat the dispersion they use as known rather than propagating + # its estimation uncertainty: every ml result therefore carries a + # calibration caveat in its diagnostics, strongest when a group is + # small, as does a glmGamPoi result whose smaller group has at most 12 + # samples. The caveat records the limitation; it changes no p-value. + # Evidence and the replication protocol: bench/null_calibration. dispersion: ml # 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 4ea20add..926d8bfd 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -446,6 +446,18 @@ export interface DifferentialResult { } diagnostics: { n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number + /** + * nb_glm only: interpretation caveat for the fixed-theta Wald test. + * Present for every `ml` result, scaled by the smaller group, and for + * `glmGamPoi` results whose smaller group has at most 12 samples. It + * changes no statistic. + */ + calibration_caveat?: { + severity: 'info' | 'caution' | 'strong' + dispersion: 'ml' | 'glmGamPoi' + min_group_n: number + message: string + } /** clr_welch only. */ zero_replacement?: { method: string; delta: number; zeros_replaced: number; n_taxa_in_composition: number diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx index 825fb812..8ee91612 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -135,6 +135,15 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { `and at most ${(100 * z.max_imputed_fraction).toPrecision(2)}% of any, and are not measurements.` })()}

+ {result.diagnostics.calibration_caveat && ( +

+ Calibration caveat: {result.diagnostics.calibration_caveat.message} +

+ )} diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index f4cc3fbe..a28a9de2 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -152,7 +152,7 @@ export const CONFIG_DESCRIPTIONS: Record = { 'analysis.beta.depth': 'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.', 'analysis.differential.method': 'Model fitted to each taxon: nb_glm (negative-binomial GLM on read counts; effect is a log2 fold change) or clr_welch (Welch\'s t-test on centred log-ratios over the taxa passing min_prevalence, after zCompositions::cmultRepl zero replacement; effect is a CLR difference, not a fold change).', 'analysis.differential.offset': 'Library-size offset in each negative-binomial model: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none). Used by nb_glm only.', - 'analysis.differential.dispersion': 'nb_glm only: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits.', + 'analysis.differential.dispersion': 'nb_glm only: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits. Both Wald tests treat the dispersion as known, so ml results always carry a calibration caveat (strongest in small groups; see bench/null_calibration).', 'analysis.differential.min_prevalence': 'Fraction of samples, 0 to 1, in which a taxon must have reads to be tested; taxa below it are listed as filtered.', 'analysis.differential.replacement_delta': 'clr_welch only: cmultRepl\'s frac, between 0 and 1; an imputed proportion above its taxon\'s smallest observed proportion is lowered to this fraction of it. Replaced zeros are not measurements.', 'analysis.alpha.show_points': 'Show individual samples overlaid on alpha comparison boxplots.', diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts index f2143e07..d035ba8e 100644 --- a/frontend/src/components/differentialTable.test.ts +++ b/frontend/src/components/differentialTable.test.ts @@ -74,3 +74,34 @@ describe('clr_welch results', () => { expect(differentialCsv(clr).split('\r\n')[2]).toBe(`# ${configLine(clr)}`) }) }) + +describe('calibration caveat', () => { + const caveated: DifferentialResult = { + ...result, + config: { method: 'nb_glm', offset: 'tss', dispersion: 'ml', min_prevalence: 0 }, + diagnostics: { + ...result.diagnostics, + calibration_caveat: { + severity: 'strong', dispersion: 'ml', min_group_n: 2, + message: 'Wald p-values treat the per-taxon theta as known.', + }, + }, + } + + test('the settings line records the dispersion the result carries', () => { + expect(configLine(caveated)).toBe('method nb_glm, offset tss, dispersion ml, min_prevalence 0') + }) + + test('the caveat travels into the CSV header, after the settings line', () => { + const lines = differentialCsv(caveated).trimEnd().split('\r\n') + expect(lines[2]).toBe(`# ${configLine(caveated)}`) + expect(lines[3]).toBe('# Calibration caveat (strong): Wald p-values treat the per-taxon theta as known.') + expect(lines[4]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,prevalence,note') + expect(lines).toHaveLength(7) + }) + + test('a result without a caveat keeps the three-line header', () => { + const lines = differentialCsv(result).trimEnd().split('\r\n') + expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,prevalence,note') + }) +}) diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts index 90a547a8..da5f07e2 100644 --- a/frontend/src/components/differentialTable.ts +++ b/frontend/src/components/differentialTable.ts @@ -22,10 +22,14 @@ export function csvField(value: string | number | null): string { */ export function differentialCsv(result: DifferentialResult): string { const { reference, contrast } = result.groups + const caveat = result.diagnostics.calibration_caveat const header = [ `# Differential abundance: ${contrast} vs ${reference} (reference), rank ${result.rank}, table ${result.table}`, `# ${result.method}`, `# ${configLine(result)}`, + // The caveat travels with the numbers it qualifies: a downloaded table + // must not look better-calibrated than the result it came from. + ...(caveat ? [`# Calibration caveat (${caveat.severity}): ${caveat.message}`] : []), ].map(l => l.replace(/[\r\n]+/g, ' ')) const rows = result.rows.map(r => COLUMNS.map(c => csvField(r[c])).join(',')) return [...header, COLUMNS.join(','), ...rows].join('\r\n') + '\r\n' @@ -48,7 +52,9 @@ export function configLine(result: DifferentialResult): string { const c = result.config return c.method === 'clr_welch' ? `method clr_welch, replacement_delta ${c.replacement_delta}, min_prevalence ${c.min_prevalence}; estimate is a CLR difference, not a fold change` - : `method nb_glm, offset ${c.offset}, min_prevalence ${c.min_prevalence}` + // Older saved results may not record the dispersion; name it only + // when the result carries it, never guess one for the record. + : `method nb_glm, offset ${c.offset}${c.dispersion ? `, dispersion ${c.dispersion}` : ''}, min_prevalence ${c.min_prevalence}` } /** The effect size of one row, from the field the result names in `effect`. */ diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index 66cbf471..33994557 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -68,6 +68,8 @@ in each negative-binomial GLM. Both options use a Wald test of the group coefficient with fixed scale 1; `"glmGamPoi"` is refused for `clr_welch`. The trend is plugged into an NB Wald model; its estimation uncertainty is not propagated. This option does not implement glmGamPoi's quasi-likelihood test. +`differential_abundance` records the fixed-theta Wald limitation of both +dispersions as `diagnostics.calibration_caveat` (see `_calibration_caveat`). Method and offset names are stripped and lowercased; dispersion names are stripped and matched without regard to case, then stored as `"ml"` or `"glmGamPoi"`. Throws `ArgumentError` for unsupported names, incompatible @@ -733,6 +735,66 @@ _row(taxon, prevalence, note) = Dict{String,Any}( "pvalue" => nothing, "padj" => nothing, "dispersion_theta" => nothing, "prevalence" => prevalence) +""" + _calibration_caveat(dispersion, n_ref, n_con) -> Union{Dict{String,Any},Nothing} + +A result-level calibration caveat for the negative-binomial Wald tests, or +`nothing` when none is recorded. Every `"ml"` result carries one, scaled by +the smaller group size: `MASS::glm.nb` estimates theta taxon by taxon and the +fixed-scale Wald test then treats that estimate as known, ignoring its +sampling uncertainty (issue #193). The `"glmGamPoi"` dispersion trend shares +information across taxa, so its caveat is recorded only through 12 samples in +the smaller group; even there, this implementation plugs the trend into the +same fixed-theta Wald test rather than propagating trend uncertainty or using +glmGamPoi's quasi-likelihood F-test. The caveat changes no statistic: it is +provenance and interpretation guidance recorded with the result. +""" +function _calibration_caveat(dispersion::AbstractString, n_ref::Integer, n_con::Integer) + min_group_n = min(n_ref, n_con) + # "caution" runs through 20, the largest group size the bench measured: + # the ml Wald test was still anti-conservative there. Beyond the measured + # range the caveat drops to "info" but, for "ml", never disappears. + severity = min_group_n <= 4 ? "strong" : min_group_n <= 20 ? "caution" : "info" + if dispersion == "ml" + message = if severity == "strong" + "The smallest group has $min_group_n sample(s). With dispersion \"ml\", " * + "theta is estimated separately for each taxon and the fixed-scale Wald test " * + "then treats that estimate as known. Null simulations for issue #193 rejected " * + "substantially too often at 4 samples per group; treat small p-values and " * + "Benjamini-Hochberg adjusted p-values at this size as exploratory rather than " * + "confirmatory. See bench/null_calibration." + elseif severity == "caution" + "The smallest group has $min_group_n samples. With dispersion \"ml\", theta " * + "is estimated separately for each taxon and the fixed-scale Wald test then " * + "treats that estimate as known. Null simulations in bench/null_calibration " * + "(issue #193) found above-nominal rejection at every measured group size " * + "through 20 per group, so interpret borderline p-values and " * + "Benjamini-Hochberg adjusted p-values cautiously." + else + "With dispersion \"ml\", theta is estimated separately for each taxon and " * + "the fixed-scale Wald test treats that estimate as known; theta uncertainty is " * + "not propagated. Null simulations in bench/null_calibration (issue #193) " * + "found above-nominal rejection through 20 samples per group, the largest " * + "size measured; beyond that range the approximation is expected to improve, " * + "but it has not been measured there." + end + return Dict{String,Any}("severity" => severity, "dispersion" => "ml", + "min_group_n" => min_group_n, "message" => message) + elseif dispersion == "glmGamPoi" && min_group_n <= 12 + message = "The smallest group has $min_group_n sample(s). The glmGamPoi " * + "dispersion trend shares information across taxa and steadies theta " * + "relative to a separate maximum-likelihood estimate for every taxon, " * + "but this implementation still plugs the trend into a fixed-theta, " * + "fixed-scale Wald test. It does not propagate trend uncertainty and " * + "does not implement glmGamPoi's quasi-likelihood F-test, so interpret " * + "borderline p-values cautiously at this sample size. " * + "See bench/null_calibration and issue #193." + return Dict{String,Any}("severity" => severity, "dispersion" => "glmGamPoi", + "min_group_n" => min_group_n, "message" => message) + end + nothing +end + """ differential_abundance(counts, samples, taxa, groups; reference, contrast, config=DifferentialConfig()) -> Dict{String,Any} @@ -746,7 +808,12 @@ enter each negative-binomial model as `log(factor)`; the estimate is a natural-log fold change. `config.dispersion` selects per-taxon maximum-likelihood theta (`"ml"`) or fixed theta from `glmGamPoi::glm_gp`'s cross-taxon dispersion trend (`"glmGamPoi"`), fitted in one call on the taxa passing the prevalence filter. -The result's `config` includes the selected dispersion name. +The result's `config` includes the selected dispersion name. Its `diagnostics` +also records a `calibration_caveat` for every `"ml"` result (scaled by the +smaller group) and for a `"glmGamPoi"` result whose smaller group has at most +12 samples: both fixed-theta Wald paths treat the dispersion used by the test +as known rather than propagating its estimation uncertainty. The caveat alters +no estimate, p-value or adjusted p-value. `clr_welch`, in this order: (1) `min_prevalence` chooses the taxa, and taxa with no reads in any sample are always left out; (2) the zeros of those taxa only @@ -834,6 +901,8 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, config_echo = Dict{String,Any}("method" => config.method, "offset" => config.offset, "min_prevalence" => config.min_prevalence, "dispersion" => config.dispersion) + caveat = _calibration_caveat(config.dispersion, n_ref, n_con) + isnothing(caveat) || (diagnostics["calibration_caveat"] = caveat) else (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( "Welch's t-test needs at least 2 samples in each group to estimate its variance; " * diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index dcab0549..17fea750 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -188,6 +188,23 @@ end end end + @testset "Calibration caveat scales with the smaller group" begin + # Pure thresholds, no R needed: every "ml" result is caveated, the + # glmGamPoi trend only through 12 in the smaller group. + caveat = Differential._calibration_caveat + @test caveat("ml", 4, 4)["severity"] == "strong" + @test caveat("ml", 30, 3)["severity"] == "strong" + @test caveat("ml", 5, 5)["severity"] == "caution" + @test caveat("ml", 12, 12)["severity"] == "caution" + @test caveat("ml", 13, 20)["severity"] == "caution" + @test caveat("ml", 13, 20)["min_group_n"] == 13 + @test caveat("ml", 21, 25)["severity"] == "info" + @test caveat("glmGamPoi", 4, 9)["severity"] == "strong" + @test caveat("glmGamPoi", 12, 12)["severity"] == "caution" + @test isnothing(caveat("glmGamPoi", 13, 13)) + @test isnothing(caveat("glmGamPoi", 20, 13)) + end + @testset "Input validation" begin taxa = ["a", "b", "c"] groups = ["A", "A", "B"] @@ -226,6 +243,9 @@ end "TSS size-factor offset; Wald test of the group coefficient; " * "Benjamini-Hochberg adjustment over the fitted taxa" @test res["diagnostics"]["n_tested"] == 3 + caveat = res["diagnostics"]["calibration_caveat"] + @test caveat["severity"] == "caution" && caveat["dispersion"] == "ml" + @test caveat["min_group_n"] == 12 && occursin("theta", caveat["message"]) @test all(r -> r["status"] == "ok", res["rows"]) @test isapprox(byt["t_up"]["estimate"], log(4); atol=0.5) @test isapprox(byt["t_up"]["log2_fold_change"], 2; atol=0.8) @@ -239,6 +259,13 @@ end reference="B", contrast="A") sb = Dict(r["taxon"] => r for r in swapped["rows"]) @test sb["t_up"]["estimate"] ≈ -byt["t_up"]["estimate"] rtol=1e-6 + # The caveat strengthens as the smaller group shrinks; the p-values + # themselves are untouched by it (checked against MASS below). + keep = [1:4; 13:16] + small = differential_abundance(effect[keep, :], effect_samples[keep], effect_taxa, + effect_groups[keep]; reference="A", contrast="B") + @test small["diagnostics"]["calibration_caveat"]["severity"] == "strong" + @test small["diagnostics"]["calibration_caveat"]["min_group_n"] == 4 @testset "matches a direct MASS::glm.nb call" begin off = log.(tss_factors(effect, effect_samples)) @@ -276,6 +303,10 @@ end "with dispersion trend from glmGamPoi::glm_gp, log link) with a " * "TSS size-factor offset; Wald test of the group coefficient; " * "Benjamini-Hochberg adjustment over the fitted taxa" + caveat = res["diagnostics"]["calibration_caveat"] + @test caveat["severity"] == "caution" && caveat["dispersion"] == "glmGamPoi" + @test caveat["min_group_n"] == 12 && + occursin("quasi-likelihood F-test", caveat["message"]) byt = Dict(r["taxon"] => r for r in res["rows"]) off = log.(tss_factors(effect, effect_samples)) g = effect_groups @@ -669,6 +700,7 @@ end @test isempty(res["size_factors"]) @test res["config"]["method"] == "clr_welch" && !haskey(res["config"], "offset") @test res["diagnostics"]["zero_replacement"]["zeros_replaced"] == 0 + @test !haskey(res["diagnostics"], "calibration_caveat") @test occursin("Welch", res["method"]) && occursin("cmultRepl", res["method"]) # Known answer, independent of R: the estimate is mean(B) - mean(A) # of each taxon's CLR. From 1335ccd01af4ff05f56c6bf395f596d3c0ca7415 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:28:51 +0000 Subject: [PATCH 2/2] fix(integration): repair cross-PR method/config/route/type unions The squash merges of #196-#207 left main unable to build: PR-branch conflict resolutions overwrote each other's differential method registry, dropped the compositional include and cladistics routes, duplicated the clade-cumulus include and frontend panel import, and corrupted _differential_settings_error. - MetaManifold.jl: restore analysis/compositional.jl, include clade_cumulus.jl once; server.jl includes routes/cladistics.jl and routes/clade_cumulus.jl. - differential.jl: rebuild the 9-method registry and DifferentialConfig union (mn/dm optimiser + mc_samples), restore ancom_bc/aldex2/ songbird and multinomial dispatch branches alongside the two-part models; routes build the union config; config validation regains the mn/dm cross-key rules. - Frontend: union DifferentialMethod/config/diagnostics types, restore cladistics client imports, alias the cladistics panel so it no longer collides with the clade-cumulus panel, restore per-method CSV config lines and PipelineStages entries for all nine methods. - Tests that pinned pre-union method lists and the 11-column CSV are updated to the union reality. Frontend: tsc clean, bun test 87/87, production build passes. Julia: exercised by CI (no Julia runtime in the authoring sandbox). Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- README.md | 2 +- config/defaults/pipeline.yml | 4 + frontend/src/api/client.ts | 1 + frontend/src/api/types.ts | 24 +- .../src/components/PipelineStages.test.ts | 11 +- frontend/src/components/PipelineStages.tsx | 16 +- .../src/components/differentialTable.test.ts | 8 +- frontend/src/components/differentialTable.ts | 45 ++- frontend/src/views/StudyView.tsx | 4 +- src/MetaManifold.jl | 2 +- src/analysis/differential.jl | 302 +++++++++++++++++- src/server/routes/analysis.jl | 10 + src/server/routes/config.jl | 31 +- src/server/server.jl | 1 + test/unit/test_differential.jl | 7 +- 15 files changed, 429 insertions(+), 39 deletions(-) diff --git a/README.md b/README.md index 7c663947..74629a96 100644 --- a/README.md +++ b/README.md @@ -667,7 +667,7 @@ analysis: normalisation: hellinger # hellinger | none | rarefy | srs depth: 0 # rarefy/srs reads per sample; 0 = smallest non-empty sample differential: - method: nb_glm # nb_glm | clr_welch | ancom_bc | aldex2 | songbird + method: nb_glm # nb_glm | clr_welch | zinb | hurdle_nb | ancom_bc | aldex2 | songbird | multinomial | dirichlet_multinomial dispersion: ml # nb_glm: ml | glmGamPoi replacement_method: bayesian_multiplicative # clr_welch: bayesian_multiplicative | multiplicative_replacement diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index a5d234c2..ef5988f0 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -350,6 +350,10 @@ analysis: # zero-replaced CLR values), zinb (zero-inflated NB), or hurdle_nb (logistic # occurrence plus zero-truncated NB). Two-part models test both group effects # jointly; their count-process fold change is reported separately. + # Compositional alternatives: ancom_bc (bias-corrected log-linear model), + # aldex2 (Dirichlet Monte-Carlo CLR), songbird (CLR-centred multinomial + # regression), multinomial and dirichlet_multinomial (reference-taxon + # log-ratio models; see the mn_* and dm_* settings below). method: nb_glm # How library size enters each model, as a log offset: tss (total reads per # sample) or rle (median-of-ratios over the taxa present in every sample). diff --git a/frontend/src/api/client.ts b/frontend/src/api/client.ts index 98a97087..5b58aaa7 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -23,6 +23,7 @@ import type { ReferenceTreeDoc, ReferenceTreeSummary, ReferenceTreeDetail, AlignmentQC, TrimSettings, CladeCumulusSourcesResponse, CladeCumulusRequest, CladeCumulusResponse, CladeMoveValidationRequest, CladeMoveValidationResponse, + CladisticsSummary, CladisticsDetail, CladisticsCreateRequest, } from './types' import type { FigureDoc } from '../figure/types' diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 9badc9b0..bd59afca 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -423,7 +423,10 @@ export interface DifferentialRow { prevalence: number } -export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' +export type DifferentialMethod = + | 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' + | 'ancom_bc' | 'aldex2' | 'songbird' + | 'multinomial' | 'dirichlet_multinomial' /** Which row field holds the effect size, and what to call it. */ export interface DifferentialEffect { @@ -457,6 +460,16 @@ export interface DifferentialResult { mc_samples?: number /** aldex2 only. */ seed?: number + /** multinomial and dirichlet_multinomial: the taxon every ratio is against. */ + mn_reference_taxon?: string + /** multinomial and dirichlet_multinomial: L2 penalty on the log-odds. */ + mn_penalty?: number + /** dirichlet_multinomial only: how the concentration was estimated. */ + dm_overdispersion_method?: 'mle' | 'mom' + /** multinomial and dirichlet_multinomial: the optimiser settings used. */ + tolerance?: number + max_iter?: number + max_seconds?: number } diagnostics: { n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number @@ -507,6 +520,15 @@ export interface DifferentialResult { n_taxa_in_composition: number; n_taxa_unobserved: number caveat?: string } + /** multinomial and dirichlet_multinomial only. */ + model?: { + reference_taxon: string; penalty: number; n_parameters: number; iterations: number + converged: boolean; log_likelihood: number; runtime_seconds: number; n_taxa_fitted: number + } + /** dirichlet_multinomial only: c and theta = 1 / (1 + c) in (0, 1). */ + overdispersion?: { + method: 'mle' | 'mom'; concentration: number; theta: number; at_bound: boolean + } } rows: DifferentialRow[] figure: unknown diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts index 0037ce55..50d94f86 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -9,9 +9,10 @@ test('dispersion uses the enum editor with both supported models', () => { }) test('differential model editor exposes supported two-part NB methods', () => { - expect(CONFIG_TYPES['analysis.differential.method']).toEqual({ - kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'], - }) + const method = CONFIG_TYPES['analysis.differential.method'] + expect(method.kind).toBe('enum') + expect(method.options).toContain('zinb') + expect(method.options).toContain('hurdle_nb') }) test('advanced zero handling exposes both methods, bounded delta and nullable alpha', () => { @@ -24,10 +25,10 @@ test('advanced zero handling exposes both methods, bounded delta and nullable al expect(CONFIG_DESCRIPTIONS['analysis.differential.bayesian_alpha']).toContain('no posterior draws') }) -test('differential method enum lists all five compositional and count models', () => { +test('differential method enum lists every count and compositional model', () => { expect(CONFIG_TYPES['analysis.differential.method']).toEqual({ kind: 'enum', - options: ['nb_glm', 'clr_welch', 'ancom_bc', 'aldex2', 'songbird'], + options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', 'multinomial', 'dirichlet_multinomial'], }) expect(CONFIG_TYPES['analysis.differential.mc_samples']).toEqual({ kind: 'int' }) expect(CONFIG_DESCRIPTIONS['analysis.differential.mc_samples']).toContain('Dirichlet') diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index 4dade241..9ef5ab24 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -150,13 +150,20 @@ export const CONFIG_DESCRIPTIONS: Record = { 'analysis.alpha.rarefaction_iterations': 'rarefy draws averaged for each alpha diversity value.', 'analysis.beta.normalisation': 'Preparation before the Bray-Curtis dissimilarity used by NMDS and PERMANOVA: hellinger (square root of relative abundances), none, rarefy or srs.', 'analysis.beta.depth': 'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.', - 'analysis.differential.method': 'Model per taxon: nb_glm (negative-binomial GLM); clr_welch (Welch test on CLR values after configured zero replacement); zinb (zero-inflated negative binomial); hurdle_nb (logistic zero process plus zero-truncated negative binomial). Two-part NB p-values jointly test group effects in both processes; use only with evidence for excess zeros and interpret cautiously.', + 'analysis.differential.method': 'Model fitted to the data. Per taxon: nb_glm (negative-binomial GLM, effect a log2 fold change); clr_welch (Welch test on CLR values after configured zero replacement; effect a CLR difference, not a fold change); zinb (zero-inflated negative binomial); hurdle_nb (logistic zero process plus zero-truncated negative binomial). Two-part NB p-values jointly test group effects in both processes; use only with evidence for excess zeros and interpret cautiously. Compositional: ancom_bc (bias-corrected log-linear model, Lin & Peddada 2020); aldex2 (Dirichlet Monte-Carlo CLR effect sizes); songbird (CLR-centred multinomial differentials); multinomial (one multinomial logistic model for the whole composition, effect a log2 fold change against mn_reference_taxon); dirichlet_multinomial (the same plus a Dirichlet prior absorbing between-sample overdispersion).', 'analysis.differential.offset': 'Library-size offset in the count component of nb_glm, zinb and hurdle_nb: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none).', 'analysis.differential.dispersion': 'nb_glm only: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits. Both Wald tests treat the dispersion as known, so ml results always carry a calibration caveat (strongest in small groups; see bench/null_calibration).', 'analysis.differential.min_prevalence': 'Fraction of samples, 0 to 1, in which a taxon must have reads to be tested; taxa below it are listed as filtered.', 'analysis.differential.replacement_method': 'clr_welch only: bayesian_multiplicative uses cmultRepl GBM posterior estimates (Martín-Fernández et al. 2015); multiplicative_replacement uses CZM (Martín-Fernández et al. 2003). Both preserve observed-part ratios, but every replacement is biased. For sparse data consider NB on raw counts, prevalence filtering, occupancy or zero-inflated models instead.', 'analysis.differential.replacement_delta': "clr_welch only: cmultRepl frac, strictly between 0 and 1. CZM starts at delta × 0.5/sample total; either method caps high imputations at delta × the taxon's smallest observed proportion. Below 0.01 or at least 0.9 warns. Do not tune delta for significance; every replacement is biased.", 'analysis.differential.bayesian_alpha': "Bayesian replacement only: null uses GBM's estimated Dirichlet concentration; a positive number uses that concentration with the same leave-one-out prior mean. The posterior parameters are recorded; no posterior draws are generated. Larger alpha gives the prior more weight.", + 'analysis.differential.mn_reference_taxon': 'multinomial and dirichlet_multinomial only: the taxon whose log-odds are the denominator of every reported log ratio; empty picks the most abundant fitted taxon. It must have reads in at least half the samples, because every ratio is anchored on it. A log ratio of two taxa is coherent under compositionality: it does not depend on the library size.', + 'analysis.differential.mn_penalty': "multinomial and dirichlet_multinomial only: L2 (ridge) penalty on the log-odds coefficients, 0 for maximum likelihood. A penalty keeps a taxon with almost no reads in one group from dominating the joint fit, and the reported standard errors are the penalised ones. Songbird's L1 sparsity path is not implemented: a lasso estimate has no valid Wald standard error.", + 'analysis.differential.dm_overdispersion_method': 'dirichlet_multinomial only: mle alternates the profile-likelihood Dirichlet concentration with the coefficients; mom takes one method-of-moments estimate from the multinomial composition. The result reports the concentration c and theta = 1/(1+c), the intra-sample correlation; theta near 0 is the multinomial.', + 'analysis.differential.max_iter': 'multinomial and dirichlet_multinomial only: Newton iteration cap for one fit; a fit that stops at the cap is refused rather than reported.', + 'analysis.differential.tolerance': 'multinomial and dirichlet_multinomial only: convergence tolerance, as the largest coefficient update relative to the size of the coefficients.', + 'analysis.differential.max_seconds': 'multinomial and dirichlet_multinomial only: wall-clock budget in seconds for one fit; exceeding it fails with an explicit error instead of holding the request open.', + 'analysis.differential.mc_samples': 'aldex2 only: Dirichlet Monte-Carlo instances per sample, 16 to 4096; more instances steady the effect sizes at linear cost.', 'analysis.alpha.show_points': 'Show individual samples overlaid on alpha comparison boxplots.', 'analysis.alpha.annotate_significance': 'Annotate alpha comparison panels with Kruskal-Wallis significance.', 'analysis.alpha.pairwise_brackets': 'Run BH-adjusted pairwise Wilcoxon rank-sum tests between every group pair and draw brackets, including n.s. results.', @@ -253,7 +260,7 @@ export const CONFIG_TYPES: Record = { 'analysis.alpha.rarefaction_iterations': { kind: 'int' }, 'analysis.beta.normalisation': { kind: 'enum', options: ['hellinger', 'none', 'rarefy', 'srs'] }, 'analysis.beta.depth': { kind: 'int' }, - 'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'] }, + 'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', 'multinomial', 'dirichlet_multinomial'] }, 'analysis.differential.offset': { kind: 'enum', options: ['tss', 'rle'] }, 'analysis.differential.dispersion': { kind: 'enum', options: ['ml', 'glmGamPoi'] }, 'analysis.differential.min_prevalence': { kind: 'float', min: 0, max: 1, step: 0.05 }, @@ -261,6 +268,11 @@ export const CONFIG_TYPES: Record = { 'analysis.differential.replacement_delta': { kind: 'float', min: 0, max: 1, step: 0.05 }, 'analysis.differential.bayesian_alpha': { kind: 'float', nullable: true, min: 0, step: 0.1 }, 'analysis.differential.mc_samples': { kind: 'int' }, + 'analysis.differential.mn_penalty': { kind: 'float', min: 0, step: 0.05 }, + 'analysis.differential.dm_overdispersion_method': { kind: 'enum', options: ['mle', 'mom'] }, + 'analysis.differential.max_iter': { kind: 'int' }, + 'analysis.differential.tolerance': { kind: 'float' }, + 'analysis.differential.max_seconds': { kind: 'float', min: 1 }, 'analysis.alpha.show_points': { kind: 'boolean' }, 'analysis.alpha.annotate_significance': { kind: 'boolean' }, 'analysis.alpha.pairwise_brackets': { kind: 'boolean' }, diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts index 1cced072..93e7f8e1 100644 --- a/frontend/src/components/differentialTable.test.ts +++ b/frontend/src/components/differentialTable.test.ts @@ -99,7 +99,7 @@ describe('multinomial results', () => { test('the reference taxon is exported with an empty statistic set, never a substitute', () => { const lines = differentialCsv(mn).trimEnd().split('\r\n') - expect(lines[4].startsWith('Bacteroides,reference,0,0,,,,,,1,')).toBe(true) + expect(lines[4].startsWith('Bacteroides,reference,0,0,,,,,,,1,')).toBe(true) expect(effectValue(mn, mn.rows[1])).toBe(1.6) }) }) @@ -163,7 +163,7 @@ describe('advanced compositional methods (ancom_bc, aldex2, songbird)', () => { expect(effectValue(ax, ax.rows[0])).toBe(1.85) expect(configLine(ax)).toBe('method aldex2, mc_samples 128, seed 123, min_prevalence 0.05; estimate is an ALDEx2 standardized CLR effect size, not a fold change') const csvRows = differentialCsv(ax).trimEnd().split('\r\n') - expect(csvRows[4]).toBe('Bacteroides,ok,,1.85,,5.545,2.9e-8,5.8e-8,,1,') + expect(csvRows[4]).toBe('Bacteroides,ok,,1.85,,5.545,2.9e-8,5.8e-8,,,1,') }) test('songbird reports CLR-centred multinomial differential metadata', () => { @@ -200,12 +200,12 @@ describe('calibration caveat', () => { const lines = differentialCsv(caveated).trimEnd().split('\r\n') expect(lines[2]).toBe(`# ${configLine(caveated)}`) expect(lines[3]).toBe('# Calibration caveat (strong): Wald p-values treat the per-taxon theta as known.') - expect(lines[4]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,prevalence,note') + expect(lines[4]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,zero_process_log_odds,prevalence,note') expect(lines).toHaveLength(7) }) test('a result without a caveat keeps the three-line header', () => { const lines = differentialCsv(result).trimEnd().split('\r\n') - expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,prevalence,note') + expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,zero_process_log_odds,prevalence,note') }) }) diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts index 3ed1feab..f9a1275d 100644 --- a/frontend/src/components/differentialTable.ts +++ b/frontend/src/components/differentialTable.ts @@ -50,17 +50,40 @@ export function formatStat(value: number | null): string { */ export function configLine(result: DifferentialResult): string { const c = result.config - if (c.method === 'nb_glm') - // Older saved results may not record the dispersion; name it only - // when the result carries it, never guess one for the record. - return `method nb_glm, offset ${c.offset}${c.dispersion ? `, dispersion ${c.dispersion}` : ''}, min_prevalence ${c.min_prevalence}` - if (c.method !== 'clr_welch') - return `method ${c.method}, offset ${c.offset}, min_prevalence ${c.min_prevalence}` - const replacement = c.replacement_method ?? 'bayesian_multiplicative' - const alpha = c.bayesian_alpha === null || c.bayesian_alpha === undefined - ? 'estimated' : String(c.bayesian_alpha) - return `method clr_welch, replacement_method ${replacement}, replacement_delta ${c.replacement_delta}, ` + - `bayesian_alpha ${alpha}, min_prevalence ${c.min_prevalence}; estimate is a CLR difference, not a fold change` + switch (c.method) { + case 'clr_welch': { + const replacement = c.replacement_method ?? 'bayesian_multiplicative' + // The Bayesian alpha belongs to the Bayesian replacement only; naming + // it for multiplicative replacement would credit the result with a + // prior it never used. + const alpha = replacement === 'bayesian_multiplicative' + ? `, bayesian_alpha ${c.bayesian_alpha === null || c.bayesian_alpha === undefined ? 'estimated' : String(c.bayesian_alpha)}` + : '' + return `method clr_welch, replacement_method ${replacement}, replacement_delta ${c.replacement_delta}${alpha}, ` + + `min_prevalence ${c.min_prevalence}; estimate is a CLR difference, not a fold change` + } + case 'ancom_bc': + return `method ancom_bc, min_prevalence ${c.min_prevalence}; effect is bias-corrected log2 fold change (Lin & Peddada 2020)` + case 'aldex2': + return `method aldex2, mc_samples ${c.mc_samples}, seed ${c.seed}, min_prevalence ${c.min_prevalence}; estimate is an ALDEx2 standardized CLR effect size, not a fold change` + case 'songbird': + return `method songbird, min_prevalence ${c.min_prevalence}; estimate is a CLR-centred multinomial differential, not an absolute fold change` + case 'multinomial': + case 'dirichlet_multinomial': { + const overdispersion = c.method === 'dirichlet_multinomial' + ? `, overdispersion ${c.dm_overdispersion_method}` + : '' + return `method ${c.method}${overdispersion}, reference ${c.mn_reference_taxon}, ` + + `penalty ${c.mn_penalty}, min_prevalence ${c.min_prevalence}; estimates are ` + + 'log2 fold changes against the reference taxon' + } + case 'nb_glm': + // Older saved results may not record the dispersion; name it only + // when the result carries it, never guess one for the record. + return `method nb_glm, offset ${c.offset}${c.dispersion ? `, dispersion ${c.dispersion}` : ''}, min_prevalence ${c.min_prevalence}` + default: + return `method ${c.method}, offset ${c.offset}, min_prevalence ${c.min_prevalence}` + } } /** The effect size of one row, from the field the result names in `effect`. */ diff --git a/frontend/src/views/StudyView.tsx b/frontend/src/views/StudyView.tsx index 89f32231..706528e6 100644 --- a/frontend/src/views/StudyView.tsx +++ b/frontend/src/views/StudyView.tsx @@ -18,7 +18,7 @@ import { EvidenceModeToggle } from '../components/EvidenceModeToggle' import { CladeCumulusPanel } from '../components/CladeCumulusPanel' import { ReportPanel } from '../components/ReportPanel' import { useNavRefresh } from '../hooks/useNavRefresh' -import { CladeCumulusPanel } from '../cladistics/CladeCumulusPanel' +import { CladeCumulusPanel as CladisticsPanel } from '../cladistics/CladeCumulusPanel' type GroupedRun = Run & { group?: string | null } @@ -300,7 +300,7 @@ export function StudyView() { source={analysisSource} /> :

No run has a results table with {analysisSource} taxonomy yet.

)} - {view === 'cladistics' && } + {view === 'cladistics' && } {view === 'trees' && <> diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 96451b6a..0371c675 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -44,7 +44,7 @@ include("canonical.jl") include("analysis/clade_cumulus.jl") include("analysis/diversity.jl") include("analysis/analysis.jl") -include("analysis/clade_cumulus.jl") +include("analysis/compositional.jl") include("analysis/differential.jl") # Opt-in validated descriptive statistics (issue #1). Ordinary workflows above # remain unchanged; these modes must be selected explicitly. diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index d7c1741b..f0fec757 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -43,8 +43,12 @@ export DifferentialConfig, ScalingRefusal, MASSUnavailable, ZCompositionsUnavail clr_transform, fit_nb, fit_two_part_nb, fit_occupancy, fit_welch, differential_abundance, volcano_chart -# Methods accepted under `analysis.differential.method`. -const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb") +# Methods accepted under `analysis.differential.method`. The multinomial +# models (fitted in `CompositionalModels`, which conditions on each sample's +# read total rather than testing taxa one at a time) are named separately so +# the config validation and the dispatch can share one definition. +const COMPOSITIONAL_METHODS = ("multinomial", "dirichlet_multinomial") +const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", "aldex2", "songbird", COMPOSITIONAL_METHODS...) # The per-taxon dispersion `nb_glm` uses: maximum likelihood (MASS::glm.nb's # own theta) or glmGamPoi's shrinkage estimator, which shares information @@ -74,6 +78,13 @@ const MAX_MC_SAMPLES = 4096 const DEFAULT_MC_SAMPLES = 128 const DEFAULT_DA_SEED = 123 +# Optimiser defaults for the two multinomial models, repeated here so the file +# documents every setting it validates; the values must match +# CompositionalModels.DEFAULT_* (asserted in the unit tests). +const COMPOSITIONAL_MAX_ITER = 500 +const COMPOSITIONAL_TOL = 1e-8 +const COMPOSITIONAL_MAX_SECONDS = 300.0 + # Canonicalise and validate a public replacement-method name in one place, so # direct library calls and YAML configuration have identical semantics. function _replacement_method(method::AbstractString) @@ -127,6 +138,12 @@ struct DifferentialConfig dispersion :: String replacement_method :: String bayesian_alpha :: Union{Float64,Nothing} + mn_reference_taxon :: Union{Nothing,String} + mn_penalty :: Float64 + dm_overdispersion_method :: String + max_iter :: Int + tol :: Float64 + max_seconds :: Float64 mc_samples :: Int function DifferentialConfig(offset::AbstractString="tss", min_prevalence::Real=0.0; method::AbstractString="nb_glm", @@ -134,6 +151,12 @@ struct DifferentialConfig dispersion::AbstractString="ml", replacement_method::AbstractString="bayesian_multiplicative", bayesian_alpha::Union{Real,Nothing}=nothing, + mn_reference_taxon::Union{Nothing,AbstractString}=nothing, + mn_penalty::Real=0.0, + dm_overdispersion_method::AbstractString="mle", + max_iter::Integer=COMPOSITIONAL_MAX_ITER, + tol::Real=COMPOSITIONAL_TOL, + max_seconds::Real=COMPOSITIONAL_MAX_SECONDS, mc_samples::Real=DEFAULT_MC_SAMPLES) m = lowercase(strip(String(method))) m in METHODS || throw(ArgumentError( @@ -190,7 +213,9 @@ struct DifferentialConfig "analysis.differential.mc_samples must be a whole number in " * "[$MIN_MC_SAMPLES, $MAX_MC_SAMPLES], not $mc_samples")) a = isnothing(bayesian_alpha) ? nothing : Float64(bayesian_alpha) - new(m, o, Float64(min_prevalence), Float64(replacement_delta), d, r, a, Int(mc_samples)) + new(m, o, Float64(min_prevalence), Float64(replacement_delta), d, r, a, ref, + Float64(mn_penalty), dm, Int(max_iter), Float64(tol), Float64(max_seconds), + Int(mc_samples)) end end @@ -2033,7 +2058,8 @@ end Test every taxon (column) of the samples x taxa read-count matrix for a difference in abundance between the samples labelled `contrast` and those labelled `reference` in `groups`, by `config.method` (`"nb_glm"`, `"clr_welch"`, -`"ancom_bc"`, `"aldex2"` or `"songbird"`). +`"zinb"`, `"hurdle_nb"`, `"ancom_bc"`, `"aldex2"`, `"songbird"`, `"multinomial"` +or `"dirichlet_multinomial"`). `nb_glm`: size factors are computed on the full matrix by `config.offset` and enter each negative-binomial model as `log(factor)`; the estimate is a @@ -2203,7 +2229,7 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, for (i, s) in enumerate(replaced.posterior.concentration)]) end diagnostics["zero_replacement"] = zero_diagnostics - else + elseif config.method in ("zinb", "hurdle_nb") factors = size_factors(counts, samples, config.offset) tested = findall(>=(config.min_prevalence), prevalence) fits = isempty(tested) ? nothing : @@ -2234,6 +2260,272 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, diagnostics["model_warning"] = config.method == "zinb" ? "Use ZINB only when diagnostics support excess zeros beyond NB; zero inflation and NB overdispersion can be difficult to distinguish." : "Hurdle NB assumes every zero comes from the occurrence process and positive counts follow a zero-truncated NB." + elseif config.method == "ancom_bc" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "ANCOM-BC needs at least 2 samples in each group to estimate within-group variance; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 3 || throw(ArgumentError( + "ancom_bc needs at least 3 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence) to separate " * + "sampling-fraction bias from taxon effects; there are $(length(tested))")) + bc = fit_ancom_bc(counts[:, tested], groups; levels, samples) + for (k, j) in enumerate(tested) + f = bc.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "log2_fold_change" => isnothing(f.estimate) ? nothing : f.estimate / log(2), + "standard_error" => f.se, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + method_text = "ANCOM-BC bias-corrected log-linear model (Lin & Peddada 2020) over the taxa " * + "passing min_prevalence: within-group sample offsets and between-group " * + "sampling-fraction bias estimated via an E-M Gaussian mixture on log(counts + 1), " * + "with bias variance propagated into two-sided Wald tests and Benjamini-Hochberg " * + "adjustment over the fitted taxa. Assumes fewer than half of the tested taxa " * + "are differentially abundant" + effect = Dict("key" => "log2_fold_change", "label" => "Bias-corrected log2 FC") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence) + diagnostics["ancom_bc"] = merge(bc.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method == "aldex2" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "ALDEx2 needs at least 2 samples in each group to estimate within-group variance; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 2 || throw(ArgumentError( + "ALDEx2 needs at least 2 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence); there are $(length(tested))")) + ax = fit_aldex2(counts[:, tested], groups; levels, samples, + mc_samples=config.mc_samples, seed=Int(seed)) + for (k, j) in enumerate(tested) + f = ax.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "standard_error" => nothing, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + method_text = "ALDEx2 Dirichlet(count + 0.5) Monte-Carlo CLR model (Fernandes et al. 2014; " * + "$(config.mc_samples) instances, seed = $(Int(seed))) over the taxa passing " * + "min_prevalence: consistent-sign posterior predictive Welch's t-test p-values " * + "with Benjamini-Hochberg adjustment over the fitted taxa. Estimates are " * + "standardized CLR effect sizes (between-group difference over max within-group " * + "dispersion), not absolute fold changes" + effect = Dict("key" => "estimate", "label" => "ALDEx2 effect size") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence, + "mc_samples" => config.mc_samples, + "seed" => Int(seed)) + diagnostics["aldex2"] = merge(ax.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method == "songbird" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "songbird needs at least 2 samples in each group; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 2 || throw(ArgumentError( + "songbird needs at least 2 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence); there are $(length(tested))")) + sb = fit_songbird(counts[:, tested], groups; levels, samples, taxa=taxa[tested]) + for (k, j) in enumerate(tested) + f = sb.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "standard_error" => f.se, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + ref_tax = something(sb.diagnostics["reference_taxon"], "none") + method_text = "Songbird-style compositional multinomial regression (Morton et al. 2019; " * + "Zhang et al. 2017; reference taxon '$ref_tax') over the identifiable taxa " * + "passing min_prevalence: CLR-centred differentials (summing to zero across " * + "fitted taxa) with delta-method standard errors, two-sided Wald z-tests, and " * + "Benjamini-Hochberg adjustment. Differentials are relative log-ratios, not " * + "absolute fold changes, and multinomial likelihood conditions on read totals " * + "without extra-multinomial overdispersion" + effect = Dict("key" => "estimate", "label" => "Differential (CLR log-ratio)") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence) + diagnostics["songbird"] = merge(sb.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method in COMPOSITIONAL_METHODS + # The multinomial models condition on each sample's total: they need no + # size factor, and a zero count is part of the likelihood rather than a + # value to replace. Two things still have to be settled before the fit. + # + # First, a taxon with no reads at all in one group has no finite log + # fold change, and because the model is fitted jointly its coefficient + # would drag the whole fit towards -Inf. Such taxa are refused by name + # and left out of the composition, exactly like taxa below the + # prevalence threshold; the note says which and why. + tested = findall(>=(config.min_prevalence), prevalence) + separated = Int[] + for j in tested + absent = [lvl for lvl in levels + if all(i -> counts[i, j] == 0, findall(==(lvl), groups))] + if !isempty(absent) + push!(separated, j) + rows[j]["status"] = "failed" + rows[j]["note"] = "no reads in any '$(absent[1])' sample (complete " * + "separation): the log fold change against the reference is unbounded, " * + "so the taxon is left out of the composition" + end + end + tested = [j for j in tested if !(j in separated)] + isempty(tested) && throw(ArgumentError( + "no taxon passes analysis.differential.min_prevalence = $(config.min_prevalence) " * + "and has reads in both '$reference' and '$contrast' samples, so there is no " * + "composition to model; raise analysis.differential.min_prevalence or use " * + "method nb_glm")) + + # Second, the reference taxon: the denominator of every reported ratio. + # An explicit choice must be a taxon at this rank and must be fitted. + ref_name = config.mn_reference_taxon + ref_pos = 0 + if isnothing(ref_name) + totals = [sum(view(counts, :, j)) for j in tested] + ref_pos = argmax(totals) + ref_name = String(taxa[tested[ref_pos]]) + else + j = findfirst(==(ref_name), taxa) + isnothing(j) && throw(ArgumentError( + "analysis.differential.mn_reference_taxon = '$(ref_name)' is not a taxon " * + "at this rank")) + pos = findfirst(==(j), tested) + isnothing(pos) && throw(ArgumentError( + "the reference taxon '$(ref_name)' is not fitted: it has reads in " * + "$(count(>(0), view(counts, :, j))) of $n samples, below " * + "analysis.differential.min_prevalence = $(config.min_prevalence), or no " * + "reads in one group. Choose another reference taxon, or lower min_prevalence")) + ref_pos = pos + end + # A reference seen in fewer than half the samples makes every ratio it + # anchors unstable, so it is refused rather than silently used. + ref_prevalence = prevalence[tested[ref_pos]] + ref_prevalence >= 0.5 || throw(ArgumentError( + "the reference taxon '$(ref_name)' has reads in only " * + "$(round(Int, ref_prevalence * n)) of $n samples; the log ratios it anchors " * + "would be unstable. Choose a reference present in at least half the samples " * + "(analysis.differential.mn_reference_taxon), or use method nb_glm or clr_welch")) + + # The design is the group indicator with an intercept: the intercept is + # the log-odds against the reference in the reference group, the slope + # the difference in those log-odds between the groups. + design = Matrix{Float64}(undef, n, 2) + for i in 1:n + design[i, 1] = 1.0 + design[i, 2] = groups[i] == contrast ? 1.0 : 0.0 + end + sub = Matrix{Int}(counts[:, tested]) + started = time() + fit = config.method == "multinomial" ? + fit_multinomial(sub, design; reference=ref_pos, penalty=config.mn_penalty, + max_iter=config.max_iter, tol=config.tol, + max_seconds=config.max_seconds) : + fit_dirichlet_multinomial(sub, design; reference=ref_pos, penalty=config.mn_penalty, + overdispersion_method=config.dm_overdispersion_method, + max_iter=config.max_iter, tol=config.tol, + max_seconds=config.max_seconds) + elapsed = time() - started + # A joint fit that has not converged describes neither taxon: there is + # no per-taxon rescue, so it is refused rather than reported. + fit.converged || throw(CompositionalRefusal( + "the $(config.method == "multinomial" ? "multinomial" : "Dirichlet-multinomial") " * + "fit did not converge in $(config.max_iter) Newton steps on $(length(tested)) taxa; " * + "raise analysis.differential.max_iter, raise min_prevalence so fewer taxa are " * + "fitted, or analyse a coarser taxonomic rank")) + c = fit.concentration + at_bound = config.method == "dirichlet_multinomial" && + !(c > CONCENTRATION_MIN * (1 + 1e-6) && c < CONCENTRATION_MAX * (1 - 1e-6)) + boundary_note = c >= CONCENTRATION_MAX * (1 - 1e-6) ? + "the Dirichlet concentration reached its upper bound: no overdispersion is " * + "detectable, so the model is effectively multinomial" : + "the Dirichlet concentration reached its lower bound: the composition is nearly " * + "degenerate and the standard errors are not trustworthy" + for (l, j) in enumerate(tested) + if l == ref_pos + rows[j]["status"] = "reference" + rows[j]["note"] = "reference taxon: the denominator of every log ratio; its " * + "coefficient is fixed at zero, so it carries no test" + rows[j]["estimate"] = 0.0 + rows[j]["log2_fold_change"] = 0.0 + continue + end + k = l < ref_pos ? l : l - 1 + est = fit.coefficients[k, 2] + se = fit.standard_errors[k, 2] + zz = fit.statistic[k, 2] + pv = fit.pvalue[k, 2] + if isnan(est) || isnan(se) || isnan(pv) + rows[j]["status"] = "failed" + rows[j]["note"] = "the fit returned a non-finite estimate or standard error" + continue + end + merge!(rows[j], Dict{String,Any}( + "status" => at_bound ? "boundary" : "ok", + "note" => at_bound ? boundary_note : "", + "estimate" => est, + "log2_fold_change" => est / log(2), + "standard_error" => se, "statistic" => zz, "pvalue" => pv, + "dispersion_theta" => nothing)) + end + model = config.method == "multinomial" ? "multinomial logistic" : + "Dirichlet-multinomial" + method_text = "$model regression over $(length(tested)) taxa with " * + "'$(ref_name)' as the reference taxon (softmax link, each sample's read " * + "composition, no size factor and no zero replacement); " * + (config.method == "multinomial" ? + "Wald test of the group coefficient from the observed information, with a " * + "ridge penalty of $(config.mn_penalty); " : + "concentration estimated by $(uppercase(config.dm_overdispersion_method)), " * + "theta = $(round(fit.overdispersion; sigdigits=3)), Wald test of the group " * + "coefficient from the observed information; ") * + "Benjamini-Hochberg adjustment over the fitted taxa. Estimates are log fold " * + "changes in relative abundance against '$(ref_name)' (a log ratio of two taxa, " * + "so coherent under compositionality): " * + (config.method == "multinomial" ? + "Songbird's model, Morton et al. 2019" : + "La Rosa et al. 2012 for the Dirichlet-multinomial") + effect = Dict("key" => "log2_fold_change", + "label" => "log2 fold change vs $(ref_name)") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence, + "mn_reference_taxon" => ref_name, + "mn_penalty" => config.mn_penalty, + "tolerance" => config.tol, + "max_iter" => config.max_iter, + "max_seconds" => config.max_seconds) + config.method == "dirichlet_multinomial" && + (config_echo["dm_overdispersion_method"] = config.dm_overdispersion_method) + diagnostics["model"] = Dict{String,Any}( + "reference_taxon" => ref_name, "penalty" => config.mn_penalty, + "n_parameters" => fit.n_parameters, "iterations" => fit.iterations, + "converged" => fit.converged, "log_likelihood" => fit.loglik, + "runtime_seconds" => elapsed, "n_taxa_fitted" => length(tested)) + if config.method == "dirichlet_multinomial" + diagnostics["overdispersion"] = Dict{String,Any}( + "method" => config.dm_overdispersion_method, + "concentration" => c, "theta" => fit.overdispersion, + "at_bound" => at_bound) + end + else + throw(ArgumentError("unsupported differential method '$(config.method)'")) end diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 2e74e3cf..85622170 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1214,6 +1214,7 @@ function _differential_config(study::String) mc_val = mc_raw isa Real ? mc_raw : throw(ArgumentError( "analysis.differential.mc_samples must be a whole number in " * "[$(Differential.MIN_MC_SAMPLES), $(Differential.MAX_MC_SAMPLES)], not '$mc_raw'")) + reference = get(get(resolved, "analysis.differential.mn_reference_taxon", (; value=nothing)), :value, nothing) DifferentialConfig(String(val("analysis.differential.offset", "tss")), Float64(val("analysis.differential.min_prevalence", 0.0)); method=String(val("analysis.differential.method", "nb_glm")), @@ -1222,6 +1223,15 @@ function _differential_config(study::String) replacement_method=String(val("analysis.differential.replacement_method", "bayesian_multiplicative")), bayesian_alpha=alpha_val, + mn_reference_taxon=isnothing(reference) ? nothing : string(reference), + mn_penalty=Float64(val("analysis.differential.mn_penalty", 0.0)), + dm_overdispersion_method=String(val("analysis.differential.dm_overdispersion_method", "mle")), + max_iter=Int(val("analysis.differential.max_iter", + Differential.COMPOSITIONAL_MAX_ITER)), + tol=Float64(val("analysis.differential.tolerance", + Differential.COMPOSITIONAL_TOL)), + max_seconds=Float64(val("analysis.differential.max_seconds", + Differential.COMPOSITIONAL_MAX_SECONDS)), mc_samples=mc_val) end diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 85749bf4..287c5ea7 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -245,11 +245,32 @@ function _differential_settings_error(effective::AbstractDict) current = get(effective, k, _DIFFERENTIAL_SETTING_DEFAULTS[k]) current isa NamedTuple ? current.value : current end - method = value(keys_[1], "nb_glm") - dispersion = value(keys_[2], "ml") - if method isa AbstractString && dispersion isa AbstractString && - lowercase(strip(method)) != "nb_glm" && lowercase(strip(dispersion)) == "glmgampoi" - return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; use dispersion: ml with $(strip(method))" + method = value("analysis.differential.method") + dispersion = value("analysis.differential.dispersion") + reference = value("analysis.differential.mn_reference_taxon") + penalty = value("analysis.differential.mn_penalty") + overdispersion = value("analysis.differential.dm_overdispersion_method") + method isa AbstractString || return nothing + m = lowercase(strip(method)) + compositional = m in Differential.COMPOSITIONAL_METHODS + if m in Differential.METHODS && m != "nb_glm" && dispersion isa AbstractString && + lowercase(strip(dispersion)) == "glmgampoi" + return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; " * + "use dispersion: ml with $m" + end + if !isnothing(reference) && !(reference isa AbstractString && isempty(strip(reference))) && + !compositional + return "analysis.differential.mn_reference_taxon requires method: multinomial or " * + "dirichlet_multinomial; method is '$method'" + end + if penalty isa Real && !(penalty isa Bool) && penalty != 0 && !compositional + return "analysis.differential.mn_penalty requires method: multinomial or " * + "dirichlet_multinomial; method is '$method'" + end + if overdispersion isa AbstractString && lowercase(strip(overdispersion)) != "mle" && + m != "dirichlet_multinomial" + return "analysis.differential.dm_overdispersion_method requires method: " * + "dirichlet_multinomial; method is '$method'" end nothing end diff --git a/src/server/server.jl b/src/server/server.jl index 3a2728b8..eca368a0 100644 --- a/src/server/server.jl +++ b/src/server/server.jl @@ -116,6 +116,7 @@ module Server include(joinpath(@__DIR__, "routes", "doi.jl")) include(joinpath(@__DIR__, "routes", "figures.jl")) include(joinpath(@__DIR__, "routes", "phylogeny.jl")) + include(joinpath(@__DIR__, "routes", "cladistics.jl")) include(joinpath(@__DIR__, "routes", "clade_cumulus.jl")) ## Bad request middleware diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 1b00f1f4..be15cfd4 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -168,7 +168,9 @@ end end @testset "Configuration" begin - @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb") + @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb", + "ancom_bc", "aldex2", "songbird", + "multinomial", "dirichlet_multinomial") @test DifferentialConfig().offset == "tss" @test DifferentialConfig().min_prevalence == 0.0 @test DifferentialConfig("rle", 0.5).offset == "rle" @@ -1232,7 +1234,8 @@ end # Both methods are selectable names, and the two are told apart from the # per-taxon models by one tuple. @test Differential.METHODS == - ("nb_glm", "clr_welch", "multinomial", "dirichlet_multinomial") + ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", "aldex2", + "songbird", "multinomial", "dirichlet_multinomial") @test Differential.COMPOSITIONAL_METHODS == ("multinomial", "dirichlet_multinomial") c = DifferentialConfig(; method="multinomial") @test c.method == "multinomial"