diff --git a/README.md b/README.md index 429d3778..b5250406 100644 --- a/README.md +++ b/README.md @@ -305,6 +305,33 @@ not exclude samples from this analysis. multiplicative (CZM) replacement; `replacement_delta` sets the replacement fraction and optional `bayesian_alpha` overrides GBM's concentration. The offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md). +- `ancom_bc`: ANCOM-BC sampling-fraction bias-corrected log-linear model + (Lin & Peddada 2020, *Nat. Commun.* 11:3514). Estimates within-group sample + offsets and between-group sampling-fraction bias via an E-M Gaussian mixture + on `log(counts + 1)`, propagating bias-estimation variance into two-sided Wald + tests. The effect is a bias-corrected log2 fold change; recommended when + asymmetric shifts or microbial load changes would bias standard CLR centering + (requires at least 3 tested taxa; assumes fewer than half are differential). +- `aldex2`: Dirichlet(`count + 0.5`) Monte-Carlo CLR model (Fernandes et al. + 2013, 2014). Draws `mc_samples` (default 128, seeded by the study's `seed`) + posterior instances per sample and reports standardized CLR effect sizes with + consistent-sign posterior predictive Welch's t-test p-values; recommended for + sparse or variable-depth cohorts where deterministic zero imputation is + undesirable. +- `songbird`: compositional baseline-category multinomial regression with + reference-frame invariant CLR-centred differentials summing to zero across + fitted taxa (Morton et al. 2019, *Nat. Commun.* 10:2719; Zhang et al. 2017). + Conditions on read totals without extra-multinomial biological overdispersion; + completely separated single-group taxa are marked `failed` rather than + diverging. +- `multinomial` and `dirichlet_multinomial`: one model for the whole + composition rather than a fit per taxon, reported as log ratios against a + reference taxon — `mn_reference_taxon`, or the most abundant fitted taxon when + empty, which must have reads in at least half the samples. `mn_penalty` is a + ridge penalty (0 for maximum likelihood); `dm_overdispersion_method` (`mle` or + `mom`) applies to the Dirichlet-multinomial only. A taxon with no reads in one + condition has an unbounded ratio: it is reported as failed with a + `complete separation` note and left out of the fit. ZINB and hurdle models can overfit sparse data; do not select ZINB solely because counts are sparse. Compare against NB diagnostics. The methods report @@ -326,6 +353,41 @@ Implementation and reference tests: `src/analysis/differential.jl`, `src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`, and `test/unit/test_differential.jl`. +### Evidence Mode + +The **Evidence** tab turns a run's outputs into a citable epistemic record. A +claim is a sentence a reader could disagree with, a **scope** naming the part of +the data it is about (its *fibre*: the rows whose value at a rank, or whose +category in a set, is one of the named labels), and a **warrant** saying why it +follows. Claims are stored in `projects/{study}/[/]/evidence/claims.yml` +and can be downloaded as that exact YAML. + +Three views share one record: + +- **Claims** — the epistemic editor (avec fibre). A claim cannot be saved + without a scope and at least one value, and cannot be marked `supported` + without a measurement, citation or inference warrant: a definition alone + supports nothing. +- **Fibre & echo** — the fibre's share of every sample, its largest members, and + the same fibre seen through a second projection: the other classifier at the + same rank, and the parent rank. Agreement is reported with a verdict; + disagreement is a limit on how strongly the claim may be stated, and an echo + that cannot be computed says so rather than vanishing. +- **Residuals** — where the fibre departs from a null you declare, either + `proportional` (the fibre distributed across samples as the library is) or a + stated `claim_share`, with Pearson residuals, a log2 ratio and a χ² statistic + printed without a p-value. The worst departure is shown first. + +The tab's Provenance card states whether `attestation.yml` is absent, degraded +or internally inconsistent before any claim is read. Longer notes on the +definitions, the schema and the limits are kept in `docs/evidence-mode.md` +(documentation under `docs/` is kept local by policy — see `.gitignore`). + +Implementation: `src/core/evidence.jl` (schema, residuals, verdicts), +`src/server/routes/evidence.jl` (the SQL), `frontend/src/evidence/` +(editor, fibre visualiser, residual explorer) and `test/unit/test_evidence.jl` +with `test/unit/test_evidence_routes.jl`. + ### Figures and report The Figures tab lays analysis charts out on a page (A4, letter or a journal @@ -662,13 +724,22 @@ analysis: normalisation: hellinger # hellinger | none | rarefy | srs depth: 0 # rarefy/srs reads per sample; 0 = smallest non-empty sample differential: - method: nb_glm # nb_glm | clr_welch | ancom_bc | aldex2 | songbird + method: nb_glm # nb_glm | clr_welch | zinb | hurdle_nb | ancom_bc + # | aldex2 | songbird | multinomial | dirichlet_multinomial + offset: tss # count models: tss | rle dispersion: ml # nb_glm: ml | glmGamPoi + min_prevalence: 0.0 # fraction of samples a taxon must appear in replacement_method: bayesian_multiplicative # clr_welch: bayesian_multiplicative | multiplicative_replacement replacement_delta: 0.65 # clr_welch: strictly between 0 and 1 bayesian_alpha: ~ # clr_welch: null (GBM estimate) or a number > 0 mc_samples: 128 # aldex2: Dirichlet Monte-Carlo instances in [16, 4096] + mn_reference_taxon: ~ # multinomial/dirichlet: null = most abundant fitted taxon + mn_penalty: 0.0 # multinomial/dirichlet: L2 penalty, 0 = maximum likelihood + dm_overdispersion_method: mle # dirichlet_multinomial: mle | mom + max_iter: 500 # multinomial models: Newton iteration cap + tolerance: 1.0e-8 # multinomial models: convergence threshold + max_seconds: 300.0 # multinomial models: wall-clock budget for one fit nmds: max_stress: 0.2 # warn if NMDS stress exceeds this value ``` diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 4d235d2e..cfde4f9a 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -346,10 +346,26 @@ 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. + # The model fitted across the two conditions. + # - nb_glm: negative-binomial GLM on read counts with a library-size offset; + # the effect is a log2 fold change. + # - zinb / hurdle_nb: a zero/occurrence process plus a count process, tested + # jointly; the count-process fold change is reported and the zero-process + # effect stays in diagnostics. + # - clr_welch: Welch's t-test on centred log-ratios over the taxa passing + # min_prevalence after zCompositions::cmultRepl zero replacement; the + # effect is a CLR difference, not a fold change. + # - ancom_bc: bias-corrected log-linear model (Lin & Peddada 2020) that + # estimates sample-specific and between-group sampling-fraction shifts via + # an E-M Gaussian mixture; the effect is a bias-corrected log2 fold change. + # - aldex2: Dirichlet(count + 0.5) Monte-Carlo CLR instances with consistent-sign + # Welch's t-test and Fernandes et al. (2014) standardized CLR effect size. + # - songbird: compositional baseline-category multinomial regression with + # reference-frame invariant CLR-centred differentials (Morton et al. 2019; + # Zhang et al. 2017); the effect is a CLR log-ratio differential, not an + # absolute fold change. + # - multinomial / dirichlet_multinomial: one joint fit of the whole + # composition, reported as log ratios against a reference taxon. method: nb_glm # How library size enters each model, as a log offset: tss (total reads per # sample) or rle (median-of-ratios over the taxa present in every sample). diff --git a/frontend/src/api/client.ts b/frontend/src/api/client.ts index 98a97087..5b58aaa7 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -23,6 +23,7 @@ import type { ReferenceTreeDoc, ReferenceTreeSummary, ReferenceTreeDetail, AlignmentQC, TrimSettings, CladeCumulusSourcesResponse, CladeCumulusRequest, CladeCumulusResponse, CladeMoveValidationRequest, CladeMoveValidationResponse, + CladisticsSummary, CladisticsDetail, CladisticsCreateRequest, } from './types' import type { FigureDoc } from '../figure/types' diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e668807d..c02aeea2 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -423,7 +423,16 @@ export interface DifferentialRow { prevalence: number } -export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' +/** + * Methods accepted by `POST .../analysis/differential`. The first four model each + * taxon separately; `multinomial` and `dirichlet_multinomial` fit the whole + * composition at once, so they take a reference taxon and a penalty where the + * per-taxon models take a size-factor offset. + */ +export type DifferentialMethod = + | 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' + | 'ancom_bc' | 'aldex2' | 'songbird' + | 'multinomial' | 'dirichlet_multinomial' /** Which row field holds the effect size, and what to call it. */ export interface DifferentialEffect { @@ -457,6 +466,16 @@ export interface DifferentialResult { mc_samples?: number /** aldex2 only. */ seed?: number + /** multinomial and dirichlet_multinomial: taxon every log ratio is taken against. */ + mn_reference_taxon?: string | null + /** multinomial and dirichlet_multinomial: L2 penalty on the coefficients. */ + mn_penalty?: number + /** dirichlet_multinomial only: how the concentration is estimated. */ + dm_overdispersion_method?: 'mle' | 'mom' + /** Compositional solvers: convergence threshold and cost ceiling. */ + tolerance?: number + max_iter?: number + max_seconds?: number } diagnostics: { n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number @@ -487,6 +506,18 @@ export interface DifferentialResult { mc_samples: number; seed: number; dirichlet_prior: number; test: string n_taxa_in_composition: number; n_taxa_unobserved: number } + /** Compositional models only: what the fit that produced these rows cost. */ + model?: { + n_taxa_fitted: number + runtime_seconds: number + /** The taxon every log ratio was taken against. */ + reference_taxon?: string | null + penalty?: number + n_parameters?: number + iterations?: number + converged?: boolean + loglik?: number + } /** songbird only. */ songbird?: { reference_taxon: string | null; reference_prevalence?: number diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx index f7dca8cf..5fa7eaa6 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -313,7 +313,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 +326,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..f3bc2a1f 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -8,9 +8,16 @@ 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 every method the backend accepts', () => { + // Every method the backend accepts, in one list: offering a method the server + // refuses, or hiding one it accepts, are both failures here. 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 +31,15 @@ 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('compositional models expose their own settings and say what they mean', () => { + expect(CONFIG_TYPES['analysis.differential.mn_reference_taxon']).toEqual({ kind: 'string', nullable: true }) + expect(CONFIG_TYPES['analysis.differential.mn_penalty']).toMatchObject({ kind: 'float', min: 0 }) + expect(CONFIG_TYPES['analysis.differential.dm_overdispersion_method']).toEqual({ + kind: 'enum', options: ['mle', 'mom'], }) expect(CONFIG_TYPES['analysis.differential.mc_samples']).toEqual({ kind: 'int' }) expect(CONFIG_DESCRIPTIONS['analysis.differential.mc_samples']).toContain('Dirichlet') + // The penalty is the one setting a user is tempted to tune for significance, + // so the description has to warn against it. + expect(CONFIG_DESCRIPTIONS['analysis.differential.mn_penalty']).toContain('do not tune it for significance') }) diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index f26198b8..785ebb6a 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -150,9 +150,16 @@ 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': '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). Whole composition: multinomial and dirichlet_multinomial, fitted once for all taxa, reporting log2 ratios against a reference taxon. Two-part NB p-values jointly test group effects in both processes; use only with evidence for excess zeros and interpret cautiously.', 'analysis.differential.offset': 'Library-size offset in the count component of nb_glm, zinb and hurdle_nb: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none).', 'analysis.differential.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.mc_samples': 'aldex2 only: Dirichlet Monte-Carlo instances drawn to estimate the CLR sampling distribution. More instances reduce Monte-Carlo noise, not biological uncertainty.', + 'analysis.differential.mn_reference_taxon': 'multinomial and dirichlet_multinomial: the taxon every log ratio is taken against. Empty chooses the most abundant fitted taxon and records which one was used.', + 'analysis.differential.mn_penalty': 'multinomial and dirichlet_multinomial: L2 penalty on the coefficients. The default 0 is maximum likelihood; a positive value is needed when separation leaves the likelihood flat, and it shrinks standard errors, so do not tune it for significance.', + 'analysis.differential.dm_overdispersion_method': 'dirichlet_multinomial only: mle profiles the Dirichlet concentration by likelihood; mom matches it to the observed between-sample variance in one step.', + 'analysis.differential.max_iter': 'multinomial and dirichlet_multinomial: the Newton iteration cap. Reaching it reports the fit as not converged rather than raising.', + 'analysis.differential.tolerance': 'multinomial and dirichlet_multinomial: convergence threshold on the largest coefficient update, relative to the coefficient size.', + 'analysis.differential.max_seconds': 'multinomial and dirichlet_multinomial: wall-clock budget for one fit. Exceeding it fails the request instead of holding it open.', 'analysis.differential.min_prevalence': 'Fraction of samples, 0 to 1, in which a taxon must have reads to be tested; taxa below it are listed as filtered.', 'analysis.differential.replacement_method': 'clr_welch only: bayesian_multiplicative uses cmultRepl GBM posterior estimates (Martín-Fernández et al. 2015); multiplicative_replacement uses CZM (Martín-Fernández et al. 2003). Both preserve observed-part ratios, but every replacement is biased. For sparse data consider NB on raw counts, prevalence filtering, occupancy or zero-inflated models instead.', 'analysis.differential.replacement_delta': "clr_welch only: cmultRepl frac, strictly between 0 and 1. CZM starts at delta × 0.5/sample total; either method caps high imputations at delta × the taxon's smallest observed proportion. Below 0.01 or at least 0.9 warns. Do not tune delta for significance; every replacement is biased.", @@ -198,6 +205,7 @@ export type ConfigType = | { kind: 'boolean' } | { kind: 'int'; nullable?: boolean } | { kind: 'float'; nullable?: boolean; min?: number; max?: number; step?: number } + | { kind: 'string'; nullable?: boolean } | { kind: 'enum'; options: string[] } | { kind: 'multiselect'; optionsFrom: string } | { kind: 'string_list' } @@ -253,7 +261,14 @@ 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 +276,12 @@ 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_reference_taxon': { kind: 'string', nullable: true }, + 'analysis.differential.mn_penalty': { kind: 'float', min: 0, step: 0.1 }, + 'analysis.differential.dm_overdispersion_method': { kind: 'enum', options: ['mle', 'mom'] }, + 'analysis.differential.max_iter': { kind: 'int' }, + 'analysis.differential.tolerance': { kind: 'float', min: 0, step: 1e-9 }, + 'analysis.differential.max_seconds': { kind: 'float', min: 0, step: 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..0af54800 100644 --- a/frontend/src/components/differentialTable.test.ts +++ b/frontend/src/components/differentialTable.test.ts @@ -68,7 +68,7 @@ describe('multinomial results', () => { n_taxa: 2, n_tested: 1, n_failed: 0, n_boundary: 0, n_filtered: 0, model: { reference_taxon: 'Bacteroides', penalty: 0.5, n_parameters: 2, iterations: 6, - converged: true, log_likelihood: -31.2, runtime_seconds: 0.02, n_taxa_fitted: 2, + converged: true, loglik: -31.2, runtime_seconds: 0.02, n_taxa_fitted: 2, }, }, rows: [ @@ -99,7 +99,7 @@ describe('multinomial results', () => { test('the reference taxon is exported with an empty statistic set, never a substitute', () => { const lines = differentialCsv(mn).trimEnd().split('\r\n') - expect(lines[4].startsWith('Bacteroides,reference,0,0,,,,,,1,')).toBe(true) + expect(lines[4].startsWith('Bacteroides,reference,0,0,,,,,,,1,')).toBe(true) expect(effectValue(mn, mn.rows[1])).toBe(1.6) }) }) @@ -163,7 +163,7 @@ describe('advanced compositional methods (ancom_bc, aldex2, songbird)', () => { expect(effectValue(ax, ax.rows[0])).toBe(1.85) expect(configLine(ax)).toBe('method aldex2, mc_samples 128, seed 123, min_prevalence 0.05; estimate is an ALDEx2 standardized CLR effect size, not a fold change') const csvRows = differentialCsv(ax).trimEnd().split('\r\n') - expect(csvRows[4]).toBe('Bacteroides,ok,,1.85,,5.545,2.9e-8,5.8e-8,,1,') + expect(csvRows[4]).toBe('Bacteroides,ok,,1.85,,5.545,2.9e-8,5.8e-8,,,1,') }) test('songbird reports CLR-centred multinomial differential metadata', () => { diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts index 2743553e..38ba6a77 100644 --- a/frontend/src/components/differentialTable.ts +++ b/frontend/src/components/differentialTable.ts @@ -46,13 +46,48 @@ 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}` + // Historical result files name methods the current API no longer offers + // (ancom_bc, aldex2, songbird). They still have to export a truthful + // provenance line, so the formatter widens the method to a plain string + // rather than narrowing the API type to what used to be selectable. + const method: string = c.method + const prevalence = `min_prevalence ${c.min_prevalence}` + // Each model family has its own settings, and a setting that does not apply + // to the fitted method is left out rather than printed as `undefined`: this + // line is exported into the CSV, where a wrong provenance note is worse than + // a short one. + if (method === 'multinomial' || method === 'dirichlet_multinomial') { + const parts = [`method ${method}`] + if (method === 'dirichlet_multinomial' && c.dm_overdispersion_method) + parts.push(`overdispersion ${c.dm_overdispersion_method}`) + if (c.mn_reference_taxon) parts.push(`reference ${c.mn_reference_taxon}`) + if (c.mn_penalty !== undefined) parts.push(`penalty ${c.mn_penalty}`) + parts.push(prevalence) + return `${parts.join(', ')}; estimates are log2 fold changes against the reference taxon` + } + if (method === 'ancom_bc') + return `method ancom_bc, ${prevalence}; effect is bias-corrected log2 fold change (Lin & Peddada 2020)` + if (method === 'aldex2') { + const parts = ['method aldex2'] + if (c.mc_samples !== undefined) parts.push(`mc_samples ${c.mc_samples}`) + if (c.seed !== undefined) parts.push(`seed ${c.seed}`) + parts.push(prevalence) + return `${parts.join(', ')}; estimate is an ALDEx2 standardized CLR effect size, not a fold change` + } + if (method === 'songbird') + return `method songbird, ${prevalence}; estimate is a CLR-centred multinomial differential, not an absolute fold change` + if (method !== 'clr_welch') { + const offset = c.offset === undefined ? '' : `offset ${c.offset}, ` + return `method ${method}, ${offset}${prevalence}` + } const replacement = c.replacement_method ?? 'bayesian_multiplicative' + // The concentration prior belongs to the Bayesian replacement only; printing + // it for the multiplicative one would claim a setting that was never used. 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` + const alphaText = replacement === 'bayesian_multiplicative' ? `, bayesian_alpha ${alpha}` : '' + return `method clr_welch, replacement_method ${replacement}, replacement_delta ${c.replacement_delta}` + + `${alphaText}, ${prevalence}; estimate is a CLR difference, not a fold change` } /** 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..087de704 100644 --- a/frontend/src/views/StudyView.tsx +++ b/frontend/src/views/StudyView.tsx @@ -18,7 +18,11 @@ import { EvidenceModeToggle } from '../components/EvidenceModeToggle' import { CladeCumulusPanel } from '../components/CladeCumulusPanel' import { ReportPanel } from '../components/ReportPanel' import { useNavRefresh } from '../hooks/useNavRefresh' -import { CladeCumulusPanel } from '../cladistics/CladeCumulusPanel' +// The study-level Cladistics workflow and the run-level CladeCumulus explorer +// are different panels that arrived on separate branches and collide on the +// exported name. The newer one is aliased here rather than renamed at its +// definition, so anything else importing it keeps working. +import { CladeCumulusPanel as CladisticsPanel } from '../cladistics/CladeCumulusPanel' type GroupedRun = Run & { group?: string | null } @@ -300,7 +304,7 @@ export function StudyView() { source={analysisSource} /> :

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

)} - {view === 'cladistics' && } + {view === 'cladistics' && } {view === 'trees' && <> diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 96451b6a..584fff91 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -41,10 +41,12 @@ include("pipeline/phylogeny.jl") include("canonical.jl") # Analysis +# The compositional models stand alone (no R) and are consumed by the +# differential module, so they precede it. +include("analysis/compositional.jl") include("analysis/clade_cumulus.jl") include("analysis/diversity.jl") include("analysis/analysis.jl") -include("analysis/clade_cumulus.jl") include("analysis/differential.jl") # Opt-in validated descriptive statistics (issue #1). Ordinary workflows above # remain unchanged; these modes must be selected explicitly. diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index a1df0e07..553e7d58 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -35,16 +35,32 @@ using ..RRuntime: with_r_lock using ..Analysis: _palette_hex, R_WAIT_SECONDS using ..CompositionalModels: CompositionalFit, CompositionalRefusal, RuntimeRefusal, DM_METHODS, CONCENTRATION_MIN, CONCENTRATION_MAX, fit_multinomial, - fit_dirichlet_multinomial + fit_dirichlet_multinomial, DEFAULT_MAX_ITER, DEFAULT_TOL, DEFAULT_MAX_SECONDS export DifferentialConfig, ScalingRefusal, MASSUnavailable, ZCompositionsUnavailable, GlmGamPoiUnavailable, RStatsUnavailable, OFFSET_METHODS, METHODS, DISPERSIONS, + COMPOSITIONAL_METHODS, DM_METHODS, REPLACEMENT_METHODS, MIN_MC_SAMPLES, MAX_MC_SAMPLES, DEFAULT_MC_SAMPLES, tss_factors, rle_factors, size_factors, bh_adjust, validate_counts, replace_zeros, clr_transform, fit_nb, fit_two_part_nb, fit_occupancy, fit_welch, + fit_ancom_bc, fit_aldex2, fit_songbird, differential_abundance, volcano_chart -# Methods accepted under `analysis.differential.method`. -const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb") +# Methods accepted under `analysis.differential.method`. The per-taxon models +# fit each taxon on its own; the compositional models fit the whole composition +# at once, which is why they take a reference taxon rather than a size-factor +# offset. Three families share this one list: +# +# count models nb_glm, zinb, hurdle_nb (R, per taxon) +# CLR models clr_welch, ancom_bc, aldex2, songbird +# composition models multinomial, dirichlet_multinomial (Julia, joint fit) +const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb", + "ancom_bc", "aldex2", "songbird", + "multinomial", "dirichlet_multinomial") + +# The subset fitted by `CompositionalModels`. Everything that belongs to one of +# these and to nothing else keys off this tuple: the reference taxon, the ridge +# penalty, the concentration estimator. +const COMPOSITIONAL_METHODS = ("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 @@ -88,12 +104,16 @@ end DifferentialConfig(offset="tss", min_prevalence=0.0; method="nb_glm", replacement_delta=0.65, dispersion="ml", replacement_method="bayesian_multiplicative", - bayesian_alpha=nothing, mc_samples=128) + bayesian_alpha=nothing, 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, read from `analysis.differential` in the pipeline config. `method` is one of `"nb_glm"`, -`"clr_welch"`, `"ancom_bc"`, `"aldex2"` or `"songbird"`. `offset` names the -size factor method (`"tss"` or `"rle"`) and is used by `nb_glm` only. +`"clr_welch"`, `"zinb"`, `"hurdle_nb"`, `"ancom_bc"`, `"aldex2"`, +`"songbird"`, `"multinomial"` or `"dirichlet_multinomial"`. `offset` names the +size factor method (`"tss"` or `"rle"`) and is used by the count models 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`. @@ -111,6 +131,14 @@ test for the group effect in both processes. Their count-component fold change i 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. +`mn_reference_taxon`, `mn_penalty`, `dm_overdispersion_method`, `max_iter`, +`tol` and `max_seconds` belong to the two compositional models: the taxon every +reported log ratio is taken against (empty means the most abundant fitted taxon), +the L2 penalty on the coefficients, how the Dirichlet concentration is +estimated, and the optimiser's iteration cap, convergence threshold and +wall-clock budget. Setting any of them for a method that has no use for it is +an error, not a silently ignored field. + 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 @@ -126,13 +154,28 @@ struct DifferentialConfig replacement_method :: String bayesian_alpha :: Union{Float64,Nothing} mc_samples :: Int + # Compositional models only. They are carried on every config so the + # validation below can refuse them for a method that has no use for them, + # rather than silently ignoring a setting the caller believes it applied. + mn_reference_taxon :: Union{String,Nothing} + 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{AbstractString,Nothing}=nothing, + mn_penalty::Real=0.0, + dm_overdispersion_method::AbstractString="mle", + max_iter::Integer=DEFAULT_MAX_ITER, + tol::Real=DEFAULT_TOL, + max_seconds::Real=DEFAULT_MAX_SECONDS) m = lowercase(strip(String(method))) m in METHODS || throw(ArgumentError( "analysis.differential.method must be one of $(join(METHODS, ", ")), not '$method'")) @@ -188,7 +231,8 @@ 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 +2014,9 @@ end Test every taxon (column) of the samples x taxa read-count matrix for a difference in abundance between the samples labelled `contrast` and those -labelled `reference` in `groups`, by `config.method` (`"nb_glm"`, `"clr_welch"`, -`"ancom_bc"`, `"aldex2"` or `"songbird"`). +labelled `reference` in `groups`, by `config.method`: `"nb_glm"`, `"clr_welch"`, +`"zinb"`, `"hurdle_nb"`, `"ancom_bc"`, `"aldex2"`, `"songbird"`, `"multinomial"` +or `"dirichlet_multinomial"`. `nb_glm`: size factors are computed on the full matrix by `config.offset` and enter each negative-binomial model as `log(factor)`; the estimate is a @@ -2000,7 +2045,39 @@ 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 +`ancom_bc`: filters by `min_prevalence` (excluding unobserved taxa; needs at +least 3 tested taxa), estimates sample-specific and between-group +sampling-fraction shifts via Lin & Peddada (2020)'s E-M Gaussian mixture on +`log(counts + 1)`, and reports bias-corrected log2 fold changes with propagated +bias-estimation variance. + +`aldex2`: filters by `min_prevalence` (excluding unobserved taxa; needs at least +2 tested taxa), draws `config.mc_samples` Dirichlet(count + 0.5) CLR instances +seeded by `seed`, and reports Fernandes et al. (2014)'s standardized CLR effect +size alongside consistent-sign posterior predictive Welch p-values. + +`songbird`: filters by `min_prevalence` (excluding unobserved taxa; needs at +least 2 identifiable taxa present in both groups and a reference taxon with +prevalence >= 0.5), fits the baseline-category multinomial logit model (Morton +et al. 2019; Zhang et al. 2017) by damped Newton-Raphson, and reports +reference-frame invariant CLR-centred differentials (`sum == 0`) with exact +delta-method standard errors and two-sided Wald z-tests. + +`multinomial` and `dirichlet_multinomial` fit the whole composition once rather +than one taxon at a time, so their estimates are log ratios against a reference +taxon and never absolute abundances. The reference is `config.mn_reference_taxon` +when set, and otherwise the most abundant taxon that passes the filter; it must +be a taxon at this rank, must be in the model, and must have reads in at least +half the samples, because every other ratio is anchored on it. A taxon with no +reads in one of the two groups has an unbounded ratio: it is reported as +`failed` with a `complete separation` note and left out of the fit. +`config.mn_penalty` is the ridge penalty, `0` for maximum likelihood. The +Dirichlet-multinomial additionally estimates a concentration by +`config.dm_overdispersion_method` (`"mle"` or `"mom"`); when that estimate +reaches its upper bound the table is multinomial to double precision, which is +reported as `boundary` rows and as `diagnostics.overdispersion.at_bound`. + +Across every method, 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,6 +2211,270 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, for (i, s) in enumerate(replaced.posterior.concentration)]) end diagnostics["zero_replacement"] = zero_diagnostics + elseif config.method == "ancom_bc" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "ANCOM-BC needs at least 2 samples in each group to estimate within-group variance; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 3 || throw(ArgumentError( + "ancom_bc needs at least 3 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence) to separate " * + "sampling-fraction bias from taxon effects; there are $(length(tested))")) + bc = fit_ancom_bc(counts[:, tested], groups; levels, samples) + for (k, j) in enumerate(tested) + f = bc.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "log2_fold_change" => isnothing(f.estimate) ? nothing : f.estimate / log(2), + "standard_error" => f.se, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + method_text = "ANCOM-BC bias-corrected log-linear model (Lin & Peddada 2020) over the taxa " * + "passing min_prevalence: within-group sample offsets and between-group " * + "sampling-fraction bias estimated via an E-M Gaussian mixture on log(counts + 1), " * + "with bias variance propagated into two-sided Wald tests and Benjamini-Hochberg " * + "adjustment over the fitted taxa. Assumes fewer than half of the tested taxa " * + "are differentially abundant" + effect = Dict("key" => "log2_fold_change", "label" => "Bias-corrected log2 FC") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence) + diagnostics["ancom_bc"] = merge(bc.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method == "aldex2" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "ALDEx2 needs at least 2 samples in each group to estimate within-group variance; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 2 || throw(ArgumentError( + "ALDEx2 needs at least 2 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence); there are $(length(tested))")) + ax = fit_aldex2(counts[:, tested], groups; levels, samples, + mc_samples=config.mc_samples, seed=Int(seed)) + for (k, j) in enumerate(tested) + f = ax.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "standard_error" => nothing, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + method_text = "ALDEx2 Dirichlet(count + 0.5) Monte-Carlo CLR model (Fernandes et al. 2014; " * + "$(config.mc_samples) instances, seed = $(Int(seed))) over the taxa passing " * + "min_prevalence: consistent-sign posterior predictive Welch's t-test p-values " * + "with Benjamini-Hochberg adjustment over the fitted taxa. Estimates are " * + "standardized CLR effect sizes (between-group difference over max within-group " * + "dispersion), not absolute fold changes" + effect = Dict("key" => "estimate", "label" => "ALDEx2 effect size") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence, + "mc_samples" => config.mc_samples, + "seed" => Int(seed)) + diagnostics["aldex2"] = merge(ax.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method == "songbird" + (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( + "songbird needs at least 2 samples in each group; " * + "'$reference' has $n_ref and '$contrast' has $n_con")) + for j in findall(==(0), present) + rows[j]["note"] = "no reads in any sample of either group, so it is not part of the composition" + end + tested = findall(j -> present[j] > 0 && prevalence[j] >= config.min_prevalence, axes(counts, 2)) + length(tested) >= 2 || throw(ArgumentError( + "songbird needs at least 2 taxa with reads at or above " * + "analysis.differential.min_prevalence = $(config.min_prevalence); there are $(length(tested))")) + sb = fit_songbird(counts[:, tested], groups; levels, samples, taxa=taxa[tested]) + for (k, j) in enumerate(tested) + f = sb.table[k, :] + merge!(rows[j], Dict{String,Any}( + "status" => f.status, "note" => f.note, "estimate" => f.estimate, + "standard_error" => f.se, "statistic" => f.statistic, + "pvalue" => f.pvalue)) + end + ref_tax = something(sb.diagnostics["reference_taxon"], "none") + method_text = "Songbird-style compositional multinomial regression (Morton et al. 2019; " * + "Zhang et al. 2017; reference taxon '$ref_tax') over the identifiable taxa " * + "passing min_prevalence: CLR-centred differentials (summing to zero across " * + "fitted taxa) with delta-method standard errors, two-sided Wald z-tests, and " * + "Benjamini-Hochberg adjustment. Differentials are relative log-ratios, not " * + "absolute fold changes, and multinomial likelihood conditions on read totals " * + "without extra-multinomial overdispersion" + effect = Dict("key" => "estimate", "label" => "Differential (CLR log-ratio)") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence) + diagnostics["songbird"] = merge(sb.diagnostics, Dict{String,Any}( + "n_taxa_in_composition" => length(tested), + "n_taxa_unobserved" => count(==(0), present))) + elseif config.method in COMPOSITIONAL_METHODS + # The multinomial models condition on each sample's total: they need no + # size factor, and a zero count is part of the likelihood rather than a + # value to replace. Two things still have to be settled before the fit. + # + # First, a taxon with no reads at all in one group has no finite log + # fold change, and because the model is fitted jointly its coefficient + # would drag the whole fit towards -Inf. Such taxa are refused by name + # and left out of the composition, exactly like taxa below the + # prevalence threshold; the note says which and why. + tested = findall(>=(config.min_prevalence), prevalence) + separated = Int[] + for j in tested + absent = [lvl for lvl in levels + if all(i -> counts[i, j] == 0, findall(==(lvl), groups))] + if !isempty(absent) + push!(separated, j) + rows[j]["status"] = "failed" + rows[j]["note"] = "no reads in any '$(absent[1])' sample (complete " * + "separation): the log fold change against the reference is unbounded, " * + "so the taxon is left out of the composition" + end + end + tested = [j for j in tested if !(j in separated)] + isempty(tested) && throw(ArgumentError( + "no taxon passes analysis.differential.min_prevalence = $(config.min_prevalence) " * + "and has reads in both '$reference' and '$contrast' samples, so there is no " * + "composition to model; raise analysis.differential.min_prevalence or use " * + "method nb_glm")) + + # Second, the reference taxon: the denominator of every reported ratio. + # An explicit choice must be a taxon at this rank and must be fitted. + ref_name = config.mn_reference_taxon + ref_pos = 0 + if isnothing(ref_name) + totals = [sum(view(counts, :, j)) for j in tested] + ref_pos = argmax(totals) + ref_name = String(taxa[tested[ref_pos]]) + else + j = findfirst(==(ref_name), taxa) + isnothing(j) && throw(ArgumentError( + "analysis.differential.mn_reference_taxon = '$(ref_name)' is not a taxon " * + "at this rank")) + pos = findfirst(==(j), tested) + isnothing(pos) && throw(ArgumentError( + "the reference taxon '$(ref_name)' is not fitted: it has reads in " * + "$(count(>(0), view(counts, :, j))) of $n samples, below " * + "analysis.differential.min_prevalence = $(config.min_prevalence), or no " * + "reads in one group. Choose another reference taxon, or lower min_prevalence")) + ref_pos = pos + end + # A reference seen in fewer than half the samples makes every ratio it + # anchors unstable, so it is refused rather than silently used. + ref_prevalence = prevalence[tested[ref_pos]] + ref_prevalence >= 0.5 || throw(ArgumentError( + "the reference taxon '$(ref_name)' has reads in only " * + "$(round(Int, ref_prevalence * n)) of $n samples; the log ratios it anchors " * + "would be unstable. Choose a reference present in at least half the samples " * + "(analysis.differential.mn_reference_taxon), or use method nb_glm or clr_welch")) + + # The design is the group indicator with an intercept: the intercept is + # the log-odds against the reference in the reference group, the slope + # the difference in those log-odds between the groups. + design = Matrix{Float64}(undef, n, 2) + for i in 1:n + design[i, 1] = 1.0 + design[i, 2] = groups[i] == contrast ? 1.0 : 0.0 + end + sub = Matrix{Int}(counts[:, tested]) + started = time() + fit = config.method == "multinomial" ? + fit_multinomial(sub, design; reference=ref_pos, penalty=config.mn_penalty, + max_iter=config.max_iter, tol=config.tol, + max_seconds=config.max_seconds) : + fit_dirichlet_multinomial(sub, design; reference=ref_pos, penalty=config.mn_penalty, + overdispersion_method=config.dm_overdispersion_method, + max_iter=config.max_iter, tol=config.tol, + max_seconds=config.max_seconds) + elapsed = time() - started + # A joint fit that has not converged describes neither taxon: there is + # no per-taxon rescue, so it is refused rather than reported. + fit.converged || throw(CompositionalRefusal( + "the $(config.method == "multinomial" ? "multinomial" : "Dirichlet-multinomial") " * + "fit did not converge in $(config.max_iter) Newton steps on $(length(tested)) taxa; " * + "raise analysis.differential.max_iter, raise min_prevalence so fewer taxa are " * + "fitted, or analyse a coarser taxonomic rank")) + c = fit.concentration + at_bound = config.method == "dirichlet_multinomial" && + !(c > CONCENTRATION_MIN * (1 + 1e-6) && c < CONCENTRATION_MAX * (1 - 1e-6)) + boundary_note = c >= CONCENTRATION_MAX * (1 - 1e-6) ? + "the Dirichlet concentration reached its upper bound: no overdispersion is " * + "detectable, so the model is effectively multinomial" : + "the Dirichlet concentration reached its lower bound: the composition is nearly " * + "degenerate and the standard errors are not trustworthy" + for (l, j) in enumerate(tested) + if l == ref_pos + rows[j]["status"] = "reference" + rows[j]["note"] = "reference taxon: the denominator of every log ratio; its " * + "coefficient is fixed at zero, so it carries no test" + rows[j]["estimate"] = 0.0 + rows[j]["log2_fold_change"] = 0.0 + continue + end + k = l < ref_pos ? l : l - 1 + est = fit.coefficients[k, 2] + se = fit.standard_errors[k, 2] + zz = fit.statistic[k, 2] + pv = fit.pvalue[k, 2] + if isnan(est) || isnan(se) || isnan(pv) + rows[j]["status"] = "failed" + rows[j]["note"] = "the fit returned a non-finite estimate or standard error" + continue + end + merge!(rows[j], Dict{String,Any}( + "status" => at_bound ? "boundary" : "ok", + "note" => at_bound ? boundary_note : "", + "estimate" => est, + "log2_fold_change" => est / log(2), + "standard_error" => se, "statistic" => zz, "pvalue" => pv, + "dispersion_theta" => nothing)) + end + model = config.method == "multinomial" ? "multinomial logistic" : + "Dirichlet-multinomial" + method_text = "$model regression over $(length(tested)) taxa with " * + "'$(ref_name)' as the reference taxon (softmax link, each sample's read " * + "composition, no size factor and no zero replacement); " * + (config.method == "multinomial" ? + "Wald test of the group coefficient from the observed information, with a " * + "ridge penalty of $(config.mn_penalty); " : + "concentration estimated by $(uppercase(config.dm_overdispersion_method)), " * + "theta = $(round(fit.overdispersion; sigdigits=3)), Wald test of the group " * + "coefficient from the observed information; ") * + "Benjamini-Hochberg adjustment over the fitted taxa. Estimates are log fold " * + "changes in relative abundance against '$(ref_name)' (a log ratio of two taxa, " * + "so coherent under compositionality): " * + (config.method == "multinomial" ? + "Songbird's model, Morton et al. 2019" : + "La Rosa et al. 2012 for the Dirichlet-multinomial") + effect = Dict("key" => "log2_fold_change", + "label" => "log2 fold change vs $(ref_name)") + size_factor_rows = Dict{String,Any}[] + config_echo = Dict{String,Any}("method" => config.method, + "min_prevalence" => config.min_prevalence, + "mn_reference_taxon" => ref_name, + "mn_penalty" => config.mn_penalty, + "tolerance" => config.tol, + "max_iter" => config.max_iter, + "max_seconds" => config.max_seconds) + config.method == "dirichlet_multinomial" && + (config_echo["dm_overdispersion_method"] = config.dm_overdispersion_method) + diagnostics["model"] = Dict{String,Any}( + "reference_taxon" => ref_name, "penalty" => config.mn_penalty, + "n_parameters" => fit.n_parameters, "iterations" => fit.iterations, + "converged" => fit.converged, "log_likelihood" => fit.loglik, + "runtime_seconds" => elapsed, "n_taxa_fitted" => length(tested)) + if config.method == "dirichlet_multinomial" + diagnostics["overdispersion"] = Dict{String,Any}( + "method" => config.dm_overdispersion_method, + "concentration" => c, "theta" => fit.overdispersion, + "at_bound" => at_bound) + end else factors = size_factors(counts, samples, config.offset) tested = findall(>=(config.min_prevalence), prevalence) diff --git a/src/core/doi.jl b/src/core/doi.jl index 006b806d..684ecf36 100644 --- a/src/core/doi.jl +++ b/src/core/doi.jl @@ -18,8 +18,6 @@ module DOI # * Durable write-ahead journal; immutable DOI/provenance receipts # * DOI insertion never rewrites scientific objects -module DOI - using Dates, SHA, HTTP, JSON3, Logging export DOIConfig, Deposition, PublicationReceipt, diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 2e74e3cf..6d931f98 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1222,7 +1222,20 @@ 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, + # Compositional models. The settings are read for every + # method because the constructor, not this route, is what + # refuses one that belongs to a different model. + mn_reference_taxon=let r = val("analysis.differential.mn_reference_taxon", nothing) + (isnothing(r) || (r isa AbstractString && isempty(strip(r)))) ? nothing : string(r) + end, + mn_penalty=Float64(val("analysis.differential.mn_penalty", 0.0)), + dm_overdispersion_method=String(val("analysis.differential.dm_overdispersion_method", + "mle")), + max_iter=Int(val("analysis.differential.max_iter", Differential.DEFAULT_MAX_ITER)), + tol=Float64(val("analysis.differential.tolerance", Differential.DEFAULT_TOL)), + max_seconds=Float64(val("analysis.differential.max_seconds", + Differential.DEFAULT_MAX_SECONDS))) end """ diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 85749bf4..e1b8c4dd 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -245,11 +245,34 @@ 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))" + # Keys are named in full rather than by position: the defaults dictionary is + # a Dict, and a merge that reorders it must not silently swap the settings + # this validation is comparing. + 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(String(method))) + compositional = m in Differential.COMPOSITIONAL_METHODS + if dispersion isa AbstractString && m != "nb_glm" && + lowercase(strip(String(dispersion))) == "glmgampoi" + return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; " * + "use dispersion: ml with $(strip(String(method)))" + end + if reference isa AbstractString && !isempty(strip(String(reference))) && !compositional + return "analysis.differential.mn_reference_taxon: a reference taxon requires " * + "method: multinomial or dirichlet_multinomial, not $(strip(String(method)))" + end + if penalty isa Real && !(penalty isa Bool) && Float64(penalty) != 0.0 && !compositional + return "analysis.differential.mn_penalty: a ridge penalty requires " * + "method: multinomial or dirichlet_multinomial, not $(strip(String(method)))" + end + if overdispersion isa AbstractString && m != "dirichlet_multinomial" && + lowercase(strip(String(overdispersion))) != "mle" + return "analysis.differential.dm_overdispersion_method: this setting requires " * + "method: dirichlet_multinomial, not $(strip(String(method)))" end nothing end diff --git a/test/runtests.jl b/test/runtests.jl index af87ec30..a785175c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,7 +62,8 @@ const UNIT_FILES = [ "test_server_security.jl", "test_canonical.jl", "test_doi.jl", - "test_clade_cumulus.jl", + "test_evidence.jl", + "test_evidence_routes.jl", ] # The files named on the command line ("test_routes.jl", "routes" or diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index f303bf34..053ebd39 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -168,7 +168,12 @@ end end @testset "Configuration" begin - @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb") + # Six selectable methods: four per-taxon models and the two + # compositional ones, which are told apart by COMPOSITIONAL_METHODS. + @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" @@ -1200,7 +1205,9 @@ 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"