From 20a97993ef458bec173570bace536586a0d0d5c9 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:15:28 +0000 Subject: [PATCH 1/8] Repair differential method union and frontend regressions Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- README.md | 46 ++- ROADMAP.md | 3 +- config/defaults/pipeline.yml | 14 +- docs/occupancy-models.md | 6 +- frontend/src/api/client.ts | 1 + frontend/src/api/types.ts | 29 +- frontend/src/components/DifferentialPanel.tsx | 10 + .../src/components/PipelineStages.test.ts | 28 +- frontend/src/components/PipelineStages.tsx | 20 +- .../src/components/differentialTable.test.ts | 4 +- frontend/src/components/differentialTable.ts | 38 +- frontend/src/views/StudyView.tsx | 8 +- src/analysis/differential.jl | 340 +++++++++++++++++- src/core/doi.jl | 2 - src/server/routes/analysis.jl | 15 +- src/server/routes/config.jl | 31 +- test/unit/test_differential.jl | 4 +- 17 files changed, 532 insertions(+), 67 deletions(-) diff --git a/README.md b/README.md index 429d3778..964a1831 100644 --- a/README.md +++ b/README.md @@ -293,18 +293,40 @@ not exclude samples from this analysis. 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`; `dispersion` selects maximum-likelihood theta or a - cross-taxon `glmGamPoi` trend. -- `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. -- `hurdle_nb`: a logistic occurrence model plus zero-truncated NB counts for - positive observations. It uses the same joint test and offset convention. + cross-taxon `glmGamPoi` trend. The glmGamPoi path uses fixed-theta NB Wald + fits, not a quasi-likelihood test. - `clr_welch`: Welch's t-test on centred log-ratios after zero replacement with `zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change. `replacement_method` selects Bayesian-multiplicative (GBM, default) or multiplicative (CZM) replacement; `replacement_delta` sets the replacement fraction and optional `bayesian_alpha` overrides GBM's concentration. The offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md). +- `multinomial`: one joint multinomial logistic composition model against a + reference taxon. It uses no library-size offset or zero replacement; estimates + are log ratios against that reference, with explicit convergence and runtime + limits. `mn_reference_taxon`, `mn_penalty`, `max_iter`, `tolerance`, and + `max_seconds` control the fit. +- `dirichlet_multinomial`: the joint reference-taxon model with an estimated + Dirichlet concentration for extra-multinomial overdispersion. Its diagnostics + report the concentration, theta, and boundary status; + `dm_overdispersion_method` selects `mle` or `mom`. +- `ancom_bc`: ANCOM-BC sampling-fraction bias-corrected log-linear model. The + effect is a bias-corrected log2 fold change; the implementation reports its + mixture and bias diagnostics and requires at least three tested taxa. 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. `mc_samples` controls + the bounded, study-seeded instance count. 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. Differentials sum to zero + across fitted taxa; the model conditions on read totals and does not add an + extra-multinomial biological-overdispersion parameter. +- `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. +- `hurdle_nb`: a logistic occurrence model plus zero-truncated NB counts for + positive observations. It uses the same joint test and offset convention. ZINB and hurdle models can overfit sparse data; do not select ZINB solely because counts are sparse. Compare against NB diagnostics. The methods report @@ -662,12 +684,20 @@ 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 | multinomial | dirichlet_multinomial + # ancom_bc | aldex2 | songbird | zinb | hurdle_nb + offset: tss # count models: tss | rle dispersion: ml # nb_glm: ml | glmGamPoi 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 + 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: 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 bc98c01c..ff0525e6 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -6,7 +6,8 @@ occupancy work. The GitHub Wiki roadmap still needs a synchronized update; see ## Available now -- `nb_glm` and `clr_welch` remain the existing differential methods. +- `nb_glm`, `clr_welch`, `multinomial`, `dirichlet_multinomial`, `ancom_bc`, + `aldex2`, and `songbird` remain available differential methods. - `zinb` and `hurdle_nb` are available in the differential-abundance config and UI. They report the count effect and zero/occurrence effect and test group effects jointly. The UI surfaces model cautions; BH adjustment remains on. diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 4d235d2e..ccc63fd2 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -346,15 +346,19 @@ analysis: # Differential abundance: compares two conditions on the raw read counts, with # mandatory Benjamini-Hochberg adjusted p-values across the fitted family. differential: - # Per-taxon model: nb_glm (negative-binomial GLM), clr_welch (Welch test on - # 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. + # Model fitted across conditions: nb_glm (negative-binomial GLM), clr_welch + # (Welch test on zero-replaced CLR values), multinomial or + # dirichlet_multinomial (joint reference-taxon composition models), ancom_bc + # (bias-corrected log-linear model), aldex2 (Monte-Carlo CLR effects), + # songbird (CLR-centred multinomial differentials), zinb (zero-inflated NB), + # or hurdle_nb (logistic occurrence plus zero-truncated NB). Two-part models + # test both group effects jointly; their count-process effect is separate. 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). # rle is refused when no taxon is present in every sample. - # Used by count models (nb_glm, zinb, hurdle_nb); clr_welch ignores it. + # Used by count models (nb_glm, zinb, hurdle_nb); CLR and joint compositional + # models do not use a library-size offset. offset: tss # nb_glm only: where each taxon's negative-binomial theta comes from. ml # (the default) takes the maximum-likelihood estimate MASS::glm.nb fits per 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 98a97087..69555060 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -21,6 +21,7 @@ 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, } from './types' diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e668807d..f08c10a2 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -423,7 +423,11 @@ export interface DifferentialRow { prevalence: number } -export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' +export type DifferentialMethod = + | 'nb_glm' | 'clr_welch' + | 'multinomial' | 'dirichlet_multinomial' + | 'ancom_bc' | 'aldex2' | 'songbird' + | 'zinb' | 'hurdle_nb' /** Which row field holds the effect size, and what to call it. */ export interface DifferentialEffect { @@ -439,7 +443,7 @@ export interface DifferentialResult { groups: { reference: string; contrast: string } n_samples: Record effect: DifferentialEffect - /** Empty for clr_welch; count models report their size factors. */ + /** Empty for CLR and compositional models; count models report their size factors. */ size_factors: { sample: string; group: string; factor: number }[] config: { method: DifferentialMethod @@ -453,9 +457,19 @@ export interface DifferentialResult { replacement_delta?: number /** Null uses GBM's estimated concentration; positive values override it. */ bayesian_alpha?: number | null + /** 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 /** aldex2 only. */ mc_samples?: number - /** aldex2 only. */ + /** aldex2 only: the study seed used for its Dirichlet instances. */ seed?: number } diagnostics: { @@ -474,6 +488,15 @@ export interface DifferentialResult { concentration: { sample: string; value: number }[] } } + /** 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 + } /** ancom_bc only. */ ancom_bc?: { delta: number | null; delta_se: number | null diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx index f7dca8cf..98f10842 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -203,6 +203,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 0037ce55..c9d90d57 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -2,16 +2,34 @@ import { expect, test } from 'bun:test' import { CONFIG_DESCRIPTIONS, CONFIG_TYPES } from './PipelineStages' +const ALL_DIFFERENTIAL_METHODS = [ + 'nb_glm', 'clr_welch', 'multinomial', 'dirichlet_multinomial', + 'ancom_bc', 'aldex2', 'songbird', 'zinb', 'hurdle_nb', +] + test('dispersion uses the enum editor with both supported models', () => { expect(CONFIG_TYPES['analysis.differential.dispersion']).toEqual({ kind: 'enum', options: ['ml', 'glmGamPoi'], }) }) -test('differential model editor exposes supported two-part NB methods', () => { +test('differential model editor exposes every supported method', () => { expect(CONFIG_TYPES['analysis.differential.method']).toEqual({ - kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'], + 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', () => { @@ -24,11 +42,7 @@ 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', () => { - expect(CONFIG_TYPES['analysis.differential.method']).toEqual({ - kind: 'enum', - options: ['nb_glm', 'clr_welch', 'ancom_bc', 'aldex2', 'songbird'], - }) +test('advanced compositional settings expose bounded ALDEx2 sampling', () => { 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 f26198b8..b6e9858e 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.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.method': 'Model fitted across conditions: nb_glm (negative-binomial GLM on read counts); clr_welch (Welch test on centred log-ratios after cmultRepl zero replacement); multinomial and dirichlet_multinomial (joint compositional models with a reference taxon); ancom_bc (sampling-fraction bias-corrected log2 fold change); aldex2 (Dirichlet Monte-Carlo CLR effect sizes); songbird (CLR-centred multinomial differentials); zinb and hurdle_nb (two-part negative-binomial models with joint zero/count tests). Two-part NB results need caution and diagnostics for excess zeros.', + '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). CLR and joint compositional models do not use a library-size offset.', '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.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.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. Filtering happens before zero replacement and compositional subcomposition construction.', '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.', + 'analysis.differential.mn_penalty': "multinomial and dirichlet_multinomial only: non-negative L2 (ridge) penalty on the log-odds coefficients, 0 for maximum likelihood. Songbird's L1 sparsity path is not implemented.", + 'analysis.differential.dm_overdispersion_method': 'dirichlet_multinomial only: mle profiles the Dirichlet concentration with the coefficients; mom takes one method-of-moments estimate. The result reports concentration and theta = 1/(1+c).', + 'analysis.differential.max_iter': 'multinomial and dirichlet_multinomial only: Newton iteration cap; a fit that stops at the cap is refused rather than reported.', + 'analysis.differential.tolerance': 'multinomial and dirichlet_multinomial only: positive convergence tolerance for the Newton fit.', + 'analysis.differential.max_seconds': 'multinomial and dirichlet_multinomial only: positive wall-clock budget in seconds; exceeding it fails explicitly instead of holding the request open.', + 'analysis.differential.mc_samples': "aldex2 only: number of Dirichlet(count + 0.5) Monte-Carlo instances drawn per sample (whole number in [16, 4096]; seeded by the study's top-level seed).", '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,13 +260,18 @@ 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', 'multinomial', 'dirichlet_multinomial', 'ancom_bc', 'aldex2', 'songbird', 'zinb', 'hurdle_nb'] }, '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 }, '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.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.differential.mc_samples': { kind: 'int' }, 'analysis.alpha.show_points': { kind: 'boolean' }, 'analysis.alpha.annotate_significance': { kind: 'boolean' }, diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts index 2bae9b2c..0090f3d1 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', () => { diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts index 2743553e..2dffa8d1 100644 --- a/frontend/src/components/differentialTable.ts +++ b/frontend/src/components/differentialTable.ts @@ -46,13 +46,37 @@ export function formatStat(value: number | null): string { */ export function configLine(result: DifferentialResult): string { const c = result.config - 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 '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 'clr_welch': { + const replacement = c.replacement_method ?? 'bayesian_multiplicative' + // Alpha is a Bayesian-prior setting; CZM must not claim to have used it. + 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 'nb_glm': + case 'zinb': + case 'hurdle_nb': + 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..c8ab0a45 100644 --- a/frontend/src/views/StudyView.tsx +++ b/frontend/src/views/StudyView.tsx @@ -15,10 +15,10 @@ import { ConfigAccordion } from '../components/ConfigAccordion' import { TreesPanel } from '../components/TreesPanel' import { PlacementPanel } from '../components/PlacementPanel' import { EvidenceModeToggle } from '../components/EvidenceModeToggle' -import { CladeCumulusPanel } from '../components/CladeCumulusPanel' +import { CladeCumulusPanel as EvidenceCladeCumulusPanel } from '../components/CladeCumulusPanel' import { ReportPanel } from '../components/ReportPanel' import { useNavRefresh } from '../hooks/useNavRefresh' -import { CladeCumulusPanel } from '../cladistics/CladeCumulusPanel' +import { CladeCumulusPanel as AsvCladeCumulusPanel } from '../cladistics/CladeCumulusPanel' type GroupedRun = Run & { group?: string | null } @@ -300,11 +300,11 @@ export function StudyView() { source={analysisSource} /> :

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

)} - {view === 'cladistics' && } + {view === 'cladistics' && } {view === 'trees' && <> - {evidenceMode && } + {evidenceMode && } setTreesVersion(v => v + 1)} /> } diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index a1df0e07..e4e6e4e6 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -4,7 +4,7 @@ ## Differential abundance # Per-taxon and compositional models comparing two groups of samples, with # mandatory Benjamini-Hochberg correction across the taxa that produced a -# p-value. Five methods: +# p-value. Nine methods: # # - nb_glm: a negative-binomial GLM per taxon (MASS::glm.nb) on the observed # integer read counts, with a library-size offset. There are no pseudocounts. @@ -23,6 +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 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 @@ -37,14 +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, OFFSET_METHODS, METHODS, DISPERSIONS, - REPLACEMENT_METHODS, MIN_MC_SAMPLES, MAX_MC_SAMPLES, DEFAULT_MC_SAMPLES, + 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`. -const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb") +# 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", COMPOSITIONAL_METHODS..., + "ancom_bc", "aldex2", "songbird", "zinb", "hurdle_nb") # The per-taxon dispersion `nb_glm` uses: maximum likelihood (MASS::glm.nb's # own theta) or glmGamPoi's shrinkage estimator, which shares information @@ -68,6 +78,12 @@ const THETA_LOWER = 1e-8 # zCompositions::cmultRepl's default `frac`, as in Martín-Fernández et al. (2003). const DEFAULT_DELTA = 0.65 +# 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). const MIN_MC_SAMPLES = 16 const MAX_MC_SAMPLES = 4096 @@ -88,18 +104,24 @@ 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) + 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"`, `"ancom_bc"`, `"aldex2"` or `"songbird"`. `offset` names the -size factor method (`"tss"` or `"rle"`) and is used by `nb_glm` only. +`"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 `"bayesian_multiplicative"`. The multinomial settings select a reference +taxon, an optional non-negative ridge penalty, the Dirichlet overdispersion +estimator, and explicit iteration/tolerance/runtime limits. `mc_samples` is +the bounded ALDEx2 Monte-Carlo instance count. `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 @@ -125,6 +147,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", @@ -132,6 +160,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( @@ -188,7 +222,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 @@ -1971,7 +2007,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"`). +`"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 @@ -1991,6 +2028,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 @@ -2134,7 +2182,271 @@ 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 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 + 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 ("zinb", "hurdle_nb") factors = size_factors(counts, samples, config.offset) tested = findall(>=(config.min_prevalence), prevalence) fits = isempty(tested) ? nothing : @@ -2165,6 +2477,8 @@ 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." + else + throw(ArgumentError("unsupported differential method '$(config.method)'")) end diff --git a/src/core/doi.jl b/src/core/doi.jl index 006b806d..684ecf36 100644 --- a/src/core/doi.jl +++ b/src/core/doi.jl @@ -18,8 +18,6 @@ module DOI # * 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 2e74e3cf..4c2a6671 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1206,10 +1206,15 @@ raise their conversion errors. """ 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 " * @@ -1222,6 +1227,12 @@ 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", 500)), + tol=Float64(val("analysis.differential.tolerance", 1e-8)), + max_seconds=Float64(val("analysis.differential.max_seconds", 300.0)), mc_samples=mc_val) end @@ -1304,6 +1315,8 @@ function _differential(study::String, body) e isa ZCompositionsUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) e isa GlmGamPoiUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) e isa RStatsUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) + e isa CompositionalRefusal && return json_error(422, "model_refused", sprint(showerror, e)) + e isa RuntimeRefusal && return json_error(422, "runtime_budget_exceeded", sprint(showerror, e)) e isa ErrorException && return json_error(422, "no_taxon_fitted", e.msg) rethrow() end diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 85749bf4..89cae411 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 dispersion isa AbstractString && m != "nb_glm" && + lowercase(strip(dispersion)) == "glmgampoi" + return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; " * + "use dispersion: ml with $(strip(method))" + 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/test/unit/test_differential.jl b/test/unit/test_differential.jl index f303bf34..bf001e3e 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -168,7 +168,7 @@ end end @testset "Configuration" begin - @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb") + @test Differential.METHODS == ("nb_glm", "clr_welch", "multinomial", "dirichlet_multinomial", "ancom_bc", "aldex2", "songbird", "zinb", "hurdle_nb") @test DifferentialConfig().offset == "tss" @test DifferentialConfig().min_prevalence == 0.0 @test DifferentialConfig("rle", 0.5).offset == "rle" @@ -1200,7 +1200,7 @@ 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", "multinomial", "dirichlet_multinomial", "ancom_bc", "aldex2", "songbird", "zinb", "hurdle_nb") @test Differential.COMPOSITIONAL_METHODS == ("multinomial", "dirichlet_multinomial") c = DifferentialConfig(; method="multinomial") @test c.method == "multinomial" From 5beded5c981a6d1b25e9f1b7c8e67d40e0f26c1a Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:18:02 +0000 Subject: [PATCH 2/8] Restore analysis module and route registration Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/MetaManifold.jl | 2 +- src/server/server.jl | 1 + test/runtests.jl | 3 ++- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 96451b6a..3e52b693 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -41,10 +41,10 @@ include("pipeline/phylogeny.jl") include("canonical.jl") # Analysis -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/server/server.jl b/src/server/server.jl index 3a2728b8..44f916f4 100644 --- a/src/server/server.jl +++ b/src/server/server.jl @@ -117,6 +117,7 @@ module Server include(joinpath(@__DIR__, "routes", "figures.jl")) include(joinpath(@__DIR__, "routes", "phylogeny.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/runtests.jl b/test/runtests.jl index af87ec30..a785175c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,7 +62,8 @@ const UNIT_FILES = [ "test_server_security.jl", "test_canonical.jl", "test_doi.jl", - "test_clade_cumulus.jl", + "test_evidence.jl", + "test_evidence_routes.jl", ] # The files named on the command line ("test_routes.jl", "routes" or From dbc75d1ace5915b062ce68b8e025c0835ead8a81 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:20:16 +0000 Subject: [PATCH 3/8] Align high-precision fixtures with independent verifier Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- test/references/check_high_precision_references.py | 2 +- test/unit/test_high_precision.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/references/check_high_precision_references.py b/test/references/check_high_precision_references.py index 3d644c8f..c63ffa23 100644 --- a/test/references/check_high_precision_references.py +++ b/test/references/check_high_precision_references.py @@ -14,7 +14,7 @@ import re import sys -PRECISION = 220 +PRECISION = 460 def references() -> dict[str, str]: diff --git a/test/unit/test_high_precision.jl b/test/unit/test_high_precision.jl index 1b903e08..9f4c9e96 100644 --- a/test/unit/test_high_precision.jl +++ b/test/unit/test_high_precision.jl @@ -10,9 +10,9 @@ const HP_LABELS = ["A", "B"] const HP_FEATURES = ["f1", "f2", "f3"] const REF_LN2 = "0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202357581305570326707516350759619307275708283714351903070386238916734711233501153644979552391204751726815749320651555247341395258829504530070953263666426541042391578149520437404303855008019441706416715186447128399681717845469570262716310645461502572074024816377733896385506952606683411372738737229289564935470257626520988596932019650586" -const REF_LN3 = "1.098612288668109691395245236922525704647490557822749451734694333637494293218608966873615754813732088787970029065957865742368004225930519821052801870767277410603162769183381367179373698844360959903742570316795911521145592" -const REF_TWO_THIRDS = "0.6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667" -const REF_H_1234 = "1.279854225833667467182185762149429928131354186921692886392239591382930699742075074568049151649059984163575294301928203177013603961545520253317776056233522504371230806742168835499680746549033237086147384098562106321046048" +const REF_LN3 = "1.098612288668109691395245236922525704647490557822749451734694333637494293218608966873615754813732088787970029065957865742368004225930519821052801870767277410603162769183381367179373698844360959903742570316795911521145591917750671347054940166775580222203170252946897560690106521505642868138036317373298577782366991654792131818149020030103823630122248652748198225991097452490896458053467008845965085748444119018857087647494867079613085829411602166121184001409826" +const REF_TWO_THIRDS = "0.6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667" +const REF_H_1234 = "1.279854225833667467182185762149429928131354186921692886392239591382930699742075074568049151649059984163575294301928203177013603961545520253317776056233522504371230806742168835499680746549033237086147384098562106321046048112868339308649672681189135461287570216937288791515535715833594380020178552966375097993963900928200120047744717834903206362014300458629735272704036751340598445797090470318251133760850603018005627582297102594277585911107434593371009711254916" function test_bigfloat_close(actual, reference, bits, ulps=16) setprecision(BigFloat, bits) do From dcd5edb4d78dfd2e31c20af4899639521e6aab43 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:22:21 +0000 Subject: [PATCH 4/8] Correct high-precision mode error expectation Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- test/unit/test_high_precision.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/unit/test_high_precision.jl b/test/unit/test_high_precision.jl index 9f4c9e96..752c9282 100644 --- a/test/unit/test_high_precision.jl +++ b/test/unit/test_high_precision.jl @@ -66,7 +66,8 @@ end sample_labels=HP_LABELS, feature_labels=HP_FEATURES, policy=policy) @test_throws ArgumentError high_precision_summary([1 -2; 3 4]; sample_labels=HP_LABELS, feature_labels=["x", "y"], policy=policy) - @test_throws UnsupportedRepresentationError high_precision_summary(HP_FIXTURE; + # A policy for another mode is a caller error, not an unrepresentable value. + @test_throws ArgumentError high_precision_summary(HP_FIXTURE; sample_labels=HP_LABELS, feature_labels=HP_FEATURES, policy=numeric_policy(:exact_counts)) From 38f39edfa9791ed8f18668b4640a8d2ae3d8559b Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:23:47 +0000 Subject: [PATCH 5/8] Align multinomial route BH assertion with reference rows Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- test/unit/test_differential.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index bf001e3e..0c296c1f 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -1710,7 +1710,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 From 27d03e23ec855b39f7b535fc0e1d7bf8c1a26055 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:35:30 +0000 Subject: [PATCH 6/8] Expose validated statistics failure annotations Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .github/workflows/ci.yml | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4854e87d..d51aa441 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -397,7 +397,16 @@ jobs: - name: Run validated-statistics tests first env: R_LIBS_SITE: /usr/local/lib/R/site-library:/usr/lib/R/site-library:/usr/lib/R/library - run: julia --project=. -t2 --compiled-modules=no test/runtests.jl numeric_policy exact_summaries high_precision derived_analyses validated_statistics + 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/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 failed (tail)::$(tail -c 3500 /tmp/julia-validated-tests.log | enc)" + fi + exit "$status" - name: Run tests env: From a0a280f1897fa007e963a2aad98b839caf2fcb82 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:37:58 +0000 Subject: [PATCH 7/8] Reject invalid differential optimizer settings Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/server/routes/analysis.jl | 12 +++++++++++- src/server/routes/config.jl | 23 ++++++++++++++++++++++ test/unit/test_differential.jl | 35 ++++++++++++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 1 deletion(-) diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 4c2a6671..88221433 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1219,6 +1219,16 @@ 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'")) + 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 DifferentialConfig(String(val("analysis.differential.offset", "tss")), Float64(val("analysis.differential.min_prevalence", 0.0)); method=String(val("analysis.differential.method", "nb_glm")), @@ -1230,7 +1240,7 @@ function _differential_config(study::String) 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", 500)), + max_iter=max_iter_val, tol=Float64(val("analysis.differential.tolerance", 1e-8)), max_seconds=Float64(val("analysis.differential.max_seconds", 300.0)), mc_samples=mc_val) diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 89cae411..4cf92cfb 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -164,6 +164,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/test/unit/test_differential.jl b/test/unit/test_differential.jl index 0c296c1f..dae06908 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -525,6 +525,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 @@ -576,6 +595,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)))) @@ -1671,6 +1694,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 From 30871b33050f79a001f6022a8e963be9b594105a Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Fri, 23 Oct 2026 07:09:33 +0000 Subject: [PATCH 8/8] Clarify differential configuration validation docstrings --- src/analysis/differential.jl | 22 +++++++++++++++------- src/server/routes/analysis.jl | 9 ++++++--- src/server/routes/config.jl | 13 +++++++++++-- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index a2fc9b96..6e0197e5 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -118,10 +118,14 @@ 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"`. The multinomial settings select a reference -taxon, an optional non-negative ridge penalty, the Dirichlet overdispersion -estimator, and explicit iteration/tolerance/runtime limits. `mc_samples` is -the bounded ALDEx2 Monte-Carlo instance count. +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 @@ -135,9 +139,13 @@ 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, incompatible -method and dispersion choices, or non-finite or out-of-range fractions, -including settings unused by the selected method. +`"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 diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 4ec2ff34..b8284432 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1200,9 +1200,12 @@ 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) diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index f2b94f08..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