diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9bad3b7b..d51aa441 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -400,11 +400,11 @@ jobs: shell: bash run: | set +e - julia --project=. -t2 --compiled-modules=no test/runtests.jl numeric_policy exact_summaries high_precision derived_analyses validated_statistics 2>&1 | tee /tmp/validated-statistics-tests.log + julia --project=. -t2 --compiled-modules=no test/runtests.jl numeric_policy exact_summaries high_precision derived_analyses validated_statistics 2>&1 | tee /tmp/julia-validated-tests.log status="${PIPESTATUS[0]}" if [ "$status" -ne 0 ]; then enc() { head -c 3500 | sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; } - echo "::error title=Validated-statistics tests failed::$(tail -n 120 /tmp/validated-statistics-tests.log | enc)" + echo "::error title=Validated statistics failed (tail)::$(tail -c 3500 /tmp/julia-validated-tests.log | enc)" fi exit "$status" diff --git a/README.md b/README.md index 951a8459..6813264c 100644 --- a/README.md +++ b/README.md @@ -299,6 +299,7 @@ not exclude samples from this analysis. 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`. + The glmGamPoi path uses fixed-theta NB Wald fits, not a quasi-likelihood test. - `zinb`: a zero-inflated negative-binomial model with separate logistic zero and count components, a count-component library-size offset, a joint 2-df likelihood-ratio test, and BH correction. @@ -316,12 +317,15 @@ not exclude samples from this analysis. on `log(counts + 1)`, propagating bias-estimation variance into two-sided Wald tests. The effect is a bias-corrected log2 fold change; useful when asymmetric shifts or load changes would bias standard CLR centering (needs at least 3 - tested taxa; assumes fewer than half are differential). + tested taxa). The implementation reports its mixture and bias diagnostics. + Its assumption that fewer than half of tested taxa are differential must be + checked rather than treated as a guarantee. - `aldex2`: Dirichlet(`count + 0.5`) Monte-Carlo CLR model (Fernandes et al. 2013, 2014). Draws `mc_samples` (default 128, seeded by the study's `seed`) instances per sample and reports standardized CLR effect sizes with consistent-sign posterior predictive Welch p-values; suited to sparse or variable-depth cohorts where deterministic zero imputation is undesirable. + The reported estimate is a standardized CLR effect size, not an absolute fold change. - `songbird`: compositional baseline-category multinomial regression with reference-frame invariant CLR-centred differentials summing to zero across fitted taxa (Morton et al. 2019, *Nat. Commun.* 10:2719; Zhang et al. 2017). @@ -336,7 +340,8 @@ not exclude samples from this analysis. composition (La Rosa et al. 2012), which absorbs between-sample overdispersion; `dm_overdispersion_method` picks `mle` (profiled) or `mom`. Both multinomial models bound their Newton fit with `max_iter`, `tolerance` - and `max_seconds`, and a fit that does not converge is refused. + and `max_seconds`, and a fit that does not converge is refused. Dirichlet- + multinomial diagnostics report the concentration, theta, and boundary status. ZINB and hurdle models can overfit sparse data; do not select ZINB solely because counts are sparse. Compare against NB diagnostics. The methods report @@ -818,23 +823,22 @@ 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 | zinb | hurdle_nb | ancom_bc - # | aldex2 | songbird | multinomial | dirichlet_multinomial + method: nb_glm # nb_glm | clr_welch | multinomial | dirichlet_multinomial + # ancom_bc | aldex2 | songbird | zinb | hurdle_nb offset: tss # count models: tss | rle dispersion: ml # nb_glm: ml | glmGamPoi min_prevalence: 0.0 # fraction of samples a taxon must be seen in replacement_method: bayesian_multiplicative # clr_welch: bayesian_multiplicative | multiplicative_replacement replacement_delta: 0.65 # clr_welch: strictly between 0 and 1 - bayesian_alpha: ~ # clr_welch: null (GBM estimate) or a number > 0 - mc_samples: 128 # aldex2: Dirichlet Monte-Carlo instances in [16, 4096] - mn_reference_taxon: ~ # multinomial models: null picks the most abundant taxon - mn_penalty: 0.0 # multinomial models: L2 ridge penalty, >= 0 - dm_overdispersion_method: mle - # dirichlet_multinomial: mle | mom + bayesian_alpha: ~ # Bayesian replacement: null (GBM estimate) or a number > 0 + mn_reference_taxon: ~ # multinomial models: null selects the most abundant fitted taxon + mn_penalty: 0.0 # multinomial models: non-negative L2 penalty + dm_overdispersion_method: mle # dirichlet_multinomial: mle | mom max_iter: 500 # multinomial models: Newton iteration cap - tolerance: 1.0e-8 # multinomial models: convergence tolerance + tolerance: 0.00000001 # multinomial models: convergence tolerance max_seconds: 300.0 # multinomial models: wall-clock budget + mc_samples: 128 # aldex2: Dirichlet Monte-Carlo instances in [16, 4096] nmds: max_stress: 0.2 # warn if NMDS stress exceeds this value ``` diff --git a/ROADMAP.md b/ROADMAP.md index 4b4beb7f..cab2edcf 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -6,7 +6,7 @@ occupancy work. The GitHub Wiki roadmap still needs a synchronized update; see ## Available now -- The differential-abundance config offers four per-taxon count models +- The differential-abundance config offers four per-taxon methods (`nb_glm`, `clr_welch`, `zinb`, `hurdle_nb`) and five compositional models (`ancom_bc`, `aldex2`, `songbird`, `multinomial`, `dirichlet_multinomial`); see [the differential-abundance section](README.md#differential-abundance). diff --git a/docs/occupancy-models.md b/docs/occupancy-models.md index 46a28661..abcef605 100644 --- a/docs/occupancy-models.md +++ b/docs/occupancy-models.md @@ -6,9 +6,9 @@ not an AnalysisConfig method or the standard endpoint. The full workflow is tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205), a follow-up to #18. -The differential-abundance configuration supports `zinb` and `hurdle_nb` in -addition to the existing `nb_glm` and `clr_welch` models. The two new methods -are fixed-effect, per-taxon models with a two-part negative-binomial likelihood. +The differential-abundance configuration supports `zinb` and `hurdle_nb` +alongside the count, CLR, and compositional methods. These two methods are +fixed-effect, per-taxon models with a two-part negative-binomial likelihood. They use the configured `tss` or `rle` library-size offset in the count component and apply Benjamini–Hochberg correction to the joint tests. diff --git a/frontend/src/api/client.ts b/frontend/src/api/client.ts index be36a201..cb9a9ccb 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -21,9 +21,9 @@ import type { TreeFileInfo, TreeFile, RunRef, ReadFunnelData, ReportItem, ReportKind, PlacementDoc, PlacementSummary, PlacementDetail, PlacementQueries, PhyloStep, ReferenceTreeDoc, ReferenceTreeSummary, ReferenceTreeDetail, AlignmentQC, TrimSettings, + CladisticsCreateRequest, CladisticsSummary, CladisticsDetail, CladeCumulusSourcesResponse, CladeCumulusRequest, CladeCumulusResponse, CladeMoveValidationRequest, CladeMoveValidationResponse, - CladisticsSummary, CladisticsDetail, CladisticsCreateRequest, } from './types' import type { FigureDoc } from '../figure/types' diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx index 545aa875..69cd1692 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -220,6 +220,16 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { `replaced values hold ${(100 * z.median_imputed_fraction).toPrecision(2)}% of the median sample ` + `and at most ${(100 * z.max_imputed_fraction).toPrecision(2)}% of any, and are not measurements.` })()} + {result.diagnostics.model && (() => { + const m = result.diagnostics.model + return ` Reference taxon: ${m.reference_taxon}; ${m.n_parameters} coefficients, ` + + `penalty ${m.penalty}, ${m.iterations} Newton steps in ${m.runtime_seconds.toPrecision(2)} s.` + })()} + {result.diagnostics.overdispersion && (() => { + const o = result.diagnostics.overdispersion + return ` Concentration ${o.concentration.toPrecision(3)} by ${o.method}, ` + + `theta = ${o.theta.toPrecision(3)}${o.at_bound ? ' (at a bound: see the row notes)' : ''}.` + })()} {result.diagnostics.ancom_bc && (() => { const b = result.diagnostics.ancom_bc const warn = b.warnings.length > 0 ? ` Warning: ${b.warnings.join('; ')}.` : '' diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts index 5496b941..a663d0f5 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -2,6 +2,11 @@ import { expect, test } from 'bun:test' import { CONFIG_DESCRIPTIONS, CONFIG_TYPES } from './PipelineStages' +const ALL_DIFFERENTIAL_METHODS = [ + 'nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', + 'multinomial', 'dirichlet_multinomial', +] + test('dispersion uses the enum editor with both supported models', () => { expect(CONFIG_TYPES['analysis.differential.dispersion']).toEqual({ kind: 'enum', options: ['ml', 'glmGamPoi'], @@ -10,14 +15,23 @@ test('dispersion uses the enum editor with both supported models', () => { test('differential model editor exposes all nine supported methods', () => { expect(CONFIG_TYPES['analysis.differential.method']).toEqual({ - kind: 'enum', - options: [ - 'nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', - 'multinomial', 'dirichlet_multinomial', - ], + kind: 'enum', options: ALL_DIFFERENTIAL_METHODS, }) }) +test('multinomial settings retain their enums and optimiser budgets', () => { + expect(CONFIG_TYPES['analysis.differential.dm_overdispersion_method']).toEqual({ + kind: 'enum', options: ['mle', 'mom'], + }) + expect(CONFIG_TYPES['analysis.differential.mn_penalty']).toEqual({ + kind: 'float', min: 0, step: 0.05, + }) + expect(CONFIG_TYPES['analysis.differential.max_seconds']).toEqual({ + kind: 'float', min: 1, + }) + expect(CONFIG_TYPES['analysis.differential.max_iter']).toEqual({ kind: 'int' }) +}) + test('advanced zero handling exposes both methods, bounded delta and nullable alpha', () => { expect(CONFIG_TYPES['analysis.differential.replacement_method']).toEqual({ kind: 'enum', options: ['bayesian_multiplicative', 'multiplicative_replacement'], diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index 9d84a40e..ddeba734 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -270,12 +270,12 @@ export const CONFIG_TYPES: Record = { 'analysis.differential.replacement_method': { kind: 'enum', options: ['bayesian_multiplicative', 'multiplicative_replacement'] }, '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', min: 0 }, 'analysis.differential.max_seconds': { kind: 'float', min: 1 }, + 'analysis.differential.mc_samples': { kind: 'int' }, 'analysis.alpha.show_points': { kind: 'boolean' }, 'analysis.alpha.annotate_significance': { kind: 'boolean' }, 'analysis.alpha.pairwise_brackets': { kind: 'boolean' }, diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index aaec5ba7..b9bec61b 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -12,8 +12,6 @@ # log-ratios, over the taxa that pass `min_prevalence`, after their zeros are # replaced by zCompositions::cmultRepl (Bayesian GBM/user or multiplicative # CZM, for counts). The estimate is a difference in CLR units, not a fold change. -# - zinb / hurdle_nb: two-part negative-binomial models with a logistic zero or -# occurrence process and a joint 2-df test of both group effects. # - ancom_bc: bias-corrected log-linear model (Lin & Peddada 2020) estimating # within-group sample-specific offsets and between-group sampling-fraction # bias via a safeguarded SQUAREM 3-component Gaussian mixture. The effect is @@ -25,8 +23,11 @@ # reference-frame invariant CLR-centred differentials (Morton et al. 2019; # Zhang et al. 2017) fitted by damped Newton-Raphson with exact O(K) # Sherman-Morrison Fisher-information inversion. -# - multinomial / dirichlet_multinomial: joint compositional regression against -# a reference taxon, with optional Dirichlet-multinomial overdispersion. +# - multinomial and dirichlet_multinomial: joint compositional models fitted in +# Julia against a reference taxon; the latter estimates extra-multinomial +# overdispersion. Neither uses a library-size offset or zero replacement. +# - zinb and hurdle_nb: two-part negative-binomial models with a joint test of +# group effects in the zero/occurrence and count processes. # # Count matrices here are samples x taxa, as everywhere else in `Analysis`. # Anything that is not a non-negative integer count is refused rather than @@ -41,16 +42,19 @@ using ..CompositionalModels: CompositionalFit, CompositionalRefusal, RuntimeRefu DM_METHODS, CONCENTRATION_MIN, CONCENTRATION_MAX, fit_multinomial, fit_dirichlet_multinomial export DifferentialConfig, ScalingRefusal, MASSUnavailable, ZCompositionsUnavailable, - GlmGamPoiUnavailable, RStatsUnavailable, CompositionalRefusal, RuntimeRefusal, - OFFSET_METHODS, METHODS, DISPERSIONS, REPLACEMENT_METHODS, COMPOSITIONAL_METHODS, + GlmGamPoiUnavailable, CompositionalRefusal, RuntimeRefusal, RStatsUnavailable, + OFFSET_METHODS, METHODS, COMPOSITIONAL_METHODS, DISPERSIONS, REPLACEMENT_METHODS, MIN_MC_SAMPLES, MAX_MC_SAMPLES, DEFAULT_MC_SAMPLES, tss_factors, rle_factors, size_factors, bh_adjust, validate_counts, replace_zeros, clr_transform, fit_nb, fit_two_part_nb, fit_occupancy, fit_welch, - differential_abundance, volcano_chart + fit_ancom_bc, fit_aldex2, fit_songbird, differential_abundance, volcano_chart -# Methods accepted under `analysis.differential.method`. +# The two joint models use a reference taxon and the optimiser settings below. +# Keep the other method groups in the public tuple as well: a configuration +# must never silently lose a model because two feature branches were merged. +const COMPOSITIONAL_METHODS = ("multinomial", "dirichlet_multinomial") const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", - "aldex2", "songbird", "multinomial", "dirichlet_multinomial") + "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,10 +78,10 @@ const THETA_LOWER = 1e-8 # zCompositions::cmultRepl's default `frac`, as in Martín-Fernández et al. (2003). const DEFAULT_DELTA = 0.65 -# Joint multinomial/Dirichlet-multinomial models and optimiser bounds. -const COMPOSITIONAL_METHODS = ("multinomial", "dirichlet_multinomial") -const COMPOSITIONAL_MAX_ITER = 500 -const COMPOSITIONAL_TOL = 1.0e-8 +# Optimiser defaults for the two multinomial models. These 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 # Bounds and default for `analysis.differential.mc_samples` (`aldex2` Dirichlet draws). @@ -100,23 +104,48 @@ end DifferentialConfig(offset="tss", min_prevalence=0.0; method="nb_glm", replacement_delta=0.65, dispersion="ml", replacement_method="bayesian_multiplicative", - bayesian_alpha=nothing, mc_samples=128, - mn_reference_taxon=nothing, mn_penalty=0.0, - dm_overdispersion_method="mle", max_iter=500, - tol=1e-8, max_seconds=300.0) - -Settings for one differential abundance analysis. `METHODS` contains the -per-taxon count models (`nb_glm`, `zinb`, `hurdle_nb`), the zero-replaced CLR -model (`clr_welch`), the advanced compositional models (`ancom_bc`, `aldex2`, -`songbird`) and the joint multinomial models in `COMPOSITIONAL_METHODS`. - -`offset` is used by the count models. Zero-replacement settings belong to -`clr_welch`; `mc_samples` belongs to `aldex2`; and the reference, penalty, -concentration and optimiser settings belong to the two joint multinomial -models. Settings incompatible with the selected method are rejected rather -than silently ignored. Method names are stripped and lowercased, dispersion -and overdispersion names are canonicalised, and every numeric setting is -validated before a fit starts. + bayesian_alpha=nothing, mn_reference_taxon=nothing, + mn_penalty=0.0, dm_overdispersion_method="mle", + max_iter=500, tol=1e-8, max_seconds=300.0, mc_samples=128) + +Settings for one differential abundance analysis, read from +`analysis.differential` in the pipeline config. `method` is one of `"nb_glm"`, +`"clr_welch"`, `"multinomial"`, `"dirichlet_multinomial"`, `"ancom_bc"`, +`"aldex2"`, `"songbird"`, `"zinb"` or `"hurdle_nb"`. `offset` names the +size factor method (`"tss"` or `"rle"`) and is used by the count models. +`min_prevalence` is the fraction of samples, in [0, 1], in which a taxon must +have at least one read to be tested. `replacement_delta`, in (0, 1), is the +`frac` passed to `zCompositions::cmultRepl` by `clr_welch`. +`replacement_method` selects `"bayesian_multiplicative"` (GBM) or +`"multiplicative_replacement"` (CZM), with optional positive `bayesian_alpha` +for `"bayesian_multiplicative"`. For the two multinomial methods, +`mn_reference_taxon=nothing` or blank text selects the most abundant fitted +taxon; other names are stripped of surrounding whitespace. `mn_penalty` is a +finite, non-negative ridge penalty. `dm_overdispersion_method` selects `"mle"` +or `"mom"` for the Dirichlet-multinomial model. `max_iter` must be at least 1; +`tol` is a positive, finite convergence tolerance and `max_seconds` a positive, +finite runtime budget in seconds. `mc_samples` is the ALDEx2 Monte-Carlo +instance count, a whole number in [16, 4096]; booleans are refused. +`dispersion` names the per-taxon dispersion `nb_glm` uses (`"ml"` or +`"glmGamPoi"`): `"ml"` takes the maximum-likelihood theta from `MASS::glm.nb`, +`"glmGamPoi"` uses `overdispersion_shrinkage_list.dispersion_trend` from +`glmGamPoi::glm_gp` for the taxa passing the prevalence filter, with `theta = 1 / dispersion` fixed +in each negative-binomial GLM. Both options use a Wald test of the group +coefficient with fixed scale 1; `"glmGamPoi"` is refused outside `nb_glm`. `zinb` and `hurdle_nb` +fit a zero-process and negative-binomial count process per taxon, then use a 2-df likelihood-ratio +test for the group effect in both processes. Their count-component fold change is reported; the zero +process effect remains in diagnostics. +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. +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, invalid numeric +settings, or incompatible method choices: an explicit reference or non-zero +penalty requires a multinomial method, and `"mom"` requires +`"dirichlet_multinomial"`. Numeric bounds are checked even for settings unused +by the selected method. Errors converting validated values to the stored +numeric types propagate, including `InexactError` for `max_iter` outside the +range of `Int`. """ struct DifferentialConfig method :: String @@ -126,26 +155,26 @@ struct DifferentialConfig dispersion :: String replacement_method :: String bayesian_alpha :: Union{Float64,Nothing} - mc_samples :: Int 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", replacement_delta::Real=DEFAULT_DELTA, dispersion::AbstractString="ml", replacement_method::AbstractString="bayesian_multiplicative", bayesian_alpha::Union{Real,Nothing}=nothing, - mc_samples::Real=DEFAULT_MC_SAMPLES, 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) + max_seconds::Real=COMPOSITIONAL_MAX_SECONDS, + mc_samples::Real=DEFAULT_MC_SAMPLES) m = lowercase(strip(String(method))) m in METHODS || throw(ArgumentError( "analysis.differential.method must be one of $(join(METHODS, ", ")), not '$method'")) @@ -158,8 +187,8 @@ struct DifferentialConfig "analysis.differential.dispersion must be one of $(join(DISPERSIONS, ", ")), not '$dispersion'")) d = DISPERSIONS[di] (d == "ml" || m == "nb_glm") || throw(ArgumentError( - "analysis.differential.dispersion is a negative-binomial GLM setting, but method is " * - "'$m'; only nb_glm uses the glmGamPoi dispersion trend")) + "analysis.differential.dispersion is a negative-binomial setting, but method is " * + "'$m'; only nb_glm has a dispersion to choose")) (isfinite(min_prevalence) && 0 <= min_prevalence <= 1) || throw(ArgumentError( "analysis.differential.min_prevalence must be a fraction in [0, 1], not $min_prevalence")) (isfinite(replacement_delta) && 0 < replacement_delta < 1) || throw(ArgumentError( @@ -178,10 +207,10 @@ struct DifferentialConfig compositional = m in COMPOSITIONAL_METHODS (isnothing(ref) || compositional) || throw(ArgumentError( "analysis.differential.mn_reference_taxon applies to method multinomial or " * - "dirichlet_multinomial, but method is '$m'")) + "dirichlet_multinomial, but method is '$m'; only those two models have a reference taxon")) (mn_penalty == 0 || compositional) || throw(ArgumentError( "analysis.differential.mn_penalty applies to method multinomial or " * - "dirichlet_multinomial, but method is '$m'")) + "dirichlet_multinomial, but method is '$m'; only those two models are penalised")) (dm == "mle" || m == "dirichlet_multinomial") || throw(ArgumentError( "analysis.differential.dm_overdispersion_method is a Dirichlet-multinomial setting, " * "but method is '$m'; set it to mle or select method dirichlet_multinomial")) @@ -201,9 +230,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), ref, Float64(mn_penalty), dm, Int(max_iter), - Float64(tol), Float64(max_seconds)) + 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 @@ -2045,8 +2074,9 @@ 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 the method selected in `config.method` -(one of `METHODS`). +labelled `reference` in `groups`, by `config.method` (`"nb_glm"`, `"clr_welch"`, +`"multinomial"`, `"dirichlet_multinomial"`, `"ancom_bc"`, `"aldex2"`, +`"songbird"`, `"zinb"` or `"hurdle_nb"`). `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 @@ -2071,6 +2101,17 @@ reference, relative to the geometric mean of the tested taxa. Replaced zeros are not measurements; `diagnostics.zero_replacement` reports their share of each sample. Each group needs at least 2 samples. +`multinomial` and `dirichlet_multinomial` fit one joint composition against a +reference taxon, without size factors or zero replacement. Their estimates are +log ratios against that reference; the Dirichlet-multinomial additionally +reports its concentration and overdispersion diagnostics. A reference that is +not fitted or is present in fewer than half the samples is refused. + +`ancom_bc`, `aldex2` and `songbird` retain their method-specific compositional +estimators and diagnostics. `ancom_bc` reports bias-corrected log2 fold changes, +`aldex2` standardized CLR effect sizes, and `songbird` CLR-centred differentials; +these are not interchangeable with an absolute fold change. + `zinb` and `hurdle_nb` use a count-process negative binomial with a library-size offset and a group-specific logistic zero/occurrence process. Their reported count effect is the NB coefficient; the p-value is a joint 2-df likelihood-ratio @@ -2080,7 +2121,7 @@ uses a zero-truncated NB for positive counts. These models can overfit; compare against NB diagnostics and do not interpret the zero component as biological absence without validation. -Across methods, taxa present in fewer than `config.min_prevalence` of the +For both methods, taxa present in fewer than `config.min_prevalence` of the samples are reported as `filtered` and not fitted, and the Benjamini-Hochberg family is the fitted taxa that produced a p-value; failed, filtered and reference taxa carry `padj = nothing`. The result's `effect` names the row key diff --git a/src/core/doi.jl b/src/core/doi.jl index f5b0bd95..6a2f926a 100644 --- a/src/core/doi.jl +++ b/src/core/doi.jl @@ -16,8 +16,6 @@ # * Durable write-ahead journal; immutable DOI/provenance receipts # * DOI insertion never rewrites scientific objects -module DOI - using Dates, SHA, HTTP, JSON3, Logging export DOIConfig, Deposition, PublicationReceipt, diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index aabc8c0d..b8284432 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1200,22 +1200,38 @@ end The study's `analysis.differential` settings, resolved through the config cascade. Missing or `nothing` values use the defaults in `DifferentialConfig`, including `dispersion="ml"`, GBM replacement and an estimated Bayesian alpha. -Invalid settings propagate the constructor's -`ArgumentError`; values that cannot be converted to strings or numbers also -raise their conversion errors. +Invalid settings propagate the constructor's `ArgumentError`; invalid types +for alpha, the reference taxon or the Monte-Carlo count also raise +`ArgumentError`. The iteration limit must be a positive whole number +representable as `Int`; an invalid limit or an `InexactError` or `OverflowError` +converting it raises `ArgumentError`. Other string and numeric conversion +errors, and errors reading or parsing the configuration, propagate. """ function _differential_config(study::String) resolved = _resolved_run_config(study, nothing, nothing, nothing) - val(k, d) = something(get(get(resolved, k, (; value=d)), :value, d), d) + val(k, d) = begin + value = get(get(resolved, k, (; value=d)), :value, d) + isnothing(value) ? d : value + end + # Null is meaningful for alpha: it selects GBM's estimated concentration. alpha = get(get(resolved, "analysis.differential.bayesian_alpha", (; value=nothing)), :value, nothing) alpha_val = isnothing(alpha) || alpha isa Real ? alpha : throw(ArgumentError( "analysis.differential.bayesian_alpha must be empty or a finite number greater than 0, not '$alpha'")) + reference = val("analysis.differential.mn_reference_taxon", nothing) mc_raw = val("analysis.differential.mc_samples", Differential.DEFAULT_MC_SAMPLES) 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) + max_iter_raw = val("analysis.differential.max_iter", 500) + max_iter_error = _config_value_error("analysis.differential.max_iter", max_iter_raw) + isnothing(max_iter_error) || throw(ArgumentError(max_iter_error)) + max_iter_val = try + Int(max_iter_raw) + catch e + (e isa InexactError || e isa OverflowError) || rethrow() + throw(ArgumentError( + "analysis.differential.max_iter must be a positive whole number representable as Int, not '$max_iter_raw'")) + end reference_val = isnothing(reference) || reference isa AbstractString ? reference : throw(ArgumentError("analysis.differential.mn_reference_taxon must be text or empty")) DifferentialConfig(String(val("analysis.differential.offset", "tss")), @@ -1231,7 +1247,7 @@ function _differential_config(study::String) 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", 500)), + max_iter=max_iter_val, tol=Float64(val("analysis.differential.tolerance", 1e-8)), max_seconds=Float64(val("analysis.differential.max_seconds", 300.0))) end diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 287c5ea7..d83b07af 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -140,8 +140,17 @@ const _CASE_FOLDED_ENUMS = Dict( """ _config_value_error(key, value) -> Union{String,Nothing} -Why a value cannot be saved under this key, or nothing. The type follows the -factory default; a key whose default is ~ accepts a string, number or boolean. +Return a reason a value cannot be saved under the dotted config `key`, or +`nothing` if it passes the checks for that key. Types follow factory defaults +unless the key has explicit enum, nullable, numeric or naming constraints; +keys without either a constraint or a typed default pass unchecked. + +Differential optimiser settings reject booleans and non-finite values: +`max_iter` must be a positive whole number representable as `Int`, tolerance +and runtime in seconds must be positive, and the ridge penalty non-negative. +`InexactError` and `OverflowError` from the iteration-limit conversion become +validation messages. This checks individual values without saving them; +method compatibility is checked separately by `_differential_settings_error`. """ function _config_value_error(key::String, value) v = value isa JSON3.Array ? collect(value) : value @@ -164,6 +173,29 @@ function _config_value_error(key::String, value) return (v isa Real && !(v isa Bool) && isfinite(v) && v > 0) ? nothing : "$key must be empty or a finite number greater than 0" end + if key == "analysis.differential.max_iter" + valid = v isa Real && !(v isa Bool) && isfinite(v) && isinteger(v) && v >= 1 + if valid + # An integral Float64 can still be outside Int's representable + # range (and typemax(Int) itself rounds when converted to Float64). + valid = try + Int(v) + true + catch e + (e isa InexactError || e isa OverflowError) || rethrow() + false + end + end + return valid ? nothing : "$key must be a positive whole number representable as Int" + end + if key in ("analysis.differential.tolerance", "analysis.differential.max_seconds") + return (v isa Real && !(v isa Bool) && isfinite(v) && v > 0) ? nothing : + "$key must be a finite number greater than 0" + end + if key == "analysis.differential.mn_penalty" + return (v isa Real && !(v isa Bool) && isfinite(v) && v >= 0) ? nothing : + "$key must be a finite number greater than or equal to 0" + end if key == "analysis.differential.mc_samples" return (v isa Real && !(v isa Bool) && isfinite(v) && isinteger(v) && Differential.MIN_MC_SAMPLES <= v <= Differential.MAX_MC_SAMPLES) ? nothing : diff --git a/src/server/server.jl b/src/server/server.jl index eca368a0..ebf5abb6 100644 --- a/src/server/server.jl +++ b/src/server/server.jl @@ -118,6 +118,7 @@ module Server include(joinpath(@__DIR__, "routes", "phylogeny.jl")) include(joinpath(@__DIR__, "routes", "cladistics.jl")) include(joinpath(@__DIR__, "routes", "clade_cumulus.jl")) + include(joinpath(@__DIR__, "routes", "cladistics.jl")) ## Bad request middleware # Helpers deep in a route (e.g. `_req_group`) throw this to reject a request. diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index b8543eff..d0f44e47 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -558,6 +558,25 @@ end for bad_mc in (0, 15, 4097, 64.5, true, "128", nothing) @test !isnothing(SV._config_value_error("analysis.differential.mc_samples", bad_mc)) end + for (key, value) in (("analysis.differential.max_iter", 1), + ("analysis.differential.max_iter", 500), + ("analysis.differential.tolerance", 1e-8), + ("analysis.differential.max_seconds", 300.0), + ("analysis.differential.mn_penalty", 0.0)) + @test isnothing(SV._config_value_error(key, value)) + end + for (key, value) in (("analysis.differential.max_iter", 0), + ("analysis.differential.max_iter", 1.5), + ("analysis.differential.max_iter", NaN), + ("analysis.differential.max_iter", big(2)^80), + ("analysis.differential.tolerance", 0), + ("analysis.differential.tolerance", Inf), + ("analysis.differential.max_seconds", -1), + ("analysis.differential.max_seconds", NaN), + ("analysis.differential.mn_penalty", -0.1), + ("analysis.differential.mn_penalty", Inf)) + @test !isnothing(SV._config_value_error(key, value)) + end for dispersion in ("ml", "glmGamPoi", " ML ", " GLMGAMPOI ") @test isnothing(SV._config_value_error("analysis.differential.dispersion", dispersion)) end @@ -609,6 +628,10 @@ end for (key, value) in (("analysis.differential.replacement_delta", 0), ("analysis.differential.replacement_delta", 1), ("analysis.differential.bayesian_alpha", -1), + ("analysis.differential.max_iter", 1.5), + ("analysis.differential.tolerance", 0), + ("analysis.differential.max_seconds", -1), + ("analysis.differential.mn_penalty", -0.1), ("analysis.differential.replacement_method", "pseudocount")) response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", url, ["Content-Type" => "application/json"], JSON3.write(Dict(key => value)))) @@ -1726,6 +1749,18 @@ end "table" => "merged", "rank" => "Genus")) @test r.status == 400 @test JSON3.read(String(r.body)).error == "invalid_config" + + # Invalid optimiser settings are a client configuration error, + # including conversion failures, never an uncaught HTTP 500. + for (key, value) in (("max_iter", "1.5"), ("tolerance", "0"), + ("max_seconds", "-1"), ("mn_penalty", "-0.1")) + write(joinpath(root, "data", "studyD", "pipeline.yml"), + "analysis:\n differential:\n method: clr_welch\n $key: $value\n") + r = _da_post("studyD", Dict("runs" => [Dict("run" => "runA"), Dict("run" => "runB")], + "table" => "merged", "rank" => "Genus")) + @test r.status == 400 + @test JSON3.read(String(r.body)).error == "invalid_config" + end end end @@ -1769,7 +1804,9 @@ end # BH adjustment is not optional: every fitted row carries a padj, and # no row carries one without a p-value. @test all(x -> (x.padj !== nothing) == (x.pvalue !== nothing), body.rows) - @test all(x -> x.padj !== nothing, body.rows) + fitted_rows = [x for x in body.rows if x.status in ("ok", "boundary")] + @test !isempty(fitted_rows) + @test all(x -> x.padj !== nothing, fitted_rows) @test haskey(body, :figure) && length(body.figure.data) >= 1 # An explicit reference reparameterises the same fit: the log ratio # between Up and Down (a difference of two coefficients) must not