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 01/11] 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 02/11] 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 From 9f49f532dbc7a187b54af92eb35ab30ceafd4127 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:30:57 +0000 Subject: [PATCH 03/11] ci: surface validated-statistics test failures 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..9bad3b7b 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/validated-statistics-tests.log + status="${PIPESTATUS[0]}" + if [ "$status" -ne 0 ]; then + enc() { head -c 3500 | sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; } + echo "::error title=Validated-statistics tests failed::$(tail -n 120 /tmp/validated-statistics-tests.log | enc)" + fi + exit "$status" - name: Run tests env: From 1079b7713f70a95ae8a3b3a472a34154abeb9ef8 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:41:33 +0000 Subject: [PATCH 04/11] fix: repair DOI module parsing and syntax gate Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- scripts/check-julia-syntax.jl | 23 +++++++++++++++++++++-- src/core/doi.jl | 4 +--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/scripts/check-julia-syntax.jl b/scripts/check-julia-syntax.jl index 5b70d0a9..bc8a9d60 100644 --- a/scripts/check-julia-syntax.jl +++ b/scripts/check-julia-syntax.jl @@ -2,16 +2,35 @@ # Parse every tracked Julia source with the pinned runtime before CI spends time # installing scientific/frontend dependencies. This checks syntax only; package # loading and behavior remain covered by the test suite. +# +# Meta.parseall does not throw for every parser failure: an unterminated block +# can be returned as Expr(:incomplete), which a bare try/catch would accept. +# Parse each expression explicitly and reject both error and incomplete results. files = sort(filter(f -> endswith(f, ".jl"), split(read(`git ls-files -z`, String), '\0'))) failed = false for file in files + issue = nothing try - Meta.parseall(read(file, String)) + source = read(file, String) + position = 1 + while position <= ncodeunits(source) + ex, next_position = Meta.parse(source, position; filename=file) + if Meta.isexpr(ex, :error) || Meta.isexpr(ex, :incomplete) + issue = ex + break + end + next_position > position || + error("Julia parser made no progress at code unit $position") + position = next_position + end catch err + issue = err + end + if !isnothing(issue) failed = true println(stderr, "Julia parse failed in $file:") - println(stderr, sprint(showerror, err)) + println(stderr, sprint(showerror, issue)) end end failed && exit(1) diff --git a/src/core/doi.jl b/src/core/doi.jl index 006b806d..f5b0bd95 100644 --- a/src/core/doi.jl +++ b/src/core/doi.jl @@ -1,8 +1,6 @@ # © 2026 Joshua Benjamin Jewell. All rights reserved. # Licensed under the GNU Affero General Public License version 3 (AGPLv3). -module DOI - # Zenodo deposition-v1 client for automated DOI minting. # # Publication is DISABLED by default. Opt-in requires DOI_ENABLED=true and a @@ -682,4 +680,4 @@ function publish_study(config::DOIConfig, study::AbstractString, run::AbstractSt return receipt end -end +end # module DOI From 8bc3b7ac13daac5691ac23815f00c0912cc4d865 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 00:59:32 +0000 Subject: [PATCH 05/11] fix(stats): aggregate DuckDB counts exactly Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .github/workflows/ci.yml | 4 +- src/server/routes/validated_statistics.jl | 85 +++++++++++++++++------ 2 files changed, 64 insertions(+), 25 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9bad3b7b..4b42f6c8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -400,11 +400,11 @@ jobs: shell: bash run: | set +e - julia --project=. -t2 --compiled-modules=no test/runtests.jl numeric_policy exact_summaries high_precision derived_analyses validated_statistics 2>&1 | tee /tmp/validated-statistics-tests.log + julia --project=. -t2 --compiled-modules=no test/runtests.jl validated_statistics numeric_policy exact_summaries high_precision derived_analyses 2>&1 | tee /tmp/validated-statistics-tests.log status="${PIPESTATUS[0]}" if [ "$status" -ne 0 ]; then enc() { head -c 3500 | sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; } - echo "::error title=Validated-statistics tests failed::$(tail -n 120 /tmp/validated-statistics-tests.log | enc)" + echo "::error title=Validated-statistics tests failed::$(tail -n 200 /tmp/validated-statistics-tests.log | enc)" fi exit "$status" diff --git a/src/server/routes/validated_statistics.jl b/src/server/routes/validated_statistics.jl index ce214c7d..398d58c9 100644 --- a/src/server/routes/validated_statistics.jl +++ b/src/server/routes/validated_statistics.jl @@ -180,8 +180,6 @@ function _validated_integer_sample_columns(con, table::String, columns::Vector{S samples end -"""Read and aggregate an integer-safe run table. Only GROUP BY is done in SQL; -all count addition is Julia BigInt. No SQL SUM and no Float64 matrix is used.""" function _reject_float_count_filters(sample_columns::Vector{String}, params::Dict{String,String}) for sample in sample_columns, prefix in ("col_min.", "col_max.") key = prefix * sample @@ -192,6 +190,9 @@ function _reject_float_count_filters(sample_columns::Vector{String}, params::Dic nothing end +"""Read and aggregate an integer-safe run table. SQL performs only projection +and filtering; feature grouping and all count addition use Julia BigInt. No SQL +SUM, arbitrary grouped value, or Float64 matrix is used.""" function _validated_duckdb_counts(con, table::String, columns::Vector{String}; prefix=nothing, taxonomic_level=nothing, where_clause::String="", where_params=Any[]) @@ -219,29 +220,67 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String}; end count_columns = join([_sql_identifier(c) for c in samples], ", ") sql = "SELECT $feature_expression AS feature, $count_columns FROM $(_sql_identifier(table)) " * - "$where_clause GROUP BY $feature_expression ORDER BY feature LIMIT $(MAX_CELLS + 1)" - frame = DataFrame(DBInterface.execute(con, sql, where_params)) - nrow(frame) <= MAX_CELLS || throw(ResourceLimitError("features", MAX_CELLS, nrow(frame))) - nrow(frame) > 0 || throw(ArgumentError("no rows remain after filtering")) - nrow(frame) * length(samples) <= MAX_CELLS || throw(ResourceLimitError( - "table cells", MAX_CELLS, nrow(frame) * length(samples))) + where_clause - features = Vector{String}(undef, nrow(frame)) - counts = Matrix{BigInt}(undef, nrow(frame), length(samples)) - for (i, row) in enumerate(eachrow(frame)) - ismissing(row.feature) && throw(UnsupportedRepresentationError( - "a non-missing feature label", "missing", "validated run statistics row $i")) - features[i] = string(row.feature) - for (j, sample) in enumerate(samples) - value = row[Symbol(sample)] - ismissing(value) && throw(UnsupportedRepresentationError( - "a non-missing integer count", "NULL", "feature '$(features[i])', sample '$sample'")) - (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError( - "an integer count", repr(value), "feature '$(features[i])', sample '$sample'")) - value >= 0 || throw(ArgumentError( - "negative count $value at feature '$(features[i])', sample '$sample'")) - counts[i, j] = BigInt(value) # every addition later is exact BigInt + # A SQL GROUP BY would either overflow fixed-width SUM or require an + # arbitrary per-group value. Stream the raw rows instead, and bound the + # aggregated matrix as features arrive rather than materialising the raw + # table in Julia. DuckDB.Tables is used through the direct DuckDB dependency. + features = String[] + rows_by_feature = Vector{Vector{BigInt}}() + feature_rows = Dict{String,Int}() + raw_rows = 0 + statement = DBInterface.prepare(con, sql, DuckDB.StreamResult) + try + result = DBInterface.execute(statement, where_params) + try + for chunk in DuckDB.Tables.partitions(result) + chunk_columns = DuckDB.Tables.columns(chunk) + for row in DuckDB.Tables.rows(chunk_columns) + raw_rows += 1 + ismissing(row.feature) && throw(UnsupportedRepresentationError( + "a non-missing feature label", "missing", + "validated run statistics row $raw_rows")) + feature = string(row.feature) + i = get(feature_rows, feature, 0) + if iszero(i) + n_features = length(features) + 1 + n_features <= MAX_CELLS || throw(ResourceLimitError( + "features", MAX_CELLS, n_features)) + n_features * length(samples) <= MAX_CELLS || throw(ResourceLimitError( + "table cells", MAX_CELLS, n_features * length(samples))) + push!(features, feature) + push!(rows_by_feature, zeros(BigInt, length(samples))) + i = length(features) + feature_rows[feature] = i + end + for (j, sample) in enumerate(samples) + value = row[Symbol(sample)] + ismissing(value) && throw(UnsupportedRepresentationError( + "a non-missing integer count", "NULL", + "feature '$feature', sample '$sample'")) + (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError( + "an integer count", repr(value), "feature '$feature', sample '$sample'")) + value >= 0 || throw(ArgumentError( + "negative count $value at feature '$feature', sample '$sample'")) + rows_by_feature[i][j] += BigInt(value) + end + end + end + finally + DBInterface.close!(result) end + finally + DBInterface.close!(statement) + end + + isempty(features) && throw(ArgumentError("no rows remain after filtering")) + order = sortperm(features) + features = features[order] + rows_by_feature = rows_by_feature[order] + counts = Matrix{BigInt}(undef, length(features), length(samples)) + for i in eachindex(features), j in eachindex(samples) + counts[i, j] = rows_by_feature[i][j] end ExactSummaries.validate_coordinates(counts, samples, features) counts, samples, features From 9d66d47c07c46516e272d3f13ca99f7b58642b3f Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 01:13:01 +0000 Subject: [PATCH 06/11] fix(stats): preserve JSON integer token validation Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/server/routes/validated_statistics.jl | 56 ++++++++++++++++++++++- test/unit/test_validated_statistics.jl | 3 ++ 2 files changed, 57 insertions(+), 2 deletions(-) diff --git a/src/server/routes/validated_statistics.jl b/src/server/routes/validated_statistics.jl index 398d58c9..038bbea1 100644 --- a/src/server/routes/validated_statistics.jl +++ b/src/server/routes/validated_statistics.jl @@ -50,6 +50,57 @@ function _statistics_body(req) parsed end +"""Reject decimal/exponent tokens whose representation JSON3 would erase. +JSON3 1.14 normalises integral `1.0` and `1e0` tokens to Int64. Every numeric +field in a direct validated-statistics request is contractually integral, so a +lexical pass over the already-validated JSON preserves that boundary.""" +function _reject_noninteger_json_tokens(source::AbstractString) + bytes = codeunits(source) + position = 1 + while position <= length(bytes) + byte = bytes[position] + if byte == UInt8('"') + position += 1 + closed = false + while position <= length(bytes) + byte = bytes[position] + if byte == UInt8('\\') + position += 2 # JSON3 has already validated the escape. + elseif byte == UInt8('"') + closed = true + position += 1 + break + else + position += 1 + end + end + closed || throw(ArgumentError("unterminated JSON string")) + elseif 0x30 <= byte <= 0x39 || byte == UInt8('-') || byte == UInt8('+') + token_start = position + noninteger = false + while position <= length(bytes) + byte = bytes[position] + if 0x30 <= byte <= 0x39 || byte == UInt8('+') || byte == UInt8('-') || + byte == UInt8('.') || byte == UInt8('e') || byte == UInt8('E') + noninteger |= byte == UInt8('.') || byte == UInt8('e') || byte == UInt8('E') + position += 1 + else + break + end + end + if noninteger + token = String(bytes[token_start:position - 1]) + throw(UnsupportedRepresentationError( + "an integer JSON token or whole-number string", token, + "validated direct statistics; decimal and exponent JSON numbers are not exact input")) + end + else + position += 1 + end + end + nothing +end + function _policy_integer(raw, key::Symbol, default) value = get(raw, key, default) (value isa Integer && !(value isa Bool)) || throw(ArgumentError( @@ -95,8 +146,8 @@ function _json_exact_count(value, feature::String, sample::String) "feature '$feature', sample '$sample'; use a canonical whole-number string")) end end - # JSON3 represents a decimal/exponent token as a float even when its value - # is mathematically integral. Reject the representation, not just 1.5. + # Non-integral decimal/exponent tokens remain floats after JSON3 parsing. + # Integral ones are rejected by the raw-token pass before this function. (value isa Integer && !(value isa Bool)) || throw(UnsupportedRepresentationError( "an integer JSON token or whole-number string", repr(value), "feature '$feature', sample '$sample'; floating-point JSON is not exact input")) @@ -320,6 +371,7 @@ end @post "/api/v1/statistics/summaries" function(req) _statistics_guard() do body = _statistics_body(req) + _reject_noninteger_json_tokens(String(req.body)) policy = _statistics_policy(body) counts, labels, features, groups = _direct_count_input(body) result = _run_statistics(counts, labels, features, groups, policy) diff --git a/test/unit/test_validated_statistics.jl b/test/unit/test_validated_statistics.jl index c19bee2c..0b1a3dab 100644 --- a/test/unit/test_validated_statistics.jl +++ b/test/unit/test_validated_statistics.jl @@ -105,6 +105,9 @@ end negative_string = merge(_VS_EXACT_BODY, Dict("persist" => false, "counts" => [["-1", 1], [1, 0], [1, 0]])) @test _vs_post("/api/v1/statistics/summaries", negative_string).status == 422 + float_policy = merge(_VS_EXACT_BODY, Dict("persist" => false, + "numeric_policy" => Dict("mode" => "exact_counts", "round_digits" => 8.0))) + @test _vs_post("/api/v1/statistics/summaries", float_policy).status == 422 ordinary = merge(_VS_EXACT_BODY, Dict("numeric_policy" => Dict("mode" => "ordinary"))) @test _vs_post("/api/v1/statistics/summaries", ordinary).status == 400 ragged = merge(_VS_EXACT_BODY, Dict("counts" => [[1], [1, 2], [1, 0]])) From 5dad73675dda16105539c33915dcdb348cb0bbfd Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 01:24:58 +0000 Subject: [PATCH 07/11] test(stats): distinguish zero counts from zero totals Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- test/unit/test_exact_summaries.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/unit/test_exact_summaries.jl b/test/unit/test_exact_summaries.jl index e1ace22d..bda7c812 100644 --- a/test/unit/test_exact_summaries.jl +++ b/test/unit/test_exact_summaries.jl @@ -14,7 +14,9 @@ const EXACT_FEATURES = ["f1", "f2", "f3"] @test !s.approximate @test s.samples[1].total == 3 @test s.samples[1].proportions == fill(big(1)//3, 3) - @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), nothing, nothing] + # Zero features in a non-empty sample have defined zero proportions; + # only a zero-total sample makes every proportion undefined. + @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), big(0), big(0)] @test length(s.groups) == 2 @test s.groups[1].label == "g1" @test s.groups[2].counts == BigInt[2, 0, 0] @@ -28,7 +30,7 @@ const EXACT_FEATURES = ["f1", "f2", "f3"] @test storage["input_status"] == "exact" @test storage["samples"][1]["proportions"] == ["1/3", "1/3", "1/3"] @test storage["samples"][1]["total"] == 3 - @test storage["samples"][2]["proportions"][2] === nothing + @test storage["samples"][2]["proportions"][2] == "0/1" @test length(jcs_sha256(storage)) == 64 display = summary_to_display(s; digits=8) @test occursin("exact", display) From a412bee8cb919bd86eaf29e40d36a47813f5bb8d Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 01:38:48 +0000 Subject: [PATCH 08/11] fix(cladistics): repair phylogeny and route integration Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/analysis/clade_cumulus.jl | 23 +++++++++++++++-------- src/pipeline/phylogeny.jl | 1 + src/server/routes/clade_cumulus.jl | 12 ++++++------ test/unit/test_clade_cumulus.jl | 3 +++ 4 files changed, 25 insertions(+), 14 deletions(-) diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index a702fde2..edc15938 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -28,7 +28,9 @@ const DEFAULT_TAXONOMIC_RANKS = [ "Class", "Order", "Family", "Genus", "Species", ] -const PLACEMENT_WEIGHT_SUM_TOL = 1e-5 +# Placement weights are often rounded or truncated in serialized jplace output. +# Accept small omissions without treating materially incomplete evidence as valid. +const PLACEMENT_WEIGHT_SUM_TOL = 0.05 ## Phylogenetic tree structures & Newick / jplace parsers mutable struct PhyloNode @@ -202,8 +204,11 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]): end end - function read_label_and_quoted()::Tuple{String,Bool} - skip_ws_and_leading_comments!() + function read_label_and_quoted(; + capture_nhx::Union{Ref{Union{Float64,Nothing}},Nothing} = nothing, + capture_edge::Union{Ref{Union{Int,Nothing}},Nothing} = nothing, + )::Tuple{String,Bool} + skip_ws_and_leading_comments!(capture_nhx = capture_nhx, capture_edge = capture_edge) pos[] > n && return ("", false) if bytes[pos[]] == UInt8('\'') pos[] += 1 @@ -260,7 +265,8 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]): nhx_sup = Ref{Union{Float64,Nothing}}(nothing) edge_num = Ref{Union{Int,Nothing}}(nothing) - raw_lbl, was_quoted = read_label_and_quoted() + raw_lbl, was_quoted = read_label_and_quoted( + capture_nhx = nhx_sup, capture_edge = edge_num) skip_ws_and_leading_comments!(capture_nhx = nhx_sup, capture_edge = edge_num) branch_len = 0.0 @@ -991,7 +997,7 @@ function _graft_jplace_queries( for q in jd.queries q.status === :valid && !isempty(q.candidates) || continue - total_w = sum((c.weight for c in q.candidates); init = 0.0) + total_w = sum((c.like_weight for c in q.candidates); init = 0.0) total_w > 0.0 || continue multi = length(q.candidates) > 1 for (qname, mult) in q.names @@ -1000,7 +1006,7 @@ function _graft_jplace_queries( tx = get(asv_deepest_taxon, qname, nothing) tconf = get(asv_tax_conf, qname, 1.0) for (ci, cand) in enumerate(q.candidates) - share = cand.weight / total_w + share = cand.like_weight / total_w tip_id = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname if multi asv_reads[tip_id] = get(asv_reads, tip_id, 0.0) + base_cnt * share @@ -1009,9 +1015,10 @@ function _graft_jplace_queries( elseif !haskey(asv_reads, tip_id) asv_reads[tip_id] = base_cnt end - query_lwr[tip_id] = cand.weight + query_lwr[tip_id] = cand.like_weight push!(get!(by_edge, cand.edge_num, _GraftedQuerySpec[]), - _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_length), max(0.0, cand.pendant_length), cand.weight)) + _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_len), + max(0.0, cand.pendant_len), cand.like_weight)) end end end diff --git a/src/pipeline/phylogeny.jl b/src/pipeline/phylogeny.jl index b394d9dc..7b0fa2af 100644 --- a/src/pipeline/phylogeny.jl +++ b/src/pipeline/phylogeny.jl @@ -33,6 +33,7 @@ export REFERENCE_STEPS, PLACEMENT_STEPS, Step, placement_settings, run_workflow, using SHA, JSON3, Dates, Logging, OrderedCollections using ..PipelineLog using ..Validation + import ..Tools using ..Tools: tool_bin, _sq, _run_logged, _safe_optional_args, _num using ..RemoteExec: run_remote, checked_target import ..Provenance diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index 45825b61..172cef21 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/src/server/routes/clade_cumulus.jl @@ -93,8 +93,8 @@ function _resolve_tree_inputs_from_request( source_digests::OrderedDict{String,String}, ) tree_source = lowercase(strip(string(get(body, :treeSource, "auto")))) - tree_file = strip(string(get(body, :treeFile, ""))) - inline_newick = strip(string(get(body, :newick, ""))) + tree_file = String(strip(string(get(body, :treeFile, "")))) + inline_newick = String(strip(string(get(body, :newick, "")))) if tree_source == "study_tree" || (!isempty(tree_file) && tree_source == "auto") _valid_tree_file(tree_file) || return (nothing, nothing, :auto, @@ -200,7 +200,7 @@ end return json_error(400, "invalid_json", "Request body must be valid JSON") end - table = strip(string(get(body, :table, "merged"))) + table = String(strip(string(get(body, :table, "merged")))) _valid_name(table) || return json_error(400, "invalid_table", "Invalid table name '$table'") mode_str = lowercase(strip(string(get(body, :mode, "combined")))) @@ -268,13 +268,13 @@ end jd_div = MetaManifold.CladeCumulus.parse_jplace(jplace_str) for q in jd_div.queries q.status === :valid && !isempty(q.candidates) || continue - tw = sum((c.weight for c in q.candidates); init = 0.0) + tw = sum((c.like_weight for c in q.candidates); init = 0.0) tw > 0.0 || continue multi = length(q.candidates) > 1 for (qname, _) in q.names for (ci, cand) in enumerate(q.candidates) tid = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname - tip_weights[tid] = cand.weight / tw + tip_weights[tid] = cand.like_weight / tw end end end @@ -391,7 +391,7 @@ end (isempty(source_id) || isempty(target_id)) && return json_error(400, "missing_clade_ids", "Both sourceId and targetId are required") - table = strip(string(get(body, :table, "merged"))) + table = String(strip(string(get(body, :table, "merged")))) _valid_name(table) || return json_error(400, "invalid_table", "Invalid table name '$table'") mode_str = lowercase(strip(string(get(body, :mode, "combined")))) diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index da09bc18..70764837 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -310,6 +310,9 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 con = DBInterface.connect(db) DuckDB.register_data_frame(con, df, "src_df") DBInterface.execute(con, "CREATE TABLE merged AS SELECT * FROM src_df") + # Registration creates a persistent view backed by this Julia object; + # only the materialized table should remain in the on-disk fixture. + DuckDB.unregister_data_frame(con, "src_df") DBInterface.close!(con) close(db) before_db_sha = bytes2hex(open(sha256, db_path)) From d5f806054ee610de09fba5afb7b37f2b7ffd40d6 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 01:53:22 +0000 Subject: [PATCH 09/11] fix(cladistics): align grafting and subtree diversity Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/analysis/clade_cumulus.jl | 9 ++++++--- test/unit/test_clade_cumulus.jl | 14 +++++++------- test/unit/test_validated_statistics.jl | 12 ++++++------ 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index edc15938..79575071 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -1028,8 +1028,8 @@ function _graft_jplace_queries( function walk(u::PhyloNode)::PhyloNode rebuilt_kids = PhyloNode[walk(ch) for ch in u.children] - cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_num, rebuilt_kids) - specs = u.edge_num === nothing ? nothing : get(by_edge, u.edge_num, nothing) + cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_number, rebuilt_kids) + specs = u.edge_number === nothing ? nothing : get(by_edge, u.edge_number, nothing) (specs === nothing || isempty(specs)) && return cur # Order grafts from most distal (closest to `cur`, smallest `distal`) to most proximal. @@ -1047,7 +1047,7 @@ function _graft_jplace_queries( sub = PhyloNode(string("jplace_clade#", seq_auto[]), "", 0.0, 1.0, nothing, PhyloNode[sub, qtip]) end rem_len = max(0.0, total_len - prev_distal) - PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_num, sub.children) + PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_number, sub.children) end walk(base_root) @@ -1822,6 +1822,9 @@ function _reaggregate_after_move( merged_taxa[k] = get(merged_taxa, k, 0) + v end end + # Match clade_faith_pd: an active non-root clade includes the branch + # connecting it to its parent as well as its active descendant branches. + any_active && !is_root && (sub_pd += max(0.0, u.branch_length)) dom_tx = nothing dom_n = 0 diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index 70764837..d4d0ee25 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -14,12 +14,12 @@ const EP = MetaManifold.Epistemic const CC = MetaManifold.CladeCumulus const DM = MetaManifold.DiversityMetrics const PH = MetaManifold.Phylogeny -const SV = MetaManifold.Server +const CLADE_SV = MetaManifold.Server _cc_req(method::String, path::String, body::String = "") = - SV.Oxygen.internalrequest( - SV.HTTP.Request(method, path, ["Host" => "127.0.0.1:8080", "Content-Type" => "application/json"], body); - middleware = SV._MIDDLEWARE, + CLADE_SV.Oxygen.internalrequest( + CLADE_SV.HTTP.Request(method, path, ["Host" => "127.0.0.1:8080", "Content-Type" => "application/json"], body); + middleware = CLADE_SV._MIDDLEWARE, ) @testset "CladeCumulus & Phylogenetic Integration" begin @@ -277,8 +277,8 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 @testset "Dual-hierarchy CladeCumulus, polyphyly detection, drag-and-drop, and HTTP routes" begin tmp = mktempdir() - old_root = SV.ServerState._root[] - SV.ServerState.set_root!(tmp) + old_root = CLADE_SV.ServerState._root[] + CLADE_SV.ServerState.set_root!(tmp) try data_dir = joinpath(tmp, "data", "StudyCC", "run1") mkpath(data_dir) @@ -424,7 +424,7 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 "treeFile" => "../secret.nwk", ))).status == 400 finally - SV.ServerState._root[] = old_root + CLADE_SV.ServerState._root[] = old_root rm(tmp; recursive = true, force = true) end end diff --git a/test/unit/test_validated_statistics.jl b/test/unit/test_validated_statistics.jl index 0b1a3dab..d9c381d8 100644 --- a/test/unit/test_validated_statistics.jl +++ b/test/unit/test_validated_statistics.jl @@ -1,25 +1,25 @@ using Test using JSON3 using MetaManifold -SV = MetaManifold.Server +const VALIDATED_SV = MetaManifold.Server function _vs_request(method::String, path::String, body; raw::Bool=false) payload = raw ? body : JSON3.write(body) - req = SV.HTTP.Request(method, path, ["Content-Type" => "application/json"], payload) - SV.Oxygen.internalrequest(req) + req = VALIDATED_SV.HTTP.Request(method, path, ["Content-Type" => "application/json"], payload) + VALIDATED_SV.Oxygen.internalrequest(req) end _vs_post(path, body; raw=false) = _vs_request("POST", path, body; raw) _vs_get(path) = _vs_request("GET", path, "") function _vs_with_root(f::Function) root = mktempdir() - old = SV.ServerState._root[] + old = VALIDATED_SV.ServerState._root[] try - SV.ServerState.set_root!(root) + VALIDATED_SV.ServerState.set_root!(root) mkpath(joinpath(root, "projects")) f(root) finally - SV.ServerState._root[] = old + VALIDATED_SV.ServerState._root[] = old rm(root; recursive=true, force=true) end end From 307c9c69e513f83b3a00a20e1fae8996c3090f8e Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 07:59:03 +0100 Subject: [PATCH 10/11] Update .github/workflows/ci.yml Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4b42f6c8..946b4145 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -404,7 +404,7 @@ jobs: status="${PIPESTATUS[0]}" if [ "$status" -ne 0 ]; then enc() { head -c 3500 | sed -e 's/%/%25/g' -e 's/\r/%0D/g' | awk '{printf "%s%%0A", $0}'; } - echo "::error title=Validated-statistics tests failed::$(tail -n 200 /tmp/validated-statistics-tests.log | enc)" + echo "::error title=Validated-statistics tests failed::$(tail -c 3500 /tmp/validated-statistics-tests.log | enc)" fi exit "$status" From 1268cb031cd77192c2e89766ba254933278031b5 Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Fri, 23 Oct 2026 07:07:52 +0000 Subject: [PATCH 11/11] fix(stats): preserve byte-backed request bodies for integer token validation --- src/server/routes/validated_statistics.jl | 12 +++++++----- test/unit/test_validated_statistics.jl | 17 +++++++++++++++++ 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/src/server/routes/validated_statistics.jl b/src/server/routes/validated_statistics.jl index 038bbea1..047302ba 100644 --- a/src/server/routes/validated_statistics.jl +++ b/src/server/routes/validated_statistics.jl @@ -38,15 +38,17 @@ function _statistics_guard(f::Function) end end -function _statistics_body(req) +function _statistics_body(req; reject_noninteger_tokens::Bool=false) length(req.body) <= _STATISTICS_REQUEST_LIMIT || throw(ResourceLimitError( "request bytes", _STATISTICS_REQUEST_LIMIT, length(req.body))) + source = String(req.body) parsed = try - JSON3.read(String(req.body)) + JSON3.read(source) catch throw(ArgumentError("request body is not valid JSON")) end parsed isa JSON3.Object || throw(ArgumentError("request body must be a JSON object")) + reject_noninteger_tokens && _reject_noninteger_json_tokens(source) parsed end @@ -252,6 +254,7 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String}; samples = filter(c -> startswith(c, prefix), samples) end isempty(samples) && throw(ArgumentError("no sample columns remain after the prefix filter")) + sample_symbols = Symbol.(samples) feature_expression = if isnothing(taxonomic_level) identity = Analysis.sequence_column_name(columns) @@ -306,7 +309,7 @@ function _validated_duckdb_counts(con, table::String, columns::Vector{String}; feature_rows[feature] = i end for (j, sample) in enumerate(samples) - value = row[Symbol(sample)] + value = row[sample_symbols[j]] ismissing(value) && throw(UnsupportedRepresentationError( "a non-missing integer count", "NULL", "feature '$feature', sample '$sample'")) @@ -370,8 +373,7 @@ end @post "/api/v1/statistics/summaries" function(req) _statistics_guard() do - body = _statistics_body(req) - _reject_noninteger_json_tokens(String(req.body)) + body = _statistics_body(req; reject_noninteger_tokens=true) policy = _statistics_policy(body) counts, labels, features, groups = _direct_count_input(body) result = _run_statistics(counts, labels, features, groups, policy) diff --git a/test/unit/test_validated_statistics.jl b/test/unit/test_validated_statistics.jl index d9c381d8..8e1330e0 100644 --- a/test/unit/test_validated_statistics.jl +++ b/test/unit/test_validated_statistics.jl @@ -102,6 +102,13 @@ end # Exponent notation is a float token even though its value is integral. raw = """{"numeric_policy":{"mode":"exact_counts"},"sample_labels":["A"],"feature_labels":["f"],"counts":[[1e0]]}""" @test _vs_post("/api/v1/statistics/summaries", raw; raw=true).status == 422 + # String(::Vector{UInt8}) can consume the request buffer. Parsing and + # lexical validation must use the same string for byte-backed bodies. + @test _vs_post("/api/v1/statistics/summaries", + Vector{UInt8}(codeunits(raw)); raw=true).status == 422 + integer_raw = replace(raw, "1e0" => "1") + @test _vs_post("/api/v1/statistics/summaries", + Vector{UInt8}(codeunits(integer_raw)); raw=true).status == 200 negative_string = merge(_VS_EXACT_BODY, Dict("persist" => false, "counts" => [["-1", 1], [1, 0], [1, 0]])) @test _vs_post("/api/v1/statistics/summaries", negative_string).status == 422 @@ -209,6 +216,16 @@ end @test loaded.config.numeric_policy.mode == "exact_counts" @test !haskey(loaded.result, :analysis_id) # response-only field is not in the stored result + # Run filters retain decimal support, and prefix selection preserves + # the correspondence between sample labels and streamed columns. + filtered = _vs_post(path, merge(request, Dict("prefix" => "s2", + "colFilters" => Dict("Pident" => Dict("min" => 99.0))))) + @test filtered.status == 200 + filtered_result = JSON3.read(String(filtered.body)) + @test length(filtered_result.samples) == 1 + @test filtered_result.samples[1].label == "s2" + @test collect(filtered_result.samples[1].counts) == [string(big(2)^126)] + @test _vs_post(path, merge(request, Dict("table" => "merged_float"))).status == 422 @test _vs_post(path, merge(request, Dict("table" => "merged_negative"))).status == 400 @test _vs_post(path, merge(request, Dict("table" => "merged_null"))).status == 422