From eab3a5117de87f336305c92d7e27a9d770b86673 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:15:31 +0000 Subject: [PATCH 1/2] fix(analysis): restore differential method integrations Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- config/defaults/pipeline.yml | 23 +- frontend/src/api/client.ts | 1 + frontend/src/api/types.ts | 35 +- frontend/src/components/DifferentialPanel.tsx | 31 +- .../src/components/PipelineStages.test.ts | 20 +- frontend/src/components/PipelineStages.tsx | 21 +- .../src/components/differentialTable.test.ts | 14 +- frontend/src/components/differentialTable.ts | 40 +- frontend/src/views/StudyView.tsx | 8 +- src/analysis/differential.jl | 365 ++++++++++++++++-- src/server/routes/analysis.jl | 15 +- src/server/routes/config.jl | 31 +- test/unit/test_differential.jl | 35 +- test/unit/test_numeric_policy.jl | 8 +- 14 files changed, 546 insertions(+), 101 deletions(-) diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 4d235d2e..0be643e4 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -346,15 +346,17 @@ 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. + # Differential model: nb_glm (negative-binomial GLM), clr_welch (Welch test + # on zero-replaced CLR values), zinb or hurdle_nb (two-part count models), + # ancom_bc (bias-corrected log-linear model), aldex2 (Monte-Carlo CLR), + # songbird (CLR-centred multinomial regression), multinomial, or + # dirichlet_multinomial (joint compositional models). Two-part models test + # both group effects jointly and report the count-process fold change. 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. + # How library size enters each count model, as a log offset: tss (total reads + # per sample) or rle (median-of-ratios over taxa present in every sample). + # rle is refused when no taxon is present in every sample. Used by nb_glm, + # zinb, and hurdle_nb; compositional methods do not use a size-factor 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 @@ -366,9 +368,8 @@ analysis: dispersion: ml # Fraction of samples, 0 to 1, in which a taxon must have reads to be # tested. Taxa below it are listed as filtered, without a p-value. Applied - # before compositional subcomposition construction in clr_welch, ancom_bc, - # aldex2, and songbird. Bayesian replacement refuses a taxon seen in fewer - # than 2 samples; CZM requires at least one. + # before compositional subcomposition construction. Bayesian replacement + # refuses a taxon seen in fewer than 2 samples; CZM requires at least one. min_prevalence: 0.0 # clr_welch only: how cmultRepl replaces zeros before the centred log-ratio. # bayesian_multiplicative uses the geometric Bayesian-multiplicative (GBM) diff --git a/frontend/src/api/client.ts b/frontend/src/api/client.ts index 98a97087..5b58aaa7 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -23,6 +23,7 @@ import type { ReferenceTreeDoc, ReferenceTreeSummary, ReferenceTreeDetail, AlignmentQC, TrimSettings, CladeCumulusSourcesResponse, CladeCumulusRequest, CladeCumulusResponse, CladeMoveValidationRequest, CladeMoveValidationResponse, + CladisticsSummary, CladisticsDetail, CladisticsCreateRequest, } from './types' import type { FigureDoc } from '../figure/types' diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e668807d..51f943d3 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -403,27 +403,31 @@ export interface DifferentialRow { status: DifferentialStatus note: string /** - * nb_glm / ancom_bc: natural-log fold change, contrast over reference. + * nb_glm, two-part NB, ancom_bc and joint multinomial models: natural-log + * fold change, contrast over reference (or over the reference taxon). * clr_welch: mean CLR difference. aldex2: standardized CLR effect size. * songbird: CLR-centred multinomial differential. */ estimate: number | null - /** nb_glm and ancom_bc; null for CLR/differential methods. */ + /** Count models, ancom_bc and joint multinomial models; null for CLR effects. */ log2_fold_change: number | null standard_error: number | null - /** Wald z / W (nb_glm, ancom_bc, songbird) or Welch's t (clr_welch, aldex2). */ + /** Wald z / W, Welch's t, or the two-part model's likelihood-ratio statistic. */ statistic: number | null pvalue: number | null /** Benjamini-Hochberg adjusted; null outside the tested family. */ padj: number | null - /** nb_glm only. */ + /** nb_glm and two-part NB count-process dispersion. */ dispersion_theta: number | null /** Group coefficient in the logistic zero/occurrence component (two-part NB only). */ zero_process_log_odds?: number | null prevalence: number } -export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' +export type DifferentialMethod = + | 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' + | 'ancom_bc' | 'aldex2' | 'songbird' + | 'multinomial' | 'dirichlet_multinomial' /** Which row field holds the effect size, and what to call it. */ export interface DifferentialEffect { @@ -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. */ + /** Count models report size factors; compositional models return an empty list. */ size_factors: { sample: string; group: string; factor: number }[] config: { method: DifferentialMethod @@ -457,6 +461,16 @@ export interface DifferentialResult { mc_samples?: number /** aldex2 only. */ seed?: number + /** multinomial and dirichlet_multinomial: the taxon every ratio is against. */ + mn_reference_taxon?: string + /** multinomial and dirichlet_multinomial: L2 penalty on the log-odds. */ + mn_penalty?: number + /** dirichlet_multinomial only: how the concentration was estimated. */ + dm_overdispersion_method?: 'mle' | 'mom' + /** multinomial and dirichlet_multinomial: the optimiser settings used. */ + tolerance?: number + max_iter?: number + max_seconds?: number } diagnostics: { n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number @@ -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..0f61f6e9 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -21,11 +21,10 @@ function csvBlob(result: DifferentialResult): Blob { } /** - * Differential abundance between exactly two runs or sub-groups: a model per - * taxon at the chosen rank (negative-binomial GLM, CLR Welch test after zero - * replacement, or one multinomial / Dirichlet-multinomial model over the whole - * composition, as configured), with BH-adjusted p-values, shown as a volcano - * plot and a table + * Differential abundance between exactly two runs or sub-groups: the configured + * per-taxon count model, zero-replaced CLR test, advanced compositional model, + * or joint multinomial / Dirichlet-multinomial model at the chosen rank, with + * BH-adjusted p-values, shown as a volcano plot and a table * that can be downloaded or added to the report. The first selected condition * is the reference. The Advanced expander reports the settings the last fit * used, lets the reference taxon of the multinomial models be changed (only @@ -168,6 +167,24 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { with a library-size offset. Tolerates overdispersion through its own dispersion parameter and suits few samples per group, but each taxon is fitted on its own. +
  • + zinb / hurdle_nb: two-part negative-binomial models with a + logistic zero or occurrence process. Their p-value jointly tests both group + effects; the reported fold change belongs to the count process. +
    • +
  • + ancom_bc: bias-corrected log-linear model for sampling-fraction + bias, reporting a bias-corrected log2 fold change. +
    • +
  • + aldex2: Monte-Carlo CLR instances with a standardized CLR + effect size, integrating over Dirichlet posterior draws rather than one SE. +
    • +
  • + songbird: compositional multinomial regression reporting + CLR-centred differentials, which are relative log-ratios rather than absolute + fold changes. +
  • multinomial: one multinomial logistic model over the whole composition, so each coefficient is a log ratio against a reference taxon and no @@ -313,7 +330,7 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { {result.rank} {result.effect.label} - {['nb_glm', 'zinb', 'hurdle_nb'].includes(result.config.method) && Estimate (ln)} + {result.effect.key === 'log2_fold_change' && Estimate (ln)} {['zinb', 'hurdle_nb'].includes(result.config.method) && Zero-process log-odds} SE p @@ -326,7 +343,7 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { {r.taxon} {formatStat(effectValue(result, r))} - {['nb_glm', 'zinb', 'hurdle_nb'].includes(result.config.method) && {formatStat(r.estimate)}} + {result.effect.key === 'log2_fold_change' && {formatStat(r.estimate)}} {['zinb', 'hurdle_nb'].includes(result.config.method) && {formatStat(r.zero_process_log_odds ?? null)}} {formatStat(r.standard_error)} {formatStat(r.pvalue)} diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts index 0037ce55..5496b941 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -8,9 +8,13 @@ test('dispersion uses the enum editor with both supported models', () => { }) }) -test('differential model editor exposes supported two-part NB methods', () => { +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'], + kind: 'enum', + options: [ + 'nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', + 'multinomial', 'dirichlet_multinomial', + ], }) }) @@ -24,11 +28,13 @@ 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 retain their editors and descriptions', () => { expect(CONFIG_TYPES['analysis.differential.mc_samples']).toEqual({ kind: 'int' }) expect(CONFIG_DESCRIPTIONS['analysis.differential.mc_samples']).toContain('Dirichlet') + expect(CONFIG_TYPES['analysis.differential.dm_overdispersion_method']).toEqual({ + kind: 'enum', options: ['mle', 'mom'], + }) + expect(CONFIG_TYPES['analysis.differential.mn_penalty']).toMatchObject({ kind: 'float', min: 0 }) + expect(CONFIG_TYPES['analysis.differential.max_iter']).toEqual({ kind: 'int' }) + expect(CONFIG_DESCRIPTIONS['analysis.differential.mn_reference_taxon']).toContain('reference') }) diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index f26198b8..86ad5d24 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -150,13 +150,20 @@ export const CONFIG_DESCRIPTIONS: Record = { 'analysis.alpha.rarefaction_iterations': 'rarefy draws averaged for each alpha diversity value.', 'analysis.beta.normalisation': 'Preparation before the Bray-Curtis dissimilarity used by NMDS and PERMANOVA: hellinger (square root of relative abundances), none, rarefy or srs.', 'analysis.beta.depth': 'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.', - 'analysis.differential.method': 'Model per taxon: nb_glm (negative-binomial GLM); clr_welch (Welch test on CLR values after configured zero replacement); zinb (zero-inflated negative binomial); hurdle_nb (logistic zero process plus zero-truncated negative binomial). Two-part NB p-values jointly test group effects in both processes; use only with evidence for excess zeros and interpret cautiously.', + 'analysis.differential.method': 'Differential model: nb_glm (negative-binomial GLM); clr_welch (Welch test on zero-replaced CLR values); zinb or hurdle_nb (two-part count models with a joint zero/count-process test); ancom_bc (bias-corrected log-linear model); aldex2 (Monte-Carlo CLR); songbird (CLR-centred multinomial regression); multinomial or dirichlet_multinomial (joint compositional models with log ratios against a reference taxon). Compositional effects are relative log-ratios, not absolute abundance changes.', 'analysis.differential.offset': 'Library-size offset in the count component of nb_glm, zinb and hurdle_nb: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none).', 'analysis.differential.dispersion': 'nb_glm only: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits.', - '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. Compositional models filter before constructing their subcomposition.', '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.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 top-level seed).', + 'analysis.differential.mn_reference_taxon': 'multinomial and dirichlet_multinomial only: the reference 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.', + 'analysis.differential.mn_penalty': 'multinomial and dirichlet_multinomial only: L2 (ridge) penalty on the log-odds coefficients, 0 for maximum likelihood. Reported standard errors are the penalised ones; an L1 sparsity path is not implemented.', + 'analysis.differential.dm_overdispersion_method': 'dirichlet_multinomial only: mle alternates profile-likelihood concentration and coefficient fits; mom takes one method-of-moments estimate. The result reports concentration c and intra-sample correlation 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: convergence tolerance for the largest coefficient update relative to coefficient size.', + 'analysis.differential.max_seconds': 'multinomial and dirichlet_multinomial only: wall-clock budget in seconds; exceeding it returns an explicit error instead of holding the request open.', 'analysis.alpha.show_points': 'Show individual samples overlaid on alpha comparison boxplots.', 'analysis.alpha.annotate_significance': 'Annotate alpha comparison panels with Kruskal-Wallis significance.', 'analysis.alpha.pairwise_brackets': 'Run BH-adjusted pairwise Wilcoxon rank-sum tests between every group pair and draw brackets, including n.s. results.', @@ -253,7 +260,10 @@ export const CONFIG_TYPES: Record = { 'analysis.alpha.rarefaction_iterations': { kind: 'int' }, 'analysis.beta.normalisation': { kind: 'enum', options: ['hellinger', 'none', 'rarefy', 'srs'] }, 'analysis.beta.depth': { kind: 'int' }, - 'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'] }, + 'analysis.differential.method': { kind: 'enum', options: [ + 'nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', + 'multinomial', 'dirichlet_multinomial', + ] }, 'analysis.differential.offset': { kind: 'enum', options: ['tss', 'rle'] }, 'analysis.differential.dispersion': { kind: 'enum', options: ['ml', 'glmGamPoi'] }, 'analysis.differential.min_prevalence': { kind: 'float', min: 0, max: 1, step: 0.05 }, @@ -261,6 +271,11 @@ export const CONFIG_TYPES: Record = { 'analysis.differential.replacement_delta': { kind: 'float', min: 0, max: 1, step: 0.05 }, 'analysis.differential.bayesian_alpha': { kind: 'float', nullable: true, min: 0, step: 0.1 }, 'analysis.differential.mc_samples': { kind: 'int' }, + 'analysis.differential.mn_penalty': { kind: 'float', min: 0, step: 0.05 }, + 'analysis.differential.dm_overdispersion_method': { kind: 'enum', options: ['mle', 'mom'] }, + 'analysis.differential.max_iter': { kind: 'int' }, + 'analysis.differential.tolerance': { kind: 'float', min: 0 }, + 'analysis.differential.max_seconds': { kind: 'float', min: 1 }, 'analysis.alpha.show_points': { kind: 'boolean' }, 'analysis.alpha.annotate_significance': { kind: 'boolean' }, 'analysis.alpha.pairwise_brackets': { kind: 'boolean' }, diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts index 2bae9b2c..e7611600 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) }) }) @@ -139,6 +139,16 @@ describe('clr_welch results', () => { }) }) +test('two-part count models keep their own method and offset in CSV metadata', () => { + for (const method of ['zinb', 'hurdle_nb'] as const) { + const twoPart: DifferentialResult = { + ...result, + config: { method, offset: 'rle', min_prevalence: 0.2 }, + } + expect(configLine(twoPart)).toBe(`method ${method}, offset rle, min_prevalence 0.2`) + } +}) + describe('advanced compositional methods (ancom_bc, aldex2, songbird)', () => { test('ancom_bc reports bias-corrected log2 fold change and clean CSV config line', () => { const bc: DifferentialResult = { @@ -163,7 +173,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..ea7a7d0f 100644 --- a/frontend/src/components/differentialTable.ts +++ b/frontend/src/components/differentialTable.ts @@ -46,13 +46,39 @@ 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' + const alpha = c.bayesian_alpha === null || c.bayesian_alpha === undefined + ? 'estimated' : String(c.bayesian_alpha) + const alphaText = replacement === 'bayesian_multiplicative' + ? `, bayesian_alpha ${alpha}` + : '' + return `method clr_welch, replacement_method ${replacement}, replacement_delta ${c.replacement_delta}` + + `${alphaText}, 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}` + default: + throw new Error(`Unsupported differential method: ${String(c.method)}`) + } } /** 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..848ec936 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 PlacementCladeCumulusPanel } from '../components/CladeCumulusPanel' import { ReportPanel } from '../components/ReportPanel' import { useNavRefresh } from '../hooks/useNavRefresh' -import { CladeCumulusPanel } from '../cladistics/CladeCumulusPanel' +import { CladeCumulusPanel as CladisticsCladeCumulusPanel } 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..be010b84 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. @@ -12,6 +12,8 @@ # 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 @@ -23,6 +25,8 @@ # 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. # # 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 +41,16 @@ 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, RStatsUnavailable, CompositionalRefusal, RuntimeRefusal, + OFFSET_METHODS, METHODS, DISPERSIONS, REPLACEMENT_METHODS, COMPOSITIONAL_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 # Methods accepted under `analysis.differential.method`. -const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb") +const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", + "aldex2", "songbird", "multinomial", "dirichlet_multinomial") # 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 +74,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 +# 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 +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,34 +100,23 @@ 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) - -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. -`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"`. -`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, incompatible -method and dispersion choices, or non-finite or out-of-range fractions, -including settings unused by the selected method. + 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. """ struct DifferentialConfig method :: String @@ -126,13 +127,25 @@ struct DifferentialConfig 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 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) + 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) m = lowercase(strip(String(method))) m in METHODS || throw(ArgumentError( "analysis.differential.method must be one of $(join(METHODS, ", ")), not '$method'")) @@ -145,8 +158,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 setting, but method is " * - "'$m'; only nb_glm has a dispersion to choose")) + "analysis.differential.dispersion is a negative-binomial GLM setting, but method is " * + "'$m'; only nb_glm uses the glmGamPoi dispersion trend")) (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( @@ -165,10 +178,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'; only those two models have a reference taxon")) + "dirichlet_multinomial, but method is '$m'")) (mn_penalty == 0 || compositional) || throw(ArgumentError( "analysis.differential.mn_penalty applies to method multinomial or " * - "dirichlet_multinomial, but method is '$m'; only those two models are penalised")) + "dirichlet_multinomial, but method is '$m'")) (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")) @@ -188,7 +201,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, + Int(mc_samples), ref, Float64(mn_penalty), dm, Int(max_iter), + Float64(tol), Float64(max_seconds)) end end @@ -1970,8 +1985,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"`). +labelled `reference` in `groups`, by the method selected in `config.method` +(one of `METHODS`). `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 @@ -2000,7 +2015,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. -For both methods, taxa present in fewer than `config.min_prevalence` of the +Across 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 @@ -2134,7 +2149,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 +2444,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/server/routes/analysis.jl b/src/server/routes/analysis.jl index 2e74e3cf..aabc8c0d 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1214,6 +1214,10 @@ function _differential_config(study::String) mc_val = mc_raw isa Real ? mc_raw : throw(ArgumentError( "analysis.differential.mc_samples must be a whole number in " * "[$(Differential.MIN_MC_SAMPLES), $(Differential.MAX_MC_SAMPLES)], not '$mc_raw'")) + reference = get(get(resolved, "analysis.differential.mn_reference_taxon", + (; value=nothing)), :value, nothing) + 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")), Float64(val("analysis.differential.min_prevalence", 0.0)); method=String(val("analysis.differential.method", "nb_glm")), @@ -1222,7 +1226,14 @@ function _differential_config(study::String) replacement_method=String(val("analysis.differential.replacement_method", "bayesian_multiplicative")), bayesian_alpha=alpha_val, - mc_samples=mc_val) + mc_samples=mc_val, + mn_reference_taxon=reference_val, + 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))) 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..287c5ea7 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -245,11 +245,32 @@ function _differential_settings_error(effective::AbstractDict) current = get(effective, k, _DIFFERENTIAL_SETTING_DEFAULTS[k]) current isa NamedTuple ? current.value : current end - method = value(keys_[1], "nb_glm") - dispersion = value(keys_[2], "ml") - if method isa AbstractString && dispersion isa AbstractString && - lowercase(strip(method)) != "nb_glm" && lowercase(strip(dispersion)) == "glmgampoi" - return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; use dispersion: ml with $(strip(method))" + method = value("analysis.differential.method") + dispersion = value("analysis.differential.dispersion") + reference = value("analysis.differential.mn_reference_taxon") + penalty = value("analysis.differential.mn_penalty") + overdispersion = value("analysis.differential.dm_overdispersion_method") + method isa AbstractString || return nothing + m = lowercase(strip(method)) + compositional = m in Differential.COMPOSITIONAL_METHODS + if m in Differential.METHODS && m != "nb_glm" && dispersion isa AbstractString && + lowercase(strip(dispersion)) == "glmgampoi" + return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; " * + "use dispersion: ml with $m" + end + if !isnothing(reference) && !(reference isa AbstractString && isempty(strip(reference))) && + !compositional + return "analysis.differential.mn_reference_taxon requires method: multinomial or " * + "dirichlet_multinomial; method is '$method'" + end + if penalty isa Real && !(penalty isa Bool) && penalty != 0 && !compositional + return "analysis.differential.mn_penalty requires method: multinomial or " * + "dirichlet_multinomial; method is '$method'" + end + if overdispersion isa AbstractString && lowercase(strip(overdispersion)) != "mle" && + m != "dirichlet_multinomial" + return "analysis.differential.dm_overdispersion_method requires method: " * + "dirichlet_multinomial; method is '$method'" end nothing end diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index f303bf34..234437dc 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -168,7 +168,9 @@ end end @testset "Configuration" begin - @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb") + @test Differential.METHODS == + ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", "aldex2", + "songbird", "multinomial", "dirichlet_multinomial") @test DifferentialConfig().offset == "tss" @test DifferentialConfig().min_prevalence == 0.0 @test DifferentialConfig("rle", 0.5).offset == "rle" @@ -589,7 +591,8 @@ end @testset "Config rejects incompatible effective dispersion pairs" begin SV = MetaManifold.Server keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" - for m in ("clr_welch", "ancom_bc", "aldex2", "songbird") + for m in Differential.METHODS + m == "nb_glm" && continue @test SV._differential_pair_error(Dict(keym=>m), Dict(keyd=>"glmGamPoi")) !== nothing @test SV._differential_pair_error(Dict(keyd=>"glmGamPoi"), Dict(keym=>m)) !== nothing @test SV._differential_pair_error(Dict(keyd=>"glmGamPoi"), Dict(keym=>m, keyd=>"ml")) === nothing @@ -615,6 +618,25 @@ end end end + @testset "Config rejects incompatible joint-compositional settings" begin + SV = MetaManifold.Server + keym = "analysis.differential.method" + keyr = "analysis.differential.mn_reference_taxon" + keyp = "analysis.differential.mn_penalty" + keyo = "analysis.differential.dm_overdispersion_method" + for m in filter(m -> !(m in Differential.COMPOSITIONAL_METHODS), Differential.METHODS) + @test SV._differential_pair_error(Dict(keym=>m), Dict(keyr=>"Up")) !== nothing + @test SV._differential_pair_error(Dict(keym=>m), Dict(keyp=>0.5)) !== nothing + @test SV._differential_pair_error(Dict(keym=>m), Dict(keyo=>"mom")) !== nothing + end + @test SV._differential_pair_error(Dict(keym=>"nb_glm"), Dict(keyr=>" ")) === nothing + @test SV._differential_pair_error(Dict(keym=>"multinomial"), + Dict(keyr=>"Up", keyp=>0.5)) === nothing + @test SV._differential_pair_error(Dict(keym=>"multinomial"), Dict(keyo=>"mom")) !== nothing + @test SV._differential_pair_error(Dict(keym=>"dirichlet_multinomial"), + Dict(keyr=>"Up", keyp=>0.5, keyo=>"mom")) === nothing + end + @testset "Config PATCH publishes canonical, complete documents" begin SV = MetaManifold.Server keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" @@ -1200,7 +1222,8 @@ end # Both methods are selectable names, and the two are told apart from the # per-taxon models by one tuple. @test Differential.METHODS == - ("nb_glm", "clr_welch", "multinomial", "dirichlet_multinomial") + ("nb_glm", "clr_welch", "zinb", "hurdle_nb", "ancom_bc", "aldex2", + "songbird", "multinomial", "dirichlet_multinomial") @test Differential.COMPOSITIONAL_METHODS == ("multinomial", "dirichlet_multinomial") c = DifferentialConfig(; method="multinomial") @test c.method == "multinomial" @@ -1684,7 +1707,7 @@ end _da_fixture() do root = MetaManifold.Server.ServerState._root[] write(joinpath(root, "data", "studyD", "pipeline.yml"), - "analysis:\n differential:\n method: multinomial\n") + "analysis:\n differential:\n method: multinomial\n mn_penalty: 0.0\n max_iter: 50\n tolerance: 0.00000001\n max_seconds: 30.0\n") r = _da_post("studyD", Dict("runs" => [Dict("run" => "runA"), Dict("run" => "runB")], "table" => "merged", "rank" => "Genus")) body = JSON3.read(String(r.body)) @@ -1704,6 +1727,10 @@ end @test byt["Up"].estimate ≈ 1.972120415794 rtol=1e-6 @test byt["Down"].estimate ≈ -2.185860818556 rtol=1e-6 @test body.config.mn_reference_taxon == "Flat" + @test body.config.mn_penalty == 0.0 + @test body.config.max_iter == 50 + @test body.config.tolerance == 1e-8 + @test body.config.max_seconds == 30.0 @test body.diagnostics.model.reference_taxon == "Flat" @test body.diagnostics.model.n_taxa_fitted == 3 @test isempty(body.size_factors) diff --git a/test/unit/test_numeric_policy.jl b/test/unit/test_numeric_policy.jl index 2f64ff8c..29d700e3 100644 --- a/test/unit/test_numeric_policy.jl +++ b/test/unit/test_numeric_policy.jl @@ -51,8 +51,12 @@ end tight = numeric_policy(mode=:exact_counts, max_denominator_bits=32) @test_throws ResourceLimitError exact_fraction(1, big(2)^40; policy=tight) @test_throws ResourceLimitError exact_fraction(big(2)^40, 1; policy=tight) - # Each addend fits the 32-bit budget on its own; their sum 1/2^31 + 1/3 needs - # 33 bits, so the limit must fire on the running total, not only on inputs. + # Check every candidate before summing: the second value has a 33-bit + # denominator, although the reduced final sum would fit the 32-bit budget. + @test_throws ResourceLimitError exact_rational_sum( + [big(1)//big(2)^31, big(1)//big(2)^31 + big(1)//3]; policy=tight) + # Independently, two reduced addends can each fit while their running total + # cannot. The limit must therefore apply to the accumulated value as well. @test_throws ResourceLimitError exact_rational_sum( [big(1)//big(2)^31, big(1)//(big(2)^31 - 1)]; policy=tight) end From 875aac52ddef6bf23d5e7926d027c8d269eb7f2e Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:23:36 +0000 Subject: [PATCH 2/2] fix(ci): check references at their embedded precision Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .../check_high_precision_references.py | 43 +++++++++++-------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/test/references/check_high_precision_references.py b/test/references/check_high_precision_references.py index 3d644c8f..bab2b9fc 100644 --- a/test/references/check_high_precision_references.py +++ b/test/references/check_high_precision_references.py @@ -8,7 +8,7 @@ from __future__ import annotations import argparse -from decimal import Decimal, getcontext +from decimal import Decimal, localcontext import json from pathlib import Path import re @@ -17,30 +17,37 @@ PRECISION = 220 -def references() -> dict[str, str]: - getcontext().prec = PRECISION - counts = (1, 2, 3, 4) - total = sum(counts) - h = sum( - -(Decimal(count) / total) * (Decimal(count) / total).ln() - for count in counts - ) - return { - "REF_LN2": str(Decimal(2).ln()), - "REF_LN3": str(Decimal(3).ln()), - "REF_TWO_THIRDS": str(Decimal(2) / Decimal(3)), - "REF_H_1234": str(h), - } +def references(precision: int = PRECISION) -> dict[str, str]: + with localcontext() as context: + context.prec = precision + counts = (1, 2, 3, 4) + total = sum(counts) + h = sum( + -(Decimal(count) / total) * (Decimal(count) / total).ln() + for count in counts + ) + return { + "REF_LN2": str(Decimal(2).ln()), + "REF_LN3": str(Decimal(3).ln()), + "REF_TWO_THIRDS": str(Decimal(2) / Decimal(3)), + "REF_H_1234": str(h), + } def check(path: Path, expected: dict[str, str]) -> list[str]: text = path.read_text(encoding="utf-8") failures = [] - for name, value in expected.items(): + for name in expected: match = re.search(rf'^const {name} = "([^"]+)"$', text, re.MULTILINE) if not match: failures.append(f"{name}: constant not found") - elif match.group(1) != value: + continue + embedded = match.group(1) + # REF_LN2 deliberately carries enough digits for the 1024-bit test; + # the other references currently use 220. Verify each value at the + # precision actually embedded instead of rejecting extra correct digits. + precision = len(Decimal(embedded).as_tuple().digits) + if embedded != references(precision)[name]: failures.append(f"{name}: embedded reference differs from Python decimal") return failures @@ -58,7 +65,7 @@ def main() -> int: if failures: print("\n".join(failures), file=sys.stderr) return 1 - print(f"Independent decimal references agree ({PRECISION}-digit context).") + print("Independent decimal references agree at their embedded precisions.") return 0 print(json.dumps({"decimal_precision": PRECISION, "references": expected}, indent=2)) return 0