diff --git a/.gitignore b/.gitignore index 181a720a..7df2318b 100644 --- a/.gitignore +++ b/.gitignore @@ -247,11 +247,13 @@ vignettes/*.pdf # R Environment Variables .Renviron -## Documentation kept local except published release, cladistics and format notes +## Documentation kept local except published release, cladistics, format, statistics, and schema notes docs/* !docs/release-notes/ !docs/cladistics/ !docs/format-migration.md +!docs/statistics/ +!docs/schemas/ # translation temp files diff --git a/.machine_readable/contractiles/intend/intend.k9.ncl b/.machine_readable/contractiles/intend/intend.k9.ncl index 017b3779..a52efae2 100644 --- a/.machine_readable/contractiles/intend/intend.k9.ncl +++ b/.machine_readable/contractiles/intend/intend.k9.ncl @@ -54,6 +54,12 @@ K9! status = "verified", summary = "Keep CladeCumulus topology proposals support-warranted, immutable as sidecars, and invariant in tip-level Faith PD across reaggregation.", }, + { + id = "compositional-constrained-ordination", + area = "src/analysis/ordination.jl", + status = "verified", + summary = "Provide RDA, CCA, CAP, and dbRDA constrained ordinations with Condition(...) partial ordination, seeded vegan::anova.cca permutation inference, adjusted-R2 variance partitioning, and KYAML/JCS provenance.", + }, { id = "honest-fleet-evidence", area = ".gitbot-fleet/enrolment.json", @@ -64,6 +70,6 @@ K9! }, validation = { - commitment_count_positive = std.array.length config.commitments >= 4, + commitment_count_positive = std.array.length config.commitments >= 5, }, } diff --git a/.machine_readable/contractiles/must/must.k9.ncl b/.machine_readable/contractiles/must/must.k9.ncl index 5aa9e87c..00bf4acb 100644 --- a/.machine_readable/contractiles/must/must.k9.ncl +++ b/.machine_readable/contractiles/must/must.k9.ncl @@ -64,6 +64,12 @@ K9! check_command = "python3 test/references/check_high_precision_references.py --check test/unit/test_high_precision.jl", rationale = "Independent arbitrary-precision reference constants must match the Julia test suite.", }, + { + id = "constrained-ordination-reference-oracle", + target = "test/references/check_ordination_references.py", + check_command = "python3 test/references/check_ordination_references.py --check test/fixtures/ordination_reference.json", + rationale = "Independent Python reference constants for RDA eigenvalues, Ezekiel adjusted-R2, two-way variance partitioning, and canonical JCS schema/fixtures must match within 1e-6.", + }, { id = "upstream-license-untouched", target = "LICENSE", @@ -76,6 +82,6 @@ K9! validation = { has_upstream_license = config.upstream_license == "AGPL-3.0-only", has_fork_license = config.fork_addition_license == "MPL-2.0", - invariant_count_positive = std.array.length config.invariants >= 5, + invariant_count_positive = std.array.length config.invariants >= 6, }, } diff --git a/.machine_readable/contractiles/trust/trust.k9.ncl b/.machine_readable/contractiles/trust/trust.k9.ncl index 6fe5b1d6..ef45f6a7 100644 --- a/.machine_readable/contractiles/trust/trust.k9.ncl +++ b/.machine_readable/contractiles/trust/trust.k9.ncl @@ -57,12 +57,17 @@ K9! surface = "test/unit/test_install_pins.jl", rule = "Workflow actions and bioinformatics tool versions remain pinned and verified by unit tests.", }, + { + id = "constrained-ordination-formula-and-compatibility-gate", + surface = "src/analysis/ordination.jl", + rule = "Constrained ordination routes enforce a strict covariate-only formula AST (rejecting arbitrary R expressions), method/distance compatibility (rejecting non-Euclidean RDA combinations such as rda + bray), permutations in [99, 9999], scaling in {1, 2}, and serialized R execution under with_r_lock.", + }, ], }, validation = { request_limit_positive = config.max_statistics_request_bytes > 0, safe_integer_bound = config.max_ieee754_safe_integer == 9007199254740991, - boundary_count_positive = std.array.length config.boundaries >= 4, + boundary_count_positive = std.array.length config.boundaries >= 5, }, } diff --git a/MetaManifold-WebUI_chora.deed b/MetaManifold-WebUI_chora.deed index 06537522..df9f568d 100644 --- a/MetaManifold-WebUI_chora.deed +++ b/MetaManifold-WebUI_chora.deed @@ -16,6 +16,8 @@ :frontend "frontend/" :analysis "src/analysis/" :cladistics "docs/cladistics/" + :statistics "docs/statistics/" + :schemas "docs/schemas/" :proofs "proofs/" :contractiles ".machine_readable/contractiles/" :enforcement ".github/workflows/") diff --git a/README.md b/README.md index 9850e50b..3e0f17a2 100644 --- a/README.md +++ b/README.md @@ -81,6 +81,7 @@ Once a run completes, analysis is performed on request through the web UI, both - Pipeline stage read-count summaries - NMDS ordination (Bray-Curtis, via R/vegan) - PERMANOVA, with PERMDISP to check for differences in dispersion (via R/vegan) +- Constrained ordinations — RDA, CCA, CAP, and dbRDA — with `Condition(...)` partial ordination, seeded `vegan::anova.cca` permutation tests, adjusted-$R^2$ variance partitioning (`vegan::varpart`), and biplot scalings 1 & 2 (see [Constrained ordinations](docs/statistics/constrained-ordinations.md)) - De novo ASV phylogenies with cumulative clade abundances, Faith's phylogenetic diversity, and weighted/unweighted UniFrac (CladeCumulus) Counts may be normalised before analysis (none, rarefaction to a fixed or auto-resolved depth, or scaling with ranked subsampling (SRS), which reaches the depth by scaling counts and keeps more of the community structure), and contamination-flagged taxa may be included or excluded. Ordinations can also apply a Hellinger transform before the Bray-Curtis dissimilarity. Composition charts can be drawn as a grid faceted on any two of run, group and sub-group. All analysis charts are returned as Plotly JSON and rendered interactively in the browser. @@ -1029,7 +1030,7 @@ The server exposes a REST API under `/api/v1/`. Key endpoint groups: | Results | `GET/POST/DELETE .../results/tables/...` | List, query, filter, save, export (`.xlsx`), and delete tables; OTU member drill-down | | QC | `GET .../results/qc`, `GET .../results/dada2` | MultiQC report metadata and DADA2 figures, logs, stats | | Analysis | `POST .../runs/{run}/analysis/{alpha,chart}`, `GET .../runs/{run}/analysis/{pipeline-stats,ranks}` | Per-run charts and rank discovery | -| Cross-run | `POST /studies/{study}/analysis/{alpha,chart,chart-facet,nmds,permanova,venn,publication-tables}` | Comparison, faceted charts, NMDS, PERMANOVA, taxon overlap and publication tables across runs | +| Cross-run | `POST /studies/{study}/analysis/{alpha,chart,chart-facet,nmds,permanova,ordination,venn,publication-tables}` | Comparison, faceted charts, NMDS, PERMANOVA, constrained ordinations (RDA/CCA/CAP/dbRDA), taxon overlap and publication tables across runs | | Composition | `POST .../runs/{run}/composition/summary`, `POST .../composition/{source}/query`, `POST .../composition/{source}/distinct/{column}` | Organism-composition summaries and tables for a run | | Composition library | `GET /composition`, `POST/DELETE /composition/{filters,sets}/{name}`, `GET /category-sets`, `POST/DELETE /category-sets/{name}` | Edit the filters and category sets in `config/composition.yml` | | Evidence | `GET .../runs/{run}/evidence`, `GET .../evidence/claims.yml`, `POST .../evidence/claims`, `DELETE .../evidence/claims/{id}`, `POST .../evidence/{fiber,residuals}` | Claims, fibres, echoes and residuals for a run; the record is exportable as YAML | diff --git a/bench/constrained_ordination/README.md b/bench/constrained_ordination/README.md new file mode 100644 index 00000000..2fa77e0f --- /dev/null +++ b/bench/constrained_ordination/README.md @@ -0,0 +1,24 @@ + + + +# Constrained Ordination Benchmarks + +Opt-in scaling and baseline-regression benchmark harness for `MetaManifold.Ordination` (`RDA`, `CCA`, `CAP`, `dbRDA`, `vegan::anova.cca`, and `vegan::varpart`). + +## Running + +```bash +# Default fast check (100 samples x 200 features, 999 permutations) +julia --project=. bench/constrained_ordination/benchmark.jl + +# Full scaling sweep including 1000 samples x 200 features (999 permutations) +julia --project=. bench/constrained_ordination/benchmark.jl --include-1000 + +# Canonical RFC 8785 JCS JSON output +julia --project=. bench/constrained_ordination/benchmark.jl --include-1000 --json +``` + +## Thresholds and Invariants + +- **5-minute warning guard**: Any method exceeding `300.0` seconds at `999` permutations is flagged (`exceeds_five_minutes: true`) in alignment with `Ordination.estimate_ordination_seconds` (`warn_five_minutes`). +- **Zero overhead on unconstrained diversity paths**: When `analysis.ordination.method` is `"none"` (the default), existing `Alpha`, `NMDS`, and `PERMANOVA` code paths in `MetaManifold.DiversityMetrics` and `MetaManifold.Analysis` execute unchanged with zero added allocations or runtime regression. diff --git a/bench/constrained_ordination/benchmark.jl b/bench/constrained_ordination/benchmark.jl new file mode 100644 index 00000000..56f15efa --- /dev/null +++ b/bench/constrained_ordination/benchmark.jl @@ -0,0 +1,239 @@ +#!/usr/bin/env julia +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Scaling and baseline-regression benchmark for Issue #19 constrained ordinations +# (RDA, CCA, CAP, dbRDA, permutation ANOVA, and variance partitioning). +# +# Usage: +# julia --project=. bench/constrained_ordination/benchmark.jl +# julia --project=. bench/constrained_ordination/benchmark.jl --include-1000 +# julia --project=. bench/constrained_ordination/benchmark.jl --json + +using DataFrames +using Dates +using MetaManifold +using MetaManifold.Analysis: r_available +using MetaManifold.Canonical +using MetaManifold.DiversityMetrics: transform_counts, shannon, simpson +using MetaManifold.Ordination +using Random + +const FIVE_MINUTE_WARN_SECONDS = 300.0 + +""" +Generate a reproducible sample-by-feature count matrix and aligned group, batch +and pH metadata using a local RNG. For positive dimensions, ensure every sample +and feature has positive total abundance. +""" +function _synthetic_microbiome(n_samples::Int, n_features::Int; seed::Int = 20261023) + rng = MersenneTwister(seed) + groups = [i <= div(n_samples, 2) ? "Control" : "Treatment" for i in 1:n_samples] + batches = [mod1(i, 3) == 1 ? "BatchA" : (mod1(i, 3) == 2 ? "BatchB" : "BatchC") for i in 1:n_samples] + ph = [6.2 + 1.4 * (i / n_samples) + 0.15 * randn(rng) for i in 1:n_samples] + + counts = Matrix{Float64}(undef, n_samples, n_features) + for i in 1:n_samples + g_shift = groups[i] == "Treatment" ? 0.85 : -0.85 + for j in 1:n_features + base_rate = exp(3.5 - 2.4 * (j / n_features) + (j <= 25 ? g_shift : 0.0) + 0.2 * (ph[i] - 6.9)) + c = floor(Int, base_rate * (0.5 + rand(rng))) + counts[i, j] = Float64(max(0, c)) + end + # Guarantee strictly positive sample depth + if sum(@view counts[i, :]) == 0.0 + counts[i, 1] = 10.0 + end + end + for j in 1:n_features + if sum(@view counts[:, j]) == 0.0 + counts[1, j] = 1.0 + end + end + + meta = DataFrame( + sample = ["S$(lpad(i, 4, '0'))" for i in 1:n_samples], + group = groups, + batch = batches, + ph = ph, + ) + (counts, meta) +end + +"""Return the median of a non-empty timing vector without changing it; empty input raises `BoundsError`.""" +function _median(xs::Vector{Float64})::Float64 + s = sort(xs) + n = length(s) + iseven(n) ? 0.5 * (s[div(n, 2)] + s[div(n, 2) + 1]) : s[div(n + 1, 2)] +end + +""" +Return median elapsed seconds over positive `reps` for Hellinger transformation +and per-sample Shannon, Simpson and squared-norm summaries, excluding warm-up. +""" +function _bench_existing_diversity(counts::Matrix{Float64}; reps::Int = 5) + # Warm-up + _ = transform_counts(counts; method="hellinger") + _ = shannon(vec(counts[1, :])) + + times = Float64[] + for _ in 1:reps + t0 = time_ns() + h = transform_counts(counts; method="hellinger") + # Pairwise Euclidean on Hellinger + Shannon + Simpson across all samples + acc = 0.0 + for i in axes(counts, 1) + acc += shannon(@view counts[i, :]) + acc += simpson(@view counts[i, :]) + acc += sum(abs2, @view h[i, :]) + end + acc >= 0.0 || error("unreachable") + push!(times, (time_ns() - t0) / 1e9) + end + _median(times) +end + +""" +Return median fit time in seconds over positive `reps` for Hellinger RDA with +group/pH constraints and a BatchA indicator as the conditioning column. +Transformation, design preparation and warm-up are excluded; fit errors propagate. +""" +function _bench_pure_julia_rda(counts::Matrix{Float64}, meta::DataFrame; reps::Int = 5) + hel = transform_counts(counts; method="hellinger") + n = size(counts, 1) + X = hcat([g == "Treatment" ? 1.0 : 0.0 for g in meta.group], Float64.(meta.ph)) + Z = reshape([b == "BatchA" ? 1.0 : 0.0 for b in meta.batch], n, 1) + _ = pure_linear_rda(hel, X; Z=Z) + + times = Float64[] + for _ in 1:reps + t0 = time_ns() + _ = pure_linear_rda(hel, X; Z=Z) + push!(times, (time_ns() - t0) / 1e9) + end + _median(times) +end + +""" +Run the 100-sample benchmark and optionally the 1,000-sample tier when +`--include-1000` is present or `BENCH_INCLUDE_1000=1`. Print canonical JSON with +`--json`, otherwise a text summary. Skip R fits when R/vegan is unavailable; +errors during selected benchmark fits propagate. Timings above 300 seconds +produce warnings in the summary without stopping the benchmark. +""" +function main(args::Vector{String}) + include_1000 = "--include-1000" in args || get(ENV, "BENCH_INCLUDE_1000", "0") == "1" + emit_json = "--json" in args + sample_tiers = include_1000 ? [100, 1000] : [100] + n_features = 200 + n_perm = 999 + + results = Dict{String,Any}[] + warnings = String[] + + for n_samples in sample_tiers + counts, meta = _synthetic_microbiome(n_samples, n_features) + div_sec = _bench_existing_diversity(counts; reps = n_samples >= 1000 ? 2 : 5) + pure_rda_sec = _bench_pure_julia_rda(counts, meta; reps = n_samples >= 1000 ? 2 : 5) + + push!(results, Dict{String,Any}( + "tier_samples" => n_samples, + "n_features" => n_features, + "engine" => "julia_baseline_diversity", + "method" => "hellinger_bray_shannon_simpson", + "permutations" => 0, + "elapsed_seconds" => round(div_sec; digits=6), + "exceeds_five_minutes" => div_sec > FIVE_MINUTE_WARN_SECONDS, + )) + push!(results, Dict{String,Any}( + "tier_samples" => n_samples, + "n_features" => n_features, + "engine" => "pure_julia_oracle", + "method" => "pure_linear_rda_partial", + "permutations" => 0, + "elapsed_seconds" => round(pure_rda_sec; digits=6), + "exceeds_five_minutes" => pure_rda_sec > FIVE_MINUTE_WARN_SECONDS, + )) + + if r_available() && vegan_available() + specs = [ + (:rda, :aitchison, "~ group + ph", true), + (:cca, :chisq, "~ group + ph", false), + (:cap, :bray, "~ group + Condition(batch)", false), + (:dbrda, :bray, "~ group + ph", true), + ] + for (m, d, form, vp) in specs + cfg = OrdinationConfig(; + ordination_method = m, + ordination_distance = d, + ordination_formula = form, + ordination_scaling = 2, + permutations = n_perm, + variance_partitioning = vp, + seed = 42, + ) + est = estimate_ordination_seconds(; + method = m, + distance = d, + n_samples = n_samples, + n_features = n_features, + permutations = n_perm, + n_terms = 2, + variance_partitioning = vp, + ) + t0 = time_ns() + res = run_constrained_ordination(counts, meta, cfg; beta_transform = (d === :euclidean ? "hellinger" : "none")) + elapsed = (time_ns() - t0) / 1e9 + slow = elapsed > FIVE_MINUTE_WARN_SECONDS + if slow + push!(warnings, "$(m)($(d)) at N=$(n_samples), perms=$(n_perm) took $(round(elapsed; digits=1))s (> 300s threshold)") + end + push!(results, Dict{String,Any}( + "tier_samples" => n_samples, + "n_features" => n_features, + "engine" => "r_vegan", + "method" => String(m), + "distance" => String(d), + "formula" => form, + "variance_partitioning" => vp, + "permutations" => n_perm, + "r2" => res.r2, + "adj_r2" => res.adj_r2, + "estimated_seconds" => est.estimated_seconds, + "elapsed_seconds" => round(elapsed; digits=4), + "exceeds_five_minutes" => slow, + )) + end + end + end + + summary_doc = Dict{String,Any}( + "benchmark" => "constrained_ordination", + "timestamp_utc" => Dates.format(now(UTC), dateformat"yyyy-mm-ddTHH:MM:SSZ"), + "r_vegan_available" => (r_available() && vegan_available()), + "five_minute_threshold_seconds" => FIVE_MINUTE_WARN_SECONDS, + "warnings" => warnings, + "runs" => results, + ) + + if emit_json + println(String(jcs_bytes(summary_doc))) + else + println("=== Constrained Ordination Benchmark ===") + println("R + vegan available: ", summary_doc["r_vegan_available"]) + for row in results + println(" N=", row["tier_samples"], + " | ", row["engine"], + " | ", row["method"], + " | ", row["elapsed_seconds"], " s", + row["exceeds_five_minutes"] ? " [WARN > 5 min]" : "") + end + for w in warnings + println("WARNING: ", w) + end + end +end + +if abspath(PROGRAM_FILE) == @__FILE__ + main(ARGS) +end diff --git a/config/defaults/ordination.kyaml b/config/defaults/ordination.kyaml new file mode 100644 index 00000000..84f82833 --- /dev/null +++ b/config/defaults/ordination.kyaml @@ -0,0 +1,151 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Canonical constrained-ordination capability and compatibility contract (KYAML). +# Loaded by MetaManifold.Ordination to keep method/distance compatibility, +# scaling semantics, permutation bounds, and primary literature citations in one +# machine-readable source of truth. +--- +{ + schema_version: 1, + defaults: { + method: "none", + distance: "bray", + scaling: 2, + permutations: 999, + variance_partitioning: false, + replacement_method: "bayesian_multiplicative", + replacement_delta: 0.65, + }, + permutation_bounds: { + min: 99, + max: 9999, + default: 999, + }, + scalings: { + "1": { + id: 1, + name: "distance", + summary: "Preserves Euclidean/principal-coordinate distances among samples; project samples onto covariate arrows at right angles.", + }, + "2": { + id: 2, + name: "correlation", + summary: "Preserves approximate angles/correlations among features and constraining covariates; distances among sample points are not approximate dissimilarities.", + }, + }, + methods: { + none: { + vegan_function: null, + geometry: "unconstrained_only", + compatible_distances: [ + "bray", + "jaccard", + "euclidean", + "aitchison", + "chisq", + "unweighted_unifrac", + "weighted_unifrac", + ], + default_distance: "bray", + summary: "Constrained ordination disabled; use unconstrained NMDS and PERMANOVA.", + }, + rda: { + vegan_function: "rda", + geometry: "linear_euclidean", + compatible_distances: [ + "aitchison", + "euclidean", + ], + default_distance: "aitchison", + summary: "Linear Redundancy Analysis on zero-replaced CLR (Aitchison) or Euclidean/Hellinger space.", + }, + cca: { + vegan_function: "cca", + geometry: "unimodal_chisq", + compatible_distances: [ + "chisq", + ], + default_distance: "chisq", + summary: "Unimodal Canonical Correspondence Analysis weighted by chi-square distance on non-negative abundances.", + }, + cap: { + vegan_function: "capscale", + geometry: "metric_pcoa", + compatible_distances: [ + "bray", + "jaccard", + "euclidean", + "aitchison", + "unweighted_unifrac", + "weighted_unifrac", + ], + default_distance: "bray", + summary: "Constrained Analysis of Principal Coordinates (vegan::capscale) on an ecological or phylogenetic dissimilarity matrix.", + }, + dbrda: { + vegan_function: "dbrda", + geometry: "gower_distance_rda", + compatible_distances: [ + "bray", + "jaccard", + "euclidean", + "aitchison", + "unweighted_unifrac", + "weighted_unifrac", + ], + default_distance: "bray", + summary: "Distance-based Redundancy Analysis (Legendre & Anderson 1999; McArdle & Anderson 2001) directly decomposing Gower-centred dissimilarities.", + }, + }, + citations: [ + { + id: "legendre_anderson_1999", + authors: "Legendre, P. & Anderson, M. J.", + year: 1999, + title: "Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments", + journal: "Ecological Monographs 69(1), 1-24", + doi: "10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2", + }, + { + id: "mcardle_anderson_2001", + authors: "McArdle, B. H. & Anderson, M. J.", + year: 2001, + title: "Fitting multivariate models to community data: a comment on distance-based redundancy analysis", + journal: "Ecology 82(1), 290-297", + doi: "10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2", + }, + { + id: "anderson_willis_2003", + authors: "Anderson, M. J. & Willis, T. J.", + year: 2003, + title: "Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology", + journal: "Ecology 84(2), 511-525", + doi: "10.1890/0012-9658(2003)084[0511:CAOPCA]2.0.CO;2", + }, + { + id: "peres_neto_2006", + authors: "Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D.", + year: 2006, + title: "Variation partitioning of species data matrices: estimation and comparison of fractions", + journal: "Ecology 87(10), 2614-2625", + doi: "10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2", + }, + { + id: "gloor_2017", + authors: "Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J.", + year: 2017, + title: "Microbiome datasets are compositional: and this is not optional", + journal: "Frontiers in Microbiology 8, 2224", + doi: "10.3389/fmicb.2017.02224", + }, + { + id: "oksanen_vegan", + authors: "Oksanen, J. et al.", + year: 2025, + title: "vegan: Community Ecology Package", + journal: "CRAN R package", + doi: "10.32614/CRAN.package.vegan", + }, + ], +} diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index d9930ae9..3c75777c 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -343,6 +343,36 @@ analysis: normalisation: hellinger # Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample. depth: 0 + # Constrained ordination: projects multivariate community variation onto + # explicit study covariates via R/vegan (rda, cca, capscale, dbrda) with + # seeded permutation inference and optional variance partitioning (varpart). + ordination: + # Constrained ordination method: none (unconstrained NMDS + PERMANOVA only), + # rda, cca, cap or dbrda. + method: none + # RHS formula over sample covariates (e.g. "~ group" or "~ group + Condition(run)"); + # ~ uses all non-constant study covariates automatically. + formula: ~ + # Dissimilarity or geometry: bray, jaccard, euclidean, aitchison, chisq, + # unweighted_unifrac or weighted_unifrac. Method compatibility is enforced: + # rda accepts aitchison or euclidean; cca accepts chisq; cap and dbrda + # accept bray, jaccard, euclidean, aitchison, unweighted_unifrac or + # weighted_unifrac. + distance: bray + # Biplot scaling: 1 preserves inter-sample distances; 2 preserves approximate + # feature and covariate correlations. + scaling: 2 + # Permutations for anova.cca inference (whole number in [99, 9999]; seeded + # by the study's top-level seed). + permutations: 999 + # Compute vegan::varpart unique and shared adjusted-R2 fractions when 2 to 4 + # constraining covariates vary across the included samples. + variance_partitioning: false + # Zero-replacement policy when distance is aitchison (CLR): + # bayesian_multiplicative or multiplicative_replacement. + replacement_method: bayesian_multiplicative + # cmultRepl frac strictly in (0, 1) when distance is aitchison. + replacement_delta: 0.65 # Differential abundance: compares two conditions on the raw read counts, with # mandatory Benjamini-Hochberg adjusted p-values across the fitted family. differential: diff --git a/docs/schemas/constrained-ordination.schema.json b/docs/schemas/constrained-ordination.schema.json new file mode 100644 index 00000000..f8c2fe0c --- /dev/null +++ b/docs/schemas/constrained-ordination.schema.json @@ -0,0 +1 @@ +{"$defs":{"eigenvalueArray":{"items":{"additionalProperties":false,"properties":{"axis":{"type":"string"},"eigenvalue":{"minimum":0,"type":"number"},"proportion_constrained":{"minimum":0,"type":"number"},"proportion_total":{"minimum":0,"type":"number"}},"required":["axis","eigenvalue","proportion_total"],"type":"object"},"type":"array"},"termTestArray":{"items":{"additionalProperties":false,"properties":{"df":{"minimum":0,"type":"integer"},"f_statistic":{"type":["number","null"]},"p_value":{"type":["number","null"]},"term":{"type":"string"},"untestable_reason":{"type":["string","null"]},"variance":{"type":["number","null"]}},"required":["term","df","variance","f_statistic","p_value","untestable_reason"],"type":"object"},"type":"array"}},"$id":"https://github.com/hyperpolymath/MetaManifold-WebUI/blob/main/docs/schemas/constrained-ordination.schema.json","$schema":"https://json-schema.org/draft/2020-12/schema","additionalProperties":false,"description":"Canonical RFC 7493 I-JSON / RFC 8785 JCS response schema for MetaManifold constrained ordinations (RDA, CCA, CAP, dbRDA).","properties":{"adj_r2":{"type":["number","null"]},"artifact":{"additionalProperties":false,"properties":{"analysis_id":{"pattern":"^[0-9a-f]{64}$","type":"string"},"created":{"type":"boolean"}},"required":["analysis_id","created"],"type":"object"},"available_covariates":{"items":{"type":"string"},"type":"array"},"axis_labels":{"items":{"type":"string"},"maxItems":2,"minItems":2,"type":"array"},"axis_names":{"items":{"type":"string"},"maxItems":2,"minItems":2,"type":"array"},"axis_tests":{"$ref":"#/$defs/termTestArray"},"biplot":{"additionalProperties":false,"properties":{"data":{"items":{"type":"object"},"type":"array"},"layout":{"type":"object"}},"required":["data","layout"],"type":"object"},"biplot_arrows":{"items":{"additionalProperties":false,"properties":{"axis1":{"type":"number"},"axis2":{"type":"number"},"label":{"type":"string"},"length":{"minimum":0,"type":"number"}},"required":["label","axis1","axis2","length"],"type":"object"},"type":"array"},"cautions":{"items":{"type":"string"},"type":"array"},"centroids":{"items":{"additionalProperties":false,"properties":{"axis1":{"type":"number"},"axis2":{"type":"number"},"label":{"type":"string"}},"required":["label","axis1","axis2"],"type":"object"},"type":"array"},"condition_terms":{"items":{"type":"string"},"type":"array"},"conditioned_inertia":{"minimum":0,"type":"number"},"conditioned_proportion":{"minimum":0,"type":"number"},"constrained_eigenvalues":{"$ref":"#/$defs/eigenvalueArray"},"constrained_inertia":{"minimum":0,"type":"number"},"constrained_proportion":{"minimum":0,"type":"number"},"constrained_terms":{"items":{"type":"string"},"minItems":1,"type":"array"},"distance":{"enum":["bray","jaccard","euclidean","aitchison","chisq","unweighted_unifrac","weighted_unifrac"],"type":"string"},"formula":{"minLength":2,"type":"string"},"imaginary_inertia":{"minimum":0,"type":"number"},"marginal_tests":{"$ref":"#/$defs/termTestArray"},"method":{"enum":["rda","cca","cap","dbrda"],"type":"string"},"n_features":{"minimum":0,"type":"integer"},"n_samples":{"minimum":3,"type":"integer"},"omnibus_test":{"additionalProperties":false,"properties":{"df_model":{"minimum":1,"type":"integer"},"df_residual":{"minimum":1,"type":"integer"},"f_statistic":{"type":["number","null"]},"min_p_value":{"exclusiveMinimum":0,"maximum":1,"type":"number"},"p_value":{"type":["number","null"]},"permutations":{"maximum":9999,"minimum":1,"type":"integer"},"text":{"type":"string"},"untestable_reason":{"type":["string","null"]},"variance_model":{"type":["number","null"]},"variance_residual":{"type":["number","null"]}},"required":["df_model","df_residual","variance_model","variance_residual","f_statistic","p_value","untestable_reason","permutations","min_p_value","text"],"type":"object"},"provenance":{"additionalProperties":true,"properties":{"beta_normalisation":{"type":"string"},"blocked":{"type":"boolean"},"condition_terms":{"items":{"type":"string"},"type":"array"},"constrained_terms":{"items":{"type":"string"},"type":"array"},"distance":{"type":"string"},"dropped_zero_features":{"minimum":0,"type":"integer"},"formula":{"type":"string"},"method":{"enum":["rda","cca","cap","dbrda"],"type":"string"},"multiple_testing_policy":{"type":"string"},"n_features":{"minimum":0,"type":"integer"},"n_samples":{"minimum":3,"type":"integer"},"permutations_performed":{"minimum":1,"type":"integer"},"permutations_requested":{"maximum":9999,"minimum":99,"type":"integer"},"rng_kind":{"type":"string"},"scaling":{"enum":[1,2],"type":"integer"},"scaling_meaning":{"enum":["distance","correlation"],"type":"string"},"seed":{"type":"integer"},"software":{"additionalProperties":false,"properties":{"r_package":{"const":"vegan","type":"string"},"vegan_version":{"type":"string"}},"required":["r_package","vegan_version"],"type":"object"},"zero_replacement":{"type":["object","null"]}},"required":["method","distance","scaling","scaling_meaning","formula","permutations_requested","permutations_performed","seed","multiple_testing_policy","software"],"type":"object"},"r2":{"type":["number","null"]},"runtime_estimate":{"additionalProperties":false,"properties":{"estimated_seconds":{"minimum":0,"type":"number"},"tier":{"enum":["fast","moderate","heavy"],"type":"string"},"warn_five_minutes":{"type":"boolean"}},"required":["estimated_seconds","tier","warn_five_minutes"],"type":"object"},"scaling":{"enum":[1,2],"type":"integer"},"sites":{"items":{"additionalProperties":false,"properties":{"axis1":{"type":"number"},"axis2":{"type":"number"},"sample":{"type":"string"}},"required":["sample","axis1","axis2"],"type":"object"},"minItems":3,"type":"array"},"term_tests":{"$ref":"#/$defs/termTestArray"},"top_feature_scores":{"items":{"type":"object"},"type":"array"},"total_inertia":{"minimum":0,"type":"number"},"unconstrained_eigenvalues":{"$ref":"#/$defs/eigenvalueArray"},"unconstrained_inertia":{"minimum":0,"type":"number"},"unconstrained_proportion":{"minimum":0,"type":"number"},"variance_partitioning":{"type":["object","null"]}},"required":["method","distance","scaling","formula","constrained_terms","condition_terms","n_samples","n_features","total_inertia","constrained_inertia","unconstrained_inertia","r2","adj_r2","sites","biplot_arrows","omnibus_test","term_tests","provenance","biplot"],"title":"ConstrainedOrdinationResult","type":"object"} diff --git a/docs/statistics/constrained-ordinations.md b/docs/statistics/constrained-ordinations.md new file mode 100644 index 00000000..dc8895cb --- /dev/null +++ b/docs/statistics/constrained-ordinations.md @@ -0,0 +1,99 @@ + + + +# Constrained Ordinations (RDA, CCA, CAP, dbRDA) + +MetaManifold provides four constrained (canonical) ordination methods via pinned R/`vegan` (`MetaManifold.Ordination`), exposed in the **Diversity -> Advanced Analysis** panel and through `POST /api/v1/studies/{study}/analysis/ordination`. + +--- + +## 1. Constrained vs. Unconstrained Ordination + +| Property | Unconstrained Ordination (`NMDS`, `PCoA`, `PCA`) | Constrained Ordination (`RDA`, `CCA`, `CAP`, `dbRDA`) | +|---|---|---| +| Objective | Extract axes capturing the largest **total** variation or rank-order dissimilarity across samples, irrespective of sample metadata. | Project community variation onto the linear subspace spanned by **explicit study covariates** (`~ group + depth + ...`), decomposing total inertia into conditioned, constrained, and residual unconstrained components. | +| Role of covariates | Overlaid *post hoc* via colour/shape or tested separately via `PERMANOVA`. | Enters the eigendecomposition directly; canonical axes are linear combinations of predictors. | +| Pitfall to avoid | Subtle treatment effects can be masked by dominant unmodelled gradients (e.g. batch or inter-individual variation). | Samples often appear separated on canonical axes **even when $R^2_{\text{adj}}$ is tiny**; always inspect $R^2_{\text{adj}}$, residual unconstrained inertia, and permutation $p$-values before making biological claims. | + +Constrained ordination in MetaManifold is **explanatory / exploratory multivariate modelling**, not standalone confirmatory causal proof. + +--- + +## 2. Methods and Compositional Geometry + +Method/distance compatibility is enforced in `config/defaults/ordination.kyaml`, `MetaManifold.Validation`, `MetaManifold.Ordination`, config cascade `PATCH` checks, and the React UI: + +| Method | `vegan` function | Allowed `distance` values | Geometry and appropriate microbiome use | +|---|---|---|---| +| `dbrda` | `vegan::dbrda` | `bray`, `jaccard`, `euclidean`, `aitchison`, `unweighted_unifrac`, `weighted_unifrac` | Direct McArdle & Anderson (2001) Gower-centred decomposition of squared dissimilarities without complex imaginary coordinates. Recommended general-purpose constrained ordination for ecological or phylogenetic dissimilarities. | +| `cap` | `vegan::capscale` | `bray`, `jaccard`, `euclidean`, `aitchison`, `unweighted_unifrac`, `weighted_unifrac` | Metric PCoA followed by linear RDA on real principal coordinates (Anderson & Willis 2003), supporting weighted-average taxon projections. | +| `rda` | `vegan::rda` | `aitchison`, `euclidean` | Linear Euclidean ordination. Accepts zero-replaced centred log-ratio (`aitchison`) coordinates (Gloor et al. 2017) or Euclidean/Hellinger-transformed abundances. **Refuses raw non-Euclidean dissimilarities such as Bray-Curtis.** | +| `cca` | `vegan::cca` | `chisq` | Unimodal ordination weighted by intrinsic chi-square distance on non-negative abundances across broad environmental niche gradients. | + +### Zero Replacement for Aitchison (`clr`) Geometry + +When `distance = "aitchison"` and the count matrix contains zeros, counts are closed after zero replacement via `MetaManifold.Compositional.replace_zeros` (`bayesian_multiplicative` GBM or `multiplicative_replacement` CZM, `replacement_delta` strictly in `(0, 1)`) before `clr` transformation. Provenance records the exact method, `delta`, and number of imputed cells, alongside an explicit bias caution (see [Zero-handling notes](zero-handling.md)). + +### Phylogenetic UniFrac Distances + +When `distance` is `unweighted_unifrac` or `weighted_unifrac`, `POST /api/v1/studies/{study}/analysis/ordination` resolves an ASV Newick phylogeny from `cladistics_id`, the study's latest completed Cladistics build, or an explicit `tree_newick`, and computes exact branch-weighted distances via `MetaManifold.CladeCumulus.unifrac_matrix`. + +--- + +## 3. Covariate Formulas and Partial Ordination (`Condition(...)`) + +Formulas are validated by a strict covariate-only AST parser (`validate_formula_syntax` and `parse_ordination_formula`) that **never passes arbitrary R code to `eval`**: + +- Main effects: `~ group + log_library_size` +- Two-way interactions: `~ group * depth` or `~ group + depth + group:depth` +- Partial constrained ordination: `~ group + Condition(run)` or `~ group + ph + Condition(run + log_library_size)` + +In a partial ordination, covariates inside `Condition(...)` are partialled out via QR projection *before* fitting constrained canonical axes, isolating treatment or environmental variance after removing batch, run, or library-size effects. Degrees of freedom are checked before calling R (`df_constrained + df_conditioned < n_samples - 1`). + +--- + +## 4. Biplot Scalings (Scaling 1 vs. Scaling 2) + +- **Scaling 1 (`distance`)**: Preserves Euclidean / principal-coordinate distances among sample points. Project sample points at right angles onto quantitative covariate arrows to approximate sample ordering along that covariate. +- **Scaling 2 (`correlation`, default)**: Preserves cosines of angles among covariate arrows (and feature loading vectors) as approximate correlations; inter-sample point distances on the plot do not approximate exact ecological distances. + +--- + +## 5. Permutation Inference and Multiple-Testing Discipline + +- Significance is evaluated via `vegan::anova.cca` (`permutations` in `[99, 9999]`, default `999`), seeded deterministically (`RNGkind("Mersenne-Twister", "Inversion", "Rejection")` + study `seed`) under `with_r_lock` with `parallel = 1L`. +- Every run reports: + 1. **Omnibus model test** (`anova.cca(fit)`) + 2. **Sequential term tests** (`by = "terms"`, Type I — order-dependent) + 3. **Marginal term tests** (`by = "margin"`, Type III — order-independent unique variance when $\ge 2$ terms) + 4. **Optional canonical axis tests** (`by = "axis"`, when `test_axes = true`) +- When repeated-measures `individual` blocks are present, permutations are restricted within individuals (`permute::Plots(strata = blocks)`). Any term constant within every block cannot be tested by within-block permutation; its $p$-value is withheld (`null`) with an explicit `untestable_reason` rather than fabricated. +- **Multiple-testing boundary**: Omnibus and term-level `anova.cca` tests are multivariate permutation tests (not adjusted by Benjamini-Hochberg FDR). Per-taxon differential abundance tests in `MetaManifold.Differential` continue to enforce mandatory Benjamini-Hochberg FDR adjustment. + +--- + +## 6. Variance Partitioning (`vegan::varpart`) + +When `variance_partitioning = true` and 2 to 4 constraining covariates vary: + +- `vegan::varpart` decomposes community variation into unique testable fractions (e.g. `[a] = X1 | X2`, `[b] = X2 | X1`), shared overlap fractions (`[c]`, `0 df`, not individually permutation-testable), and residual unexplained variation (`[d]`) using Ezekiel's unbiased adjusted $R^2$ (Peres-Neto et al. 2006). +- Small **negative adjusted-$R^2$ fractions** can arise legitimately when predictors explain less variance than random normal variates or under suppressor relationships; MetaManifold reports them unclamped with an explanatory note so fractions continue to sum to `1.0`. + +--- + +## 7. Wire Contracts and Verification + +- Capability & compatibility contract (KYAML): `config/defaults/ordination.kyaml` (`GET /api/v1/ordination/contract`) +- Canonical response schema (RFC 7493 I-JSON + RFC 8785 JCS): `docs/schemas/constrained-ordination.schema.json` +- Independent Python stdlib reference oracle & canonical fixture: `test/references/check_ordination_references.py` and `test/fixtures/ordination_reference.json` + +--- + +## References + +1. Legendre, P. & Anderson, M. J. (1999). Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. *Ecological Monographs* 69(1), 1–24. `doi:10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2` +2. McArdle, B. H. & Anderson, M. J. (2001). Fitting multivariate models to community data: a comment on distance-based redundancy analysis. *Ecology* 82(1), 290–297. `doi:10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2` +3. Anderson, M. J. & Willis, T. J. (2003). Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology. *Ecology* 84(2), 511–525. `doi:10.1890/0012-9658(2003)084[0511:CAOPCA]2.0.CO;2` +4. Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D. (2006). Variation partitioning of species data matrices: estimation and comparison of fractions. *Ecology* 87(10), 2614–2625. `doi:10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2` +5. Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. (2017). Microbiome datasets are compositional: and this is not optional. *Frontiers in Microbiology* 8, 2224. `doi:10.3389/fmicb.2017.02224` +6. Oksanen, J. et al. (2025). *vegan: Community Ecology Package* (CRAN). `doi:10.32614/CRAN.package.vegan` diff --git a/frontend/src/api/client.ts b/frontend/src/api/client.ts index f1325275..4a89d575 100644 --- a/frontend/src/api/client.ts +++ b/frontend/src/api/client.ts @@ -24,6 +24,7 @@ import type { CladisticsCreateRequest, CladisticsSummary, CladisticsDetail, CladeCumulusSourcesResponse, CladeCumulusRequest, CladeCumulusResponse, CladeMoveValidationRequest, CladeMoveValidationResponse, + OrdinationRequest, OrdinationResult, } from './types' import type { FigureDoc } from '../figure/types' @@ -255,6 +256,12 @@ export const api = { post(`/api/v1/studies/${study}/analysis/nmds`, body), permanova: (study: string, body: ComparisonRequest) => post(`/api/v1/studies/${study}/analysis/permanova`, body), + /** Request a study ordination; reject on serialisation, network, HTTP or response-parsing errors. */ + ordination: (study: string, body: OrdinationRequest) => + post(`/api/v1/studies/${study}/analysis/ordination`, body), + /** Fetch the ordination capability/citation contract; request and JSON parsing errors propagate. */ + ordinationContract: () => + get>('/api/v1/ordination/contract'), capabilities: () => get<{ r_available: boolean }>('/api/v1/capabilities'), venn: (study: string, body: VennRequest) => post(`/api/v1/studies/${study}/analysis/venn`, body), diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 9b72935c..7ab20e58 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -1069,3 +1069,186 @@ export interface CladeMoveValidationResponse { proposed_root: CladeNodeData | null proposed_plotly?: { data: unknown[]; layout: Record } } + +export type OrdinationMethod = 'none' | 'rda' | 'cca' | 'cap' | 'dbrda' +export type ActiveOrdinationMethod = Exclude +export type OrdinationDistance = + | 'bray' + | 'jaccard' + | 'euclidean' + | 'aitchison' + | 'chisq' + | 'unweighted_unifrac' + | 'weighted_unifrac' +export type OrdinationScaling = 1 | 2 + +export interface OrdinationEigenvalue { + axis: string + eigenvalue: number + proportion_total: number + proportion_constrained?: number +} + +export interface OrdinationSiteScore { + sample: string + axis1: number + axis2: number +} + +export interface OrdinationBiplotArrow { + label: string + axis1: number + axis2: number + length: number +} + +export interface OrdinationCentroid { + label: string + axis1: number + axis2: number +} + +export interface OrdinationFeatureScore { + feature: string + axis1: number + axis2: number + length: number +} + +export interface OrdinationTermTest { + term: string + df: number + variance: number | null + f_statistic: number | null + p_value: number | null + untestable_reason: string | null +} + +export interface OrdinationOmnibusTest { + df_model: number + df_residual: number + variance_model: number | null + variance_residual: number | null + f_statistic: number | null + p_value: number | null + untestable_reason: string | null + permutations: number + min_p_value: number + text: string +} + +export interface VariancePartitionFraction { + label: string + df: number + r2: number | null + adj_r2: number | null + testable?: boolean +} + +export interface VariancePartitionResult { + status: 'computed' | 'skipped' + reason?: string + covariates: string[] + conditioned_on: string[] + individual_fractions: VariancePartitionFraction[] + gross_fractions: VariancePartitionFraction[] + has_negative_fraction: boolean + note: string | null +} + +export interface OrdinationProvenance { + method: ActiveOrdinationMethod + distance: OrdinationDistance + scaling: OrdinationScaling + scaling_meaning: 'distance' | 'correlation' + formula: string + constrained_terms: string[] + condition_terms: string[] + permutations_requested: number + permutations_performed: number + seed: number + rng_kind: string + blocked: boolean + beta_normalisation: string + n_samples: number + n_features: number + dropped_zero_features: number + zero_replacement: { + method: string + delta: number + zeros_replaced: number + } | null + phylogeny_source?: string + multiple_testing_policy: string + software: { + r_package: string + vegan_version: string + } +} + +export interface OrdinationRequest extends ComparisonRequest { + method: ActiveOrdinationMethod + distance?: OrdinationDistance + scaling?: OrdinationScaling + permutations?: number + formula?: string | null + variance_partitioning?: boolean + replacement_method?: 'bayesian_multiplicative' | 'multiplicative_replacement' + replacement_delta?: number + seed?: number + test_axes?: boolean + show_centroids?: boolean + show_features?: boolean + cladistics_id?: string | null + tree_newick?: string | null + persist?: boolean + covariates?: Record> +} + +export interface OrdinationResult { + method: ActiveOrdinationMethod + distance: OrdinationDistance + scaling: OrdinationScaling + formula: string + constrained_terms: string[] + condition_terms: string[] + available_covariates: string[] + n_samples: number + n_features: number + total_inertia: number + conditioned_inertia: number + constrained_inertia: number + unconstrained_inertia: number + imaginary_inertia: number + conditioned_proportion: number + constrained_proportion: number + unconstrained_proportion: number + r2: number | null + adj_r2: number | null + axis_names: string[] + axis_labels: string[] + constrained_eigenvalues: OrdinationEigenvalue[] + unconstrained_eigenvalues: OrdinationEigenvalue[] + sites: OrdinationSiteScore[] + biplot_arrows: OrdinationBiplotArrow[] + centroids: OrdinationCentroid[] + top_feature_scores: OrdinationFeatureScore[] + omnibus_test: OrdinationOmnibusTest + term_tests: OrdinationTermTest[] + marginal_tests: OrdinationTermTest[] + axis_tests: OrdinationTermTest[] + variance_partitioning: VariancePartitionResult | null + cautions: string[] + runtime_estimate: { + estimated_seconds: number + tier: 'fast' | 'moderate' | 'heavy' + warn_five_minutes: boolean + } + provenance: OrdinationProvenance + biplot: { data: Record[]; layout: Record } + artifact?: { + analysis_id: string + created: boolean + } +} + diff --git a/frontend/src/components/ConstrainedOrdinationPanel.tsx b/frontend/src/components/ConstrainedOrdinationPanel.tsx new file mode 100644 index 00000000..4baf9fa2 --- /dev/null +++ b/frontend/src/components/ConstrainedOrdinationPanel.tsx @@ -0,0 +1,656 @@ +// SPDX-License-Identifier: MPL-2.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +import { useEffect, useMemo, useRef, useState } from 'react' +import { api } from '../api/client' +import { errorMessage } from '../api/errorMessage' +import { saveBlob } from '../utils/download' +import { AnalysisChart } from './AnalysisChart' +import { useToast } from './Toast' +import type { + ActiveOrdinationMethod, + CladisticsSummary, + ColFilter, + ComparisonRunSpec, + OrdinationDistance, + OrdinationResult, + OrdinationScaling, + OrdinationTermTest, +} from '../api/types' +import { + ORDINATION_CITATIONS, + ORDINATION_DISTANCE_LABELS, + ORDINATION_METHODS_INFO, + estimateOrdinationRuntime, + formatOrdinationCsv, + getMethodInfo, + isDistanceCompatible, + resolveCompatibleDistance, + scalingDescription, + validateOrdinationFormulaSyntax, +} from './ordinationModel' + +const ALL_DISTANCES: readonly OrdinationDistance[] = [ + 'bray', + 'jaccard', + 'aitchison', + 'euclidean', + 'chisq', + 'weighted_unifrac', + 'unweighted_unifrac', +] + +/** Format a finite value to fixed decimal places; show an em dash for missing/non-finite values. */ +function fmtNum(v: number | null | undefined, digits = 4): string { + if (v === null || v === undefined || !Number.isFinite(v)) return '—' + return v.toFixed(digits) +} + +/** Display a proportion as a percentage to one decimal place, or an em dash when unavailable. */ +function fmtPct(v: number | null | undefined): string { + if (v === null || v === undefined || !Number.isFinite(v)) return '—' + return `${(100 * v).toFixed(1)}%` +} + +/** Render term tests with reasons for withheld p-values; return null when there are no rows. */ +function TermAnovaTable({ title, subtitle, rows }: { + title: string + subtitle: string + rows: OrdinationTermTest[] +}) { + if (rows.length === 0) return null + return ( +
+
{title}
+
{subtitle}
+
+ + + + + + + + + + + + {rows.map(r => ( + + + + + + + + ))} + +
TermdfVariance / InertiaPseudo-FPermutation p
{r.term}{r.df}{fmtNum(r.variance, 4)}{fmtNum(r.f_statistic, 3)} + {r.untestable_reason + ? withheld (blocked) + : fmtNum(r.p_value, 4)} +
+
+
+ ) +} + +/** Render expandable method/interpretation guidance and links to primary literature. */ +function OrdinationHelpDetails() { + return ( +
+ + Method selection, interpretation guide & primary literature + +
+

+ Constrained vs. unconstrained ordination: Unconstrained ordination (NMDS, PCoA, PCA) + finds axes of largest total community variation first and overlays groups afterwards. Constrained + ordination (dbRDA, CAP, RDA, CCA) projects community variation directly onto the + subspace spanned by explicit study covariates, separating explained inertia from residual unconstrained + variation. Always inspect adjusted R², residual unconstrained inertia, and permutation p-values: + visual separation on constrained axes alone is not confirmatory proof when R² is small. +

+

+ Compositional compatibility: Linear rda operates in Euclidean space and + therefore accepts zero-replaced CLR (aitchison; Gloor et al. 2017) or Euclidean/Hellinger + coordinates, refusing raw non-Euclidean dissimilarities such as Bray-Curtis. Use dbrda or{' '} + cap for Bray-Curtis, Jaccard, and UniFrac dissimilarities, and cca for + unimodal chi-square niche gradients. +

+

+ Multiple-testing discipline: Overall model and term tests use seeded{' '} + vegan::anova.cca permutation inference (not Benjamini-Hochberg FDR). Per-taxon differential + abundance tests in the Differential tab continue to enforce Benjamini-Hochberg adjustment. +

+
References
+
    + {ORDINATION_CITATIONS.map(c => ( +
  • + {c.authors}. {c.title}{' '} + + doi:{c.doi} + +
    • + ))} +
+
+
+ ) +} + +/** Render the biplot, inference and variance partitions, with a CSV/provenance download named for the study. */ +function OrdinationResultsView({ study, result }: { study: string; result: OrdinationResult }) { + /** Download the displayed result and provenance as a UTF-8 CSV file. */ + const downloadCsv = () => { + const csv = formatOrdinationCsv(result) + saveBlob( + new Blob([csv], { type: 'text/csv;charset=utf-8' }), + `${study}_${result.method}_${result.distance}_ordination.csv`, + ) + } + + const vp = result.variance_partitioning + const om = result.omnibus_test + + return ( +
+ + +
+
+
+ Constrained Ordination Summary ({result.method.toUpperCase()} · {result.distance} · scaling {result.scaling}) +
+ +
+ +
+
+
Constrained R² / Adj. R²
+
+ {fmtNum(result.r2, 3)} / {fmtNum(result.adj_r2, 3)} +
+
+
+
Omnibus Permutation Test
+
+ F({om.df_model}, {om.df_residual}) = {fmtNum(om.f_statistic, 2)}, p = {om.untestable_reason ? 'withheld' : fmtNum(om.p_value, 4)} +
+
+
+
Inertia Partition
+
+ Constrained {fmtPct(result.constrained_proportion)} + {result.conditioned_proportion > 0 ? ` · Conditioned ${fmtPct(result.conditioned_proportion)}` : ''} + {' · '}Residual {fmtPct(result.unconstrained_proportion)} +
+
+
+
Provenance
+
+ vegan {result.provenance.software.vegan_version} · {om.permutations} perms · seed {result.provenance.seed} +
+
+
+ + {result.cautions.length > 0 && ( +
+ {result.cautions.map((c, idx) => ( +
0 ? 4 : 0 }}>{c}
+ ))} +
+ )} + + + + + + {vp && ( +
+
Variance Partitioning (vegan::varpart)
+ {vp.status === 'skipped' ? ( +

{vp.reason}

+ ) : ( + <> + {vp.note && ( +
{vp.note}
+ )} +
+ + + + + + + + + + + + {vp.individual_fractions.map(f => ( + + + + + + + + ))} + +
FractiondfUnadjusted R²Adjusted R²Testability
{f.label}{f.df}{fmtNum(f.r2, 4)}{fmtNum(f.adj_r2, 4)}{f.testable ? 'Testable (partial ordination)' : 'Shared / residual (not permutation-testable)'}
+
+ + )} +
+ )} +
+
+ ) +} + +/** + * Render expandable constrained-ordination controls for the selected study runs + * and table. Opening the panel loads study defaults; users can save edits back + * to study config or request an analysis. `aggregate` is forwarded to the server. + * A false `rAvailable` disables running; null allows a request. Results are + * discarded when study, runs, table or aggregation changes. + */ +export function ConstrainedOrdinationPanel({ + study, + runs, + table, + aggregate, + rAvailable, +}: { + study: string + runs: ComparisonRunSpec[] + table: string + aggregate: boolean + rAvailable: boolean | null +}) { + const toast = useToast() + const [open, setOpen] = useState(false) + const [method, setMethod] = useState('dbrda') + const [distance, setDistance] = useState('bray') + const [scaling, setScaling] = useState(2) + const [permutations, setPermutations] = useState(999) + const [formula, setFormula] = useState('~ group') + const [variancePartitioning, setVariancePartitioning] = useState(false) + const [testAxes, setTestAxes] = useState(false) + const [showFeatures, setShowFeatures] = useState(false) + const [cladisticsId, setCladisticsId] = useState('') + const [cladisticsTrees, setCladisticsTrees] = useState([]) + const [savingDefaults, setSavingDefaults] = useState(false) + const [computing, setComputing] = useState(false) + const [result, setResult] = useState(null) + const reqCounter = useRef(0) + const loadedStudyCfg = useRef(null) + + useEffect(() => { + reqCounter.current++ + setResult(null) + setComputing(false) + }, [study, runs, table, aggregate]) + + useEffect(() => { + if (!open || loadedStudyCfg.current === study) return + loadedStudyCfg.current = study + api.config.getStudy(study).then(cfg => { + const mRaw = String(cfg['analysis.ordination.method']?.value ?? 'none') + const activeM: ActiveOrdinationMethod = + (mRaw === 'rda' || mRaw === 'cca' || mRaw === 'cap' || mRaw === 'dbrda') ? mRaw : 'dbrda' + const dRaw = String(cfg['analysis.ordination.distance']?.value ?? 'bray') as OrdinationDistance + const compatD = resolveCompatibleDistance(activeM, ALL_DISTANCES.includes(dRaw) ? dRaw : 'bray') + setMethod(activeM) + setDistance(compatD) + const scRaw = Number(cfg['analysis.ordination.scaling']?.value ?? 2) + setScaling(scRaw === 1 ? 1 : 2) + const permRaw = Number(cfg['analysis.ordination.permutations']?.value ?? 999) + if (Number.isInteger(permRaw) && permRaw >= 99 && permRaw <= 9999) setPermutations(permRaw) + const vpRaw = Boolean(cfg['analysis.ordination.variance_partitioning']?.value ?? false) + setVariancePartitioning(vpRaw) + const fRaw = cfg['analysis.ordination.formula']?.value + if (typeof fRaw === 'string' && fRaw.trim() !== '') setFormula(fRaw) + }).catch(() => { + if (loadedStudyCfg.current === study) loadedStudyCfg.current = null + }) + }, [open, study]) + + const needsPhylo = distance === 'weighted_unifrac' || distance === 'unweighted_unifrac' + useEffect(() => { + if (!open || !needsPhylo) return + api.cladistics.list(study) + .then(items => setCladisticsTrees(items.filter(t => t.built))) + .catch(() => setCladisticsTrees([])) + }, [open, needsPhylo, study]) + + /** Select a method and replace the distance only if it is incompatible. */ + const handleMethodChange = (next: ActiveOrdinationMethod) => { + setMethod(next) + setDistance(prev => resolveCompatibleDistance(next, prev)) + } + + const formulaCheck = useMemo(() => validateOrdinationFormulaSyntax(formula), [formula]) + const permValid = Number.isInteger(permutations) && permutations >= 99 && permutations <= 9999 + + const runtimeEst = useMemo(() => estimateOrdinationRuntime({ + method, + distance, + nSamples: Math.max(6, runs.length * 6), + nFeatures: 300, + permutations: permValid ? permutations : 999, + nTerms: Math.max(1, formulaCheck.constrainedCount), + variancePartitioning, + }), [method, distance, runs.length, permValid, permutations, formulaCheck.constrainedCount, variancePartitioning]) + + const availableCovariateChips = result?.available_covariates ?? ['group', 'run', 'log_library_size'] + + /** Start a blank RHS or append an additive term; syntax is checked separately. */ + const appendFormulaToken = (token: string) => { + setFormula(prev => { + const trimmed = prev.trim() + if (!trimmed || trimmed === '~') return `~ ${token}` + return `${trimmed} + ${token}` + }) + } + + /** Save valid form values as study defaults; report request failures through a toast. */ + const saveStudyDefaults = async () => { + if (!formulaCheck.valid || !permValid) return + setSavingDefaults(true) + try { + await api.config.patchStudy(study, { + 'analysis.ordination.method': method, + 'analysis.ordination.distance': distance, + 'analysis.ordination.scaling': scaling, + 'analysis.ordination.permutations': permutations, + 'analysis.ordination.variance_partitioning': variancePartitioning, + 'analysis.ordination.formula': formula.trim() || null, + }) + toast.success('Saved constrained ordination defaults to study config') + } catch (err) { + toast.error(`Could not save ordination defaults: ${errorMessage(err)}`) + } finally { + setSavingDefaults(false) + } + } + + /** + * Request an analysis for valid form values, showing only the latest response + * for the current inputs. Request failures clear the result and show a toast. + */ + const runOrdination = async () => { + if (!formulaCheck.valid || !permValid) return + const reqId = ++reqCounter.current + setComputing(true) + try { + const res = await api.analysis.ordination(study, { + table, + runs, + colFilters: {} as Record, + aggregate, + method, + distance, + scaling, + permutations, + formula: formula.trim() || null, + variance_partitioning: variancePartitioning, + test_axes: testAxes, + show_features: showFeatures, + cladistics_id: needsPhylo && cladisticsId ? cladisticsId : null, + }) + if (reqId === reqCounter.current) setResult(res) + } catch (err) { + if (reqId !== reqCounter.current) return + setResult(null) + toast.error(`Constrained ordination failed: ${errorMessage(err)}`) + } finally { + if (reqId === reqCounter.current) setComputing(false) + } + } + + const methodInfo = getMethodInfo(method) + + return ( +
setOpen((e.currentTarget as HTMLDetailsElement).open)} + > + + Advanced Analysis: Constrained Ordinations (RDA, CCA, CAP, dbRDA) + + +
+

+ {methodInfo.summary} {methodInfo.whenToUse} +

+ +
+ + + + + + + +
+ +
+ +
+ Insert term: + {availableCovariateChips.map(cov => ( + + ))} + + +
+ {!formulaCheck.valid && ( +
{formulaCheck.error}
+ )} + {!permValid && ( +
+ Permutations must be a whole number between 99 and 9999. +
+ )} +
+ + {needsPhylo && ( +
+ + {cladisticsTrees.length === 0 && ( +
+ No completed study Cladistics tree found yet; build an ASV phylogeny in the study Cladistics tab before running UniFrac ordination. +
+ )} +
+ )} + +
+ + + +
+ +
+ {scalingDescription(scaling)} +
+ +
+ + + + Estimated runtime: {runtimeEst.label} + {runtimeEst.warnFiveMinutes ? ' (> 5 min warning: consider fewer permutations for exploration)' : ''} + +
+ + {result && } + +
+
+ ) +} diff --git a/frontend/src/components/DiversityPanel.tsx b/frontend/src/components/DiversityPanel.tsx index ae84d450..54985df1 100644 --- a/frontend/src/components/DiversityPanel.tsx +++ b/frontend/src/components/DiversityPanel.tsx @@ -3,13 +3,14 @@ import { useEffect, useRef, useState } from 'react' import { api } from '../api/client' import { AnalysisChart } from './AnalysisChart' +import { ConstrainedOrdinationPanel } from './ConstrainedOrdinationPanel' import { useAlphaMetricSelection, AlphaMetricToggles, ALPHA_METRICS, extractAlphaPanel } from './alphaMetrics' import { useToast } from './Toast' import { errorMessage } from '../api/errorMessage' import type { AnalysisOption } from './annotationShared' import type { ColFilter, ComparisonRunSpec, PermanovaResult } from '../api/types' -/** Alpha comparison, NMDS and PERMANOVA across the runs/sub-groups in `runs`. */ +/** Alpha comparison, NMDS, PERMANOVA and constrained ordination across the runs/sub-groups in `runs`. */ export function DiversityPanel({ study, runs, option, aggregate }: { study: string runs: ComparisonRunSpec[] @@ -170,6 +171,16 @@ export function DiversityPanel({ study, runs, option, aggregate }: { )} )} + + {option && ( + + )} ) } diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts index a663d0f5..a6779b5c 100644 --- a/frontend/src/components/PipelineStages.test.ts +++ b/frontend/src/components/PipelineStages.test.ts @@ -52,3 +52,18 @@ test('advanced compositional settings retain their editors and descriptions', () expect(CONFIG_TYPES['analysis.differential.max_iter']).toEqual({ kind: 'int' }) expect(CONFIG_DESCRIPTIONS['analysis.differential.mn_reference_taxon']).toContain('reference') }) + +test('constrained ordination settings expose methods, distances, scaling, permutations and variance partitioning', () => { + expect(CONFIG_TYPES['analysis.ordination.method']).toEqual({ + kind: 'enum', options: ['none', 'dbrda', 'cap', 'rda', 'cca'], + }) + expect(CONFIG_TYPES['analysis.ordination.distance']).toEqual({ + kind: 'enum', + options: ['bray', 'jaccard', 'euclidean', 'aitchison', 'chisq', 'unweighted_unifrac', 'weighted_unifrac'], + }) + expect(CONFIG_TYPES['analysis.ordination.scaling']).toEqual({ kind: 'int' }) + expect(CONFIG_TYPES['analysis.ordination.permutations']).toEqual({ kind: 'int' }) + expect(CONFIG_TYPES['analysis.ordination.variance_partitioning']).toEqual({ kind: 'boolean' }) + expect(CONFIG_DESCRIPTIONS['analysis.ordination.distance']).toContain('rda + bray') +}) + diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index ddeba734..d9f61818 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -39,7 +39,7 @@ export const STAGE_CONFIG_PREFIXES: Record = { dada2_classify: ['dada2.taxonomy.', 'dada2.output.'], swarm: ['swarm.'], vsearch: ['vsearch.'], - analysis: ['analysis.alpha.', 'analysis.beta.', 'analysis.differential.', 'analysis.taxa_bar.'], + analysis: ['analysis.alpha.', 'analysis.beta.', 'analysis.ordination.', 'analysis.differential.', 'analysis.taxa_bar.'], phylogeny: ['phylogeny.threads', 'phylogeny.reference.align.', 'phylogeny.reference.trim.', 'phylogeny.reference.tree.', 'phylogeny.placement.align.', 'phylogeny.placement.trim.', 'phylogeny.placement.place.', 'phylogeny.placement.accumulate.'], @@ -150,6 +150,14 @@ 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.ordination.method': 'Constrained ordination method via R/vegan: none (unconstrained NMDS + PERMANOVA only), dbrda (distance-based RDA; Legendre & Anderson 1999), cap (capscale; Anderson & Willis 2003), rda (linear Redundancy Analysis on CLR/Aitchison or Euclidean space), or cca (unimodal Canonical Correspondence Analysis on chi-square geometry).', + 'analysis.ordination.formula': 'Right-hand-side covariate formula (e.g. "~ group" or "~ group + Condition(run)" for partial ordination); null uses all varying study covariates.', + 'analysis.ordination.distance': 'Dissimilarity or geometry for constrained ordination: bray, jaccard, euclidean, aitchison (zero-replaced CLR), chisq (CCA only), unweighted_unifrac, or weighted_unifrac. Incompatible method/distance pairs (such as rda + bray) are refused.', + 'analysis.ordination.scaling': 'Biplot scaling: 1 (distance scaling, preserving inter-sample distances) or 2 (correlation scaling, preserving approximate feature and covariate correlations).', + 'analysis.ordination.permutations': 'Permutations for vegan::anova.cca significance tests (whole number in [99, 9999]; seeded by the study top-level seed). Omnibus and term tests are permutation-based; BH remains mandatory for per-taxon tests.', + 'analysis.ordination.variance_partitioning': 'Run vegan::varpart adjusted-R2 variance partitioning to separate unique and shared contributions when 2 to 4 constraining covariates vary.', + 'analysis.ordination.replacement_method': 'Zero-replacement policy when distance is aitchison (CLR): bayesian_multiplicative (GBM) or multiplicative_replacement (CZM). Every zero replacement is biased.', + 'analysis.ordination.replacement_delta': 'cmultRepl frac strictly in (0, 1) when distance is aitchison (CLR).', '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. Both Wald tests treat the dispersion as known, so ml results always carry a calibration caveat (strongest in small groups; see bench/null_calibration).', @@ -260,6 +268,13 @@ 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.ordination.method': { kind: 'enum', options: ['none', 'dbrda', 'cap', 'rda', 'cca'] }, + 'analysis.ordination.distance': { kind: 'enum', options: ['bray', 'jaccard', 'euclidean', 'aitchison', 'chisq', 'unweighted_unifrac', 'weighted_unifrac'] }, + 'analysis.ordination.scaling': { kind: 'int' }, + 'analysis.ordination.permutations': { kind: 'int' }, + 'analysis.ordination.variance_partitioning': { kind: 'boolean' }, + 'analysis.ordination.replacement_method': { kind: 'enum', options: ['bayesian_multiplicative', 'multiplicative_replacement'] }, + 'analysis.ordination.replacement_delta': { kind: 'float', min: 0.01, max: 0.99, step: 0.05 }, 'analysis.differential.method': { kind: 'enum', options: [ 'nb_glm', 'clr_welch', 'zinb', 'hurdle_nb', 'ancom_bc', 'aldex2', 'songbird', 'multinomial', 'dirichlet_multinomial', diff --git a/frontend/src/components/ordinationModel.test.ts b/frontend/src/components/ordinationModel.test.ts new file mode 100644 index 00000000..f07d7aec --- /dev/null +++ b/frontend/src/components/ordinationModel.test.ts @@ -0,0 +1,185 @@ +// SPDX-License-Identifier: MPL-2.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +import { describe, expect, test } from 'bun:test' +import type { OrdinationResult } from '../api/types' +import { + ORDINATION_CITATIONS, + estimateOrdinationRuntime, + formatOrdinationCsv, + isDistanceCompatible, + resolveCompatibleDistance, + scalingDescription, + validateOrdinationFormulaSyntax, +} from './ordinationModel' + +describe('method and distance compatibility', () => { + test('RDA requires linear Euclidean or Aitchison CLR geometry and rejects Bray-Curtis', () => { + expect(isDistanceCompatible('rda', 'aitchison')).toBe(true) + expect(isDistanceCompatible('rda', 'euclidean')).toBe(true) + expect(isDistanceCompatible('rda', 'bray')).toBe(false) + expect(isDistanceCompatible('rda', 'weighted_unifrac')).toBe(false) + expect(resolveCompatibleDistance('rda', 'bray')).toBe('aitchison') + expect(resolveCompatibleDistance('rda', 'euclidean')).toBe('euclidean') + }) + + test('CCA requires chi-square geometry', () => { + expect(isDistanceCompatible('cca', 'chisq')).toBe(true) + expect(isDistanceCompatible('cca', 'bray')).toBe(false) + expect(resolveCompatibleDistance('cca', 'bray')).toBe('chisq') + }) + + test('CAP and dbRDA accept Bray-Curtis, Jaccard, Euclidean, Aitchison, and UniFrac', () => { + for (const m of ['cap', 'dbrda'] as const) { + expect(isDistanceCompatible(m, 'bray')).toBe(true) + expect(isDistanceCompatible(m, 'jaccard')).toBe(true) + expect(isDistanceCompatible(m, 'euclidean')).toBe(true) + expect(isDistanceCompatible(m, 'aitchison')).toBe(true) + expect(isDistanceCompatible(m, 'unweighted_unifrac')).toBe(true) + expect(isDistanceCompatible(m, 'weighted_unifrac')).toBe(true) + expect(isDistanceCompatible(m, 'chisq')).toBe(false) + } + }) +}) + +describe('formula syntax validation', () => { + test('accepts empty formula, main effects, interactions, and Condition(...) partial ordination', () => { + expect(validateOrdinationFormulaSyntax('')).toEqual({ + valid: true, error: null, constrainedCount: 1, hasCondition: false, + }) + expect(validateOrdinationFormulaSyntax('~ group + run')).toEqual({ + valid: true, error: null, constrainedCount: 2, hasCondition: false, + }) + expect(validateOrdinationFormulaSyntax('~ group * depth + Condition(run + batch)')).toEqual({ + valid: true, error: null, constrainedCount: 3, hasCondition: true, + }) + }) + + test('rejects R code injection, Condition-only formulas, and overlapping Condition covariates', () => { + expect(validateOrdinationFormulaSyntax('~ system(ls)').valid).toBe(false) + expect(validateOrdinationFormulaSyntax('y ~ group').valid).toBe(false) + expect(validateOrdinationFormulaSyntax('~ Condition(run)').valid).toBe(false) + expect(validateOrdinationFormulaSyntax('~ group + Condition(group)').valid).toBe(false) + expect(validateOrdinationFormulaSyntax('~ group +').valid).toBe(false) + }) +}) + +describe('runtime estimation, scaling, citations, and CSV export', () => { + test('estimates fast runtime on small studies and flags > 5 minute heavy runs', () => { + const small = estimateOrdinationRuntime({ + method: 'dbrda', + distance: 'bray', + nSamples: 24, + nFeatures: 120, + permutations: 999, + nTerms: 2, + variancePartitioning: true, + }) + expect(small.tier).toBe('fast') + expect(small.warnFiveMinutes).toBe(false) + + const huge = estimateOrdinationRuntime({ + method: 'cap', + distance: 'weighted_unifrac', + nSamples: 2600, + nFeatures: 4000, + permutations: 9999, + nTerms: 3, + variancePartitioning: true, + }) + expect(huge.warnFiveMinutes).toBe(true) + expect(huge.tier).toBe('heavy') + }) + + test('scaling descriptions and primary literature citations are complete', () => { + expect(scalingDescription(1)).toContain('Scaling 1') + expect(scalingDescription(2)).toContain('Scaling 2') + expect(ORDINATION_CITATIONS.map(c => c.doi)).toContain('10.3389/fmicb.2017.02224') + expect(ORDINATION_CITATIONS.map(c => c.doi)).toContain('10.1890/0012-9658(2003)084[0511:CAOPCA]2.0.CO;2') + }) + + test('formats full provenance and result tables into CSV', () => { + const mock: OrdinationResult = { + method: 'dbrda', + distance: 'bray', + scaling: 2, + formula: '~ group + Condition(run)', + constrained_terms: ['group'], + condition_terms: ['run'], + available_covariates: ['group', 'run'], + n_samples: 4, + n_features: 12, + total_inertia: 1.25, + conditioned_inertia: 0.25, + constrained_inertia: 0.55, + unconstrained_inertia: 0.45, + imaginary_inertia: 0, + conditioned_proportion: 0.2, + constrained_proportion: 0.44, + unconstrained_proportion: 0.36, + r2: 0.55, + adj_r2: 0.41, + axis_names: ['dbRDA1', 'MDS1'], + axis_labels: ['dbRDA1 (44.0% total, 100.0% constrained)', 'MDS1 (36.0% total)'], + constrained_eigenvalues: [{ axis: 'dbRDA1', eigenvalue: 0.55, proportion_total: 0.44, proportion_constrained: 1 }], + unconstrained_eigenvalues: [{ axis: 'MDS1', eigenvalue: 0.45, proportion_total: 0.36 }], + sites: [{ sample: 's1', axis1: -0.4, axis2: 0.1 }], + biplot_arrows: [{ label: 'groupTreatment', axis1: 0.85, axis2: 0.0, length: 0.85 }], + centroids: [{ label: 'group:Control', axis1: -0.5, axis2: 0.0 }], + top_feature_scores: [], + omnibus_test: { + df_model: 1, + df_residual: 1, + variance_model: 0.55, + variance_residual: 0.45, + f_statistic: 4.12, + p_value: 0.012, + untestable_reason: null, + permutations: 999, + min_p_value: 0.001, + text: 'Permutation test for dbrda', + }, + term_tests: [{ term: 'group', df: 1, variance: 0.55, f_statistic: 4.12, p_value: 0.012, untestable_reason: null }], + marginal_tests: [{ term: 'group', df: 1, variance: 0.55, f_statistic: 4.12, p_value: 0.012, untestable_reason: null }], + axis_tests: [], + variance_partitioning: { + status: 'computed', + covariates: ['group', 'depth'], + conditioned_on: ['run'], + individual_fractions: [{ label: '[a] = X1|X2', df: 1, r2: 0.32, adj_r2: 0.28, testable: true }], + gross_fractions: [], + has_negative_fraction: false, + note: null, + }, + cautions: ['Constrained ordination is exploratory/explanatory.'], + runtime_estimate: { estimated_seconds: 0.12, tier: 'fast', warn_five_minutes: false }, + provenance: { + method: 'dbrda', + distance: 'bray', + scaling: 2, + scaling_meaning: 'correlation', + formula: '~ group + Condition(run)', + constrained_terms: ['group'], + condition_terms: ['run'], + permutations_requested: 999, + permutations_performed: 999, + seed: 20260409, + rng_kind: 'Mersenne-Twister / Inversion / Rejection', + blocked: false, + beta_normalisation: 'hellinger', + n_samples: 4, + n_features: 12, + dropped_zero_features: 0, + zero_replacement: null, + multiple_testing_policy: 'omnibus_and_term_permutation_tests_not_bh; per_taxon_tests_require_bh', + software: { r_package: 'vegan', vegan_version: '2.7-3' }, + }, + biplot: { data: [], layout: {} }, + } + const csv = formatOrdinationCsv(mock) + expect(csv).toContain('# method=dbrda; distance=bray; scaling=2(correlation); formula=~ group + Condition(run)') + expect(csv).toContain('site_score,s1,-0.4,0.1,,,') + expect(csv).toContain('term_anova_sequential,group,1,0.55,4.12,0.012,') + expect(csv).toContain('varpart_fraction,[a] = X1|X2,1,,0.32,0.28,testable') + }) +}) diff --git a/frontend/src/components/ordinationModel.ts b/frontend/src/components/ordinationModel.ts new file mode 100644 index 00000000..9ed3b438 --- /dev/null +++ b/frontend/src/components/ordinationModel.ts @@ -0,0 +1,432 @@ +// SPDX-License-Identifier: MPL-2.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +import type { + ActiveOrdinationMethod, + OrdinationDistance, + OrdinationResult, + OrdinationScaling, +} from '../api/types' + +export interface OrdinationMethodInfo { + id: ActiveOrdinationMethod + label: string + veganFn: string + defaultDistance: OrdinationDistance + compatibleDistances: readonly OrdinationDistance[] + summary: string + whenToUse: string +} + +export const ORDINATION_METHODS_INFO: readonly OrdinationMethodInfo[] = [ + { + id: 'dbrda', + label: 'dbRDA (Distance-based RDA)', + veganFn: 'vegan::dbrda', + defaultDistance: 'bray', + compatibleDistances: ['bray', 'jaccard', 'euclidean', 'aitchison', 'unweighted_unifrac', 'weighted_unifrac'], + summary: 'Direct constrained decomposition of Gower-centred ecological or phylogenetic dissimilarities (Legendre & Anderson 1999; McArdle & Anderson 2001).', + whenToUse: 'Recommended general-purpose constrained ordination for Bray-Curtis, Jaccard, or UniFrac dissimilarities without complex imaginary axes.', + }, + { + id: 'cap', + label: 'CAP (Constrained Principal Coordinates)', + veganFn: 'vegan::capscale', + defaultDistance: 'bray', + compatibleDistances: ['bray', 'jaccard', 'euclidean', 'aitchison', 'unweighted_unifrac', 'weighted_unifrac'], + summary: 'Metric PCoA followed by linear RDA on real principal coordinates (Anderson & Willis 2003).', + whenToUse: 'Use when classical PCoA coordinate scores and weighted-average taxon projections on arbitrary dissimilarities are desired.', + }, + { + id: 'rda', + label: 'RDA (Redundancy Analysis)', + veganFn: 'vegan::rda', + defaultDistance: 'aitchison', + compatibleDistances: ['aitchison', 'euclidean'], + summary: 'Linear constrained ordination in Euclidean geometry (compositional CLR / Aitchison space or Hellinger-transformed Euclidean space).', + whenToUse: 'Best for compositional log-ratio (Aitchison) analysis (Gloor et al. 2017) or Hellinger transformation-based RDA with direct taxon loading vectors.', + }, + { + id: 'cca', + label: 'CCA (Canonical Correspondence Analysis)', + veganFn: 'vegan::cca', + defaultDistance: 'chisq', + compatibleDistances: ['chisq'], + summary: 'Unimodal constrained ordination weighted by chi-square distance on non-negative abundances.', + whenToUse: 'Use across broad environmental gradients where taxa exhibit bell-shaped unimodal niche optima rather than monotonic linear responses.', + }, +] as const + +export const ORDINATION_DISTANCE_LABELS: Record = { + bray: 'Bray-Curtis dissimilarity', + jaccard: 'Jaccard dissimilarity', + euclidean: 'Euclidean (uses study beta transform)', + aitchison: 'Aitchison (zero-replaced CLR)', + chisq: 'Chi-square (intrinsic to CCA)', + unweighted_unifrac: 'Unweighted UniFrac (phylogenetic)', + weighted_unifrac: 'Weighted UniFrac (phylogenetic)', +} + +export const ORDINATION_CITATIONS = [ + { + authors: 'Legendre, P. & Anderson, M. J. (1999)', + title: 'Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecological Monographs 69(1), 1–24.', + doi: '10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2', + }, + { + authors: 'McArdle, B. H. & Anderson, M. J. (2001)', + title: 'Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82(1), 290–297.', + doi: '10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2', + }, + { + authors: 'Anderson, M. J. & Willis, T. J. (2003)', + title: 'Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology. Ecology 84(2), 511–525.', + doi: '10.1890/0012-9658(2003)084[0511:CAOPCA]2.0.CO;2', + }, + { + authors: 'Peres-Neto, P. R., Legendre, P., Dray, S. & Borcard, D. (2006)', + title: 'Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87(10), 2614–2625.', + doi: '10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2', + }, + { + authors: 'Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. (2017)', + title: 'Microbiome datasets are compositional: and this is not optional. Frontiers in Microbiology 8, 2224.', + doi: '10.3389/fmicb.2017.02224', + }, + { + authors: 'Oksanen, J. et al. (2025)', + title: 'vegan: Community Ecology Package (CRAN).', + doi: '10.32614/CRAN.package.vegan', + }, +] as const + +const RESERVED_R_WORDS = new Set([ + 'if', 'else', 'repeat', 'while', 'function', 'for', 'in', 'next', 'break', + 'TRUE', 'FALSE', 'NULL', 'Inf', 'NaN', 'NA', 'NA_integer_', 'NA_real_', + 'NA_complex_', 'NA_character_', '.', '...', 'Condition', +]) + +/** Return method metadata, falling back to dbRDA for an unrecognised runtime value. */ +export function getMethodInfo(method: ActiveOrdinationMethod): OrdinationMethodInfo { + const found = ORDINATION_METHODS_INFO.find(m => m.id === method) + return found ?? ORDINATION_METHODS_INFO[0] +} + +/** Check the method's allowed geometries; an unrecognised method uses the dbRDA list. */ +export function isDistanceCompatible( + method: ActiveOrdinationMethod, + distance: OrdinationDistance, +): boolean { + return getMethodInfo(method).compatibleDistances.includes(distance) +} + +/** Keep the current distance when compatible; otherwise use the method's default. */ +export function resolveCompatibleDistance( + method: ActiveOrdinationMethod, + current: OrdinationDistance, +): OrdinationDistance { + const info = getMethodInfo(method) + return info.compatibleDistances.includes(current) ? current : info.defaultDistance +} + +/** Return interpretation guidance for distance scaling (1) or correlation scaling (2). */ +export function scalingDescription(scaling: OrdinationScaling): string { + return scaling === 1 + ? 'Scaling 1 (distance biplot): distances between sample points approximate their ecological/Euclidean distances; project sample points onto covariate arrows at right angles.' + : 'Scaling 2 (correlation biplot): cosines of angles between covariate arrows (and feature vectors) approximate correlations; distances between sample points are not approximate dissimilarities.' +} + +/** Split additive terms outside parentheses, returning an error for empty terms or unbalanced parentheses. */ +function splitTopLevelPlus(rhs: string): string[] | { error: string } { + const parts: string[] = [] + let depth = 0 + let start = 0 + for (let i = 0; i < rhs.length; i++) { + const ch = rhs[i] + if (ch === '(') depth++ + else if (ch === ')') { + depth-- + if (depth < 0) return { error: "Unbalanced ')' in formula" } + } else if (ch === '+' && depth === 0) { + const seg = rhs.slice(start, i).trim() + if (!seg) return { error: "Empty '+' term in formula" } + parts.push(seg) + start = i + 1 + } + } + if (depth !== 0) return { error: "Unbalanced '(' in formula" } + const tail = rhs.slice(start).trim() + if (!tail) return { error: "Trailing '+' in formula" } + parts.push(tail) + return parts +} + +/** Return a trimmed ASCII covariate identifier, or null for invalid or reserved names. */ +function validIdent(s: string): string | null { + const t = s.trim() + if (!/^[A-Za-z][A-Za-z0-9_.]*$/.test(t)) return null + if (RESERVED_R_WORDS.has(t) || /^\.\.[0-9]+$/.test(t)) return null + return t +} + +/** + * Client-side preflight validation matching `Ordination.validate_formula_syntax` + * on the backend. Empty/null formula is valid (uses study default covariates). + * Only two-way interactions are supported; `a*b` counts as `a`, `b` and `a:b`. + * On success, return the number of distinct expanded constraints and whether + * conditioning is present; invalid syntax returns an error message. Empty input + * uses a count of 1 for runtime estimation. This checks syntax only; metadata + * and model degrees of freedom are checked by the server. + */ +export function validateOrdinationFormulaSyntax(formula: string | null | undefined): { + valid: boolean + error: string | null + constrainedCount: number + hasCondition: boolean +} { + if (formula == null || formula.trim() === '') { + return { valid: true, error: null, constrainedCount: 1, hasCondition: false } + } + if (!/^[A-Za-z0-9_.+*:()~ \t]+$/.test(formula)) { + return { + valid: false, + error: "Formula may contain only covariate identifiers, '+', ':', '*', and Condition(...)", + constrainedCount: 0, + hasCondition: false, + } + } + let rhs = formula.trim() + if (rhs.includes('~')) { + const pieces = rhs.split('~') + if (pieces.length !== 2) { + return { valid: false, error: "Formula may contain at most one '~'", constrainedCount: 0, hasCondition: false } + } + if (pieces[0].trim() !== '') { + return { valid: false, error: "Formula must be right-hand-side only (start with '~')", constrainedCount: 0, hasCondition: false } + } + rhs = pieces[1].trim() + if (!rhs) { + return { valid: false, error: 'Formula right-hand side cannot be empty', constrainedCount: 0, hasCondition: false } + } + } + + const splitRes = splitTopLevelPlus(rhs) + if (!Array.isArray(splitRes)) { + return { valid: false, error: splitRes.error, constrainedCount: 0, hasCondition: false } + } + + const constrained: string[] = [] + const conditioned = new Set() + for (const term of splitRes) { + const condMatch = /^Condition\s*\((.*)\)$/.exec(term) + if (condMatch) { + const inner = condMatch[1].trim() + if (!inner || inner.includes('(') || inner.includes(')')) { + return { valid: false, error: 'Invalid Condition(...) clause', constrainedCount: 0, hasCondition: false } + } + for (const sub of inner.split('+')) { + const id = validIdent(sub) + if (!id) return { valid: false, error: `Invalid covariate '${sub.trim()}' in Condition(...)`, constrainedCount: 0, hasCondition: false } + conditioned.add(id) + } + } else if (term.includes('(') || term.includes(')')) { + return { valid: false, error: `Only Condition(...) calls are allowed in formulas; rejected '${term}'`, constrainedCount: 0, hasCondition: false } + } else if (term.includes('*') || term.includes(':')) { + const op = term.includes('*') ? '*' : ':' + const other = op === '*' ? ':' : '*' + if (term.includes(other)) { + return { valid: false, error: `Do not mix '*' and ':' in '${term}'`, constrainedCount: 0, hasCondition: false } + } + const pair = term.split(op) + if (pair.length !== 2) { + return { valid: false, error: `Only two-way interactions are supported; rejected '${term}'`, constrainedCount: 0, hasCondition: false } + } + const a = validIdent(pair[0]) + const b = validIdent(pair[1]) + if (!a || !b || a === b) { + return { valid: false, error: `Invalid interaction term '${term}'`, constrainedCount: 0, hasCondition: false } + } + if (op === '*') { + if (!constrained.includes(a)) constrained.push(a) + if (!constrained.includes(b)) constrained.push(b) + } + const inter = `${a}:${b}` + if (!constrained.includes(inter)) constrained.push(inter) + } else { + const id = validIdent(term) + if (!id) { + return { valid: false, error: `Invalid covariate identifier '${term}'`, constrainedCount: 0, hasCondition: false } + } + if (!constrained.includes(id)) constrained.push(id) + } + } + + if (constrained.length === 0) { + return { + valid: false, + error: 'Include at least one constraining covariate outside Condition(...)', + constrainedCount: 0, + hasCondition: conditioned.size > 0, + } + } + + for (const t of constrained) { + for (const part of t.split(':')) { + if (conditioned.has(part)) { + return { + valid: false, + error: `Covariate '${part}' cannot appear both as a constraining term and inside Condition(...)`, + constrainedCount: 0, + hasCondition: true, + } + } + } + } + + return { + valid: true, + error: null, + constrainedCount: constrained.length, + hasCondition: conditioned.size > 0, + } +} + +/** + * Estimate wall-clock seconds for UI guidance, with a display label and cost tier. + * Missing feature/term counts use 250/1; counts are truncated and clamped to at + * least 3 samples, 10 features, 99 permutations and 1 term. Tiers change at 5 and + * 60 seconds; the five-minute warning is set only above 300 seconds. This is a + * heuristic, not a measured duration or a validation of the request. + */ +export function estimateOrdinationRuntime(opts: { + method: ActiveOrdinationMethod + distance: OrdinationDistance + nSamples: number + nFeatures?: number + permutations: number + nTerms?: number + variancePartitioning?: boolean +}): { + estimatedSeconds: number + tier: 'fast' | 'moderate' | 'heavy' + warnFiveMinutes: boolean + label: string +} { + const n = Math.max(3, Math.trunc(opts.nSamples)) + const p = Math.max(10, Math.trunc(opts.nFeatures ?? 250)) + const b = Math.max(99, Math.trunc(opts.permutations)) + const k = Math.max(1, Math.trunc(opts.nTerms ?? 1)) + + const baseS = 0.06 + const prepS = (opts.distance === 'unweighted_unifrac' || opts.distance === 'weighted_unifrac') + ? 4.5e-7 * n * n * p + : opts.distance === 'aitchison' + ? 2.5e-7 * n * p + : 1.2e-7 * n * n * p + const isDistModel = opts.method === 'cap' || opts.method === 'dbrda' + const decompS = isDistModel ? 2.2e-9 * n * n * n : 1.5e-9 * n * p * Math.min(n, p) + const passes = k >= 2 ? 3 : 2 + const permUnit = isDistModel ? 2.8e-8 * n * n : 1.8e-8 * n * Math.min(p, 250) + const permS = passes * b * permUnit + const vpS = opts.variancePartitioning && k >= 2 + ? 0.35 * (decompS + (opts.method === 'cca' ? b * permUnit : 0.02)) + : 0 + + const est = Math.round((baseS + prepS + decompS + permS + vpS) * 100) / 100 + const tier = est >= 60 ? 'heavy' : est >= 5 ? 'moderate' : 'fast' + const label = est < 1 ? '< 1 s' : est < 60 ? `~${est.toFixed(1)} s` : `~${(est / 60).toFixed(1)} min` + return { + estimatedSeconds: est, + tier, + warnFiveMinutes: est > 300, + label, + } +} + +/** Escape commas, quotes and line breaks for one CSV cell; null/undefined become empty cells. */ +function csvCell(val: string | number | boolean | null | undefined): string { + if (val === null || val === undefined) return '' + const s = String(val) + if (/[",\n\r]/.test(s)) return `"${s.replace(/"/g, '""')}"` + return s +} + +/** + * Format a constrained ordination result into a self-documenting CSV carrying + * provenance headers, sample site scores, covariate biplot vectors, permutation + * ANOVA results, and optional variance partitioning fractions. + */ +export function formatOrdinationCsv(res: OrdinationResult): string { + const lines: string[] = [] + const prov = res.provenance + lines.push( + `# MetaManifold constrained ordination (${prov.software.r_package} ${prov.software.vegan_version})`, + ) + const zr = prov.zero_replacement + ? `; zero_replacement=${prov.zero_replacement.method}(delta=${prov.zero_replacement.delta},replaced=${prov.zero_replacement.zeros_replaced})` + : '' + lines.push( + `# method=${res.method}; distance=${res.distance}; scaling=${res.scaling}(${prov.scaling_meaning}); formula=${res.formula}; permutations=${prov.permutations_performed}; seed=${prov.seed}; beta_normalisation=${prov.beta_normalisation}${zr}`, + ) + lines.push( + `# total_inertia=${res.total_inertia}; constrained_inertia=${res.constrained_inertia}; r2=${res.r2 ?? ''}; adj_r2=${res.adj_r2 ?? ''}; omnibus_F=${res.omnibus_test.f_statistic ?? ''}; omnibus_p=${res.omnibus_test.p_value ?? ''}`, + ) + for (const c of res.cautions) { + lines.push(`# caution: ${c}`) + } + lines.push('section,name,axis1_or_df,axis2_or_variance,f_or_r2,p_or_adj_r2,note') + for (const s of res.sites) { + lines.push(['site_score', s.sample, s.axis1, s.axis2, '', '', ''].map(csvCell).join(',')) + } + for (const b of res.biplot_arrows) { + lines.push(['covariate_arrow', b.label, b.axis1, b.axis2, b.length, '', ''].map(csvCell).join(',')) + } + for (const c of res.centroids) { + lines.push(['centroid', c.label, c.axis1, c.axis2, '', '', ''].map(csvCell).join(',')) + } + lines.push([ + 'omnibus_anova', + 'Model', + res.omnibus_test.df_model, + res.omnibus_test.variance_model, + res.omnibus_test.f_statistic, + res.omnibus_test.p_value, + res.omnibus_test.untestable_reason ?? `residual_df=${res.omnibus_test.df_residual};permutations=${res.omnibus_test.permutations}`, + ].map(csvCell).join(',')) + for (const t of res.term_tests) { + lines.push([ + 'term_anova_sequential', + t.term, + t.df, + t.variance, + t.f_statistic, + t.p_value, + t.untestable_reason ?? '', + ].map(csvCell).join(',')) + } + for (const t of res.marginal_tests) { + lines.push([ + 'term_anova_marginal', + t.term, + t.df, + t.variance, + t.f_statistic, + t.p_value, + t.untestable_reason ?? '', + ].map(csvCell).join(',')) + } + if (res.variance_partitioning?.status === 'computed') { + for (const f of res.variance_partitioning.individual_fractions) { + lines.push([ + 'varpart_fraction', + f.label, + f.df, + '', + f.r2, + f.adj_r2, + f.testable ? 'testable' : 'untestable_shared_or_residual', + ].map(csvCell).join(',')) + } + } + return lines.join('\n') + '\n' +} diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 49547ae9..c2681b57 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -55,6 +55,7 @@ include("analysis/method_catalogue.jl") include("analysis/exact_summaries.jl") include("analysis/high_precision.jl") include("analysis/derived_analyses.jl") +include("analysis/ordination.jl") # Web server include("server/server.jl") diff --git a/src/analysis/ordination.jl b/src/analysis/ordination.jl new file mode 100644 index 00000000..c6e122c9 --- /dev/null +++ b/src/analysis/ordination.jl @@ -0,0 +1,2042 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +module Ordination + +using DataFrames, OrderedCollections, YAML, RCall +using ..RRuntime: with_r_lock, RBusyError +using ..Analysis: R_WAIT_SECONDS, r_available +using ..CladeCumulus: PhyloNode, parse_newick, unifrac_matrix +using ..DiversityMetrics: transform_counts +using ..Differential: clr_transform, replace_zeros + +export OrdinationConfig, ParsedFormula, VeganUnavailable, vegan_available, + ORDINATION_METHODS, ORDINATION_DISTANCES, ORDINATION_SCALINGS, + MIN_ORDINATION_PERMUTATIONS, MAX_ORDINATION_PERMUTATIONS, DEFAULT_ORDINATION_PERMUTATIONS, + compatible_distances, default_distance_for_method, validate_method_distance!, + validate_formula_syntax, parse_ordination_formula, ordination_config_from_dict, + load_ordination_contract, estimate_ordination_seconds, + pure_linear_rda, pure_ezekiel_adj_r2, pure_varpart_two_way, + run_constrained_ordination, ordination_biplot_chart + +const ORDINATION_METHODS = (:none, :rda, :cca, :cap, :dbrda) +const ACTIVE_ORDINATION_METHODS = (:rda, :cca, :cap, :dbrda) +const ORDINATION_DISTANCES = ( + :bray, :jaccard, :euclidean, :aitchison, :chisq, + :unweighted_unifrac, :weighted_unifrac, +) +const ORDINATION_SCALINGS = (1, 2) +const MIN_ORDINATION_PERMUTATIONS = 99 +const MAX_ORDINATION_PERMUTATIONS = 9999 +const DEFAULT_ORDINATION_PERMUTATIONS = 999 +const REPLACEMENT_METHODS = ("bayesian_multiplicative", "multiplicative_replacement") +const DEFAULT_REPLACEMENT_DELTA = 0.65 + +const _METHOD_DISTANCES = Dict{Symbol,Tuple{Vararg{Symbol}}}( + :none => ORDINATION_DISTANCES, + :rda => (:aitchison, :euclidean), + :cca => (:chisq,), + :cap => (:bray, :jaccard, :euclidean, :aitchison, :unweighted_unifrac, :weighted_unifrac), + :dbrda => (:bray, :jaccard, :euclidean, :aitchison, :unweighted_unifrac, :weighted_unifrac), +) + +const _DEFAULT_DISTANCE = Dict{Symbol,Symbol}( + :none => :bray, + :rda => :aitchison, + :cca => :chisq, + :cap => :bray, + :dbrda => :bray, +) + +const _R_FORMULA_RESERVED_NAMES = Set([ + "if", "else", "repeat", "while", "function", "for", "in", "next", "break", + "TRUE", "FALSE", "NULL", "Inf", "NaN", "NA", "NA_integer_", "NA_real_", + "NA_complex_", "NA_character_", ".", "...", "Condition", +]) + +""" +Thrown when R's `vegan` package cannot be loaded inside the shared R runtime. +""" +struct VeganUnavailable <: Exception end + +Base.showerror(io::IO, ::VeganUnavailable) = print(io, + "the R package vegan is not available, so constrained ordinations cannot be fitted; " * + "run install.sh or renv::restore() to install the pinned vegan release") + +_string_keyed(d::AbstractDict) = Dict{String,Any}(string(k) => _string_keyed(v) for (k, v) in d) +_string_keyed(v::AbstractVector) = Any[_string_keyed(x) for x in v] +_string_keyed(x) = x + +""" + load_ordination_contract() -> Dict{String,Any} + +Load the canonical KYAML capability and citation contract from +`config/defaults/ordination.kyaml`. File access and YAML parsing errors propagate. +""" +function load_ordination_contract()::Dict{String,Any} + contract_path = normpath(joinpath(@__DIR__, "..", "..", "config", "defaults", "ordination.kyaml")) + _string_keyed(YAML.load_file(contract_path)) +end + +"""Trim and lowercase a setting; throw `ArgumentError` labelled with `label` if empty.""" +function _norm_symbol(val::Union{Symbol,AbstractString}, label::AbstractString)::Symbol + s = lowercase(strip(String(val))) + isempty(s) && throw(ArgumentError("$label cannot be empty")) + Symbol(s) +end + +""" + compatible_distances(method) -> Tuple{Vararg{Symbol}} + +Return the allowed dissimilarity/geometry symbols for `method`, ignoring case +and surrounding whitespace. Throw `ArgumentError` for an unknown or empty method. +""" +function compatible_distances(method::Union{Symbol,AbstractString}) + m = _norm_symbol(method, "ordination_method") + haskey(_METHOD_DISTANCES, m) || throw(ArgumentError( + "ordination_method must be one of $(join(string.(ORDINATION_METHODS), ", ")), not '$method'")) + _METHOD_DISTANCES[m] +end + +""" + default_distance_for_method(method) -> Symbol + +Return the canonical default distance for `method`, ignoring case and surrounding +whitespace. Throw `ArgumentError` for an unknown or empty method. +""" +function default_distance_for_method(method::Union{Symbol,AbstractString})::Symbol + m = _norm_symbol(method, "ordination_method") + haskey(_DEFAULT_DISTANCE, m) || throw(ArgumentError( + "ordination_method must be one of $(join(string.(ORDINATION_METHODS), ", ")), not '$method'")) + _DEFAULT_DISTANCE[m] +end + +""" + validate_method_distance!(method, distance) -> (Symbol, Symbol) + +Enforce method/distance geometric compatibility, throwing `ArgumentError` on +unknown settings or incompatible combinations (such as linear RDA with raw +Bray-Curtis or CCA with non-chi-square distances). Return the normalised method +and distance symbols without mutating either argument. +""" +function validate_method_distance!(method::Union{Symbol,AbstractString}, + distance::Union{Symbol,AbstractString}) + m = _norm_symbol(method, "ordination_method") + m in ORDINATION_METHODS || throw(ArgumentError( + "analysis.ordination.method must be one of $(join(string.(ORDINATION_METHODS), ", ")), not '$method'")) + d = _norm_symbol(distance, "ordination_distance") + d in ORDINATION_DISTANCES || throw(ArgumentError( + "analysis.ordination.distance must be one of $(join(string.(ORDINATION_DISTANCES), ", ")), not '$distance'")) + allowed = _METHOD_DISTANCES[m] + if !(d in allowed) + hint = m === :rda ? + "RDA operates in linear Euclidean space (use 'aitchison' for CLR-transformed compositions or 'euclidean', or switch method to 'dbrda'/'cap' for '$d')" : + m === :cca ? + "CCA uses unimodal chi-square geometry exclusively (set distance to 'chisq')" : + "method '$m' accepts $(join(string.(allowed), ", "))" + throw(ArgumentError( + "analysis.ordination.distance='$d' is incompatible with analysis.ordination.method='$m': $hint")) + end + (m, d) +end + +""" + OrdinationConfig(; kwargs...) + +Validated configuration for constrained ordination (`rda`, `cca`, `cap`, `dbrda`, +or `:none` when disabled). + +The constructor rejects incompatible method/distance pairs, invalid formula +syntax, scaling outside `1:2`, permutations outside `99:9999`, an unknown +replacement method or a replacement fraction that is non-finite or outside +`(0, 1)` with `ArgumentError`. Scaling 1 is distance scaling; scaling 2 is +correlation scaling. An empty formula becomes `nothing` for automatic selection +later; metadata are not checked here. The distance defaults to `:bray` even when +another method is selected. +""" +struct OrdinationConfig + ordination_method :: Symbol + ordination_formula :: Union{String,Nothing} + ordination_distance :: Symbol + ordination_scaling :: Int + permutations :: Int + variance_partitioning :: Bool + replacement_method :: String + replacement_delta :: Float64 + seed :: Int + + function OrdinationConfig(; + ordination_method::Union{Symbol,AbstractString} = :none, + ordination_formula::Union{Nothing,AbstractString} = nothing, + ordination_distance::Union{Symbol,AbstractString} = :bray, + ordination_scaling::Integer = 2, + permutations::Integer = DEFAULT_ORDINATION_PERMUTATIONS, + variance_partitioning::Bool = false, + replacement_method::AbstractString = "bayesian_multiplicative", + replacement_delta::Real = DEFAULT_REPLACEMENT_DELTA, + seed::Integer = 42, + ) + ordination_scaling isa Bool && throw(ArgumentError( + "analysis.ordination.scaling must be integer 1 or 2, not Bool")) + permutations isa Bool && throw(ArgumentError( + "analysis.ordination.permutations must be an integer in [$MIN_ORDINATION_PERMUTATIONS, $MAX_ORDINATION_PERMUTATIONS], not Bool")) + seed isa Bool && throw(ArgumentError("seed must be an integer, not Bool")) + + m, d = validate_method_distance!(ordination_method, ordination_distance) + + sc = Int(ordination_scaling) + sc in ORDINATION_SCALINGS || throw(ArgumentError( + "analysis.ordination.scaling must be 1 (distance) or 2 (correlation), not $ordination_scaling")) + + nperm = Int(permutations) + (MIN_ORDINATION_PERMUTATIONS <= nperm <= MAX_ORDINATION_PERMUTATIONS) || throw(ArgumentError( + "analysis.ordination.permutations must be in [$MIN_ORDINATION_PERMUTATIONS, $MAX_ORDINATION_PERMUTATIONS], not $permutations")) + + rmeth = lowercase(strip(String(replacement_method))) + rmeth in REPLACEMENT_METHODS || throw(ArgumentError( + "analysis.ordination.replacement_method must be one of $(join(REPLACEMENT_METHODS, ", ")), not '$replacement_method'")) + + (isfinite(replacement_delta) && 0.0 < replacement_delta < 1.0) || throw(ArgumentError( + "analysis.ordination.replacement_delta must be finite and strictly in (0, 1), not $replacement_delta")) + + form_clean = if isnothing(ordination_formula) || isempty(strip(String(ordination_formula))) + nothing + else + s = String(strip(String(ordination_formula))) + validate_formula_syntax(s) + s + end + + new(m, form_clean, d, sc, nperm, variance_partitioning, rmeth, Float64(replacement_delta), Int(seed)) + end +end + +## Formula grammar parser and R-injection safety boundary +# +# R's `as.formula` + `model.frame` evaluates arbitrary call expressions embedded +# in a formula string. User-supplied formulas are therefore parsed strictly in +# Julia into an explicit identifier/interaction/Condition(...) AST, and the +# string handed to R is reconstructed exclusively from validated tokens. + +"""Return a trimmed ASCII covariate name, or throw `ArgumentError` for invalid or reserved names.""" +function _check_identifier(ident::AbstractString) + s = strip(String(ident)) + occursin(r"\A[A-Za-z][A-Za-z0-9_.]*\z", s) || throw(ArgumentError( + "Invalid covariate identifier '$ident' in ordination formula; " * + "only plain ASCII R identifiers starting with a letter are permitted")) + if s in _R_FORMULA_RESERVED_NAMES || occursin(r"\A\.\.[0-9]+\z", s) + throw(ArgumentError( + "Reserved identifier '$s' cannot be used as a covariate in an ordination formula")) + end + String(s) +end + +"""Split additive terms outside parentheses; throw `ArgumentError` for empty terms or unbalanced parentheses.""" +function _split_top_level_plus(rhs::AbstractString)::Vector{String} + parts = String[] + depth = 0 + start_idx = firstindex(rhs) + i = firstindex(rhs) + while i <= lastindex(rhs) + c = rhs[i] + if c == '(' + depth += 1 + elseif c == ')' + depth -= 1 + depth < 0 && throw(ArgumentError("Unbalanced ')' in ordination formula '$rhs'")) + elseif c == '+' && depth == 0 + seg = strip(rhs[start_idx:prevind(rhs, i)]) + isempty(seg) && throw(ArgumentError("Empty '+' term in ordination formula '$rhs'")) + push!(parts, String(seg)) + start_idx = nextind(rhs, i) + end + i = nextind(rhs, i) + end + depth == 0 || throw(ArgumentError("Unbalanced '(' in ordination formula '$rhs'")) + tail = strip(rhs[start_idx:lastindex(rhs)]) + isempty(tail) && throw(ArgumentError("Trailing '+' in ordination formula '$rhs'")) + push!(parts, String(tail)) + parts +end + +""" + ParsedFormula + +Validated AST representation of a constrained-ordination RHS formula. +""" +struct ParsedFormula + canonical_rhs :: String + constrained_terms :: Vector{String} + condition_terms :: Vector{String} + referenced_columns :: Vector{String} + df_constrained :: Int + df_conditioned :: Int +end + +""" + validate_formula_syntax(formula) -> NamedTuple + +Validate an ordination formula string against the safe identifier, interaction +(`:` / `*`), and `Condition(...)` grammar without requiring sample metadata. +Throws `ArgumentError` on malformed syntax, disallowed characters, LHS +expressions, arbitrary function calls, or missing constraining terms. +Return `canonical_rhs` without a leading `~`, deduplicated constrained and +condition terms, and referenced column names. Expand `a*b` to `a + b + a:b`; +only two-way interactions are accepted. A column cannot be both constrained +and conditioned. +""" +function validate_formula_syntax(formula::AbstractString) + raw_str = String(formula) + occursin(r"\A[A-Za-z0-9_.+*:()~ \t]+\z", raw_str) || throw(ArgumentError( + "Ordination formula contains disallowed characters; only covariate names, '+', ':', '*', and Condition(...) are allowed")) + s = strip(raw_str) + isempty(s) && throw(ArgumentError("Ordination formula cannot be empty")) + + if occursin('~', s) + tilde_parts = split(s, '~'; keepempty=true) + length(tilde_parts) == 2 || throw(ArgumentError( + "Ordination formula may contain at most one '~'")) + isempty(strip(tilde_parts[1])) || throw(ArgumentError( + "Ordination formula must be right-hand-side only (e.g. '~ group + Condition(run)'); left-hand side '$(strip(tilde_parts[1]))' is not permitted")) + s = strip(tilde_parts[2]) + isempty(s) && throw(ArgumentError("Ordination formula right-hand side cannot be empty")) + end + + raw_terms = _split_top_level_plus(s) + constrained = String[] + conditioned = String[] + referenced = OrderedSet{String}() + + for raw_t in raw_terms + m_cond = match(r"\ACondition\s*\((.*)\)\z", raw_t) + if !isnothing(m_cond) + inner = strip(m_cond.captures[1]) + isempty(inner) && throw(ArgumentError("Condition() in ordination formula cannot be empty")) + if occursin('(', inner) || occursin(')', inner) + throw(ArgumentError("Nested parentheses inside Condition(...) are not permitted")) + end + for sub in split(inner, '+') + id = _check_identifier(sub) + id in conditioned || push!(conditioned, id) + push!(referenced, id) + end + elseif occursin('(', raw_t) || occursin(')', raw_t) + throw(ArgumentError( + "Only Condition(...) calls are permitted in ordination formulas; rejected '$raw_t'")) + elseif occursin('*', raw_t) || occursin(':', raw_t) + op = occursin('*', raw_t) ? '*' : ':' + other = op == '*' ? ':' : '*' + occursin(other, raw_t) && throw(ArgumentError( + "Mixed '*' and ':' in single term '$raw_t' is not supported; expand terms explicitly")) + pieces = split(raw_t, op) + length(pieces) == 2 || throw(ArgumentError( + "Only two-way interactions ('a:b' or 'a*b') are supported in ordination formulas; got '$raw_t'")) + a = _check_identifier(pieces[1]) + b = _check_identifier(pieces[2]) + a == b && throw(ArgumentError("Self-interaction '$raw_t' is invalid")) + for id in (a, b) + push!(referenced, id) + end + if op == '*' + a in constrained || push!(constrained, a) + b in constrained || push!(constrained, b) + inter = "$a:$b" + inter in constrained || push!(constrained, inter) + else + inter = "$a:$b" + inter in constrained || push!(constrained, inter) + end + else + id = _check_identifier(raw_t) + id in constrained || push!(constrained, id) + push!(referenced, id) + end + end + + isempty(constrained) && throw(ArgumentError( + "Constrained ordination formula must include at least one constraining covariate outside Condition(...)")) + + # Reject covariates appearing both as constrained predictors and inside Condition(...) + cond_set = Set(conditioned) + for t in constrained + for part in split(t, ':') + if part in cond_set + throw(ArgumentError( + "Covariate '$part' appears both as a constraining predictor and inside Condition(...), making the partial model aliased")) + end + end + end + + rhs_pieces = copy(constrained) + if !isempty(conditioned) + push!(rhs_pieces, "Condition(" * join(conditioned, " + ") * ")") + end + (; canonical_rhs = join(rhs_pieces, " + "), + constrained_terms = constrained, + condition_terms = conditioned, + referenced_columns = collect(referenced)) +end + +""" +Return 0 for a constant column, 1 for a varying numeric column, or levels minus +1 for categorical values. Throw `ArgumentError` for missing +values, non-finite numeric entries or empty categorical labels. +""" +function _column_df(col::AbstractVector)::Int + if eltype(col) <: Union{Real,Missing} && !(eltype(col) <: Bool) + vals = Float64[] + for v in col + ismissing(v) && throw(ArgumentError("Missing values are not allowed in ordination metadata")) + isfinite(v) || throw(ArgumentError("Non-finite numeric values are not allowed in ordination metadata")) + push!(vals, Float64(v)) + end + length(unique(vals)) >= 2 || return 0 + return 1 + else + vals = String[] + for v in col + (ismissing(v) || isnothing(v)) && throw(ArgumentError("Missing values are not allowed in ordination metadata")) + s = strip(string(v)) + isempty(s) && throw(ArgumentError("Empty categorical labels are not allowed in ordination metadata")) + push!(vals, s) + end + k = length(unique(vals)) + return max(0, k - 1) + end +end + +"""Estimate a term's degrees of freedom, multiplying column degrees of freedom for an interaction.""" +function _term_df(term::String, metadata::DataFrame)::Int + if occursin(':', term) + parts = split(term, ':') + return _column_df(metadata[!, parts[1]]) * _column_df(metadata[!, parts[2]]) + end + _column_df(metadata[!, term]) +end + +""" + parse_ordination_formula(formula, metadata) -> ParsedFormula + +Resolve and validate `formula` (or auto-construct a default formula over varying +metadata columns except `sample` when the formula is `nothing` or blank, +deduplicating identical numeric columns and equivalent categorical partitions) +against `metadata` (`DataFrame` with one row per sample). Verifies column +existence, non-constancy, absence of missing/non-finite entries, and positive +residual degrees of freedom (`df_total < n_samples - 1`). +Throw `ArgumentError` when fewer than three samples are supplied or any formula +or metadata check fails. Degrees of freedom are estimated from individual +columns and interactions; this does not establish full design-matrix rank. +""" +function parse_ordination_formula(formula::Union{Nothing,AbstractString}, + metadata::DataFrame)::ParsedFormula + n_samples = nrow(metadata) + n_samples >= 3 || throw(ArgumentError( + "Constrained ordination requires at least 3 samples (got $n_samples)")) + + avail = Set(names(metadata)) + if isnothing(formula) || isempty(strip(String(formula))) + auto_cols = String[] + seen_signatures = Set{Any}() + for col_name in names(metadata) + lowercase(col_name) == "sample" && continue + _check_identifier(col_name) + col_vec = metadata[!, col_name] + _column_df(col_vec) >= 1 || continue + sig = if eltype(col_vec) <: Union{Real,Missing} && !(eltype(col_vec) <: Bool) + (:num, Float64[Float64(v) for v in col_vec]) + else + # Canonical partition signature (relabelled 1..k in order of first appearance) + lbl_map = Dict{String,Int}() + (:cat, Int[get!(lbl_map, strip(string(v)), length(lbl_map) + 1) for v in col_vec]) + end + sig in seen_signatures && continue + push!(seen_signatures, sig) + push!(auto_cols, col_name) + end + isempty(auto_cols) && throw(ArgumentError( + "No metadata covariate has at least 2 distinct values across the $n_samples samples")) + formula = "~ " * join(auto_cols, " + ") + end + + syn = validate_formula_syntax(String(formula)) + for col_name in syn.referenced_columns + col_name in avail || throw(ArgumentError( + "Formula covariate '$col_name' is not in sample metadata (available: $(join(sort(collect(avail)), ", ")))")) + df_c = _column_df(metadata[!, col_name]) + df_c >= 1 || throw(ArgumentError( + "Covariate '$col_name' is constant across all $n_samples samples and has 0 degrees of freedom")) + end + + df_con = sum(_term_df(t, metadata) for t in syn.constrained_terms; init=0) + df_cnd = sum(_term_df(t, metadata) for t in syn.condition_terms; init=0) + if df_con + df_cnd >= n_samples - 1 + throw(ArgumentError( + "Ordination model is saturated (constrained df=$df_con + conditioned df=$df_cnd >= n_samples - 1 = $(n_samples - 1)), leaving 0 residual degrees of freedom for permutation testing")) + end + + ParsedFormula( + syn.canonical_rhs, + syn.constrained_terms, + syn.condition_terms, + syn.referenced_columns, + df_con, + df_cnd, + ) +end + +""" + ordination_config_from_dict(cfg; overrides=Dict()) -> OrdinationConfig + +Construct a validated `OrdinationConfig` from a merged `pipeline.yml` config +dictionary, applying optional request-level `overrides`. +Both dictionaries use string keys; override keys are names such as `"method"`. +Null overrides inherit configured values, except a null `"formula"` explicitly +requests automatic selection. The seed comes from the top-level config unless +overridden. An inherited distance is retained and checked for compatibility; +the method's default is used only when no distance is supplied in either scope. +Configuration validation and conversion errors propagate. +""" +function ordination_config_from_dict(cfg::AbstractDict; + overrides::AbstractDict = Dict{String,Any}())::OrdinationConfig + analysis_blk = get(cfg, "analysis", Dict{String,Any}()) + ord_blk = analysis_blk isa AbstractDict ? get(analysis_blk, "ordination", Dict{String,Any}()) : Dict{String,Any}() + ord_blk = ord_blk isa AbstractDict ? ord_blk : Dict{String,Any}() + + pick(k, default) = haskey(overrides, k) && !isnothing(overrides[k]) ? overrides[k] : + haskey(ord_blk, k) && !isnothing(ord_blk[k]) ? ord_blk[k] : default + + raw_formula = haskey(overrides, "formula") ? overrides["formula"] : get(ord_blk, "formula", nothing) + top_seed = haskey(overrides, "seed") && !isnothing(overrides["seed"]) ? overrides["seed"] : get(cfg, "seed", 42) + + method_val = pick("method", "none") + # If the caller overrides `method` without explicitly overriding `distance`, + # and the inherited config distance was left on a different method's default, + # only auto-select when `distance` was not explicitly supplied in `overrides`. + dist_val = if haskey(overrides, "distance") && !isnothing(overrides["distance"]) + overrides["distance"] + elseif haskey(ord_blk, "distance") && !isnothing(ord_blk["distance"]) + ord_blk["distance"] + else + string(default_distance_for_method(string(method_val))) + end + + OrdinationConfig(; + ordination_method = string(method_val), + ordination_formula = isnothing(raw_formula) ? nothing : string(raw_formula), + ordination_distance = string(dist_val), + ordination_scaling = pick("scaling", 2), + permutations = pick("permutations", DEFAULT_ORDINATION_PERMUTATIONS), + variance_partitioning = pick("variance_partitioning", false), + replacement_method = string(pick("replacement_method", "bayesian_multiplicative")), + replacement_delta = pick("replacement_delta", DEFAULT_REPLACEMENT_DELTA), + seed = top_seed, + ) +end + +""" + estimate_ordination_seconds(; method, distance, n_samples, n_features, permutations, n_terms=1, variance_partitioning=false) -> NamedTuple + +Analytical wall-clock cost heuristic for progressive-disclosure UI feedback and +preflight warnings. Returns `(; estimated_seconds, tier, warn_five_minutes)`. +Clamp each count to at least 1. Tiers are `"fast"` below 5 seconds, `"moderate"` +below 60 seconds and `"heavy"` otherwise; warn only above 300 seconds. This does +not measure execution time or validate method/distance compatibility. +""" +function estimate_ordination_seconds(; + method::Union{Symbol,AbstractString} = :dbrda, + distance::Union{Symbol,AbstractString} = :bray, + n_samples::Integer, + n_features::Integer, + permutations::Integer = DEFAULT_ORDINATION_PERMUTATIONS, + n_terms::Integer = 1, + variance_partitioning::Bool = false, +) + m = _norm_symbol(method, "method") + d = _norm_symbol(distance, "distance") + n = max(1, Int(n_samples)) + p = max(1, Int(n_features)) + b = max(1, Int(permutations)) + k = max(1, Int(n_terms)) + + # Base R/RCall dispatch + matrix marshalling overhead (~0.06s) + base_s = 0.06 + # Distance / CLR preparation cost + prep_s = if d in (:unweighted_unifrac, :weighted_unifrac) + 4.5e-7 * n * n * p + elseif d === :aitchison + 2.5e-7 * n * p + else + 1.2e-7 * n * n * p + end + # Eigendecomposition / SVD cost (Gower N x N for cap/dbrda vs N x min(N, P) for rda/cca) + decomp_s = (m in (:cap, :dbrda)) ? (2.2e-9 * n^3) : (1.5e-9 * n * p * min(n, p)) + # Permutation passes (omnibus + by="terms" + optional by="margin") + passes = k >= 2 ? 3 : 2 + perm_unit = (m in (:cap, :dbrda)) ? (2.8e-8 * n * n) : (1.8e-8 * n * min(p, 250)) + perm_s = passes * b * perm_unit + vp_s = (variance_partitioning && k >= 2) ? (0.35 * (decomp_s + (m === :cca ? b * perm_unit : 0.02))) : 0.0 + + est = round(base_s + prep_s + decomp_s + perm_s + vp_s; digits=2) + tier = est >= 60.0 ? "heavy" : est >= 5.0 ? "moderate" : "fast" + (; estimated_seconds = est, tier, warn_five_minutes = est > 300.0) +end + +## Pure-Julia mathematical building blocks (zero extra stdlibs; cross-checkable without R) + +""" + pure_ezekiel_adj_r2(r2, n_samples, df_model; df_conditioned=0) -> Float64 + +Analytical Ezekiel (1930) / Peres-Neto et al. (2006) adjusted R²: +`1 - (1 - r2) * (n - 1 - df_cond) / (n - 1 - df_cond - df_model)` when +conditioning removes `df_cond` degrees of freedom. With `df_cond == 0`, this +matches `vegan::RsquareAdj` for unconditioned linear RDA, CAP and dbRDA: +`1 - (1 - r2) * (n - 1) / (n - 1 - df_model)`. +For `df_conditioned > 0`, `r2` is the fraction of inertia remaining after +conditioning; this is not vegan's partial-model R² relative to total inertia. +Throw `ArgumentError` when residual degrees of freedom are not positive. +""" +function pure_ezekiel_adj_r2(r2::Real, n_samples::Integer, df_model::Integer; + df_conditioned::Integer = 0)::Float64 + n = Int(n_samples) + k = Int(df_model) + q = Int(df_conditioned) + df_res = n - 1 - q - k + df_res > 0 || throw(ArgumentError("Residual degrees of freedom must be positive (got $df_res)")) + df_tot = n - 1 - q + 1.0 - (1.0 - Float64(r2)) * (Float64(df_tot) / Float64(df_res)) +end + +""" + pure_varpart_two_way(adj_r2_a, adj_r2_b, adj_r2_ab) -> NamedTuple + +Exact two-table Peres-Neto et al. (2006) adjusted-R² variance partitioning into +unique `[a]` (`X1 | X2`), shared `[b]`, unique `[c]` (`X2 | X1`), and `[d]` +(unexplained residual `1 - adj_r2_ab`). +Return `unique_x1`, `shared`, `unique_x2` and `residual` without clamping negative +fractions. +""" +function pure_varpart_two_way(adj_r2_a::Real, adj_r2_b::Real, adj_r2_ab::Real) + a_tot = Float64(adj_r2_a) + b_tot = Float64(adj_r2_b) + ab_tot = Float64(adj_r2_ab) + uniq_a = ab_tot - b_tot + uniq_b = ab_tot - a_tot + shared = a_tot + b_tot - ab_tot + resid = 1.0 - ab_tot + (; unique_x1 = uniq_a, shared = shared, unique_x2 = uniq_b, residual = resid) +end + +""" +Approximate eigenvalues of the symmetric part of square `S` in descending order, +retaining only values above the absolute/relative cutoff. Stop after +`max_sweeps` even if off-diagonal entries still exceed `tol`; leave `S` unchanged. +""" +function _sym_eigenvalues_desc(S::Matrix{Float64}; tol::Float64 = 1e-13, max_sweeps::Int = 80)::Vector{Float64} + m = size(S, 1) + m == 0 && return Float64[] + m == 1 && return S[1, 1] > tol ? [S[1, 1]] : Float64[] + A = 0.5 .* (S .+ permutedims(S)) + for _ in 1:max_sweeps + max_off = 0.0 + for q in 2:m, p in 1:(q - 1) + v = abs(A[p, q]) + v > max_off && (max_off = v) + end + max_off <= tol && break + for p in 1:(m - 1), q in (p + 1):m + apq = A[p, q] + abs(apq) <= tol && continue + tau = (A[q, q] - A[p, p]) / (2.0 * apq) + t = tau >= 0.0 ? (1.0 / (tau + sqrt(1.0 + tau * tau))) : (-1.0 / (-tau + sqrt(1.0 + tau * tau))) + c = 1.0 / sqrt(1.0 + t * t) + s = t * c + app = A[p, p] + aqq = A[q, q] + A[p, p] = app - t * apq + A[q, q] = aqq + t * apq + A[p, q] = 0.0 + A[q, p] = 0.0 + for r in 1:m + (r == p || r == q) && continue + arp = A[r, p] + arq = A[r, q] + A[r, p] = c * arp - s * arq + A[p, r] = A[r, p] + A[r, q] = s * arp + c * arq + A[q, r] = A[r, q] + end + end + end + vals = [A[i, i] for i in 1:m] + sort!(vals; rev=true) + cutoff = isempty(vals) ? tol : max(tol, 1e-12 * max(1.0, vals[1])) + [v for v in vals if v > cutoff] +end + +"""Return retained positive Gram eigenvalues of `M`, divided by `scale_dof`, in descending order.""" +function _gram_eigenvalues(M::Matrix{Float64}, scale_dof::Float64)::Vector{Float64} + n, p = size(M) + G = n <= p ? ((M * permutedims(M)) ./ scale_dof) : ((permutedims(M) * M) ./ scale_dof) + _sym_eigenvalues_desc(G) +end + +"""Return an orthonormal column basis without changing `X`, dropping residual columns of norm at most `tol`.""" +function _modified_gram_schmidt(X::Matrix{Float64}; tol::Float64 = 1e-11)::Matrix{Float64} + n, m = size(X) + cols = Vector{Float64}[] + for j in 1:m + v = copy(X[:, j]) + # Two-pass MGS for numerical orthogonality + for _ in 1:2 + for q in cols + proj = 0.0 + for i in 1:n + proj += q[i] * v[i] + end + for i in 1:n + v[i] -= proj * q[i] + end + end + end + nrm = sqrt(sum(x * x for x in v; init=0.0)) + if nrm > tol + push!(cols, v ./ nrm) + end + end + isempty(cols) && return zeros(Float64, n, 0) + Q = zeros(Float64, n, length(cols)) + for (j, col) in enumerate(cols) + Q[:, j] = col + end + Q +end + +"""Return a floating-point copy with each column's mean removed.""" +function _center_columns(M::AbstractMatrix{<:Real})::Matrix{Float64} + n, p = size(M) + out = zeros(Float64, n, p) + for j in 1:p + mu = sum(Float64(M[i, j]) for i in 1:n; init=0.0) / n + for i in 1:n + out[i, j] = Float64(M[i, j]) - mu + end + end + out +end + +""" + pure_linear_rda(Y, X; Z=nothing) -> NamedTuple + +Pure-Julia linear Redundancy Analysis (with optional partial conditioning matrix +`Z`) using two-pass modified Gram-Schmidt projection and Ezekiel adjusted R². +Used as an independent in-Julia algebraic cross-check of total, conditioned, +constrained, and residual inertia (`tot.chi`, `pCCA\$tot.chi`, `CCA\$tot.chi`, +`CA\$tot.chi`), `R2`, `adj_r2`, and pseudo-F without invoking R. + +Rows are samples in the same order in `Y` (responses), `X` (constraints) and +`Z` (conditioning covariates); columns are centred internally. Return these +statistics, positive eigenvalues in descending order and model/condition/residual +degrees of freedom. For partial models, R² uses total inertia and adjusted R² +subtracts the adjusted conditioning-only fit from the joint fit. Pseudo-F is +`Inf` when residual inertia is zero. Throw `ArgumentError` for mismatched row +counts, fewer than three samples, rank-zero constraints after conditioning or +non-positive residual degrees of freedom. +""" +function pure_linear_rda(Y::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}; + Z::Union{Nothing,AbstractMatrix{<:Real}} = nothing) + n, p = size(Y) + size(X, 1) == n || throw(ArgumentError("X row count ($(size(X, 1))) != Y row count ($n)")) + n >= 3 || throw(ArgumentError("pure_linear_rda requires at least 3 samples")) + + Yc = _center_columns(Y) + scale_dof = Float64(n - 1) + ss_total = sum(x * x for x in Yc; init=0.0) / scale_dof + + ss_cond = 0.0 + q_rank = 0 + Y_work = copy(Yc) + X_work = _center_columns(X) + + if !isnothing(Z) + size(Z, 1) == n || throw(ArgumentError("Z row count ($(size(Z, 1))) != Y row count ($n)")) + Zc = _center_columns(Z) + Qz = _modified_gram_schmidt(Zc) + q_rank = size(Qz, 2) + if q_rank > 0 + # Project Y_work and X_work onto orthogonal complement of Qz + Y_hat_z = zeros(Float64, n, p) + for k in 1:q_rank, j in 1:p + dot_val = sum(Qz[i, k] * Y_work[i, j] for i in 1:n; init=0.0) + for i in 1:n + Y_hat_z[i, j] += Qz[i, k] * dot_val + end + end + ss_cond = sum(x * x for x in Y_hat_z; init=0.0) / scale_dof + Y_work .-= Y_hat_z + + for k in 1:q_rank, j in 1:size(X_work, 2) + dot_val = sum(Qz[i, k] * X_work[i, j] for i in 1:n; init=0.0) + for i in 1:n + X_work[i, j] -= Qz[i, k] * dot_val + end + end + end + end + + Qx = _modified_gram_schmidt(X_work) + k_rank = size(Qx, 2) + k_rank > 0 || throw(ArgumentError("Constraining design matrix has rank 0 after centring/conditioning")) + df_res = n - 1 - q_rank - k_rank + df_res > 0 || throw(ArgumentError("Model is saturated (residual df = $df_res)")) + + Y_hat = zeros(Float64, n, p) + for k in 1:k_rank, j in 1:p + dot_val = sum(Qx[i, k] * Y_work[i, j] for i in 1:n; init=0.0) + for i in 1:n + Y_hat[i, j] += Qx[i, k] * dot_val + end + end + Y_res = Y_work .- Y_hat + + ss_constrained = sum(x * x for x in Y_hat; init=0.0) / scale_dof + ss_residual = sum(x * x for x in Y_res; init=0.0) / scale_dof + + # vegan defines RsquareAdj(rda)$r.squared = CCA$tot.chi / tot.chi even in + # partial RDA, and adjusts against the unconditioned proportion via Ezekiel + r2 = ss_total > 0 ? (ss_constrained / ss_total) : 0.0 + adj_r2 = if q_rank == 0 + pure_ezekiel_adj_r2(r2, n, k_rank) + else + # vegan:::RsquareAdj.rda for partial RDA: + # RsquareAdj(tot_XZ) - RsquareAdj(cond_Z) + r2_xz = ss_total > 0 ? ((ss_constrained + ss_cond) / ss_total) : 0.0 + r2_z = ss_total > 0 ? (ss_cond / ss_total) : 0.0 + pure_ezekiel_adj_r2(r2_xz, n, k_rank + q_rank) - pure_ezekiel_adj_r2(r2_z, n, q_rank) + end + f_stat = ss_residual > 0 ? ((ss_constrained / k_rank) / (ss_residual / df_res)) : Inf + con_eigs = _gram_eigenvalues(Y_hat, scale_dof) + unc_eigs = _gram_eigenvalues(Y_res, scale_dof) + + (; total_inertia = ss_total, + conditioned_inertia = ss_cond, + constrained_inertia = ss_constrained, + unconstrained_inertia = ss_residual, + constrained_eigenvalues = con_eigs, + unconstrained_eigenvalues = unc_eigs, + r2 = r2, + adj_r2 = adj_r2, + f_statistic = f_stat, + df_model = k_rank, + df_conditioned = q_rank, + df_residual = df_res) +end + +## Block-permutation testability helper (mirrors Analysis._varies_within_any_block) + +"""Return whether any block contains differing values; paired inputs are traversed with `zip`.""" +function _varies_within_any_block(values::AbstractVector, blocks::AbstractVector)::Bool + first_seen = Dict{Any,Any}() + for (v, b) in zip(values, blocks) + isequal(get!(first_seen, b, v), v) || return true + end + false +end + +"""Return whether every constituent column of `term` is constant within every block.""" +function _term_untestable_in_blocks(term::String, metadata::DataFrame, + blocks::AbstractVector)::Bool + for part in split(term, ':') + _varies_within_any_block(metadata[!, part], blocks) && return false + end + true +end + +_num_or_nothing(x) = (ismissing(x) || isnothing(x) || !(x isa Real) || !isfinite(x)) ? nothing : Float64(x) + +_as_str_vec(x::AbstractVector) = String[String(v) for v in x] +_as_str_vec(x::AbstractString) = String[String(x)] +_as_str_vec(::Nothing) = String[] +_as_str_vec(::Missing) = String[] + +_as_f64_vec(x::AbstractVector) = Float64[Float64(v) for v in x] +_as_f64_vec(x::Real) = Float64[Float64(x)] +_as_f64_vec(::Nothing) = Float64[] +_as_f64_vec(::Missing) = Float64[] + +_as_int_vec(x::AbstractVector) = Int[Int(v) for v in x] +_as_int_vec(x::Integer) = Int[Int(x)] +_as_int_vec(::Nothing) = Int[] +_as_int_vec(::Missing) = Int[] + +_as_bool_vec(x::AbstractVector) = Bool[Bool(v) for v in x] +_as_bool_vec(x::Bool) = Bool[x] +_as_bool_vec(::Nothing) = Bool[] +_as_bool_vec(::Missing) = Bool[] + +_as_any_vec(x::AbstractVector) = collect(x) +_as_any_vec(x) = (ismissing(x) || isnothing(x)) ? Any[] : Any[x] + +""" + vegan_available() -> Bool + +Return `true` when R can be initialised and `vegan` can be loaded under `with_r_lock`. +Return `false` if R is unavailable or checking the package raises an error, +including a timeout waiting for the shared runtime. +""" +function vegan_available()::Bool + r_available() || return false + try + with_r_lock(; timeout=R_WAIT_SECONDS[]) do + Bool(RCall.rcopy(RCall.reval("requireNamespace('vegan', quietly = TRUE)"))) + end + catch + false + end +end + +""" +Return `nothing` for an `n × n` finite distance matrix with diagonal magnitudes +at most `1e-10`, upper-triangle values at least `-1e-12` and symmetry within +`1e-8`; otherwise throw `ArgumentError`. The matrix is not changed. +""" +function _validate_distance_matrix!(D::AbstractMatrix{<:Real}, n::Int) + size(D) == (n, n) || throw(ArgumentError( + "Distance matrix must have shape ($n, $n) matching sample count, got $(size(D))")) + for i in 1:n + d_ii = Float64(D[i, i]) + isfinite(d_ii) || throw(ArgumentError("Distance matrix contains non-finite diagonal at ($i, $i)")) + abs(d_ii) <= 1e-10 || throw(ArgumentError( + "Distance matrix diagonal must be 0 (got $d_ii at ($i, $i))")) + for j in (i + 1):n + d_ij = Float64(D[i, j]) + d_ji = Float64(D[j, i]) + (isfinite(d_ij) && isfinite(d_ji)) || throw(ArgumentError( + "Distance matrix contains non-finite entry at ($i, $j)")) + d_ij >= -1e-12 || throw(ArgumentError( + "Distance matrix contains negative dissimilarity $d_ij at ($i, $j)")) + abs(d_ij - d_ji) <= 1e-8 || throw(ArgumentError( + "Distance matrix is not symmetric at ($i, $j): $d_ij != $d_ji")) + end + end + nothing +end + +""" +Copy selected, already validated columns for R and return `(clean, factor_cols)`. +Numeric columns become floats; categorical columns become +trimmed strings and are listed for conversion to R factors. +""" +function _prepare_metadata_for_r(metadata::DataFrame, cols::Vector{String}) + clean = DataFrame() + factor_cols = String[] + for c in cols + raw_col = metadata[!, c] + if eltype(raw_col) <: Union{Real,Missing} && !(eltype(raw_col) <: Bool) + clean[!, c] = Float64[Float64(v) for v in raw_col] + else + clean[!, c] = String[String(strip(string(v))) for v in raw_col] + push!(factor_cols, c) + end + end + (clean, factor_cols) +end + +const _ORD_R = raw""" +ord_vegan_ver <- as.character(utils::packageDescription("vegan", fields = "Version")) +rownames(ord_mat) <- ord_sample_ids +if (ncol(ord_mat) == length(ord_feature_ids)) { + colnames(ord_mat) <- ord_feature_ids +} +for (col_nm in ord_factor_cols) { + ord_meta[[col_nm]] <- factor(ord_meta[[col_nm]], levels = unique(ord_meta[[col_nm]])) +} +rownames(ord_meta) <- ord_sample_ids + +ord_seed_reset <- function() { + set.seed(ord_seed, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection") +} + +ord_perm_ctrl <- if (length(ord_blocks) > 0L) { + permute::how(blocks = factor(ord_blocks, levels = unique(ord_blocks)), nperm = ord_nperm) +} else { + ord_nperm +} + +ord_vegan_dist <- if (ord_distance %in% c("aitchison", "euclidean")) "euclidean" else ord_distance + +ord_seed_reset() +if (ord_use_dist) { + rownames(ord_dist_mat) <- ord_sample_ids + colnames(ord_dist_mat) <- ord_sample_ids + ord_d <- stats::as.dist(ord_dist_mat) + ord_form <- stats::as.formula(paste("ord_d ~", ord_formula_rhs), env = globalenv()) + if (ord_method == "cap") { + ord_fit <- if (ncol(ord_mat) > 0L) { + vegan::capscale(ord_form, data = ord_meta, comm = ord_mat) + } else { + vegan::capscale(ord_form, data = ord_meta) + } + } else { + ord_fit <- vegan::dbrda(ord_form, data = ord_meta) + } +} else { + ord_form <- stats::as.formula(paste("ord_mat ~", ord_formula_rhs), env = globalenv()) + if (ord_method == "rda") { + ord_fit <- vegan::rda(ord_form, data = ord_meta) + } else if (ord_method == "cca") { + ord_fit <- vegan::cca(ord_form, data = ord_meta) + } else if (ord_method == "cap") { + ord_fit <- vegan::capscale(ord_form, data = ord_meta, distance = ord_vegan_dist) + } else { + ord_fit <- vegan::dbrda(ord_form, data = ord_meta, distance = ord_vegan_dist) + } +} + +ord_tot_chi <- as.numeric(ord_fit$tot.chi) +ord_cond_chi <- if (!is.null(ord_fit$pCCA) && !is.null(ord_fit$pCCA$tot.chi)) as.numeric(ord_fit$pCCA$tot.chi) else 0.0 +ord_cca_chi <- if (!is.null(ord_fit$CCA) && !is.null(ord_fit$CCA$tot.chi)) as.numeric(ord_fit$CCA$tot.chi) else 0.0 +ord_ca_chi <- if (!is.null(ord_fit$CA) && !is.null(ord_fit$CA$tot.chi)) as.numeric(ord_fit$CA$tot.chi) else 0.0 + +ord_cca_raw_eig <- if (!is.null(ord_fit$CCA) && !is.null(ord_fit$CCA$eig)) as.numeric(ord_fit$CCA$eig) else numeric(0) +ord_cca_raw_names <- if (!is.null(ord_fit$CCA) && !is.null(ord_fit$CCA$eig)) as.character(names(ord_fit$CCA$eig)) else character(0) +ord_ca_raw_eig <- if (!is.null(ord_fit$CA) && !is.null(ord_fit$CA$eig)) as.numeric(ord_fit$CA$eig) else numeric(0) +ord_ca_raw_names <- if (!is.null(ord_fit$CA) && !is.null(ord_fit$CA$eig)) as.character(names(ord_fit$CA$eig)) else character(0) + +ord_cca_pos <- ord_cca_raw_eig > 0 +ord_ca_pos <- ord_ca_raw_eig > 0 +ord_cca_eig <- ord_cca_raw_eig[ord_cca_pos] +ord_cca_names <- ord_cca_raw_names[ord_cca_pos] +ord_ca_eig <- ord_ca_raw_eig[ord_ca_pos] +ord_ca_names <- ord_ca_raw_names[ord_ca_pos] + +ord_imag_chi <- if (!is.null(ord_fit$CA) && !is.null(ord_fit$CA$imaginary.chi)) { + abs(as.numeric(ord_fit$CA$imaginary.chi)) +} else { + abs(sum(c(ord_cca_raw_eig[!ord_cca_pos], ord_ca_raw_eig[!ord_ca_pos]))) +} + +ord_seed_reset() +ord_r2_obj <- if (ord_method == "cca") { + tryCatch(vegan::RsquareAdj(ord_fit, permutations = ord_perm_ctrl), + error = function(e) vegan::RsquareAdj(ord_fit, permutations = ord_nperm)) +} else { + vegan::RsquareAdj(ord_fit) +} +ord_r2 <- if (is.null(ord_r2_obj$r.squared) || length(ord_r2_obj$r.squared) == 0L) NA_real_ else as.numeric(ord_r2_obj$r.squared[1]) +ord_adj_r2 <- if (is.null(ord_r2_obj$adj.r.squared) || length(ord_r2_obj$adj.r.squared) == 0L) NA_real_ else as.numeric(ord_r2_obj$adj.r.squared[1]) + +ord_n_cca_u <- if (!is.null(ord_fit$CCA) && !is.null(ord_fit$CCA$u)) ncol(as.matrix(ord_fit$CCA$u)) else 0L +ord_n_ca_u <- if (!is.null(ord_fit$CA) && !is.null(ord_fit$CA$u)) ncol(as.matrix(ord_fit$CA$u)) else 0L +ord_n_avail <- ord_n_cca_u + ord_n_ca_u +ord_choices <- seq_len(min(2L, max(1L, ord_n_avail))) + +ord_pad2 <- function(m) { + m <- as.matrix(m) + if (ncol(m) >= 2L) { + matrix(as.numeric(m[, 1:2, drop = FALSE]), nrow = nrow(m), ncol = 2L) + } else if (ncol(m) == 1L) { + matrix(c(as.numeric(m[, 1]), rep(0.0, nrow(m))), nrow = nrow(m), ncol = 2L) + } else { + matrix(0.0, nrow = nrow(m), ncol = 2L) + } +} + +ord_sites_raw <- as.matrix(vegan::scores(ord_fit, choices = ord_choices, display = "sites", scaling = ord_scaling)) +ord_site_scores <- ord_pad2(ord_sites_raw) +ord_axis_names <- colnames(ord_sites_raw) +if (is.null(ord_axis_names)) ord_axis_names <- paste0("Axis", seq_len(ncol(ord_sites_raw))) +if (length(ord_axis_names) < 2L) ord_axis_names <- c(ord_axis_names, "Axis2") +ord_axis_names <- as.character(ord_axis_names[1:2]) + +if (ord_n_cca_u >= 1L && !is.null(ord_fit$CCA$biplot) && nrow(as.matrix(ord_fit$CCA$biplot)) > 0L) { + ord_bp_choices <- seq_len(min(2L, ord_n_cca_u)) + ord_bp_raw <- as.matrix(vegan::scores(ord_fit, choices = ord_bp_choices, display = "bp", scaling = ord_scaling)) + ord_bp_scores <- ord_pad2(ord_bp_raw) + ord_bp_names <- as.character(rownames(ord_bp_raw)) +} else { + ord_bp_scores <- matrix(0.0, nrow = 0L, ncol = 2L) + ord_bp_names <- character(0) +} + +if (ord_n_cca_u >= 1L && !is.null(ord_fit$CCA$centroids) && nrow(as.matrix(ord_fit$CCA$centroids)) > 0L) { + ord_cn_choices <- seq_len(min(2L, ord_n_cca_u)) + ord_cn_raw <- tryCatch( + as.matrix(vegan::scores(ord_fit, choices = ord_cn_choices, display = "cn", scaling = ord_scaling)), + error = function(e) NULL + ) + if (!is.null(ord_cn_raw) && nrow(ord_cn_raw) > 0L) { + ord_cn_scores <- ord_pad2(ord_cn_raw) + ord_cn_names <- as.character(rownames(ord_cn_raw)) + } else { + ord_cn_scores <- matrix(0.0, nrow = 0L, ncol = 2L) + ord_cn_names <- character(0) + } +} else { + ord_cn_scores <- matrix(0.0, nrow = 0L, ncol = 2L) + ord_cn_names <- character(0) +} + +ord_sp_raw <- tryCatch( + vegan::scores(ord_fit, choices = ord_choices, display = "species", scaling = ord_scaling), + error = function(e) NULL +) +if (!is.null(ord_sp_raw) && length(ord_sp_raw) > 0L && !all(is.na(ord_sp_raw))) { + ord_sp_mat <- as.matrix(ord_sp_raw) + ord_sp_keep <- stats::complete.cases(ord_sp_mat) + ord_sp_mat <- ord_sp_mat[ord_sp_keep, , drop = FALSE] + if (nrow(ord_sp_mat) > 0L) { + ord_sp_scores <- ord_pad2(ord_sp_mat) + ord_sp_names <- as.character(rownames(ord_sp_mat)) + if (length(ord_sp_names) == 0L) ord_sp_names <- paste0("feature_", seq_len(nrow(ord_sp_mat))) + } else { + ord_sp_scores <- matrix(0.0, nrow = 0L, ncol = 2L) + ord_sp_names <- character(0) + } +} else { + ord_sp_scores <- matrix(0.0, nrow = 0L, ncol = 2L) + ord_sp_names <- character(0) +} + +ord_extract_aov <- function(aov_obj) { + if (is.null(aov_obj)) { + return(list(terms = character(0), df = integer(0), var = numeric(0), + f = numeric(0), p = numeric(0), nperm = ord_nperm, text = "")) + } + rn <- as.character(rownames(aov_obj)) + df_col <- as.integer(aov_obj[["Df"]]) + var_col <- as.numeric(aov_obj[[2]]) + f_col <- as.numeric(aov_obj[["F"]]) + p_col <- as.numeric(aov_obj[["Pr(>F)"]]) + ctrl <- attr(aov_obj, "control") + np <- if (!is.null(ctrl) && !is.null(ctrl[["nperm"]])) as.integer(ctrl[["nperm"]]) else as.integer(ord_nperm) + txt <- paste(utils::capture.output(print(aov_obj)), collapse = "\n") + list(terms = rn, df = df_col, var = var_col, f = f_col, p = p_col, nperm = np, text = txt) +} + +ord_seed_reset() +ord_aov_overall <- ord_extract_aov(vegan::anova.cca(ord_fit, permutations = ord_perm_ctrl, parallel = 1L)) + +ord_seed_reset() +ord_aov_terms <- ord_extract_aov(vegan::anova.cca(ord_fit, permutations = ord_perm_ctrl, by = "terms", parallel = 1L)) + +ord_aov_margin <- if (ord_test_margin) { + ord_seed_reset() + ord_extract_aov(tryCatch(vegan::anova.cca(ord_fit, permutations = ord_perm_ctrl, by = "margin", parallel = 1L), + error = function(e) NULL)) +} else { + ord_extract_aov(NULL) +} + +ord_aov_axes <- if (ord_test_axes) { + ord_seed_reset() + ord_extract_aov(tryCatch(vegan::anova.cca(ord_fit, permutations = ord_perm_ctrl, by = "axis", parallel = 1L), + error = function(e) NULL)) +} else { + ord_extract_aov(NULL) +} + +ord_vp_ok <- FALSE +ord_vp_ind_labels <- character(0) +ord_vp_ind_df <- integer(0) +ord_vp_ind_r2 <- numeric(0) +ord_vp_ind_adj <- numeric(0) +ord_vp_ind_testable <- logical(0) +ord_vp_gross_labels <- character(0) +ord_vp_gross_df <- integer(0) +ord_vp_gross_r2 <- numeric(0) +ord_vp_gross_adj <- numeric(0) +ord_vp_err <- "" + +if (ord_do_varpart) { + tryCatch({ + vp_formulas <- lapply(ord_vp_terms, function(tm) stats::as.formula(paste("~", tm), env = globalenv())) + vp_y <- if (ord_use_dist) { + ord_d + } else if (ord_method %in% c("cap", "dbrda")) { + vegan::vegdist(ord_mat, method = ord_vegan_dist) + } else { + ord_mat + } + ord_seed_reset() + vp_args <- c(list(vp_y), vp_formulas, list(data = ord_meta)) + if (ord_method == "cca") { + vp_args$chisquare <- TRUE + vp_args$permutations <- ord_perm_ctrl + } + vp_obj <- do.call(vegan::varpart, vp_args) + vp_ind <- vp_obj$part$indfract + vp_gr <- vp_obj$part$fract + r2_col <- function(df) { + if ("R.squared" %in% colnames(df)) as.numeric(df[["R.squared"]]) + else if ("R.square" %in% colnames(df)) as.numeric(df[["R.square"]]) + else rep(NA_real_, nrow(df)) + } + adj_col <- function(df) { + if ("Adj.R.squared" %in% colnames(df)) as.numeric(df[["Adj.R.squared"]]) + else if ("Adj.R.square" %in% colnames(df)) as.numeric(df[["Adj.R.square"]]) + else rep(NA_real_, nrow(df)) + } + ord_vp_ind_labels <- as.character(rownames(vp_ind)) + ord_vp_ind_df <- ifelse(is.na(vp_ind[["Df"]]), 0L, as.integer(vp_ind[["Df"]])) + ord_vp_ind_r2 <- r2_col(vp_ind) + ord_vp_ind_adj <- adj_col(vp_ind) + ord_vp_ind_testable <- ifelse(is.na(vp_ind[["Testable"]]), FALSE, as.logical(vp_ind[["Testable"]])) + ord_vp_gross_labels <- as.character(rownames(vp_gr)) + ord_vp_gross_df <- ifelse(is.na(vp_gr[["Df"]]), 0L, as.integer(vp_gr[["Df"]])) + ord_vp_gross_r2 <- r2_col(vp_gr) + ord_vp_gross_adj <- adj_col(vp_gr) + ord_vp_ok <- TRUE + }, error = function(e) { + ord_vp_err <<- conditionMessage(e) + }) +} +""" + +const _ORD_CLEANUP_R = + "rm(list = intersect(ls(), c(" * + "'ord_method', 'ord_distance', 'ord_use_dist', 'ord_mat', 'ord_dist_mat', " * + "'ord_sample_ids', 'ord_feature_ids', 'ord_meta', 'ord_factor_cols', " * + "'ord_formula_rhs', 'ord_scaling', 'ord_nperm', 'ord_seed', 'ord_blocks', " * + "'ord_test_margin', 'ord_test_axes', 'ord_do_varpart', 'ord_vp_terms', " * + "'ord_vegan_ver', 'col_nm', 'ord_seed_reset', 'ord_perm_ctrl', 'ord_vegan_dist', " * + "'ord_d', 'ord_form', 'ord_fit', 'ord_tot_chi', 'ord_cond_chi', 'ord_cca_chi', " * + "'ord_ca_chi', 'ord_cca_raw_eig', 'ord_cca_raw_names', 'ord_ca_raw_eig', 'ord_ca_raw_names', " * + "'ord_cca_pos', 'ord_ca_pos', 'ord_cca_eig', 'ord_cca_names', 'ord_ca_eig', 'ord_ca_names', " * + "'ord_imag_chi', 'ord_r2_obj', 'ord_r2', 'ord_adj_r2', 'ord_n_cca_u', " * + "'ord_n_ca_u', 'ord_n_avail', 'ord_choices', 'ord_pad2', 'ord_sites_raw', " * + "'ord_site_scores', 'ord_axis_names', 'ord_bp_choices', 'ord_bp_raw', " * + "'ord_bp_scores', 'ord_bp_names', 'ord_cn_choices', 'ord_cn_raw', " * + "'ord_cn_scores', 'ord_cn_names', 'ord_sp_raw', 'ord_sp_mat', 'ord_sp_keep', 'ord_sp_scores', " * + "'ord_sp_names', 'ord_extract_aov', 'ord_aov_overall', 'ord_aov_terms', " * + "'ord_aov_margin', 'ord_aov_axes', 'ord_vp_ok', 'ord_vp_ind_labels', " * + "'ord_vp_ind_df', 'ord_vp_ind_r2', 'ord_vp_ind_adj', 'ord_vp_ind_testable', " * + "'ord_vp_gross_labels', 'ord_vp_gross_df', 'ord_vp_gross_r2', " * + "'ord_vp_gross_adj', 'ord_vp_err', 'vp_formulas', 'vp_y', 'vp_args', " * + "'vp_obj', 'vp_ind', 'vp_gr', 'r2_col', 'adj_col'))); invisible(gc())" + +""" +Convert an extracted vegan ANOVA table into term rows plus residual statistics, +permutation count and printed text. Omit total/residual rows from the term list, +replace non-finite statistics with `nothing`, and withhold p-values for terms +in `untestable`, retaining their F statistics and an explanatory reason. +""" +function _parse_aov_list(d, untestable::Set{String}) + terms = _as_str_vec(d[:terms]) + dfs = _as_any_vec(d[:df]) + vars = _as_any_vec(d[:var]) + fs = _as_any_vec(d[:f]) + ps = _as_any_vec(d[:p]) + nperm = first(_as_int_vec(d[:nperm])) + txt = first(_as_str_vec(d[:text])) + + rows = NamedTuple[] + res_df = 0 + res_var = 0.0 + for k in eachindex(terms) + t = strip(terms[k]) + df_val = ismissing(dfs[k]) ? 0 : Int(dfs[k]) + v_val = _num_or_nothing(vars[k]) + if t in ("Residual", "Residuals") + res_df = df_val + res_var = isnothing(v_val) ? 0.0 : v_val + continue + end + t == "Total" && continue + is_untest = t in untestable + f_val = _num_or_nothing(fs[k]) + p_raw = _num_or_nothing(ps[k]) + p_val = is_untest ? nothing : p_raw + reason = is_untest ? + "'$t' is constant within every permutation block, so within-block permutations never exchange its labels" : + nothing + push!(rows, (; + term = t, + df = df_val, + variance = v_val, + f_statistic = f_val, + p_value = p_val, + untestable_reason = reason, + )) + end + (; rows, residual_df = res_df, residual_variance = res_var, nperm, text = txt) +end + +""" + run_constrained_ordination(counts, metadata, cfg; kwargs...) -> NamedTuple + +Fit a constrained ordination (`:rda`, `:cca`, `:cap`, or `:dbrda`) via R's +`vegan` package under the shared R runtime lock, run seeded permutation tests +(`anova.cca` omnibus, sequential `by = "terms"`, marginal `by = "margin"`, and +optional `by = "axis"`), and optionally compute `vegan::varpart` variance +partitioning. + +`counts` is `n_samples × n_features` (rows = samples, columns = ASVs/OTUs/taxa). +For `:unweighted_unifrac` or `:weighted_unifrac`, pass either a precomputed +`distance_matrix` (`n_samples × n_samples`) or the `tree` and `feature_ids` keywords so +`CladeCumulus.unifrac_matrix` computes the phylogenetic distances. + +Metadata and sample IDs follow count-row order; feature IDs follow column order. +At least three samples and positive residual model degrees of freedom are +required. Aitchison input must be non-negative: with zeros, retained features +must occur in at least two samples for Bayesian replacement or one for CZM, +and retained counts must be integer-valued. `beta_transform` accepts `"none"` +or `"hellinger"` where applied; Aitchison and tree-derived UniFrac ignore it. + +`blocks` restricts permutations only when its length matches the sample count +and at least one label repeats; otherwise permutations are unrestricted. +Terms constant within every block have their p-values withheld. `r_timeout` +limits lock acquisition in seconds, not model execution. Fitting resets and +advances the shared R random-number generator using `cfg.seed`. + +Return inertia and R² statistics, eigenvalues, two-axis sample/feature/covariate +scores, permutation tests, cautions and provenance. `max_feature_scores` limits +the feature scores selected by vector length, while all scores are also returned. +Variance partitioning is `nothing` when disabled, otherwise a computed or +skipped result with a reason. Failures of optional marginal/axis tests or +centroid/species scoring yield empty collections; caught `vegan::varpart` +failures yield a skipped result. + +Throw `ArgumentError` for invalid data or formulae, `VeganUnavailable` if vegan +cannot load, and `RBusyError` if the runtime lock times out. Zero replacement +can propagate `ZCompositionsUnavailable`; uncaught R fitting, omnibus/sequential +testing and conversion errors also propagate. Unknown beta transforms raise +`ErrorException`. +""" +function run_constrained_ordination( + counts::AbstractMatrix{<:Real}, + metadata::DataFrame, + cfg::OrdinationConfig; + sample_ids::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + feature_ids::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + beta_transform::AbstractString = "none", + distance_matrix::Union{Nothing,AbstractMatrix{<:Real}} = nothing, + tree::Union{Nothing,PhyloNode,AbstractString} = nothing, + blocks::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + test_axes::Bool = false, + max_feature_scores::Integer = 25, + r_timeout::Real = R_WAIT_SECONDS[], +) + cfg.ordination_method === :none && throw(ArgumentError( + "Constrained ordination method is ':none'; select 'rda', 'cca', 'cap', or 'dbrda'")) + + n_samples, n_features = size(counts) + nrow(metadata) == n_samples || throw(ArgumentError( + "metadata has $(nrow(metadata)) rows for $n_samples samples in counts")) + + s_ids = isnothing(sample_ids) ? ["sample_$i" for i in 1:n_samples] : String.(sample_ids) + length(s_ids) == n_samples || throw(ArgumentError( + "sample_ids has length $(length(s_ids)) for $n_samples samples")) + length(unique(s_ids)) == n_samples || throw(ArgumentError( + "sample_ids must be unique across all $n_samples samples")) + + f_ids = isnothing(feature_ids) ? ["feature_$j" for j in 1:n_features] : String.(feature_ids) + length(f_ids) == n_features || throw(ArgumentError( + "feature_ids has length $(length(f_ids)) for $n_features features")) + + parsed = parse_ordination_formula(cfg.ordination_formula, metadata) + r_meta, factor_cols = _prepare_metadata_for_r(metadata, parsed.referenced_columns) + + blocked = !isnothing(blocks) && length(blocks) == n_samples && !allunique(blocks) + untestable_set = if blocked + Set([t for t in parsed.constrained_terms if _term_untestable_in_blocks(t, metadata, blocks)]) + else + Set{String}() + end + + cautions = String[] + zero_replacement_info = nothing + dropped_zero_features = 0 + + use_dist = false + dist_mat_r = zeros(Float64, 0, 0) + ord_mat = Matrix{Float64}(counts) + kept_feature_ids = copy(f_ids) + + for x in ord_mat + isfinite(x) || throw(ArgumentError("Input community matrix contains non-finite values")) + end + + if !isnothing(distance_matrix) + cfg.ordination_method in (:cap, :dbrda) || throw(ArgumentError( + "Precomputed distance_matrix is only supported for 'cap' and 'dbrda', not '$(cfg.ordination_method)'")) + _validate_distance_matrix!(distance_matrix, n_samples) + use_dist = true + dist_mat_r = Matrix{Float64}(distance_matrix) + if lowercase(strip(String(beta_transform))) != "none" && n_features > 0 && all(>=(0.0), ord_mat) + ord_mat = transform_counts(ord_mat; method=String(beta_transform)) + end + elseif cfg.ordination_distance in (:unweighted_unifrac, :weighted_unifrac) + isnothing(tree) && throw(ArgumentError( + "Phylogenetic distance '$(cfg.ordination_distance)' requires a phylogenetic `tree` or precomputed `distance_matrix`")) + n_features >= 2 || throw(ArgumentError( + "UniFrac ordination requires at least 2 features (got $n_features)")) + all(>=(0.0), ord_mat) || throw(ArgumentError( + "UniFrac ordination requires non-negative abundances")) + for i in 1:n_samples + sum(@view ord_mat[i, :]) > 0.0 || throw(ArgumentError( + "Sample '$(s_ids[i])' has zero total abundance; empty samples have undefined UniFrac distances")) + end + metric_sym = cfg.ordination_distance === :unweighted_unifrac ? :unweighted : :weighted + phylo_node = tree isa AbstractString ? parse_newick(tree) : tree + dist_mat_r = unifrac_matrix(permutedims(ord_mat), kept_feature_ids, s_ids, phylo_node; metric=metric_sym) + _validate_distance_matrix!(dist_mat_r, n_samples) + use_dist = true + elseif cfg.ordination_distance === :aitchison + n_features >= 2 || throw(ArgumentError( + "Aitchison (CLR) ordination requires at least 2 features (got $n_features)")) + all(>=(0.0), ord_mat) || throw(ArgumentError( + "Aitchison (CLR) ordination expects non-negative counts or proportions; pass distance='euclidean' if data are already CLR-transformed")) + for i in 1:n_samples + sum(@view ord_mat[i, :]) > 0.0 || throw(ArgumentError( + "Sample '$(s_ids[i])' has zero total reads; CLR is undefined on empty samples")) + end + has_zeros = any(==(0.0), ord_mat) + if has_zeros + min_pres = cfg.replacement_method == "bayesian_multiplicative" ? 2 : 1 + keep_cols = [j for j in 1:n_features if count(>(0.0), @view ord_mat[:, j]) >= min_pres] + dropped_zero_features = n_features - length(keep_cols) + length(keep_cols) >= 2 || throw(ArgumentError( + "Aitchison ($(cfg.replacement_method)) ordination requires at least 2 features present in >= $min_pres sample(s)")) + ord_mat = ord_mat[:, keep_cols] + kept_feature_ids = kept_feature_ids[keep_cols] + for i in 1:n_samples + sum(@view ord_mat[i, :]) > 0.0 || throw(ArgumentError( + "Sample '$(s_ids[i])' has zero reads across retained features after sparsity filtering")) + end + int_counts = round.(Int, ord_mat) + all(isapprox.(ord_mat, int_counts; atol=1e-9)) || throw(ArgumentError( + "Aitchison zero replacement expects raw integer counts before normalisation (got non-integer values with zeros)")) + repl = replace_zeros( + int_counts; + method = cfg.replacement_method, + frac = cfg.replacement_delta, + samples = s_ids, + taxa = kept_feature_ids, + ) + ord_mat = clr_transform(repl.proportions) + zero_replacement_info = Dict{String,Any}( + "method" => cfg.replacement_method, + "delta" => cfg.replacement_delta, + "zeros_replaced" => Int(repl.zeros_replaced), + ) + push!(cautions, + "Zeros were imputed via $(cfg.replacement_method) (delta=$(cfg.replacement_delta), $(repl.zeros_replaced) cell(s)) before CLR; imputed zeros are biased adjustments, not observed abundances.") + else + # Close each positive row to a composition and compute exact CLR + n_kept = size(ord_mat, 2) + props = similar(ord_mat) + for i in 1:n_samples + rs = sum(@view ord_mat[i, :]) + for j in 1:n_kept + props[i, j] = ord_mat[i, j] / rs + end + end + ord_mat = clr_transform(props) + zero_replacement_info = Dict{String,Any}( + "method" => cfg.replacement_method, + "delta" => cfg.replacement_delta, + "zeros_replaced" => 0, + ) + end + elseif cfg.ordination_method === :cca + n_features >= 2 || throw(ArgumentError( + "CCA requires at least 2 features (got $n_features)")) + all(>=(0.0), ord_mat) || throw(ArgumentError( + "CCA chi-square ordination requires non-negative abundances; negative entries (such as CLR values) are incompatible with CCA")) + btrans = lowercase(strip(String(beta_transform))) + if btrans != "none" + ord_mat = transform_counts(ord_mat; method=btrans) + end + keep_cols = [j for j in 1:n_features if sum(@view ord_mat[:, j]) > 0.0] + dropped_zero_features = n_features - length(keep_cols) + length(keep_cols) >= 2 || throw(ArgumentError( + "CCA requires at least 2 features with positive total abundance across samples")) + ord_mat = ord_mat[:, keep_cols] + kept_feature_ids = kept_feature_ids[keep_cols] + for i in 1:n_samples + sum(@view ord_mat[i, :]) > 0.0 || throw(ArgumentError( + "Sample '$(s_ids[i])' has zero total abundance; CCA chi-square weights are undefined on empty samples")) + end + else + # rda(euclidean) or cap/dbrda(bray, jaccard, euclidean) + n_features >= 2 || throw(ArgumentError( + "Constrained ordination requires at least 2 features (got $n_features)")) + btrans = lowercase(strip(String(beta_transform))) + if btrans != "none" + all(>=(0.0), ord_mat) || throw(ArgumentError( + "Beta transform '$btrans' requires non-negative counts")) + ord_mat = transform_counts(ord_mat; method=btrans) + elseif cfg.ordination_method === :rda && all(>=(0.0), ord_mat) && any(==(0.0), ord_mat) + push!(cautions, + "Linear Euclidean RDA on untransformed sparse counts is sensitive to library-size differences and shared zeros (Legendre & Gallagher 2001; Gloor et al. 2017); consider distance='aitchison' (CLR) or Hellinger beta normalisation.") + end + if cfg.ordination_distance in (:bray, :jaccard) + all(>=(0.0), ord_mat) || throw(ArgumentError( + "Dissimilarity '$(cfg.ordination_distance)' requires non-negative abundances")) + for i in 1:n_samples + sum(@view ord_mat[i, :]) > 0.0 || throw(ArgumentError( + "Sample '$(s_ids[i])' has zero total abundance; '$(cfg.ordination_distance)' is undefined on empty samples")) + end + end + end + + can_varpart_r = cfg.variance_partitioning && + (2 <= length(parsed.constrained_terms) <= 4) && + isempty(parsed.condition_terms) + + r_out = with_r_lock(; timeout=Float64(r_timeout)) do + RCall.rcopy(RCall.reval("requireNamespace('vegan', quietly = TRUE)")) || throw(VeganUnavailable()) + RCall.reval("suppressPackageStartupMessages(library(vegan))") + + RCall.globalEnv[:ord_method] = string(cfg.ordination_method) + RCall.globalEnv[:ord_distance] = string(cfg.ordination_distance) + RCall.globalEnv[:ord_use_dist] = use_dist + RCall.globalEnv[:ord_mat] = ord_mat + RCall.globalEnv[:ord_dist_mat] = dist_mat_r + RCall.globalEnv[:ord_sample_ids] = s_ids + RCall.globalEnv[:ord_feature_ids] = kept_feature_ids + RCall.globalEnv[:ord_meta] = r_meta + RCall.globalEnv[:ord_factor_cols] = factor_cols + RCall.globalEnv[:ord_formula_rhs] = parsed.canonical_rhs + RCall.globalEnv[:ord_scaling] = cfg.ordination_scaling + RCall.globalEnv[:ord_nperm] = cfg.permutations + RCall.globalEnv[:ord_seed] = cfg.seed + RCall.globalEnv[:ord_blocks] = blocked ? String.(blocks) : String[] + RCall.globalEnv[:ord_test_margin] = length(parsed.constrained_terms) >= 2 + RCall.globalEnv[:ord_test_axes] = test_axes + RCall.globalEnv[:ord_do_varpart] = can_varpart_r + RCall.globalEnv[:ord_vp_terms] = parsed.constrained_terms + try + RCall.reval(_ORD_R) + (; + vegan_version = String(RCall.rcopy(RCall.reval("ord_vegan_ver"))), + tot_chi = Float64(RCall.rcopy(RCall.reval("ord_tot_chi"))), + cond_chi = Float64(RCall.rcopy(RCall.reval("ord_cond_chi"))), + cca_chi = Float64(RCall.rcopy(RCall.reval("ord_cca_chi"))), + ca_chi = Float64(RCall.rcopy(RCall.reval("ord_ca_chi"))), + imag_chi = Float64(RCall.rcopy(RCall.reval("ord_imag_chi"))), + cca_eig = _as_f64_vec(RCall.rcopy(RCall.reval("ord_cca_eig"))), + cca_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_cca_names"))), + ca_eig = _as_f64_vec(RCall.rcopy(RCall.reval("ord_ca_eig"))), + ca_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_ca_names"))), + r2 = _num_or_nothing(RCall.rcopy(RCall.reval("ord_r2"))), + adj_r2 = _num_or_nothing(RCall.rcopy(RCall.reval("ord_adj_r2"))), + site_scores = Matrix{Float64}(RCall.rcopy(RCall.reval("ord_site_scores"))), + axis_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_axis_names"))), + bp_scores = Matrix{Float64}(RCall.rcopy(RCall.reval("ord_bp_scores"))), + bp_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_bp_names"))), + cn_scores = Matrix{Float64}(RCall.rcopy(RCall.reval("ord_cn_scores"))), + cn_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_cn_names"))), + sp_scores = Matrix{Float64}(RCall.rcopy(RCall.reval("ord_sp_scores"))), + sp_names = _as_str_vec(RCall.rcopy(RCall.reval("ord_sp_names"))), + aov_overall = RCall.rcopy(RCall.reval("ord_aov_overall")), + aov_terms = RCall.rcopy(RCall.reval("ord_aov_terms")), + aov_margin = RCall.rcopy(RCall.reval("ord_aov_margin")), + aov_axes = RCall.rcopy(RCall.reval("ord_aov_axes")), + vp_ok = Bool(RCall.rcopy(RCall.reval("ord_vp_ok"))), + vp_ind_labels = _as_str_vec(RCall.rcopy(RCall.reval("ord_vp_ind_labels"))), + vp_ind_df = _as_int_vec(RCall.rcopy(RCall.reval("ord_vp_ind_df"))), + vp_ind_r2 = _as_any_vec(RCall.rcopy(RCall.reval("ord_vp_ind_r2"))), + vp_ind_adj = _as_any_vec(RCall.rcopy(RCall.reval("ord_vp_ind_adj"))), + vp_ind_test = _as_bool_vec(RCall.rcopy(RCall.reval("ord_vp_ind_testable"))), + vp_gr_labels = _as_str_vec(RCall.rcopy(RCall.reval("ord_vp_gross_labels"))), + vp_gr_df = _as_int_vec(RCall.rcopy(RCall.reval("ord_vp_gross_df"))), + vp_gr_r2 = _as_any_vec(RCall.rcopy(RCall.reval("ord_vp_gross_r2"))), + vp_gr_adj = _as_any_vec(RCall.rcopy(RCall.reval("ord_vp_gross_adj"))), + vp_err = String(RCall.rcopy(RCall.reval("ord_vp_err"))), + ) + finally + RCall.reval(_ORD_CLEANUP_R) + end + end + + if abs(r_out.imag_chi) > 0.0 && r_out.tot_chi > 0.0 && (abs(r_out.imag_chi) / r_out.tot_chi) > 0.05 + pct = round(100.0 * abs(r_out.imag_chi) / r_out.tot_chi; digits=1) + push!(cautions, + "Non-Euclidean dissimilarity produced negative eigenvalues ($(pct)% of total inertia); biplot axes display the real Euclidean subspace.") + end + + push!(cautions, + "Constrained ordination axes maximise variance aligned with the supplied covariates; interpret biplots alongside permutation F/p-values and residual unconstrained inertia rather than reading clustering alone as confirmatory proof.") + + tot = r_out.tot_chi + prop_of_tot(x::Float64) = tot > 0.0 ? (x / tot) : 0.0 + + constrained_eigs = [ + (; axis = r_out.cca_names[i], + eigenvalue = r_out.cca_eig[i], + proportion_total = prop_of_tot(r_out.cca_eig[i]), + proportion_constrained = r_out.cca_chi > 0.0 ? (r_out.cca_eig[i] / r_out.cca_chi) : 0.0) + for i in eachindex(r_out.cca_eig) + ] + unconstrained_eigs = [ + (; axis = r_out.ca_names[i], + eigenvalue = r_out.ca_eig[i], + proportion_total = prop_of_tot(r_out.ca_eig[i])) + for i in eachindex(r_out.ca_eig) + ] + + eig_lookup = Dict{String,Float64}() + for e in constrained_eigs + eig_lookup[e.axis] = e.proportion_total + end + for e in unconstrained_eigs + eig_lookup[e.axis] = e.proportion_total + end + + axis_labels = [ + let ax = r_out.axis_names[k] + if haskey(eig_lookup, ax) + pct = round(100.0 * max(0.0, eig_lookup[ax]); digits=1) + kind = any(e -> e.axis == ax, constrained_eigs) ? "constrained" : "unconstrained" + "$ax ($pct%, $kind)" + else + ax + end + end + for k in 1:2 + ] + + overall_parsed = _parse_aov_list(r_out.aov_overall, Set{String}()) + terms_parsed = _parse_aov_list(r_out.aov_terms, untestable_set) + margin_parsed = _parse_aov_list(r_out.aov_margin, untestable_set) + axes_parsed = _parse_aov_list(r_out.aov_axes, Set{String}()) + + overall_untestable = blocked && !isempty(parsed.constrained_terms) && + all(t -> t in untestable_set, parsed.constrained_terms) + overall_row = isempty(overall_parsed.rows) ? nothing : overall_parsed.rows[1] + actual_nperm = overall_parsed.nperm + min_p = 1.0 / (actual_nperm + 1) + + omnibus = (; + df_model = isnothing(overall_row) ? parsed.df_constrained : overall_row.df, + df_residual = overall_parsed.residual_df, + variance_model = isnothing(overall_row) ? r_out.cca_chi : overall_row.variance, + variance_residual = overall_parsed.residual_variance, + f_statistic = isnothing(overall_row) ? nothing : overall_row.f_statistic, + p_value = (isnothing(overall_row) || overall_untestable) ? nothing : overall_row.p_value, + untestable_reason = overall_untestable ? + "All constraining covariates are constant within every permutation block, so within-block permutations cannot test the model" : + nothing, + permutations = actual_nperm, + min_p_value = min_p, + text = overall_parsed.text, + ) + + site_rows = [ + (; sample = s_ids[i], + axis1 = r_out.site_scores[i, 1], + axis2 = r_out.site_scores[i, 2]) + for i in 1:n_samples + ] + + biplot_arrows = [ + let x = r_out.bp_scores[i, 1], y = r_out.bp_scores[i, 2] + (; label = r_out.bp_names[i], + axis1 = x, + axis2 = y, + length = sqrt(x * x + y * y)) + end + for i in eachindex(r_out.bp_names) + ] + + centroids = [ + (; label = r_out.cn_names[i], + axis1 = r_out.cn_scores[i, 1], + axis2 = r_out.cn_scores[i, 2]) + for i in eachindex(r_out.cn_names) + ] + + all_feature_rows = [ + let x = r_out.sp_scores[i, 1], y = r_out.sp_scores[i, 2] + (; feature = r_out.sp_names[i], + axis1 = x, + axis2 = y, + length = sqrt(x * x + y * y)) + end + for i in eachindex(r_out.sp_names) + ] + sort!(all_feature_rows; by = r -> (-r.length, r.feature)) + top_features = all_feature_rows[1:min(Int(max_feature_scores), length(all_feature_rows))] + + varpart_result = if !cfg.variance_partitioning + nothing + elseif can_varpart_r && r_out.vp_ok + ind_rows = [ + (; label = strip(r_out.vp_ind_labels[i]), + df = r_out.vp_ind_df[i], + r2 = _num_or_nothing(r_out.vp_ind_r2[i]), + adj_r2 = _num_or_nothing(r_out.vp_ind_adj[i]), + testable = r_out.vp_ind_test[i]) + for i in eachindex(r_out.vp_ind_labels) + ] + gr_rows = [ + (; label = strip(r_out.vp_gr_labels[i]), + df = r_out.vp_gr_df[i], + r2 = _num_or_nothing(r_out.vp_gr_r2[i]), + adj_r2 = _num_or_nothing(r_out.vp_gr_adj[i])) + for i in eachindex(r_out.vp_gr_labels) + ] + has_neg = any(r -> !isnothing(r.adj_r2) && r.adj_r2 < -1e-6, ind_rows) + vp_note = has_neg ? + "Adjusted-R2 subtraction produced small negative fraction(s); per Peres-Neto et al. (2006), negative fractions indicate predictors explaining less variance than random normal variates (or suppressor effects) and are reported unclamped so fractions sum to 1." : + "Adjusted-R2 variance partitioning (vegan::varpart; Peres-Neto et al. 2006). Shared overlap fractions have 0 df and cannot be tested individually by permutation." + (; status = "computed", + covariates = copy(parsed.constrained_terms), + conditioned_on = String[], + individual_fractions = ind_rows, + gross_fractions = gr_rows, + has_negative_fraction = has_neg, + note = vp_note) + elseif !isempty(parsed.condition_terms) && length(parsed.constrained_terms) == 2 + # Partial 2-way decomposition conditioned on Condition(...) + t1, t2 = parsed.constrained_terms[1], parsed.constrained_terms[2] + cond_str = "Condition(" * join(parsed.condition_terms, " + ") * ")" + cfg_sub1 = OrdinationConfig(; + ordination_method = cfg.ordination_method, + ordination_formula = "~ $t1 + $cond_str", + ordination_distance = cfg.ordination_distance, + ordination_scaling = cfg.ordination_scaling, + permutations = cfg.permutations, + variance_partitioning = false, + replacement_method = cfg.replacement_method, + replacement_delta = cfg.replacement_delta, + seed = cfg.seed, + ) + cfg_sub2 = OrdinationConfig(; + ordination_method = cfg.ordination_method, + ordination_formula = "~ $t2 + $cond_str", + ordination_distance = cfg.ordination_distance, + ordination_scaling = cfg.ordination_scaling, + permutations = cfg.permutations, + variance_partitioning = false, + replacement_method = cfg.replacement_method, + replacement_delta = cfg.replacement_delta, + seed = cfg.seed, + ) + sub1 = run_constrained_ordination(counts, metadata, cfg_sub1; + sample_ids=s_ids, feature_ids=f_ids, beta_transform, distance_matrix, tree, blocks, r_timeout) + sub2 = run_constrained_ordination(counts, metadata, cfg_sub2; + sample_ids=s_ids, feature_ids=f_ids, beta_transform, distance_matrix, tree, blocks, r_timeout) + if !isnothing(sub1.adj_r2) && !isnothing(sub2.adj_r2) && !isnothing(r_out.adj_r2) + vp2 = pure_varpart_two_way(sub1.adj_r2, sub2.adj_r2, r_out.adj_r2) + ind_rows = [ + (; label = "[a] = $t1 | $t2", df = _term_df(t1, metadata), r2 = nothing, adj_r2 = vp2.unique_x1, testable = true), + (; label = "[b] = Shared ($t1 & $t2)", df = 0, r2 = nothing, adj_r2 = vp2.shared, testable = false), + (; label = "[c] = $t2 | $t1", df = _term_df(t2, metadata), r2 = nothing, adj_r2 = vp2.unique_x2, testable = true), + (; label = "[d] = Residuals incl. conditioned share", df = omnibus.df_residual, r2 = nothing, adj_r2 = vp2.residual, testable = false), + ] + has_neg = any(r -> !isnothing(r.adj_r2) && r.adj_r2 < -1e-6, ind_rows) + (; status = "computed", + covariates = copy(parsed.constrained_terms), + conditioned_on = copy(parsed.condition_terms), + individual_fractions = ind_rows, + gross_fractions = [ + (; label = "[a+b] = $t1", df = _term_df(t1, metadata), r2 = sub1.r2, adj_r2 = sub1.adj_r2), + (; label = "[b+c] = $t2", df = _term_df(t2, metadata), r2 = sub2.r2, adj_r2 = sub2.adj_r2), + (; label = "[a+b+c] = $t1 + $t2", df = omnibus.df_model, r2 = r_out.r2, adj_r2 = r_out.adj_r2), + ], + has_negative_fraction = has_neg, + note = "Partial adjusted-R2 variance partitioning conditioned on $(join(parsed.condition_terms, ", ")).") + else + (; status = "skipped", + reason = "Adjusted R2 could not be computed for partial sub-models", + covariates = copy(parsed.constrained_terms), + conditioned_on = copy(parsed.condition_terms), + individual_fractions = NamedTuple[], + gross_fractions = NamedTuple[], + has_negative_fraction = false, + note = nothing) + end + else + why = if length(parsed.constrained_terms) < 2 + "Variance partitioning requires 2 to 4 constraining terms in the formula (got $(length(parsed.constrained_terms)): $(join(parsed.constrained_terms, ", ")))" + elseif length(parsed.constrained_terms) > 4 + "vegan::varpart supports up to 4 constraining tables (got $(length(parsed.constrained_terms)))" + elseif !isempty(r_out.vp_err) + "vegan::varpart could not partition this model: $(r_out.vp_err)" + else + "Partial variance partitioning with Condition(...) currently supports 2 constraining terms (got $(length(parsed.constrained_terms)))" + end + (; status = "skipped", + reason = why, + covariates = copy(parsed.constrained_terms), + conditioned_on = copy(parsed.condition_terms), + individual_fractions = NamedTuple[], + gross_fractions = NamedTuple[], + has_negative_fraction = false, + note = nothing) + end + + provenance = Dict{String,Any}( + "method" => string(cfg.ordination_method), + "distance" => string(cfg.ordination_distance), + "scaling" => cfg.ordination_scaling, + "scaling_meaning" => cfg.ordination_scaling == 1 ? "distance" : "correlation", + "formula" => "~ " * parsed.canonical_rhs, + "constrained_terms" => copy(parsed.constrained_terms), + "condition_terms" => copy(parsed.condition_terms), + "permutations_requested" => cfg.permutations, + "permutations_performed" => actual_nperm, + "seed" => cfg.seed, + "rng_kind" => "R Mersenne-Twister / Inversion / Rejection", + "blocked" => blocked, + "beta_normalisation" => lowercase(strip(String(beta_transform))), + "n_samples" => n_samples, + "n_features" => length(kept_feature_ids), + "dropped_zero_features" => dropped_zero_features, + "zero_replacement" => zero_replacement_info, + "multiple_testing_policy" => "omnibus_and_term_permutation_unadjusted; per_taxon_tests_require_bh", + "software" => Dict{String,Any}( + "r_package" => "vegan", + "vegan_version" => r_out.vegan_version, + ), + ) + + (; + method = string(cfg.ordination_method), + distance = string(cfg.ordination_distance), + scaling = cfg.ordination_scaling, + formula = "~ " * parsed.canonical_rhs, + constrained_terms = copy(parsed.constrained_terms), + condition_terms = copy(parsed.condition_terms), + n_samples = n_samples, + n_features = length(kept_feature_ids), + total_inertia = r_out.tot_chi, + conditioned_inertia = r_out.cond_chi, + constrained_inertia = r_out.cca_chi, + unconstrained_inertia = r_out.ca_chi, + imaginary_inertia = r_out.imag_chi, + conditioned_proportion = prop_of_tot(r_out.cond_chi), + constrained_proportion = prop_of_tot(r_out.cca_chi), + unconstrained_proportion = prop_of_tot(r_out.ca_chi), + r2 = r_out.r2, + adj_r2 = r_out.adj_r2, + axis_names = copy(r_out.axis_names), + axis_labels = axis_labels, + constrained_eigenvalues = constrained_eigs, + unconstrained_eigenvalues = unconstrained_eigs, + site_scores = r_out.site_scores, + sites = site_rows, + biplot_arrows = biplot_arrows, + centroids = centroids, + top_feature_scores = top_features, + all_feature_scores = all_feature_rows, + omnibus_test = omnibus, + term_tests = terms_parsed.rows, + marginal_tests = margin_parsed.rows, + axis_tests = axes_parsed.rows, + variance_partitioning = varpart_result, + cautions = cautions, + provenance = provenance, + ) +end + +## Plotly Biplot Builder + +const _BIPLOT_PALETTE = [ + "#2563eb", "#dc2626", "#16a34a", "#d97706", "#7c3aed", + "#0891b2", "#db2777", "#65a30d", "#ea580c", "#4f46e5", +] + +const _BIPLOT_SYMBOLS = [ + "circle", "square", "diamond", "triangle-up", "triangle-down", + "pentagon", "hexagon", "star", "cross", "x", +] + +"""Scale the longest arrow to `fill_frac` of the largest site radius; return 1 for empty arrows or zero radii.""" +function _arrow_multiplier(site_scores::Matrix{Float64}, arrows::AbstractVector; fill_frac::Float64 = 0.78)::Float64 + isempty(arrows) && return 1.0 + max_site = 0.0 + for i in axes(site_scores, 1) + r = hypot(site_scores[i, 1], site_scores[i, 2]) + r > max_site && (max_site = r) + end + max_arrow = 0.0 + for a in arrows + r = hypot(Float64(a.axis1), Float64(a.axis2)) + r > max_arrow && (max_arrow = r) + end + (max_site <= 0.0 || max_arrow <= 0.0) && return 1.0 + fill_frac * (max_site / max_arrow) +end + +""" + ordination_biplot_chart(ord_result; colour_by=nothing, shape_by=nothing, show_centroids=true, show_features=false, max_features=8) -> Dict{String,Any} + +Build an interactive, I-JSON-compliant Plotly biplot dictionary showing sample +points (coloured by `colour_by`, shaped by `shape_by`), covariate arrows scaled +to the sample radii, optional factor centroids, and optional top taxon vectors. +Use a result from `run_constrained_ordination`; grouping labels follow site +order and mismatched lengths raise `ArgumentError`. Feature overlays use up to +`max_features` entries from `top_feature_scores`. Return `data` and `layout` +without changing the result or rendering a chart. +""" +function ordination_biplot_chart( + ord_result; + colour_by::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + shape_by::Union{Nothing,AbstractVector{<:AbstractString}} = nothing, + show_centroids::Bool = true, + show_features::Bool = false, + max_features::Integer = 8, +)::Dict{String,Any} + coords = ord_result.site_scores + labels = [s.sample for s in ord_result.sites] + n = size(coords, 1) + + cb = isnothing(colour_by) ? fill("Samples", n) : String.(colour_by) + sb = isnothing(shape_by) ? fill("all", n) : String.(shape_by) + length(cb) == n || throw(ArgumentError("colour_by length ($(length(cb))) != sample count ($n)")) + length(sb) == n || throw(ArgumentError("shape_by length ($(length(sb))) != sample count ($n)")) + + cgroups = unique(cb) + sgroups = unique(sb) + + traces = Dict{String,Any}[] + for (ci, cg) in enumerate(cgroups) + col_hex = _BIPLOT_PALETTE[mod1(ci, length(_BIPLOT_PALETTE))] + for (si, sg) in enumerate(sgroups) + idxs = [i for i in 1:n if cb[i] == cg && sb[i] == sg] + isempty(idxs) && continue + same_label = cg == sg + trace_name = length(sgroups) <= 1 ? cg : (same_label ? cg : "$cg / $sg") + sym = _BIPLOT_SYMBOLS[mod1(si, length(_BIPLOT_SYMBOLS))] + push!(traces, Dict{String,Any}( + "type" => "scatter", + "mode" => "markers", + "name" => trace_name, + "legendgroup" => trace_name, + "x" => Float64[coords[i, 1] for i in idxs], + "y" => Float64[coords[i, 2] for i in idxs], + "text" => String[labels[i] for i in idxs], + "hovertemplate" => "Sample: %{text}
Group: " * trace_name * + "
%{xaxis.title.text}: %{x:.4f}
%{yaxis.title.text}: %{y:.4f}", + "marker" => Dict{String,Any}( + "color" => col_hex, + "size" => 10, + "symbol" => sym, + "line" => Dict{String,Any}("color" => "#ffffff", "width" => 1), + ), + )) + end + end + + if show_centroids && !isempty(ord_result.centroids) + push!(traces, Dict{String,Any}( + "type" => "scatter", + "mode" => "markers+text", + "name" => "Factor centroids", + "legendgroup" => "__centroids__", + "x" => Float64[c.axis1 for c in ord_result.centroids], + "y" => Float64[c.axis2 for c in ord_result.centroids], + "text" => String[c.label for c in ord_result.centroids], + "textposition" => "top center", + "hovertemplate" => "Centroid: %{text}
(%{x:.4f}, %{y:.4f})", + "marker" => Dict{String,Any}( + "color" => "#111827", + "size" => 11, + "symbol" => "diamond-open", + "line" => Dict{String,Any}("color" => "#111827", "width" => 2), + ), + )) + end + + annotations = Dict{String,Any}[] + mult = _arrow_multiplier(coords, ord_result.biplot_arrows) + for arr in ord_result.biplot_arrows + ax = arr.axis1 * mult + ay = arr.axis2 * mult + # Arrow annotation from origin (0, 0) to scaled tip (ax, ay) + push!(annotations, Dict{String,Any}( + "x" => ax, + "y" => ay, + "ax" => 0.0, + "ay" => 0.0, + "xref" => "x", + "yref" => "y", + "axref" => "x", + "ayref" => "y", + "text" => "", + "showarrow" => true, + "arrowhead" => 2, + "arrowsize" => 1.1, + "arrowwidth" => 2.0, + "arrowcolor" => "#991b1b", + )) + # Label slightly beyond tip + push!(annotations, Dict{String,Any}( + "x" => ax * 1.08, + "y" => ay * 1.08, + "xref" => "x", + "yref" => "y", + "text" => "$(arr.label)", + "showarrow" => false, + "font" => Dict{String,Any}("size" => 11, "color" => "#991b1b"), + )) + end + + if !isempty(ord_result.biplot_arrows) + # Invisible hover markers on arrow tips so exact unscaled biplot cosines remain inspectable + push!(traces, Dict{String,Any}( + "type" => "scatter", + "mode" => "markers", + "name" => "Covariate vectors", + "legendgroup" => "__covariates__", + "x" => Float64[a.axis1 * mult for a in ord_result.biplot_arrows], + "y" => Float64[a.axis2 * mult for a in ord_result.biplot_arrows], + "text" => String[a.label for a in ord_result.biplot_arrows], + "customdata" => Vector{Float64}[[a.axis1, a.axis2] for a in ord_result.biplot_arrows], + "hovertemplate" => "Covariate: %{text}
Biplot score: (%{customdata[0]:.4f}, %{customdata[1]:.4f})", + "marker" => Dict{String,Any}("color" => "#991b1b", "size" => 6, "symbol" => "circle-open"), + )) + end + + if show_features && !isempty(ord_result.top_feature_scores) + feats = ord_result.top_feature_scores[1:min(Int(max_features), length(ord_result.top_feature_scores))] + f_mult = _arrow_multiplier(coords, feats; fill_frac=0.65) + push!(traces, Dict{String,Any}( + "type" => "scatter", + "mode" => "markers+text", + "name" => "Top feature loadings", + "legendgroup" => "__features__", + "x" => Float64[f.axis1 * f_mult for f in feats], + "y" => Float64[f.axis2 * f_mult for f in feats], + "text" => String[f.feature for f in feats], + "textposition" => "bottom center", + "customdata" => Vector{Float64}[[f.axis1, f.axis2] for f in feats], + "hovertemplate" => "Feature: %{text}
Score: (%{customdata[0]:.4f}, %{customdata[1]:.4f})", + "marker" => Dict{String,Any}("color" => "#475569", "size" => 6, "symbol" => "triangle-up-open"), + "textfont" => Dict{String,Any}("size" => 10, "color" => "#475569"), + )) + end + + r2_txt = isnothing(ord_result.r2) ? "n/a" : string(round(ord_result.r2; digits=3)) + adj_txt = isnothing(ord_result.adj_r2) ? "n/a" : string(round(ord_result.adj_r2; digits=3)) + p_txt = isnothing(ord_result.omnibus_test.p_value) ? "withheld" : + string(round(ord_result.omnibus_test.p_value; digits=4)) + summary_badge = "$(uppercase(ord_result.method)) ($(ord_result.distance), scaling $(ord_result.scaling)) · R²=$r2_txt, adj. R²=$adj_txt, p=$p_txt" + + push!(annotations, Dict{String,Any}( + "text" => summary_badge, + "xref" => "paper", + "yref" => "paper", + "x" => 0.01, + "y" => 0.99, + "xanchor" => "left", + "yanchor" => "top", + "showarrow" => false, + "font" => Dict{String,Any}("size" => 11, "color" => "#334155"), + )) + + method_upper = ord_result.method == "dbrda" ? "dbRDA" : uppercase(ord_result.method) + layout = Dict{String,Any}( + "title" => Dict{String,Any}("text" => "$method_upper constrained ordination ($(ord_result.formula))"), + "xaxis" => Dict{String,Any}("title" => ord_result.axis_labels[1], "zeroline" => true, "zerolinecolor" => "#cbd5e1"), + "yaxis" => Dict{String,Any}("title" => ord_result.axis_labels[2], "zeroline" => true, "zerolinecolor" => "#cbd5e1"), + "annotations" => annotations, + ) + + Dict{String,Any}("data" => traces, "layout" => layout) +end + +end # module Ordination diff --git a/src/core/validate.jl b/src/core/validate.jl index b98d55d3..dc9b0b76 100644 --- a/src/core/validate.jl +++ b/src/core/validate.jl @@ -540,6 +540,47 @@ export validate_environment, validate_project, ValidationError, isnothing(lt) || (_is_number(lt) && 0 < lt <= 1) || _err(errors, ctx, "clade_cumulus.lwr_threshold must be in (0, 1] (got: $(repr(lt)))") end + + an = get(cfg, "analysis", nothing) + if an isa Dict + ord = get(an, "ordination", nothing) + if ord isa Dict + om = get(ord, "method", "none") + valid_methods = ("none", "rda", "cca", "cap", "dbrda") + if !(om isa AbstractString && lowercase(strip(om)) in valid_methods) + _err(errors, ctx, "analysis.ordination.method must be one of $(join(valid_methods, ", ")) (got: $(repr(om)))") + end + od = get(ord, "distance", "bray") + valid_dists = ("bray", "jaccard", "euclidean", "aitchison", "chisq", "unweighted_unifrac", "weighted_unifrac") + if !(od isa AbstractString && lowercase(strip(od)) in valid_dists) + _err(errors, ctx, "analysis.ordination.distance must be one of $(join(valid_dists, ", ")) (got: $(repr(od)))") + elseif om isa AbstractString && lowercase(strip(om)) in valid_methods + m_norm = lowercase(strip(om)) + d_norm = lowercase(strip(od)) + compat = m_norm == "rda" ? ("aitchison", "euclidean") : + m_norm == "cca" ? ("chisq",) : + m_norm in ("cap", "dbrda") ? ("bray", "jaccard", "euclidean", "aitchison", "unweighted_unifrac", "weighted_unifrac") : + valid_dists + d_norm in compat || _err(errors, ctx, + "analysis.ordination.distance '$d_norm' is incompatible with analysis.ordination.method '$m_norm'") + end + osc = get(ord, "scaling", 2) + (_whole(osc) && osc in (1, 2)) || + _err(errors, ctx, "analysis.ordination.scaling must be integer 1 or 2 (got: $(repr(osc)))") + operm = get(ord, "permutations", 999) + (_whole(operm) && 99 <= operm <= 9999) || + _err(errors, ctx, "analysis.ordination.permutations must be an integer between 99 and 9999 (got: $(repr(operm)))") + ovp = get(ord, "variance_partitioning", false) + ovp isa Bool || + _err(errors, ctx, "analysis.ordination.variance_partitioning must be Bool (got: $(repr(ovp)))") + orm = get(ord, "replacement_method", "bayesian_multiplicative") + (orm isa AbstractString && lowercase(strip(orm)) in ("bayesian_multiplicative", "multiplicative_replacement")) || + _err(errors, ctx, "analysis.ordination.replacement_method must be bayesian_multiplicative or multiplicative_replacement (got: $(repr(orm)))") + ord_delta = get(ord, "replacement_delta", 0.65) + (_is_number(ord_delta) && 0 < ord_delta < 1) || + _err(errors, ctx, "analysis.ordination.replacement_delta must be in (0, 1) (got: $(repr(ord_delta)))") + end + end end ## Per-project Validation diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 3d8073d1..a745d615 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -111,6 +111,9 @@ const _ENUM_VALUES = Dict( "analysis.differential.dispersion" => Differential.DISPERSIONS, "analysis.differential.dm_overdispersion_method" => Differential.DM_METHODS, "analysis.differential.replacement_method" => Differential.REPLACEMENT_METHODS, + "analysis.ordination.method" => string.(Ordination.ORDINATION_METHODS), + "analysis.ordination.distance" => string.(Ordination.ORDINATION_DISTANCES), + "analysis.ordination.replacement_method" => Ordination.REPLACEMENT_METHODS, "phylogeny.reference.align.strategy" => Validation.PHYLOGENY_ALIGN_STRATEGIES, "phylogeny.placement.align.strategy" => Validation.PHYLOGENY_ALIGN_STRATEGIES, "phylogeny.reference.trim.method" => Validation.PHYLOGENY_TRIM_METHODS, @@ -127,7 +130,8 @@ const _NULLABLE_NUMBERS = ("phylogeny.placement.place.heuristic", "analysis.differential.bayesian_alpha") # Text keys whose factory default is ~: an empty value clears the setting. -const _NULLABLE_STRINGS = ("analysis.differential.mn_reference_taxon",) +const _NULLABLE_STRINGS = ("analysis.differential.mn_reference_taxon", + "analysis.ordination.formula") # Enum keys whose submitted value is matched without regard to case and then # stored canonically, so " ML " and " glmGAMPOI " reach the pipeline as written. @@ -135,6 +139,9 @@ const _CASE_FOLDED_ENUMS = Dict( "analysis.differential.dispersion" => Differential.DISPERSIONS, "analysis.differential.dm_overdispersion_method" => Differential.DM_METHODS, "analysis.differential.replacement_method" => Differential.REPLACEMENT_METHODS, + "analysis.ordination.method" => string.(Ordination.ORDINATION_METHODS), + "analysis.ordination.distance" => string.(Ordination.ORDINATION_DISTANCES), + "analysis.ordination.replacement_method" => Ordination.REPLACEMENT_METHODS, ) """ @@ -150,7 +157,10 @@ Differential optimiser settings reject booleans and non-finite values: and runtime in seconds must be positive, and the ridge penalty non-negative. `InexactError` and `OverflowError` from the iteration-limit conversion become validation messages. This checks individual values without saving them; -method compatibility is checked separately by `_differential_settings_error`. +method compatibility is checked separately by `_differential_settings_error` +and `_ordination_settings_error`. Non-blank ordination formulae undergo syntax +validation; formula `ArgumentError`s become validation messages. Ordination +scaling must be 1 or 2 and permutation counts must lie in `99:9999`. """ function _config_value_error(key::String, value) v = value isa JSON3.Array ? collect(value) : value @@ -162,12 +172,30 @@ function _config_value_error(key::String, value) end if key in _NULLABLE_STRINGS (isnothing(v) || v isa AbstractString) || return "$key must be text or empty" + if key == "analysis.ordination.formula" && v isa AbstractString && !isempty(strip(v)) + try + Ordination.validate_formula_syntax(v) + catch e + e isa ArgumentError && return "$key is invalid: $(e.msg)" + rethrow() + end + end return nothing end - if key == "analysis.differential.replacement_delta" + if key in ("analysis.differential.replacement_delta", "analysis.ordination.replacement_delta") return (v isa Real && !(v isa Bool) && isfinite(v) && 0 < v < 1) ? nothing : "$key must be a finite number in (0, 1)" end + if key == "analysis.ordination.scaling" + return (v isa Real && !(v isa Bool) && isfinite(v) && isinteger(v) && + 1 <= v <= 2 && Int(v) in Ordination.ORDINATION_SCALINGS) ? nothing : + "$key must be 1 (distance) or 2 (correlation)" + end + if key == "analysis.ordination.permutations" + return (v isa Real && !(v isa Bool) && isfinite(v) && isinteger(v) && + Ordination.MIN_ORDINATION_PERMUTATIONS <= v <= Ordination.MAX_ORDINATION_PERMUTATIONS) ? nothing : + "$key must be a whole number in [$(Ordination.MIN_ORDINATION_PERMUTATIONS), $(Ordination.MAX_ORDINATION_PERMUTATIONS)]" + end if key == "analysis.differential.bayesian_alpha" isnothing(v) && return nothing return (v isa Real && !(v isa Bool) && isfinite(v) && v > 0) ? nothing : @@ -260,6 +288,107 @@ const _DIFFERENTIAL_SETTING_DEFAULTS = Dict{String,Any}( const _DIFFERENTIAL_SETTING_KEYS = collect(keys(_DIFFERENTIAL_SETTING_DEFAULTS)) +const _ORDINATION_SETTING_DEFAULTS = Dict{String,Any}( + "analysis.ordination.method" => "none", + "analysis.ordination.distance" => "bray", +) + +const _ORDINATION_SETTING_KEYS = collect(keys(_ORDINATION_SETTING_DEFAULTS)) + +""" + _ordination_settings_error(effective) -> Union{String,Nothing} + +Return a method/distance validation message or `nothing`. Dotted keys may map +to raw values or resolved `(; value, source)` entries; absent keys use factory +defaults. Non-text values are left to individual value checks. Compatibility +`ArgumentError`s become messages; other errors propagate. +""" +function _ordination_settings_error(effective::AbstractDict) + value(k) = begin + current = get(effective, k, _ORDINATION_SETTING_DEFAULTS[k]) + current isa NamedTuple ? current.value : current + end + method = value("analysis.ordination.method") + distance = value("analysis.ordination.distance") + (method isa AbstractString && distance isa AbstractString) || return nothing + try + Ordination.validate_method_distance!(method, distance) + nothing + catch e + e isa ArgumentError && return e.msg + rethrow() + end +end + +""" + _ordination_cascade_error(path, document) -> Union{String,Nothing} + +Check ordination method/distance compatibility in scopes affected by replacing +the config at `path` with `document`, without writing it. Combine factory and +global defaults with each affected study/group/run override. Return the first +invalid scope and reason, or `nothing`. Errors while loading/checking the +cascade become validation messages. +""" +function _ordination_cascade_error(path::String, document) + data = ServerState.data_dir() + root = dirname(data) + global_path = joinpath(root, "config", "pipeline.yml") + prefix = path == global_path ? String[] : splitpath(relpath(dirname(path), data)) + cache = Dict{String,Any}(path => _flatten(document)) + flat(file) = get!(cache, file) do + _flatten(_load_yml(file)) + end + function check(values, label) + effective = Dict{String,Any}(k => get(values, k, v) + for (k, v) in _ORDINATION_SETTING_DEFAULTS) + bad = _ordination_settings_error(effective) + isnothing(bad) ? nothing : "$label: $bad" + end + try + base = merge(flat(joinpath(root, "config", "defaults", "pipeline.yml")), flat(global_path)) + if isempty(prefix) + bad = check(base, "global config") + isnothing(bad) || return bad + end + targets = Set{String}() + isempty(prefix) || push!(targets, joinpath(prefix...)) + if length(prefix) >= 2 + if length(prefix) == 2 && !get(document, "pool_children", false) + for run in _group_run_names(prefix[1], prefix[2]) + push!(targets, joinpath(prefix..., run)) + end + end + else + for study in _study_names() + isempty(prefix) || first(prefix) == study || continue + push!(targets, study) + for run in _run_names(study) + push!(targets, joinpath(study, run)) + end + for group in _group_names(study) + push!(targets, joinpath(study, group)) + for run in _group_run_names(study, group) + push!(targets, joinpath(study, group, run)) + end + end + end + end + for label in sort!(collect(targets)) + parts = splitpath(label) + length(parts) >= length(prefix) && parts[1:length(prefix)] == prefix || continue + effective = copy(base) + for i in eachindex(parts) + merge!(effective, flat(joinpath(data, parts[1:i]..., "pipeline.yml"))) + end + bad = check(effective, label) + isnothing(bad) || return bad + end + nothing + catch e + "Cannot validate the affected config cascade: $(sprint(showerror, e))" + end +end + """ _differential_settings_error(effective) -> Union{String,Nothing} @@ -385,6 +514,15 @@ end # Validate, prepare the complete document, and publish once under one lock. # Readers do not take this lock, so per-key renames would expose partial PATCHes. +""" + _write_config_updates(path, changes; remove=()) -> Union{HTTP.Response,Nothing} + +Apply dotted-key updates and removals under the config lock, normalising known +case-insensitive enums. Validate affected differential and ordination cascades +before publishing the whole YAML document atomically. Return an HTTP 400 response +without writing on a cascade failure; return `nothing` after a write or no-op. +File access, update and write errors propagate. +""" function _write_config_updates(path::String, changes::AbstractDict; remove=()) lock(_config_file_lock) do native = Dict{String,Any}(string(k) => _to_plain(v) for (k, v) in changes) @@ -408,6 +546,10 @@ function _write_config_updates(path::String, changes::AbstractDict; remove=()) bad = _differential_cascade_error(path, cfg) isnothing(bad) || return json_error(400, "invalid_config_value", bad) end + if any(k -> haskey(native, k) || k in remove, _ORDINATION_SETTING_KEYS) + bad = _ordination_cascade_error(path, cfg) + isnothing(bad) || return json_error(400, "invalid_config_value", bad) + end _atomic_write_yaml(path, cfg) nothing end diff --git a/src/server/routes/ordination.jl b/src/server/routes/ordination.jl new file mode 100644 index 00000000..63bcff1a --- /dev/null +++ b/src/server/routes/ordination.jl @@ -0,0 +1,534 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Routes: +# GET /api/v1/ordination/contract +# POST /api/v1/studies/{study}/analysis/ordination +# GET /api/v1/studies/{study}/analysis/ordination/{id} + +using JSON3, DataFrames, OrderedCollections +using ..CladeCumulus +using ..DerivedAnalyses +using ..RRuntime: RBusyError + +const _ORDINATION_REQUEST_LIMIT = 16 * 1024 * 1024 + +_ordination_derived_dir(study::String) = + joinpath(ServerState.projects_dir(), study, "derived_analyses", "ordination") + +""" +Return sample labels, prefixing all labels with their runs if any repeat and +adding row indices if run/sample pairs still repeat. `run_labels` must align +with `samples`; neither input is changed. +""" +function _unique_ordination_sample_ids(samples::Vector{String}, run_labels::Vector{String})::Vector{String} + length(unique(samples)) == length(samples) && return copy(samples) + combined = ["$(run_labels[i]) / $(samples[i])" for i in eachindex(samples)] + length(unique(combined)) == length(combined) && return combined + ["$(combined[i]) #$i" for i in eachindex(combined)] +end + +""" +Read an I-JSON request object without consuming its body. Return an HTTP error +response for bodies over 16 MiB (413), parse failures or non-object values (400). +""" +function _read_ordination_ijson(req::HTTP.Request) + body_bytes = req.body isa AbstractString ? ncodeunits(req.body) : length(req.body) + body_bytes <= _ORDINATION_REQUEST_LIMIT || return json_error( + 413, "request_too_large", + "Ordination request exceeds the $(_ORDINATION_REQUEST_LIMIT)-byte limit") + raw = req.body isa AbstractString ? String(req.body) : String(copy(req.body)) + parsed = try + ijson_read(raw) + catch e + return json_error(400, "invalid_json", + "Request body is not valid RFC 7493 I-JSON: $(sprint(showerror, e))") + end + parsed isa JSON3.Object || return json_error(400, "invalid_json", + "Request body must be a JSON object") + parsed +end + +"""Return nested ordination settings and the seed from the study's resolved config, filling absent keys with defaults.""" +function _study_ordination_dict(study::String)::Dict{String,Any} + resolved = _resolve_config(study) + val(k, d) = begin + entry = get(resolved, k, (; value=d, source="default")) + entry isa NamedTuple ? entry.value : entry + end + Dict{String,Any}( + "seed" => val("seed", 42), + "analysis" => Dict{String,Any}( + "ordination" => Dict{String,Any}( + "method" => val("analysis.ordination.method", "none"), + "formula" => val("analysis.ordination.formula", nothing), + "distance" => val("analysis.ordination.distance", "bray"), + "scaling" => val("analysis.ordination.scaling", 2), + "permutations" => val("analysis.ordination.permutations", 999), + "variance_partitioning" => val("analysis.ordination.variance_partitioning", false), + "replacement_method" => val("analysis.ordination.replacement_method", "bayesian_multiplicative"), + "replacement_delta" => val("analysis.ordination.replacement_delta", 0.65), + ), + ), + ) +end + +""" +Copy recognised request settings to a plain dictionary, preserving explicit +nulls and ignoring other keys. Throw `ArgumentError` for non-null scaling, +permutation or seed values that are not integers (booleans excluded), or a +non-boolean variance-partitioning flag. Range checks occur in `OrdinationConfig`. +""" +function _extract_ordination_overrides(body::JSON3.Object)::Dict{String,Any} + overrides = Dict{String,Any}() + for key in ("method", "formula", "distance", "scaling", "permutations", + "variance_partitioning", "replacement_method", "replacement_delta", "seed") + sym = Symbol(key) + if haskey(body, sym) + v = body[sym] + # Validate integer-typed fields against float tokens with fractional parts or Bool + if key in ("scaling", "permutations", "seed") && !isnothing(v) + (v isa Integer && !(v isa Bool)) || throw(ArgumentError( + "'$key' must be an integer, not $(typeof(v))")) + elseif key == "variance_partitioning" && !isnothing(v) + v isa Bool || throw(ArgumentError("'variance_partitioning' must be true or false")) + end + overrides[key] = _to_plain(v) + end + end + overrides +end + +""" +Add custom columns to `metadata` and return their sorted names, or an empty list +for `nothing`. Arrays follow `sample_ids` order; maps must supply every sample. +All-numeric values become finite floats, while other values become non-empty +strings. Throw `ArgumentError` for invalid names, built-in column collisions, +shape mismatches or invalid values. Columns added before a later error remain. +""" +function _attach_custom_covariates!(metadata::DataFrame, raw_covs, sample_ids::Vector{String}) + isnothing(raw_covs) && return String[] + (raw_covs isa JSON3.Object || raw_covs isa AbstractDict) || throw(ArgumentError( + "'covariates' must be a JSON object mapping covariate names to per-sample arrays or sample->value maps")) + n = length(sample_ids) + added = String[] + reserved_cols = Set(["sample", "group", "run", "individual", "library_size", "log_library_size"]) + for (k_raw, v_raw) in pairs(raw_covs) + col_name = String(k_raw) + col_name in reserved_cols && throw(ArgumentError( + "Custom covariate '$col_name' collides with built-in sample metadata column")) + Ordination.validate_formula_syntax("~ $col_name") + vals = if v_raw isa AbstractVector + length(v_raw) == n || throw(ArgumentError( + "Custom covariate '$col_name' has length $(length(v_raw)), expected $n")) + collect(v_raw) + elseif v_raw isa JSON3.Object || v_raw isa AbstractDict + [haskey(v_raw, s) ? v_raw[s] : + haskey(v_raw, Symbol(s)) ? v_raw[Symbol(s)] : + throw(ArgumentError("Custom covariate '$col_name' is missing value for sample '$s'")) + for s in sample_ids] + else + throw(ArgumentError("Custom covariate '$col_name' must be an array or sample->value object")) + end + if all(x -> x isa Real && !(x isa Bool), vals) + num_vec = Float64[] + for x in vals + fx = Float64(x) + isfinite(fx) || throw(ArgumentError("Custom covariate '$col_name' contains non-finite value")) + push!(num_vec, fx) + end + metadata[!, col_name] = num_vec + else + str_vec = String[] + for x in vals + (isnothing(x) || ismissing(x)) && throw(ArgumentError( + "Custom covariate '$col_name' contains null/missing value")) + sx = String(strip(string(x))) + isempty(sx) && throw(ArgumentError( + "Custom covariate '$col_name' contains empty string value")) + push!(str_vec, sx) + end + metadata[!, col_name] = str_vec + end + push!(added, col_name) + end + sort!(added) + added +end + +""" +Return `(tree, tips, source)` from inline `tree_newick`, an explicit +`cladistics_id`, or the latest completed study tree, in that precedence order. +Blank IDs are ignored; failures resolving a non-blank ID or supplied Newick do +not fall back. Throw `ArgumentError` for empty Newick, invalid IDs or missing +trees; tree parsing and file access errors propagate. +""" +function _resolve_study_unifrac_tree(study::String, body::JSON3.Object) + if haskey(body, :tree_newick) && !isnothing(body[:tree_newick]) + nwk = String(strip(string(body[:tree_newick]))) + isempty(nwk) && throw(ArgumentError("'tree_newick' cannot be empty")) + idx_tree = CladeCumulus.parse_phylogenetic_tree(nwk) + return (CladeCumulus.parse_newick(nwk), idx_tree.tips, "inline_newick") + end + req_id = haskey(body, :cladistics_id) && !isnothing(body[:cladistics_id]) ? + strip(string(body[:cladistics_id])) : nothing + if !isnothing(req_id) && !isempty(req_id) + _valid_phylo_id(req_id) || throw(ArgumentError( + "cladistics_id must be an 8-character lowercase hex identifier")) + tree_file = _clad_files(study, req_id)["reference.treefile"] + isfile(tree_file) || throw(ArgumentError( + "Cladistics analysis '$req_id' has no completed ASV tree yet")) + nwk = read(tree_file, String) + idx_tree = CladeCumulus.parse_phylogenetic_tree(nwk) + return (CladeCumulus.parse_newick(nwk), idx_tree.tips, "cladistics:$req_id") + end + root = _clad_root(study) + if isdir(root) + candidates = String[] + for id in readdir(root) + _valid_phylo_id(id) || continue + isfile(_clad_files(study, id)["reference.treefile"]) && push!(candidates, id) + end + if !isempty(candidates) + sort!(candidates; by = id -> _clad_summary(study, id).created, rev=true) + best_id = first(candidates) + tree_file = _clad_files(study, best_id)["reference.treefile"] + nwk = read(tree_file, String) + idx_tree = CladeCumulus.parse_phylogenetic_tree(nwk) + return (CladeCumulus.parse_newick(nwk), idx_tree.tips, "cladistics:$best_id") + end + end + throw(ArgumentError( + "UniFrac constrained ordination requires an ASV phylogenetic tree; build one in the study's Cladistics tab first or pass 'cladistics_id' / 'tree_newick'")) +end + +""" +Select public result fields for the HTTP payload and attach available covariates +and the runtime estimate. The biplot and persistence receipt are added by the caller; +internal score matrices and the full feature-score list are omitted. +""" +function _serializable_ordination_payload(res, available_covariates::Vector{String}, + runtime_estimate) + Dict{String,Any}( + "method" => res.method, + "distance" => res.distance, + "scaling" => res.scaling, + "formula" => res.formula, + "constrained_terms" => res.constrained_terms, + "condition_terms" => res.condition_terms, + "available_covariates" => available_covariates, + "n_samples" => res.n_samples, + "n_features" => res.n_features, + "total_inertia" => res.total_inertia, + "conditioned_inertia" => res.conditioned_inertia, + "constrained_inertia" => res.constrained_inertia, + "unconstrained_inertia" => res.unconstrained_inertia, + "imaginary_inertia" => res.imaginary_inertia, + "conditioned_proportion" => res.conditioned_proportion, + "constrained_proportion" => res.constrained_proportion, + "unconstrained_proportion" => res.unconstrained_proportion, + "r2" => res.r2, + "adj_r2" => res.adj_r2, + "axis_names" => res.axis_names, + "axis_labels" => res.axis_labels, + "constrained_eigenvalues" => res.constrained_eigenvalues, + "unconstrained_eigenvalues" => res.unconstrained_eigenvalues, + "sites" => res.sites, + "biplot_arrows" => res.biplot_arrows, + "centroids" => res.centroids, + "top_feature_scores" => res.top_feature_scores, + "omnibus_test" => res.omnibus_test, + "term_tests" => res.term_tests, + "marginal_tests" => res.marginal_tests, + "axis_tests" => res.axis_tests, + "variance_partitioning" => res.variance_partitioning, + "cautions" => res.cautions, + "runtime_estimate" => runtime_estimate, + "provenance" => res.provenance, + ) +end + +"""Read a boolean request option, using `default` for absent/null values and throwing `ArgumentError` for other types.""" +function _opt_ordination_bool(body::JSON3.Object, key::Symbol, default::Bool)::Bool + (!haskey(body, key) || isnothing(body[key])) && return default + v = body[key] + v isa Bool || throw(ArgumentError("'$key' must be true or false")) + v +end + +@get "/api/v1/ordination/contract" function(req) + jcs_json(Ordination.load_ordination_contract()) +end + +@get "/api/v1/studies/{study}/analysis/ordination/{id}" function(req, study::String, id::String) + _study_exists(study) || return json_error(404, "study_not_found", "Study '$study' not found") + try + rec = DerivedAnalyses.load_derived_analysis(_ordination_derived_dir(study), id) + return jcs_json(rec) + catch err + err isa DerivedAnalyses.DerivedAnalysisError && return json_error( + err.code == "not_found" ? 404 : err.code == "conflict" ? 409 : 400, + "derived_analysis_$(err.code)", err.reason) + rethrow() + end +end + +@post "/api/v1/studies/{study}/analysis/ordination" function(req, study::String) + _study_exists(study) || return json_error(404, "study_not_found", "Study '$study' not found") + body = _read_ordination_ijson(req) + body isa HTTP.Response && return body + + cfg = try + overrides = _extract_ordination_overrides(body) + Ordination.ordination_config_from_dict(_study_ordination_dict(study); overrides) + catch err + err isa ArgumentError && return json_error(400, "invalid_ordination_config", err.msg) + rethrow() + end + + if cfg.ordination_method === :none + return json_error(400, "invalid_ordination_config", + "Constrained ordination method is 'none'; choose 'rda', 'cca', 'cap', or 'dbrda'") + end + + runs_spec = get(body, :runs, []) + (runs_spec isa AbstractVector && !isempty(runs_spec)) || + return json_error(400, "no_runs", "Provide at least one run in 'runs'") + + table = String(strip(string(get(body, :table, "merged")))) + isempty(table) && return json_error(400, "invalid_table", "Table name cannot be empty") + (aggregate, test_axes, show_centroids, show_features, persist_record) = try + ( + _opt_ordination_bool(body, :aggregate, false), + _opt_ordination_bool(body, :test_axes, false), + _opt_ordination_bool(body, :show_centroids, true), + _opt_ordination_bool(body, :show_features, false), + _opt_ordination_bool(body, :persist, false), + ) + catch err + err isa ArgumentError && return json_error(400, "invalid_request", err.msg) + rethrow() + end + params = try + _body_filter_params(body) + catch err + err isa BadRequest && return json_error(400, "bad_request", err.message) + rethrow() + end + + run_data, group_labels, run_labels, individual_ids, missing_seq = try + _collect_cross_run_asvs(study, runs_spec, table, params; aggregate) + catch err + err isa BadRequest && return json_error(400, "bad_request", err.message) + rethrow() + end + missing_seq && return json_error(400, "no_sequence_data", + "Table '$table' does not contain a DNA sequence column required for cross-run ordination") + isempty(run_data) && return json_error(400, "no_data", "No count data found for the selected runs") + + raw_mat, raw_samples, sequences, _ = combined_asv_counts_across_runs(run_data) + (filtered, kept, _) = _drop_empty_samples(raw_mat, raw_samples, group_labels, run_labels, individual_ids) + mat, raw_samples, group_labels, run_labels, individual_ids = filtered + kept < 3 && return json_error(400, "too_few_samples", + "Constrained ordination requires at least 3 non-empty samples (found $kept)") + + (norm_method, norm_depth, default_beta_transform) = _beta_normalisation(study) + use_raw_counts = cfg.ordination_distance in (:aitchison, :unweighted_unifrac, :weighted_unifrac) + if !use_raw_counts && norm_method != "none" + eff_depth = _resolve_cross_run_depth(norm_method, norm_depth, (mat,)) + (filtered_norm, nkept) = _apply_normalisation( + mat, norm_method, eff_depth, cfg.seed, + raw_samples, group_labels, run_labels, individual_ids) + mat, raw_samples, group_labels, run_labels, individual_ids = filtered_norm + nkept < 3 && return json_error(400, "too_few_samples", + "Constrained ordination requires at least 3 samples after '$norm_method' depth normalisation (found $nkept)") + end + + beta_transform = if haskey(body, :beta_transform) && !isnothing(body[:beta_transform]) + requested = lowercase(String(strip(string(body[:beta_transform])))) + if cfg.ordination_distance === :aitchison && requested != "none" + return json_error(400, "incompatible_beta_transform", + "beta_transform '$requested' is not applied to distance 'aitchison'; omit beta_transform or set it to 'none'") + end + requested + elseif cfg.ordination_distance in (:aitchison, :chisq, :unweighted_unifrac, :weighted_unifrac) + "none" + else + default_beta_transform + end + + sample_ids = _unique_ordination_sample_ids(raw_samples, run_labels) + lib_sizes = Float64[sum(@view mat[i, :]) for i in axes(mat, 1)] + log_lib_sizes = Float64[log10(max(v, 1.0)) for v in lib_sizes] + + metadata = DataFrame( + sample = sample_ids, + group = group_labels, + run = run_labels, + individual = individual_ids, + library_size = lib_sizes, + log_library_size = log_lib_sizes, + ) + custom_cols = try + _attach_custom_covariates!(metadata, get(body, :covariates, nothing), sample_ids) + catch err + err isa ArgumentError && return json_error(400, "invalid_covariates", err.msg) + rethrow() + end + + # Determine explicit or sensible study default formula before invoking R + effective_cfg = if !isnothing(cfg.ordination_formula) + cfg + else + auto_rhs = if length(unique(group_labels)) >= 2 + "~ group" + elseif length(unique(run_labels)) >= 2 + "~ run" + elseif !isempty(custom_cols) + "~ " * join(custom_cols, " + ") + else + nothing + end + isnothing(auto_rhs) ? cfg : Ordination.OrdinationConfig(; + ordination_method = cfg.ordination_method, + ordination_formula = auto_rhs, + ordination_distance = cfg.ordination_distance, + ordination_scaling = cfg.ordination_scaling, + permutations = cfg.permutations, + variance_partitioning = cfg.variance_partitioning, + replacement_method = cfg.replacement_method, + replacement_delta = cfg.replacement_delta, + seed = cfg.seed, + ) + end + + # Pre-validate formula against metadata before touching R + parsed_formula = try + Ordination.parse_ordination_formula(effective_cfg.ordination_formula, metadata) + catch err + err isa ArgumentError && return json_error(400, "invalid_ordination_formula", err.msg) + rethrow() + end + + available_covariates = String[] + for c in names(metadata) + c == "sample" && continue + length(unique(metadata[!, c])) >= 2 && push!(available_covariates, c) + end + + # Resolve UniFrac distance matrix if requested + dist_mat = nothing + tree_source = nothing + feature_ids = CladeCumulus.asv_tip_name.(sequences) + if effective_cfg.ordination_distance in (:unweighted_unifrac, :weighted_unifrac) + phylo_node, tip_map, tree_src = try + _resolve_study_unifrac_tree(study, body) + catch err + err isa ArgumentError && return json_error(400, "missing_phylogeny", err.msg) + rethrow() + end + tree_source = tree_src + on_tree = [j for j in eachindex(feature_ids) if haskey(tip_map, feature_ids[j])] + length(on_tree) >= 2 || return json_error(400, "insufficient_tree_overlap", + "Fewer than 2 ASVs from table '$table' match tips on phylogeny '$tree_src' (matched $(length(on_tree)) of $(length(feature_ids)))") + metric_sym = effective_cfg.ordination_distance === :unweighted_unifrac ? :unweighted : :weighted + dist_mat = try + CladeCumulus.unifrac_matrix( + permutedims(mat[:, on_tree]), + feature_ids[on_tree], + sample_ids, + phylo_node; + metric = metric_sym, + totals = lib_sizes, + ) + catch err + err isa ArgumentError && return json_error(400, "unifrac_failed", err.msg) + rethrow() + end + end + + !r_available() && return json_error(503, "r_unavailable", + "R/vegan is not available - install R and the pinned vegan package") + + runtime_est = Ordination.estimate_ordination_seconds(; + method = effective_cfg.ordination_method, + distance = effective_cfg.ordination_distance, + n_samples = size(mat, 1), + n_features = size(mat, 2), + permutations = effective_cfg.permutations, + n_terms = length(parsed_formula.constrained_terms), + variance_partitioning = effective_cfg.variance_partitioning, + ) + + ord_res = try + Ordination.run_constrained_ordination( + mat, + metadata, + effective_cfg; + sample_ids = sample_ids, + feature_ids = feature_ids, + beta_transform = beta_transform, + distance_matrix = dist_mat, + blocks = individual_ids, + test_axes = test_axes, + ) + catch err + err isa Ordination.VeganUnavailable && return json_error(503, "vegan_unavailable", sprint(showerror, err)) + err isa Differential.ZCompositionsUnavailable && return json_error(503, "zcompositions_unavailable", sprint(showerror, err)) + err isa RBusyError && return json_error(503, "r_busy", "R runtime is busy; retry shortly") + err isa ArgumentError && return json_error(400, "ordination_refused", err.msg) + rethrow() + end + + if !isnothing(tree_source) + ord_res.provenance["phylogeny_source"] = tree_source + end + + unique_groups = unique(group_labels) + unique_runs = unique(run_labels) + colour_by = length(unique_groups) >= 2 ? group_labels : run_labels + shape_by = (length(unique_groups) >= 2 && length(unique_runs) >= 2 && group_labels != run_labels) ? + run_labels : nothing + + biplot = Ordination.ordination_biplot_chart( + ord_res; + colour_by = colour_by, + shape_by = shape_by, + show_centroids = show_centroids, + show_features = show_features, + ) + + payload = _serializable_ordination_payload(ord_res, available_covariates, runtime_est) + payload["biplot"] = biplot + + if persist_record + src_doc = Dict{String,Any}( + "study" => study, + "table" => table, + "aggregate" => aggregate, + "n_samples" => ord_res.n_samples, + "n_features" => ord_res.n_features, + ) + cfg_doc = ord_res.provenance + result_doc = copy(payload) + saved = try + DerivedAnalyses.persist_derived_analysis( + _ordination_derived_dir(study); + source = src_doc, + config = cfg_doc, + result = result_doc, + ) + catch err + err isa DerivedAnalyses.DerivedAnalysisError && return json_error( + err.code == "conflict" ? 409 : err.code == "resource_limit" ? 413 : 500, + "derived_analysis_$(err.code)", err.reason) + rethrow() + end + payload["artifact"] = Dict{String,Any}( + "analysis_id" => saved.id, + "created" => saved.created, + ) + end + + jcs_json(payload) +end diff --git a/src/server/server.jl b/src/server/server.jl index 83abc4ff..2a543ec9 100644 --- a/src/server/server.jl +++ b/src/server/server.jl @@ -19,7 +19,7 @@ module Server using ..Databases, ..DuckDBStore, ..Validation using ..Tools, ..TaxonomyTableTools, ..ProjectSetup using ..DADA2, ..OTUPipeline - using ..DiversityMetrics, ..Analysis, ..Differential + using ..DiversityMetrics, ..Analysis, ..Differential, ..Ordination ## EPIPE log filter # HTTP.jl logs every broken-pipe error from SSE streams as @error @@ -155,6 +155,7 @@ module Server include(joinpath(@__DIR__, "routes", "phylogeny.jl")) include(joinpath(@__DIR__, "routes", "cladistics.jl")) include(joinpath(@__DIR__, "routes", "clade_cumulus.jl")) + include(joinpath(@__DIR__, "routes", "ordination.jl")) ## Bad request middleware # Helpers deep in a route (e.g. `_req_group`) throw this to reject a request. diff --git a/test/fixtures/ordination_reference.json b/test/fixtures/ordination_reference.json new file mode 100644 index 00000000..6848c2d2 --- /dev/null +++ b/test/fixtures/ordination_reference.json @@ -0,0 +1 @@ +{"datasets":[{"description":"vegan dune meadow 10-site x 8-species benchmark with A1 and Moisture","expected_rda":{"adj_r2":0.084661352984,"conditioned_inertia":0,"constrained_eigenvalues":[5.979195015336,0.761641608973],"constrained_inertia":6.740836624309,"df_conditioned":0,"df_model":2,"df_residual":7,"f_statistic":1.416213266718,"r2":0.28806994121,"total_inertia":23.4,"unconstrained_eigenvalues":[8.929798982119,4.107583292904,2.641722238559,0.807784030477,0.081896851769,0.065990315981,0.024387663883],"unconstrained_inertia":16.659163375691},"expected_varpart":{"gross_x12_adj_r2":0.084661352984,"gross_x1_adj_r2":-0.084152684191,"gross_x2_adj_r2":0.15960351377,"residual":0.915338647016,"shared":-0.009210523405,"unique_x1":-0.074942160786,"unique_x2":0.168814037175},"id":"dune_meadow_subset","x":[[2.8,1],[3.5,1],[4.3,2],[4.2,2],[6.3,1],[4.3,1],[2.8,1],[4.2,5],[3.7,4],[3.3,2]],"x1":[[2.8],[3.5],[4.3],[4.2],[6.3],[4.3],[2.8],[4.2],[3.7],[3.3]],"x2":[[1],[1],[2],[2],[1],[1],[1],[5],[4],[2]],"y":[[1,0,0,0,0,0,0,0],[3,0,0,2,0,3,4,0],[0,4,0,7,0,2,0,0],[0,8,0,2,0,2,3,0],[2,0,0,0,4,2,2,0],[2,0,0,0,3,0,0,0],[2,0,0,0,2,0,2,0],[0,4,0,5,0,0,0,0],[0,3,0,3,0,0,0,1],[4,0,0,0,4,2,4,0]],"z":null},{"description":"Gut microbiome 12-sample x 6-taxon CLR (Aitchison) ordination across antibiotic exposure and fibre","expected_rda":{"adj_r2":0.977265213493,"conditioned_inertia":0,"constrained_eigenvalues":[1.97525299637,0.043732724422],"constrained_inertia":2.018985720793,"df_conditioned":0,"df_model":2,"df_residual":9,"f_statistic":237.420019717123,"r2":0.98139881104,"total_inertia":2.05725307396,"unconstrained_eigenvalues":[0.02619735965,0.010187160782,0.001189928426,0.000525828444,0.000167075866],"unconstrained_inertia":0.038267353168},"expected_varpart":{"gross_x12_adj_r2":0.977265213493,"gross_x1_adj_r2":0.851325852215,"gross_x2_adj_r2":0.913738511672,"residual":0.022734786507,"shared":0.787799150394,"unique_x1":0.063526701821,"unique_x2":0.125939361278},"id":"microbiome_gut_clr_antibiotic","raw_counts":[[420,310,120,85,45,20],[395,330,115,90,48,22],[450,290,110,80,50,20],[410,305,130,82,49,24],[140,160,290,210,130,70],[125,175,310,195,125,70],[155,150,280,215,135,65],[135,165,305,205,120,70],[75,95,220,270,210,130],[80,90,210,280,215,125],[70,100,230,265,205,130],[85,88,215,275,212,125]],"transform":"clr","x":[[0,28],[0,26],[0,31],[0,29],[1,18],[1,19],[1,17],[1,20],[1,9],[1,8],[1,10],[1,9.5]],"x1":[[0],[0],[0],[0],[1],[1],[1],[1],[1],[1],[1],[1]],"x2":[[28],[26],[31],[29],[18],[19],[17],[20],[9],[8],[10],[9.5]],"y":[[1.405360582719,1.10167816892,0.152597614223,-0.192242872069,-0.828231638789,-1.639161855005],[1.314725151181,1.13493204074,0.080771514643,-0.16435094339,-0.792959602812,-1.573118160362],[1.481015785232,1.041649125448,0.07224856826,-0.246205162858,-0.716208792104,-1.632499523978],[1.336057699284,1.040212316193,0.187434990042,-0.27338021315,-0.788279162303,-1.502045630066],[-0.083329974732,0.050201417893,0.64490852564,0.322135133376,-0.157437946885,-0.776477155292],[-0.184933020468,0.151539216153,0.723325539709,0.259752800794,-0.184933020468,-0.764751515721],[0.020233435145,-0.012556387678,0.611597921395,0.347446346353,-0.117916903336,-0.848804411879],[-0.109813387048,0.090857308415,0.705223611121,0.307921813652,-0.227596422704,-0.766592923437],[-0.69552130174,-0.459132523676,0.380618131076,0.585412543722,0.334098115441,-0.145474964821],[-0.62842089938,-0.510637863723,0.336659996664,0.624342069116,0.360190494074,-0.182133796751],[-0.761841256707,-0.405166312769,0.427742810166,0.569393327229,0.312673480382,-0.142802048301],[-0.572731611633,-0.538046053645,0.355255160005,0.601388229544,0.341203406549,-0.187069130821]],"z":null},{"description":"Oral microbiome 12-sample x 6-taxon partial CLR RDA with Condition(batch)","expected_rda":{"adj_r2":0.965412795307,"conditioned_inertia":0.121033850097,"constrained_eigenvalues":[0.952730999572],"constrained_inertia":0.952730999572,"df_conditioned":1,"df_model":1,"df_residual":9,"f_statistic":819.669292153208,"r2":0.878720037893,"total_inertia":1.084225872277,"unconstrained_eigenvalues":[0.005322092311,0.003402710738,0.00120542606,0.000331262987,0.000199530511],"unconstrained_inertia":0.010461022608},"expected_varpart":null,"id":"microbiome_oral_partial_batch","raw_counts":[[500,250,120,70,40,20],[480,270,115,72,42,21],[520,240,110,68,41,21],[190,180,260,190,110,70],[205,175,250,195,108,67],[185,185,268,188,106,68],[310,380,150,85,50,25],[295,395,155,80,52,23],[320,370,148,86,49,27],[110,240,310,185,95,60],[118,232,305,190,97,58],[105,248,315,180,92,60]],"transform":"clr","x":[[0],[0],[0],[1],[1],[1],[0],[0],[0],[1],[1],[1]],"y":[[1.638496810292,0.945349629732,0.211380454651,-0.327616046081,-0.887231834017,-1.580379014577],[1.577786376188,1.002422231284,0.148932400649,-0.319333608698,-0.858330109431,-1.551477289991],[1.685070426523,0.91188053829,0.13172198074,-0.349250679876,-0.855186318348,-1.524235947329],[0.214247033149,0.160179811879,0.527904592004,0.214247033149,-0.332296673219,-0.784281796962],[0.294829933138,0.136605927923,0.493280871862,0.244819512563,-0.346048818876,-0.82348742661],[0.195174592947,0.195174592947,0.56580574838,0.211260730699,-0.361742138019,-0.805673526955],[1.026417475632,1.230016430873,0.300480472249,-0.267503565357,-0.798131816419,-1.491278996979],[0.990633555116,1.282543963677,0.347083315695,-0.31431516655,-0.745098082643,-1.560847585295],[1.048147539866,1.193329549711,0.277038817837,-0.265826159674,-0.828353157817,-1.424336589923],[-0.26389777846,0.51626077909,0.772194153227,0.255977680826,-0.410501252652,-0.87003358203],[-0.199300681193,0.476752066008,0.750326470949,0.277038766502,-0.395274327155,-0.909542295112],[-0.300881604061,0.558586791947,0.797730684607,0.238114896672,-0.433053377169,-0.860497391996]],"z":[[0],[0],[0],[0],[0],[0],[1],[1],[1],[1],[1],[1]]},{"description":"Marine microbiome 10-sample x 7-taxon Hellinger RDA with salinity, depth, and interaction","expected_rda":{"adj_r2":0.986681729333,"conditioned_inertia":0,"constrained_eigenvalues":[0.086941578763,0.005267434306,0.000064634526],"constrained_inertia":0.092273647595,"df_conditioned":0,"df_model":3,"df_residual":6,"f_statistic":223.254469968811,"r2":0.991121152889,"total_inertia":0.093100270664,"unconstrained_eigenvalues":[0.000768244521,0.000027123432,0.000014620449,0.000013467428,0.000002651096,5.16142e-7],"unconstrained_inertia":0.000826623069},"expected_varpart":{"gross_x12_adj_r2":0.983054634187,"gross_x1_adj_r2":0.907216758143,"gross_x2_adj_r2":0.820716537415,"residual":0.016945365813,"shared":0.744878661371,"unique_x1":0.162338096772,"unique_x2":0.075837876044},"id":"microbiome_marine_salinity_depth","raw_counts":[[610,220,90,40,20,12,8],[580,245,95,38,22,11,9],[430,310,140,65,30,15,10],[390,330,155,70,28,16,11],[240,290,230,125,65,30,20],[210,275,250,140,70,33,22],[130,195,270,195,115,58,37],[115,180,280,205,120,62,38],[65,110,240,255,175,95,60],[55,102,232,265,182,100,64]],"transform":"hellinger","x":[[30.1,5,150.5],[30.4,8,243.2],[31.8,15,477],[32.1,18,577.8],[33.5,35,1172.5],[33.9,42,1423.8],[34.8,75,2610],[35,85,2975],[35.7,140,4998],[36,160,5760]],"x1":[[30.1],[30.4],[31.8],[32.1],[33.5],[33.9],[34.8],[35],[35.7],[36]],"x2":[[5],[8],[15],[18],[35],[42],[75],[85],[140],[160]],"y":[[0.781024967591,0.469041575982,0.3,0.2,0.141421356237,0.109544511501,0.0894427191],[0.761577310586,0.494974746831,0.308220700148,0.194935886896,0.148323969742,0.104880884817,0.094868329805],[0.65574385243,0.556776436283,0.374165738677,0.25495097568,0.173205080757,0.122474487139,0.1],[0.62449979984,0.574456264654,0.393700393701,0.264575131106,0.167332005307,0.126491106407,0.104880884817],[0.489897948557,0.538516480713,0.479583152331,0.353553390593,0.25495097568,0.173205080757,0.141421356237],[0.458257569496,0.524404424085,0.5,0.374165738677,0.264575131106,0.181659021246,0.148323969742],[0.360555127546,0.441588043316,0.519615242271,0.441588043316,0.339116499156,0.240831891576,0.192353840617],[0.339116499156,0.424264068712,0.529150262213,0.452769256907,0.346410161514,0.24899799196,0.194935886896],[0.25495097568,0.331662479036,0.489897948557,0.504975246918,0.418330013267,0.308220700148,0.244948974278],[0.234520787991,0.319374388453,0.481663783152,0.514781507049,0.426614580154,0.316227766017,0.252982212813]],"z":null}],"schema_version":"1.0.0","tolerance":0.000001} diff --git a/test/references/check_ordination_references.py b/test/references/check_ordination_references.py new file mode 100755 index 00000000..8cdedaae --- /dev/null +++ b/test/references/check_ordination_references.py @@ -0,0 +1,868 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +""" +Independent Python stdlib reference oracle for Issue #19 constrained ordinations +(RDA, partial RDA with Condition(...), CLR/Aitchison and Hellinger geometries, +Ezekiel adjusted-R2, Jacobi eigenvalues, and Peres-Neto two-way variance +partitioning), plus RFC 7493 I-JSON / RFC 8785 JCS canonical verification for +the constrained-ordination schema and reference fixture. + +Zero third-party dependencies: uses only Python 3 standard library modules. +""" + +from __future__ import annotations + +import argparse +import json +import math +import pathlib +import re +import sys +from typing import Any + +SAFE_INT_MAX = (1 << 53) - 1 +REPO_ROOT = pathlib.Path(__file__).resolve().parents[2] +DEFAULT_FIXTURE = REPO_ROOT / "test" / "fixtures" / "ordination_reference.json" +DEFAULT_SCHEMA = REPO_ROOT / "docs" / "schemas" / "constrained-ordination.schema.json" +DEFAULT_KYAML = REPO_ROOT / "config" / "defaults" / "ordination.kyaml" + + +def _utf16_units(s: str) -> tuple[int, ...]: + raw = s.encode("utf-16-be") + return tuple((raw[i] << 8) | raw[i + 1] for i in range(0, len(raw), 2)) + + +def _check_ijson_str(s: str) -> None: + for ch in s: + cp = ord(ch) + if 0xD800 <= cp <= 0xDFFF: + raise ValueError(f"Lone surrogate U+{cp:04X} forbidden in I-JSON") + if 0xFDD0 <= cp <= 0xFDEF or (cp & 0xFFFE) == 0xFFFE: + raise ValueError(f"Unicode noncharacter U+{cp:04X} forbidden in I-JSON") + + +def _jcs_escape_str(s: str) -> str: + _check_ijson_str(s) + out = ['"'] + for ch in s: + cp = ord(ch) + if ch == '"': + out.append('\\"') + elif ch == "\\": + out.append("\\\\") + elif cp == 0x08: + out.append("\\b") + elif cp == 0x09: + out.append("\\t") + elif cp == 0x0A: + out.append("\\n") + elif cp == 0x0C: + out.append("\\f") + elif cp == 0x0D: + out.append("\\r") + elif cp < 0x20: + out.append(f"\\u{cp:04x}") + else: + out.append(ch) + out.append('"') + return "".join(out) + + +def _jcs_format_float(x: float) -> str: + if not math.isfinite(x): + raise ValueError("Non-finite float forbidden in I-JSON/JCS") + if x == 0.0: + return "0" + neg = x < 0.0 + ax = -x if neg else x + # Python's repr(float) uses David Gay / Ryū shortest round-trip digits + s = repr(ax) + if "e" in s or "E" in s: + mant, exp_s = s.lower().split("e") + digits = mant.replace(".", "") + exp10 = int(exp_s) + else: + if "." in s: + int_p, frac_p = s.split(".") + if int_p == "0": + stripped = frac_p.lstrip("0") + zeros = len(frac_p) - len(stripped) + digits = stripped.rstrip("0") + exp10 = -(zeros + 1) + else: + digits = (int_p + frac_p).rstrip("0") + exp10 = len(int_p) - 1 + else: + digits = s.rstrip("0") + exp10 = len(s) - 1 + + k = len(digits) + n = exp10 + 1 # ECMAScript Number::toString n (10^(n-1) <= ax < 10^n) + if k <= n <= 21: + body = digits + ("0" * (n - k)) + elif 0 < n <= 21: + body = digits[:n] + "." + digits[n:] + elif -5 <= n <= 0: + body = "0." + ("0" * (-n)) + digits + else: + sign_e = "+" if (n - 1) >= 0 else "-" + mant_e = digits if k == 1 else (digits[0] + "." + digits[1:]) + body = f"{mant_e}e{sign_e}{abs(n - 1)}" + return ("-" + body) if neg else body + + +def jcs_dumps(val: Any) -> str: + if val is None: + return "null" + if val is True: + return "true" + if val is False: + return "false" + if isinstance(val, int): + if not (-SAFE_INT_MAX <= val <= SAFE_INT_MAX): + raise ValueError(f"Integer {val} outside I-JSON safe range") + return str(val) + if isinstance(val, float): + return _jcs_format_float(val) + if isinstance(val, str): + return _jcs_escape_str(val) + if isinstance(val, (list, tuple)): + return "[" + ",".join(jcs_dumps(item) for item in val) + "]" + if isinstance(val, dict): + items: list[tuple[tuple[int, ...], str, Any]] = [] + for k, v in val.items(): + if not isinstance(k, str): + raise ValueError("JSON object keys must be strings") + items.append((_utf16_units(k), k, v)) + items.sort(key=lambda t: t[0]) + return "{" + ",".join(f"{_jcs_escape_str(k)}:{jcs_dumps(v)}" for _, k, v in items) + "}" + raise TypeError(f"Unsupported type for JCS: {type(val)!r}") + + +# --- Pure-Python statistical reference engine --- + +def _centre_cols(mat: list[list[float]]) -> list[list[float]]: + n = len(mat) + p = len(mat[0]) + means = [sum(mat[i][j] for i in range(n)) / n for j in range(p)] + return [[mat[i][j] - means[j] for j in range(p)] for i in range(n)] + + +def _mgs_orthonormal(mat: list[list[float]], tol: float = 1e-11) -> list[list[float]]: + """Return orthonormal basis columns (as list of column vectors of length n).""" + n = len(mat) + m = len(mat[0]) + basis: list[list[float]] = [] + for j in range(m): + v = [mat[i][j] for i in range(n)] + for _ in range(2): + for q in basis: + proj = sum(q[i] * v[i] for i in range(n)) + for i in range(n): + v[i] -= proj * q[i] + norm = math.sqrt(sum(x * x for x in v)) + if norm > tol: + basis.append([x / norm for x in v]) + return basis + + +def _project(basis: list[list[float]], mat: list[list[float]]) -> list[list[float]]: + n = len(mat) + p = len(mat[0]) + out = [[0.0] * p for _ in range(n)] + for q in basis: + for j in range(p): + dot = sum(q[i] * mat[i][j] for i in range(n)) + for i in range(n): + out[i][j] += q[i] * dot + return out + + +def _sub_mat(a: list[list[float]], b: list[list[float]]) -> list[list[float]]: + return [[a[i][j] - b[i][j] for j in range(len(a[0]))] for i in range(len(a))] + + +def _frob_sq(mat: list[list[float]]) -> float: + return sum(x * x for row in mat for x in row) + + +def _jacobi_sym_eigenvalues(sym: list[list[float]], tol: float = 1e-13) -> list[float]: + m = len(sym) + if m == 0: + return [] + if m == 1: + return [sym[0][0]] if sym[0][0] > tol else [] + a = [[0.5 * (sym[i][j] + sym[j][i]) for j in range(m)] for i in range(m)] + for _ in range(90): + max_off = 0.0 + for p in range(m - 1): + for q in range(p + 1, m): + if abs(a[p][q]) > max_off: + max_off = abs(a[p][q]) + if max_off <= tol: + break + for p in range(m - 1): + for q in range(p + 1, m): + apq = a[p][q] + if abs(apq) <= tol: + continue + tau = (a[q][q] - a[p][p]) / (2.0 * apq) + t = (1.0 / (tau + math.sqrt(1.0 + tau * tau))) if tau >= 0.0 else (-1.0 / (-tau + math.sqrt(1.0 + tau * tau))) + c = 1.0 / math.sqrt(1.0 + t * t) + s = t * c + app = a[p][p] + aqq = a[q][q] + a[p][p] = app - t * apq + a[q][q] = aqq + t * apq + a[p][q] = 0.0 + a[q][p] = 0.0 + for r in range(m): + if r == p or r == q: + continue + arp = a[r][p] + arq = a[r][q] + a[r][p] = c * arp - s * arq + a[p][r] = a[r][p] + a[r][q] = s * arp + c * arq + a[q][r] = a[r][q] + vals = sorted((a[i][i] for i in range(m)), reverse=True) + cutoff = max(tol, 1e-12 * max(1.0, vals[0])) if vals else tol + return [v for v in vals if v > cutoff] + + +def _gram_eigenvalues(mat: list[list[float]], dof: float) -> list[float]: + n = len(mat) + p = len(mat[0]) + if n <= p: + g = [[sum(mat[i][k] * mat[j][k] for k in range(p)) / dof for j in range(n)] for i in range(n)] + else: + g = [[sum(mat[k][i] * mat[k][j] for k in range(n)) / dof for j in range(p)] for i in range(p)] + return _jacobi_sym_eigenvalues(g) + + +def ezekiel_adj_r2(r2: float, n: int, m: int) -> float: + return 1.0 - (1.0 - r2) * (n - 1) / (n - 1 - m) + + +def clr_transform(counts: list[list[float]]) -> list[list[float]]: + out: list[list[float]] = [] + for row in counts: + logs = [math.log(v) for v in row] + mean_log = sum(logs) / len(logs) + out.append([lv - mean_log for lv in logs]) + return out + + +def hellinger_transform(counts: list[list[float]]) -> list[list[float]]: + out: list[list[float]] = [] + for row in counts: + tot = sum(row) + out.append([math.sqrt(v / tot) for v in row]) + return out + + +def reference_linear_rda( + y: list[list[float]], + x: list[list[float]], + z: list[list[float]] | None = None, +) -> dict[str, Any]: + n = len(y) + dof = float(n - 1) + yc = _centre_cols(y) + xc = _centre_cols(x) + tot_inertia = _frob_sq(yc) / dof + + cond_inertia = 0.0 + q_rank = 0 + y_work = [row[:] for row in yc] + x_work = [row[:] for row in xc] + if z is not None: + zc = _centre_cols(z) + qz = _mgs_orthonormal(zc) + q_rank = len(qz) + if q_rank > 0: + yz = _project(qz, y_work) + cond_inertia = _frob_sq(yz) / dof + y_work = _sub_mat(y_work, yz) + xz = _project(qz, x_work) + x_work = _sub_mat(x_work, xz) + + qx = _mgs_orthonormal(x_work) + k_rank = len(qx) + df_res = n - 1 - q_rank - k_rank + y_hat = _project(qx, y_work) + y_res = _sub_mat(y_work, y_hat) + con_inertia = _frob_sq(y_hat) / dof + unc_inertia = _frob_sq(y_res) / dof + r2 = con_inertia / tot_inertia + if q_rank == 0: + adj_r2 = ezekiel_adj_r2(r2, n, k_rank) + else: + r2_xz = (con_inertia + cond_inertia) / tot_inertia + r2_z = cond_inertia / tot_inertia + adj_r2 = ezekiel_adj_r2(r2_xz, n, k_rank + q_rank) - ezekiel_adj_r2(r2_z, n, q_rank) + f_stat = (con_inertia / k_rank) / (unc_inertia / df_res) + con_eigs = [round(v, 12) for v in _gram_eigenvalues(y_hat, dof)] + unc_eigs = [round(v, 12) for v in _gram_eigenvalues(y_res, dof)] + + return { + "adj_r2": round(adj_r2, 12), + "conditioned_inertia": round(cond_inertia, 12), + "constrained_eigenvalues": con_eigs, + "constrained_inertia": round(con_inertia, 12), + "df_conditioned": q_rank, + "df_model": k_rank, + "df_residual": df_res, + "f_statistic": round(f_stat, 12), + "r2": round(r2, 12), + "total_inertia": round(tot_inertia, 12), + "unconstrained_eigenvalues": unc_eigs, + "unconstrained_inertia": round(unc_inertia, 12), + } + + +def reference_varpart_twoway( + y: list[list[float]], + x1: list[list[float]], + x2: list[list[float]], +) -> dict[str, float]: + r_a = reference_linear_rda(y, x1)["adj_r2"] + r_b = reference_linear_rda(y, x2)["adj_r2"] + x12 = [r1 + r2 for r1, r2 in zip(x1, x2)] + r_ab = reference_linear_rda(y, x12)["adj_r2"] + uniq_x1 = r_ab - r_b + uniq_x2 = r_ab - r_a + shared = r_a + r_b - r_ab + resid = 1.0 - r_ab + return { + "gross_x1_adj_r2": round(r_a, 12), + "gross_x12_adj_r2": round(r_ab, 12), + "gross_x2_adj_r2": round(r_b, 12), + "residual": round(resid, 12), + "shared": round(shared, 12), + "unique_x1": round(uniq_x1, 12), + "unique_x2": round(uniq_x2, 12), + } + + +def build_reference_document() -> dict[str, Any]: + # 1. Ecological benchmark: first 10 sites x 8 species of vegan's dune meadow dataset + dune_y = [ + [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [3.0, 0.0, 0.0, 2.0, 0.0, 3.0, 4.0, 0.0], + [0.0, 4.0, 0.0, 7.0, 0.0, 2.0, 0.0, 0.0], + [0.0, 8.0, 0.0, 2.0, 0.0, 2.0, 3.0, 0.0], + [2.0, 0.0, 0.0, 0.0, 4.0, 2.0, 2.0, 0.0], + [2.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0], + [2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 0.0], + [0.0, 4.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 3.0, 0.0, 3.0, 0.0, 0.0, 0.0, 1.0], + [4.0, 0.0, 0.0, 0.0, 4.0, 2.0, 4.0, 0.0], + ] + dune_a1 = [[2.8], [3.5], [4.3], [4.2], [6.3], [4.3], [2.8], [4.2], [3.7], [3.3]] + dune_moist = [[1.0], [1.0], [2.0], [2.0], [1.0], [1.0], [1.0], [5.0], [4.0], [2.0]] + dune_x = [[a[0], m[0]] for a, m in zip(dune_a1, dune_moist)] + + # 2. Microbiome dataset 1: Gut amplicon counts -> CLR (Aitchison geometry) + gut_counts = [ + [420.0, 310.0, 120.0, 85.0, 45.0, 20.0], + [395.0, 330.0, 115.0, 90.0, 48.0, 22.0], + [450.0, 290.0, 110.0, 80.0, 50.0, 20.0], + [410.0, 305.0, 130.0, 82.0, 49.0, 24.0], + [140.0, 160.0, 290.0, 210.0, 130.0, 70.0], + [125.0, 175.0, 310.0, 195.0, 125.0, 70.0], + [155.0, 150.0, 280.0, 215.0, 135.0, 65.0], + [135.0, 165.0, 305.0, 205.0, 120.0, 70.0], + [75.0, 95.0, 220.0, 270.0, 210.0, 130.0], + [80.0, 90.0, 210.0, 280.0, 215.0, 125.0], + [70.0, 100.0, 230.0, 265.0, 205.0, 130.0], + [85.0, 88.0, 215.0, 275.0, 212.0, 125.0], + ] + gut_y = clr_transform(gut_counts) + gut_abx = [[0.0], [0.0], [0.0], [0.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0]] + gut_fibre = [[28.0], [26.0], [31.0], [29.0], [18.0], [19.0], [17.0], [20.0], [9.0], [8.0], [10.0], [9.5]] + gut_x = [[a[0], f[0]] for a, f in zip(gut_abx, gut_fibre)] + + # 3. Microbiome dataset 2: Oral microbiome CLR with Condition(batch) partial RDA + oral_counts = [ + [500.0, 250.0, 120.0, 70.0, 40.0, 20.0], + [480.0, 270.0, 115.0, 72.0, 42.0, 21.0], + [520.0, 240.0, 110.0, 68.0, 41.0, 21.0], + [190.0, 180.0, 260.0, 190.0, 110.0, 70.0], + [205.0, 175.0, 250.0, 195.0, 108.0, 67.0], + [185.0, 185.0, 268.0, 188.0, 106.0, 68.0], + [310.0, 380.0, 150.0, 85.0, 50.0, 25.0], + [295.0, 395.0, 155.0, 80.0, 52.0, 23.0], + [320.0, 370.0, 148.0, 86.0, 49.0, 27.0], + [110.0, 240.0, 310.0, 185.0, 95.0, 60.0], + [118.0, 232.0, 305.0, 190.0, 97.0, 58.0], + [105.0, 248.0, 315.0, 180.0, 92.0, 60.0], + ] + oral_y = clr_transform(oral_counts) + oral_status = [[0.0], [0.0], [0.0], [1.0], [1.0], [1.0], [0.0], [0.0], [0.0], [1.0], [1.0], [1.0]] + oral_batch = [[0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [1.0], [1.0], [1.0], [1.0], [1.0], [1.0]] + + # 4. Microbiome dataset 3: Marine Hellinger-transformed community across salinity, depth & interaction + marine_counts = [ + [610.0, 220.0, 90.0, 40.0, 20.0, 12.0, 8.0], + [580.0, 245.0, 95.0, 38.0, 22.0, 11.0, 9.0], + [430.0, 310.0, 140.0, 65.0, 30.0, 15.0, 10.0], + [390.0, 330.0, 155.0, 70.0, 28.0, 16.0, 11.0], + [240.0, 290.0, 230.0, 125.0, 65.0, 30.0, 20.0], + [210.0, 275.0, 250.0, 140.0, 70.0, 33.0, 22.0], + [130.0, 195.0, 270.0, 195.0, 115.0, 58.0, 37.0], + [115.0, 180.0, 280.0, 205.0, 120.0, 62.0, 38.0], + [65.0, 110.0, 240.0, 255.0, 175.0, 95.0, 60.0], + [55.0, 102.0, 232.0, 265.0, 182.0, 100.0, 64.0], + ] + marine_y = hellinger_transform(marine_counts) + marine_sal = [[30.1], [30.4], [31.8], [32.1], [33.5], [33.9], [34.8], [35.0], [35.7], [36.0]] + marine_dep = [[5.0], [8.0], [15.0], [18.0], [35.0], [42.0], [75.0], [85.0], [140.0], [160.0]] + marine_x = [[s[0], d[0], s[0] * d[0]] for s, d in zip(marine_sal, marine_dep)] + + return { + "datasets": [ + { + "description": "vegan dune meadow 10-site x 8-species benchmark with A1 and Moisture", + "expected_rda": reference_linear_rda(dune_y, dune_x), + "expected_varpart": reference_varpart_twoway(dune_y, dune_a1, dune_moist), + "id": "dune_meadow_subset", + "x": [[round(v, 10) for v in row] for row in dune_x], + "x1": dune_a1, + "x2": dune_moist, + "y": dune_y, + "z": None, + }, + { + "description": "Gut microbiome 12-sample x 6-taxon CLR (Aitchison) ordination across antibiotic exposure and fibre", + "expected_rda": reference_linear_rda(gut_y, gut_x), + "expected_varpart": reference_varpart_twoway(gut_y, gut_abx, gut_fibre), + "id": "microbiome_gut_clr_antibiotic", + "raw_counts": gut_counts, + "transform": "clr", + "x": gut_x, + "x1": gut_abx, + "x2": gut_fibre, + "y": [[round(v, 12) for v in row] for row in gut_y], + "z": None, + }, + { + "description": "Oral microbiome 12-sample x 6-taxon partial CLR RDA with Condition(batch)", + "expected_rda": reference_linear_rda(oral_y, oral_status, oral_batch), + "expected_varpart": None, + "id": "microbiome_oral_partial_batch", + "raw_counts": oral_counts, + "transform": "clr", + "x": oral_status, + "y": [[round(v, 12) for v in row] for row in oral_y], + "z": oral_batch, + }, + { + "description": "Marine microbiome 10-sample x 7-taxon Hellinger RDA with salinity, depth, and interaction", + "expected_rda": reference_linear_rda(marine_y, marine_x), + "expected_varpart": reference_varpart_twoway(marine_y, marine_sal, marine_dep), + "id": "microbiome_marine_salinity_depth", + "raw_counts": marine_counts, + "transform": "hellinger", + "x": [[round(v, 10) for v in row] for row in marine_x], + "x1": marine_sal, + "x2": marine_dep, + "y": [[round(v, 12) for v in row] for row in marine_y], + "z": None, + }, + ], + "schema_version": "1.0.0", + "tolerance": 0.000001, + } + + +def build_schema_document() -> dict[str, Any]: + return { + "$id": "https://github.com/hyperpolymath/MetaManifold-WebUI/blob/main/docs/schemas/constrained-ordination.schema.json", + "$schema": "https://json-schema.org/draft/2020-12/schema", + "additionalProperties": False, + "description": "Canonical RFC 7493 I-JSON / RFC 8785 JCS response schema for MetaManifold constrained ordinations (RDA, CCA, CAP, dbRDA).", + "properties": { + "adj_r2": {"type": ["number", "null"]}, + "artifact": { + "additionalProperties": False, + "properties": { + "analysis_id": {"pattern": "^[0-9a-f]{64}$", "type": "string"}, + "created": {"type": "boolean"}, + }, + "required": ["analysis_id", "created"], + "type": "object", + }, + "available_covariates": {"items": {"type": "string"}, "type": "array"}, + "axis_labels": {"items": {"type": "string"}, "maxItems": 2, "minItems": 2, "type": "array"}, + "axis_names": {"items": {"type": "string"}, "maxItems": 2, "minItems": 2, "type": "array"}, + "axis_tests": {"$ref": "#/$defs/termTestArray"}, + "biplot": { + "additionalProperties": False, + "properties": { + "data": {"items": {"type": "object"}, "type": "array"}, + "layout": {"type": "object"}, + }, + "required": ["data", "layout"], + "type": "object", + }, + "biplot_arrows": { + "items": { + "additionalProperties": False, + "properties": { + "axis1": {"type": "number"}, + "axis2": {"type": "number"}, + "label": {"type": "string"}, + "length": {"minimum": 0, "type": "number"}, + }, + "required": ["label", "axis1", "axis2", "length"], + "type": "object", + }, + "type": "array", + }, + "cautions": {"items": {"type": "string"}, "type": "array"}, + "centroids": { + "items": { + "additionalProperties": False, + "properties": { + "axis1": {"type": "number"}, + "axis2": {"type": "number"}, + "label": {"type": "string"}, + }, + "required": ["label", "axis1", "axis2"], + "type": "object", + }, + "type": "array", + }, + "condition_terms": {"items": {"type": "string"}, "type": "array"}, + "conditioned_inertia": {"minimum": 0, "type": "number"}, + "conditioned_proportion": {"minimum": 0, "type": "number"}, + "constrained_eigenvalues": {"$ref": "#/$defs/eigenvalueArray"}, + "constrained_inertia": {"minimum": 0, "type": "number"}, + "constrained_proportion": {"minimum": 0, "type": "number"}, + "constrained_terms": {"items": {"type": "string"}, "minItems": 1, "type": "array"}, + "distance": { + "enum": [ + "bray", + "jaccard", + "euclidean", + "aitchison", + "chisq", + "unweighted_unifrac", + "weighted_unifrac", + ], + "type": "string", + }, + "formula": {"minLength": 2, "type": "string"}, + "imaginary_inertia": {"minimum": 0, "type": "number"}, + "marginal_tests": {"$ref": "#/$defs/termTestArray"}, + "method": {"enum": ["rda", "cca", "cap", "dbrda"], "type": "string"}, + "n_features": {"minimum": 0, "type": "integer"}, + "n_samples": {"minimum": 3, "type": "integer"}, + "omnibus_test": { + "additionalProperties": False, + "properties": { + "df_model": {"minimum": 1, "type": "integer"}, + "df_residual": {"minimum": 1, "type": "integer"}, + "f_statistic": {"type": ["number", "null"]}, + "min_p_value": {"exclusiveMinimum": 0, "maximum": 1, "type": "number"}, + "p_value": {"type": ["number", "null"]}, + "permutations": {"maximum": 9999, "minimum": 1, "type": "integer"}, + "text": {"type": "string"}, + "untestable_reason": {"type": ["string", "null"]}, + "variance_model": {"type": ["number", "null"]}, + "variance_residual": {"type": ["number", "null"]}, + }, + "required": [ + "df_model", + "df_residual", + "variance_model", + "variance_residual", + "f_statistic", + "p_value", + "untestable_reason", + "permutations", + "min_p_value", + "text", + ], + "type": "object", + }, + "provenance": { + "additionalProperties": True, + "properties": { + "beta_normalisation": {"type": "string"}, + "blocked": {"type": "boolean"}, + "condition_terms": {"items": {"type": "string"}, "type": "array"}, + "constrained_terms": {"items": {"type": "string"}, "type": "array"}, + "distance": {"type": "string"}, + "dropped_zero_features": {"minimum": 0, "type": "integer"}, + "formula": {"type": "string"}, + "method": {"enum": ["rda", "cca", "cap", "dbrda"], "type": "string"}, + "multiple_testing_policy": {"type": "string"}, + "n_features": {"minimum": 0, "type": "integer"}, + "n_samples": {"minimum": 3, "type": "integer"}, + "permutations_performed": {"minimum": 1, "type": "integer"}, + "permutations_requested": {"maximum": 9999, "minimum": 99, "type": "integer"}, + "rng_kind": {"type": "string"}, + "scaling": {"enum": [1, 2], "type": "integer"}, + "scaling_meaning": {"enum": ["distance", "correlation"], "type": "string"}, + "seed": {"type": "integer"}, + "software": { + "additionalProperties": False, + "properties": { + "r_package": {"const": "vegan", "type": "string"}, + "vegan_version": {"type": "string"}, + }, + "required": ["r_package", "vegan_version"], + "type": "object", + }, + "zero_replacement": {"type": ["object", "null"]}, + }, + "required": [ + "method", + "distance", + "scaling", + "scaling_meaning", + "formula", + "permutations_requested", + "permutations_performed", + "seed", + "multiple_testing_policy", + "software", + ], + "type": "object", + }, + "r2": {"type": ["number", "null"]}, + "runtime_estimate": { + "additionalProperties": False, + "properties": { + "estimated_seconds": {"minimum": 0, "type": "number"}, + "tier": {"enum": ["fast", "moderate", "heavy"], "type": "string"}, + "warn_five_minutes": {"type": "boolean"}, + }, + "required": ["estimated_seconds", "tier", "warn_five_minutes"], + "type": "object", + }, + "scaling": {"enum": [1, 2], "type": "integer"}, + "sites": { + "items": { + "additionalProperties": False, + "properties": { + "axis1": {"type": "number"}, + "axis2": {"type": "number"}, + "sample": {"type": "string"}, + }, + "required": ["sample", "axis1", "axis2"], + "type": "object", + }, + "minItems": 3, + "type": "array", + }, + "term_tests": {"$ref": "#/$defs/termTestArray"}, + "top_feature_scores": {"items": {"type": "object"}, "type": "array"}, + "total_inertia": {"minimum": 0, "type": "number"}, + "unconstrained_eigenvalues": {"$ref": "#/$defs/eigenvalueArray"}, + "unconstrained_inertia": {"minimum": 0, "type": "number"}, + "unconstrained_proportion": {"minimum": 0, "type": "number"}, + "variance_partitioning": {"type": ["object", "null"]}, + }, + "required": [ + "method", + "distance", + "scaling", + "formula", + "constrained_terms", + "condition_terms", + "n_samples", + "n_features", + "total_inertia", + "constrained_inertia", + "unconstrained_inertia", + "r2", + "adj_r2", + "sites", + "biplot_arrows", + "omnibus_test", + "term_tests", + "provenance", + "biplot", + ], + "title": "ConstrainedOrdinationResult", + "type": "object", + "$defs": { + "eigenvalueArray": { + "items": { + "additionalProperties": False, + "properties": { + "axis": {"type": "string"}, + "eigenvalue": {"minimum": 0, "type": "number"}, + "proportion_constrained": {"minimum": 0, "type": "number"}, + "proportion_total": {"minimum": 0, "type": "number"}, + }, + "required": ["axis", "eigenvalue", "proportion_total"], + "type": "object", + }, + "type": "array", + }, + "termTestArray": { + "items": { + "additionalProperties": False, + "properties": { + "df": {"minimum": 0, "type": "integer"}, + "f_statistic": {"type": ["number", "null"]}, + "p_value": {"type": ["number", "null"]}, + "term": {"type": "string"}, + "untestable_reason": {"type": ["string", "null"]}, + "variance": {"type": ["number", "null"]}, + }, + "required": ["term", "df", "variance", "f_statistic", "p_value", "untestable_reason"], + "type": "object", + }, + "type": "array", + }, + }, + } + + +def verify_jcs_file(path: pathlib.Path) -> Any: + raw = path.read_bytes() + if not raw.endswith(b"\n") or raw.endswith(b"\n\n"): + raise AssertionError(f"{path}: must end with exactly one trailing LF") + text = raw[:-1].decode("utf-8") + parsed = json.loads(text) + canonical = jcs_dumps(parsed) + if text != canonical: + raise AssertionError(f"{path}: content is not byte-for-byte RFC 8785 JCS") + return parsed + + +def verify_kyaml_contract(path: pathlib.Path) -> None: + text = path.read_text(encoding="utf-8") + lines = text.splitlines() + body = [ln for ln in lines if ln.strip() and not ln.strip().startswith("#")] + if not body or body[0].strip() != "---": + raise AssertionError(f"{path}: KYAML must open with explicit '---' after SPDX comments") + if body[1].strip() != "{": + raise AssertionError(f"{path}: KYAML top-level mapping must use explicit '{{' flow syntax") + if body[-1].strip() != "}": + raise AssertionError(f"{path}: KYAML top-level mapping must close with '}}'") + for req_token in ("rda:", "cca:", "cap:", "dbrda:", '"capscale"', '"aitchison"', '"chisq"', "99", "9999"): + if req_token not in text: + raise AssertionError(f"{path}: missing expected contract token {req_token}") + # Check there are no bare unquoted string values after ':' + for idx, ln in enumerate(lines, start=1): + stripped = ln.strip() + if not stripped or stripped.startswith("#") or stripped in ("---", "{", "}", "},", "],"): + continue + m = re.match(r"^[A-Za-z0-9_\"']+:\s*(.+?)\s*,?$", stripped) + if m: + val = m.group(1).rstrip(",").strip() + if val in ("{", "[") or val.startswith('"') or val in ("true", "false", "null"): + continue + if re.fullmatch(r"-?[0-9]+(\.[0-9]+)?", val): + continue + raise AssertionError(f"{path}:{idx}: unquoted scalar value {val!r} violates KYAML") + + +def verify_fixture_math(doc: dict[str, Any], tol: float = 1e-9) -> None: + datasets = doc.get("datasets", []) + if len(datasets) < 4: + raise AssertionError(f"Expected >= 4 reference datasets, found {len(datasets)}") + for ds in datasets: + ds_id = ds["id"] + y = ds["y"] + if ds.get("transform") == "clr": + y_recomp = clr_transform(ds["raw_counts"]) + elif ds.get("transform") == "hellinger": + y_recomp = hellinger_transform(ds["raw_counts"]) + else: + y_recomp = y + for r1, r2 in zip(y, y_recomp): + for a, b in zip(r1, r2): + if abs(a - b) > 1e-10: + raise AssertionError(f"{ds_id}: stored transformed matrix drifted from raw_counts") + + actual_rda = reference_linear_rda(y_recomp, ds["x"], ds["z"]) + exp_rda = ds["expected_rda"] + for key in ( + "total_inertia", + "conditioned_inertia", + "constrained_inertia", + "unconstrained_inertia", + "r2", + "adj_r2", + "f_statistic", + ): + if abs(actual_rda[key] - exp_rda[key]) > tol: + raise AssertionError( + f"{ds_id}.{key}: expected {exp_rda[key]}, recomputed {actual_rda[key]}" + ) + for eig_key in ("constrained_eigenvalues", "unconstrained_eigenvalues"): + if len(actual_rda[eig_key]) != len(exp_rda[eig_key]): + raise AssertionError(f"{ds_id}.{eig_key}: length mismatch") + for a, b in zip(actual_rda[eig_key], exp_rda[eig_key]): + if abs(a - b) > tol: + raise AssertionError(f"{ds_id}.{eig_key}: {a} != {b}") + + # Verify inertia conservation: total == conditioned + constrained + unconstrained + inertia_sum = ( + actual_rda["conditioned_inertia"] + + actual_rda["constrained_inertia"] + + actual_rda["unconstrained_inertia"] + ) + if abs(actual_rda["total_inertia"] - inertia_sum) > 1e-10: + raise AssertionError(f"{ds_id}: inertia partition does not sum to total_inertia") + + if ds["expected_varpart"] is not None: + actual_vp = reference_varpart_twoway(y_recomp, ds["x1"], ds["x2"]) + exp_vp = ds["expected_varpart"] + for k, v in exp_vp.items(): + if abs(actual_vp[k] - v) > tol: + raise AssertionError(f"{ds_id}.varpart.{k}: expected {v}, got {actual_vp[k]}") + # Verify Peres-Neto identity: [a] + [b] + [c] + [d] == 1 + part_sum = ( + actual_vp["unique_x1"] + + actual_vp["unique_x2"] + + actual_vp["shared"] + + actual_vp["residual"] + ) + if abs(part_sum - 1.0) > 1e-10: + raise AssertionError(f"{ds_id}: varpart fractions sum to {part_sum}, expected 1.0") + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--generate", action="store_true", help="Write canonical JCS schema and reference fixture") + ap.add_argument("--check", type=pathlib.Path, default=DEFAULT_FIXTURE, help="Fixture path to verify") + args = ap.parse_args() + + if args.generate: + DEFAULT_SCHEMA.parent.mkdir(parents=True, exist_ok=True) + DEFAULT_FIXTURE.parent.mkdir(parents=True, exist_ok=True) + DEFAULT_SCHEMA.write_bytes((jcs_dumps(build_schema_document()) + "\n").encode("utf-8")) + DEFAULT_FIXTURE.write_bytes((jcs_dumps(build_reference_document()) + "\n").encode("utf-8")) + + # Cross-check our JCS serializer against the repo's jcs_oracle.jsonl first + oracle_path = REPO_ROOT / "test" / "fixtures" / "jcs_oracle.jsonl" + if oracle_path.exists(): + for idx, line in enumerate(oracle_path.read_text(encoding="utf-8").splitlines(), start=1): + if not line.strip(): + continue + parsed = json.loads(line) + if jcs_dumps(parsed["input"]) != parsed["expected"]: + raise AssertionError(f"jcs_oracle.jsonl line {idx} mismatch") + + verify_kyaml_contract(DEFAULT_KYAML) + verify_jcs_file(DEFAULT_SCHEMA) + fixture_doc = verify_jcs_file(args.check) + verify_fixture_math(fixture_doc) + print( + f"OK: verified {len(fixture_doc['datasets'])} ordination reference datasets, " + f"JCS schema ({DEFAULT_SCHEMA.name}), and KYAML contract ({DEFAULT_KYAML.name})." + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/test/runtests.jl b/test/runtests.jl index d04c468c..4763ea0a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,7 @@ const UNIT_FILES = [ "test_evidence.jl", "test_evidence_routes.jl", "test_permanova.jl", + "test_ordination.jl", "test_write_gates.jl", "test_server_security.jl", "test_canonical.jl", diff --git a/test/unit/test_ordination.jl b/test/unit/test_ordination.jl new file mode 100644 index 00000000..823c3170 --- /dev/null +++ b/test/unit/test_ordination.jl @@ -0,0 +1,478 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +using Test +using DataFrames +using DuckDB +using DBInterface +using JSON3 +using MetaManifold +using MetaManifold.Analysis: r_available +using MetaManifold.Canonical +using MetaManifold.CladeCumulus +using MetaManifold.Ordination +using MetaManifold.RRuntime +using MetaManifold.Validation +using RCall + +const _ORD_SV = MetaManifold.Server +const _ORD_REPO_ROOT = normpath(joinpath(@__DIR__, "..", "..")) +const _ORD_REF_PATH = joinpath(_ORD_REPO_ROOT, "test", "fixtures", "ordination_reference.json") +const _ORD_PY_ORACLE = joinpath(_ORD_REPO_ROOT, "test", "references", "check_ordination_references.py") + +function _ord_http(method::String, path::String, body = nothing; raw::Bool = false) + payload = isnothing(body) ? UInt8[] : (raw ? body : jcs_bytes(body)) + headers = isnothing(body) ? Pair{String,String}[] : ["Content-Type" => "application/json"] + req = _ORD_SV.HTTP.Request(method, path, headers, payload) + _ORD_SV.Oxygen.internalrequest(req) +end + +function _ord_with_study_root(f::Function) + root = mktempdir() + saved = _ORD_SV.ServerState._root[] + try + mkpath(joinpath(root, "config")) + cp(joinpath(_ORD_REPO_ROOT, "config", "defaults"), joinpath(root, "config", "defaults")) + for (grp, rname, counts_rows) in ( + ("Control", "run1", [ + ("seq1", "ACGTACGTACGTACGT", "Bacteroides", 410, 395, 430), + ("seq2", "TGCAACGTACGTACGT", "Prevotella", 290, 310, 275), + ("seq3", "GGCCACGTACGTACGT", "Faecalibacterium", 120, 115, 125), + ("seq4", "TTAAACGTACGTACGT", "Roseburia", 80, 90, 75), + ]), + ("Treatment", "run2", [ + ("seq1", "ACGTACGTACGTACGT", "Bacteroides", 95, 110, 88), + ("seq2", "TGCAACGTACGTACGT", "Prevotella", 140, 130, 150), + ("seq3", "GGCCACGTACGTACGT", "Faecalibacterium", 320, 305, 335), + ("seq4", "TTAAACGTACGTACGT", "Roseburia", 265, 280, 255), + ]), + ) + mkpath(joinpath(root, "data", "study1", grp, rname)) + for s in ("s1", "s2", "s3"), rd in ("R1", "R2") + touch(joinpath(root, "data", "study1", grp, rname, "$(s)_L001_$(rd)_001.fastq.gz")) + end + merged_dir = joinpath(root, "projects", "study1", grp, rname, "merged") + mkpath(merged_dir) + db = DuckDB.DB(joinpath(merged_dir, "results.duckdb")) + con = DBInterface.connect(db) + DBInterface.execute(con, + "CREATE TABLE merged (SeqName VARCHAR, Sequence VARCHAR, Genus_dada2 VARCHAR, s1 BIGINT, s2 BIGINT, s3 BIGINT)") + for (sn, sq, gn, c1, c2, c3) in counts_rows + DBInterface.execute(con, + "INSERT INTO merged VALUES ('$sn', '$sq', '$gn', $c1, $c2, $c3)") + end + DBInterface.close!(con) + close(db) + end + _ORD_SV.ServerState.set_root!(root) + f(root) + finally + _ORD_SV.ServerState._root[] = saved + rm(root; recursive = true, force = true) + end +end + +function _json_matrix(rows)::Matrix{Float64} + n = length(rows) + p = length(rows[1]) + m = Matrix{Float64}(undef, n, p) + for i in 1:n, j in 1:p + m[i, j] = Float64(rows[i][j]) + end + m +end + +@testset "Constrained ordinations (Issue #19)" begin + + @testset "Independent Python reference oracle, KYAML contract, and JCS schema" begin + if !isnothing(Sys.which("python3")) + proc = run(`python3 $_ORD_PY_ORACLE --check $_ORD_REF_PATH`) + @test proc.exitcode == 0 + end + contract = load_ordination_contract() + @test contract["schema_version"] == 1 + @test haskey(contract["methods"], "rda") + @test haskey(contract["methods"], "cca") + @test haskey(contract["methods"], "cap") + @test haskey(contract["methods"], "dbrda") + @test contract["permutation_bounds"]["min"] == 99 + @test contract["permutation_bounds"]["max"] == 9999 + # Contract serialises cleanly to canonical RFC 8785 JCS + raw_jcs = jcs_bytes(contract) + @test !isempty(raw_jcs) + @test ijson_read(raw_jcs)[:schema_version] == 1 + end + + @testset "OrdinationConfig and method/distance compatibility validation" begin + cfg_default = OrdinationConfig() + @test cfg_default.ordination_method === :none + @test cfg_default.ordination_distance === :bray + @test cfg_default.ordination_scaling == 2 + @test cfg_default.permutations == 999 + @test !cfg_default.variance_partitioning + + # RDA requires linear Euclidean / Aitchison CLR space, rejecting Bray-Curtis & UniFrac + @test OrdinationConfig(; ordination_method=:rda, ordination_distance=:aitchison).ordination_distance === :aitchison + @test OrdinationConfig(; ordination_method=:rda, ordination_distance=:euclidean).ordination_distance === :euclidean + @test_throws ArgumentError OrdinationConfig(; ordination_method=:rda, ordination_distance=:bray) + @test_throws ArgumentError OrdinationConfig(; ordination_method=:rda, ordination_distance=:weighted_unifrac) + + # CCA requires chi-square geometry + @test OrdinationConfig(; ordination_method=:cca, ordination_distance=:chisq).ordination_distance === :chisq + @test_throws ArgumentError OrdinationConfig(; ordination_method=:cca, ordination_distance=:bray) + + # CAP and dbRDA accept Bray-Curtis, Jaccard, Euclidean, Aitchison, and UniFrac + for m in (:cap, :dbrda), d in (:bray, :jaccard, :euclidean, :aitchison, :unweighted_unifrac, :weighted_unifrac) + c = OrdinationConfig(; ordination_method=m, ordination_distance=d) + @test c.ordination_method === m + @test c.ordination_distance === d + end + @test_throws ArgumentError OrdinationConfig(; ordination_method=:dbrda, ordination_distance=:chisq) + + # Permutation bounds [99, 9999] and scaling {1, 2} + @test OrdinationConfig(; ordination_method=:dbrda, permutations=99).permutations == 99 + @test OrdinationConfig(; ordination_method=:dbrda, permutations=9999).permutations == 9999 + @test_throws ArgumentError OrdinationConfig(; ordination_method=:dbrda, permutations=98) + @test_throws ArgumentError OrdinationConfig(; ordination_method=:dbrda, permutations=10000) + @test OrdinationConfig(; ordination_method=:dbrda, ordination_scaling=1).ordination_scaling == 1 + @test_throws ArgumentError OrdinationConfig(; ordination_method=:dbrda, ordination_scaling=3) + @test_throws ArgumentError OrdinationConfig(; ordination_method=:rda, ordination_distance=:aitchison, replacement_delta=0.0) + @test_throws ArgumentError OrdinationConfig(; ordination_method=:rda, ordination_distance=:aitchison, replacement_delta=1.0) + end + + @testset "Formula AST parsing, Condition(...) partial terms, and collinearity guards" begin + meta = DataFrame( + sample = ["s1", "s2", "s3", "s4", "s5", "s6", "s7", "s8"], + group = ["A", "A", "A", "A", "B", "B", "B", "B"], + run = ["R1", "R1", "R2", "R2", "R1", "R1", "R2", "R2"], + ph = [6.1, 6.3, 6.5, 6.7, 7.0, 7.2, 7.4, 7.6], + const_col = fill("same", 8), + ) + + p1 = parse_ordination_formula("~ group + ph + Condition(run)", meta) + @test p1.constrained_terms == ["group", "ph"] + @test p1.condition_terms == ["run"] + @test p1.df_constrained == 2 + @test p1.df_conditioned == 1 + + p_int = parse_ordination_formula("~ group * ph", meta) + @test p_int.constrained_terms == ["group", "ph", "group:ph"] + @test p_int.df_constrained == 3 + + # Auto-formula deduplicates identical categorical partitions (e.g. group == run) + meta_collinear = DataFrame( + sample = ["s1", "s2", "s3", "s4", "s5", "s6"], + group = ["Control", "Control", "Control", "Treatment", "Treatment", "Treatment"], + run = ["run1", "run1", "run1", "run2", "run2", "run2"], + ) + p_auto = parse_ordination_formula(nothing, meta_collinear) + @test p_auto.constrained_terms == ["group"] + @test p_auto.df_constrained == 1 + + # Reject injection, LHS formulas, Condition-only, overlap, constant covariates, saturation + @test_throws ArgumentError validate_formula_syntax("~ system(ls)") + @test_throws ArgumentError validate_formula_syntax("y ~ group") + @test_throws ArgumentError validate_formula_syntax("~ Condition(run)") + @test_throws ArgumentError validate_formula_syntax("~ group + Condition(group)") + @test_throws ArgumentError parse_ordination_formula("~ const_col", meta) + @test_throws ArgumentError parse_ordination_formula("~ missing_cov", meta) + meta_tiny = DataFrame(sample=["a", "b", "c"], g=["x", "y", "z"]) + @test_throws ArgumentError parse_ordination_formula("~ g", meta_tiny) + end + + @testset "Pure-Julia RDA, Ezekiel adjusted-R2, and two-way varpart vs Python reference fixture" begin + ref_doc = ijson_read(read(_ORD_REF_PATH)) + for ds in ref_doc[:datasets] + Y = _json_matrix(ds[:y]) + X = _json_matrix(ds[:x]) + Z = isnothing(ds[:z]) ? nothing : _json_matrix(ds[:z]) + res = pure_linear_rda(Y, X; Z=Z) + exp_rda = ds[:expected_rda] + @test isapprox(res.total_inertia, Float64(exp_rda[:total_inertia]); atol=1e-9) + @test isapprox(res.conditioned_inertia, Float64(exp_rda[:conditioned_inertia]); atol=1e-9) + @test isapprox(res.constrained_inertia, Float64(exp_rda[:constrained_inertia]); atol=1e-9) + @test isapprox(res.unconstrained_inertia, Float64(exp_rda[:unconstrained_inertia]); atol=1e-9) + @test isapprox(res.r2, Float64(exp_rda[:r2]); atol=1e-9) + @test isapprox(res.adj_r2, Float64(exp_rda[:adj_r2]); atol=1e-9) + @test isapprox(res.f_statistic, Float64(exp_rda[:f_statistic]); atol=1e-9) + @test length(res.constrained_eigenvalues) == length(exp_rda[:constrained_eigenvalues]) + for (a, b) in zip(res.constrained_eigenvalues, exp_rda[:constrained_eigenvalues]) + @test isapprox(a, Float64(b); atol=1e-9) + end + for (a, b) in zip(res.unconstrained_eigenvalues, exp_rda[:unconstrained_eigenvalues]) + @test isapprox(a, Float64(b); atol=1e-9) + end + + if !isnothing(ds[:expected_varpart]) + X1 = _json_matrix(ds[:x1]) + X2 = _json_matrix(ds[:x2]) + r1 = pure_linear_rda(Y, X1).adj_r2 + r2 = pure_linear_rda(Y, X2).adj_r2 + r12 = pure_linear_rda(Y, hcat(X1, X2)).adj_r2 + vp = pure_varpart_two_way(r1, r2, r12) + exp_vp = ds[:expected_varpart] + @test isapprox(vp.unique_x1, Float64(exp_vp[:unique_x1]); atol=1e-9) + @test isapprox(vp.unique_x2, Float64(exp_vp[:unique_x2]); atol=1e-9) + @test isapprox(vp.shared, Float64(exp_vp[:shared]); atol=1e-9) + @test isapprox(vp.residual, Float64(exp_vp[:residual]); atol=1e-9) + end + end + end + + @testset "Pipeline validation, config cascade, and HTTP route pre-R validation" begin + _ord_with_study_root() do root + # Contract endpoint returns canonical JCS + res_contract = _ord_http("GET", "/api/v1/ordination/contract") + @test res_contract.status == 200 + parsed_contract = ijson_read(res_contract.body) + @test parsed_contract[:schema_version] == 1 + + # Config cascade rejects rda + inherited bray distance, accepts atomic rda + aitchison + bad_patch = _ord_http("PATCH", "/api/v1/studies/study1/config", + Dict("analysis.ordination.method" => "rda")) + @test bad_patch.status == 400 + @test occursin("incompatible", String(bad_patch.body)) + + good_patch = _ord_http("PATCH", "/api/v1/studies/study1/config", + Dict("analysis.ordination.method" => "rda", + "analysis.ordination.distance" => "aitchison", + "analysis.ordination.scaling" => 1, + "analysis.ordination.permutations" => 499)) + @test good_patch.status == 200 + + # Pre-R HTTP 400 checks on /api/v1/studies/{study}/analysis/ordination + runs_payload = [ + Dict("run" => "run1", "group" => "Control"), + Dict("run" => "run2", "group" => "Treatment"), + ] + r404 = _ord_http("POST", "/api/v1/studies/unknown_study/analysis/ordination", + Dict("method" => "dbrda", "runs" => runs_payload)) + @test r404.status == 404 + + r_bad_compat = _ord_http("POST", "/api/v1/studies/study1/analysis/ordination", + Dict("method" => "rda", "distance" => "bray", "runs" => runs_payload)) + @test r_bad_compat.status == 400 + + r_bad_formula = _ord_http("POST", "/api/v1/studies/study1/analysis/ordination", + Dict("method" => "dbrda", "distance" => "bray", "formula" => "~ system(id)", "runs" => runs_payload)) + @test r_bad_formula.status == 400 + + r_missing_tree = _ord_http("POST", "/api/v1/studies/study1/analysis/ordination", + Dict("method" => "dbrda", "distance" => "weighted_unifrac", "runs" => runs_payload)) + @test r_missing_tree.status == 400 + @test occursin("missing_phylogeny", String(r_missing_tree.body)) + end + end + + if r_available() && vegan_available() + @testset "vegan accuracy within 1e-6 (dune + 3 microbiome datasets, all 4 methods, Condition, varpart, UniFrac)" begin + # 1. Cross-check all 4 reference datasets against live R vegan::rda & vegan::varpart within 1e-6 + ref_doc = ijson_read(read(_ORD_REF_PATH)) + for ds in ref_doc[:datasets] + Y = _json_matrix(ds[:y]) + X = _json_matrix(ds[:x]) + n = size(Y, 1) + meta_df = DataFrame(sample = ["s$i" for i in 1:n]) + x_names = String[] + for j in 1:size(X, 2) + nm = "x$j" + meta_df[!, nm] = X[:, j] + push!(x_names, nm) + end + rhs = join(x_names, " + ") + if !isnothing(ds[:z]) + Z = _json_matrix(ds[:z]) + z_names = String[] + for j in 1:size(Z, 2) + nm = "z$j" + meta_df[!, nm] = Z[:, j] + push!(z_names, nm) + end + rhs *= " + Condition(" * join(z_names, " + ") * ")" + end + cfg = OrdinationConfig(; + ordination_method = :rda, + ordination_distance = :euclidean, + ordination_formula = "~ " * rhs, + ordination_scaling = 2, + permutations = 199, + variance_partitioning = !isnothing(ds[:expected_varpart]) && size(X, 2) == 2, + seed = 12345, + ) + out = run_constrained_ordination(Y, meta_df, cfg; beta_transform = "none") + exp_rda = ds[:expected_rda] + @test isapprox(out.total_inertia, Float64(exp_rda[:total_inertia]); atol=1e-6) + @test isapprox(out.conditioned_inertia, Float64(exp_rda[:conditioned_inertia]); atol=1e-6) + @test isapprox(out.constrained_inertia, Float64(exp_rda[:constrained_inertia]); atol=1e-6) + @test isapprox(out.unconstrained_inertia, Float64(exp_rda[:unconstrained_inertia]); atol=1e-6) + @test isapprox(out.r2, Float64(exp_rda[:r2]); atol=1e-6) + @test isapprox(out.adj_r2, Float64(exp_rda[:adj_r2]); atol=1e-6) + @test isapprox(out.omnibus_test.f_statistic, Float64(exp_rda[:f_statistic]); atol=1e-6) + for (e_got, e_exp) in zip(out.constrained_eigenvalues, exp_rda[:constrained_eigenvalues]) + @test isapprox(e_got.eigenvalue, Float64(e_exp); atol=1e-6) + end + if cfg.variance_partitioning + exp_vp = ds[:expected_varpart] + @test out.variance_partitioning.status == "computed" + ind = out.variance_partitioning.individual_fractions + @test isapprox(ind[1].adj_r2, Float64(exp_vp[:unique_x1]); atol=1e-6) + @test isapprox(ind[2].adj_r2, Float64(exp_vp[:shared]); atol=1e-6) + @test isapprox(ind[3].adj_r2, Float64(exp_vp[:unique_x2]); atol=1e-6) + @test isapprox(ind[4].adj_r2, Float64(exp_vp[:residual]); atol=1e-6) + end + end + + # 2. Direct vegan::dune comparison for RDA, CCA, CAP, dbRDA including seeded permutation p-values & scores + dune_mat, dune_a1, dune_mgmt, ref_vals = with_r_lock() do + RCall.reval(raw""" + data(dune, package = "vegan") + data(dune.env, package = "vegan") + d_mat <- as.matrix(dune) + d_a1 <- as.numeric(dune.env$A1) + d_mgmt <- as.character(dune.env$Management) + d_df <- data.frame(A1 = d_a1, Management = factor(d_mgmt, levels = unique(d_mgmt))) + RNGkind("Mersenne-Twister", "Inversion", "Rejection") + set.seed(777) + fit_dbrda <- vegan::dbrda(d_mat ~ A1 + Management, data = d_df, distance = "bray") + r2_dbrda <- vegan::RsquareAdj(fit_dbrda) + set.seed(777) + aov_dbrda <- vegan::anova.cca(fit_dbrda, permutations = permute::how(nperm = 199L), parallel = 1L) + sc_dbrda <- as.matrix(vegan::scores(fit_dbrda, choices = c(1, 2), display = "sites", scaling = 2L)) + """) + out_tuple = ( + Matrix{Float64}(RCall.rcopy(RCall.reval("d_mat"))), + Vector{Float64}(RCall.rcopy(RCall.reval("d_a1"))), + Vector{String}(RCall.rcopy(RCall.reval("d_mgmt"))), + (; + tot = Float64(RCall.rcopy(RCall.reval("as.numeric(fit_dbrda[['tot.chi']])"))), + con = Float64(RCall.rcopy(RCall.reval("as.numeric(fit_dbrda[['CCA']][['tot.chi']])"))), + r2 = Float64(RCall.rcopy(RCall.reval("as.numeric(r2_dbrda[['r.squared']])"))), + adj = Float64(RCall.rcopy(RCall.reval("as.numeric(r2_dbrda[['adj.r.squared']])"))), + f = Float64(RCall.rcopy(RCall.reval("as.numeric(aov_dbrda[['F']][1])"))), + p = Float64(RCall.rcopy(RCall.reval("as.numeric(aov_dbrda[['Pr(>F)']][1])"))), + sc = Matrix{Float64}(RCall.rcopy(RCall.reval("sc_dbrda"))), + ), + ) + RCall.reval("rm(list = intersect(ls(), c('d_mat', 'd_a1', 'd_mgmt', 'd_df', 'fit_dbrda', 'r2_dbrda', 'aov_dbrda', 'sc_dbrda')))") + out_tuple + end + dune_meta = DataFrame( + sample = ["site$i" for i in 1:size(dune_mat, 1)], + A1 = dune_a1, + Management = dune_mgmt, + ) + res_dbrda = run_constrained_ordination( + dune_mat, + dune_meta, + OrdinationConfig(; + ordination_method = :dbrda, + ordination_distance = :bray, + ordination_formula = "~ A1 + Management", + ordination_scaling = 2, + permutations = 199, + variance_partitioning = true, + seed = 777, + ); + beta_transform = "none", + test_axes = true, + ) + @test isapprox(res_dbrda.total_inertia, ref_vals.tot; atol=1e-6) + @test isapprox(res_dbrda.constrained_inertia, ref_vals.con; atol=1e-6) + @test isapprox(res_dbrda.r2, ref_vals.r2; atol=1e-6) + @test isapprox(res_dbrda.adj_r2, ref_vals.adj; atol=1e-6) + @test isapprox(res_dbrda.omnibus_test.f_statistic, ref_vals.f; atol=1e-6) + @test isapprox(res_dbrda.omnibus_test.p_value, ref_vals.p; atol=1e-6) + for i in 1:size(dune_mat, 1) + @test isapprox(res_dbrda.sites[i].axis1, ref_vals.sc[i, 1]; atol=1e-6) + @test isapprox(res_dbrda.sites[i].axis2, ref_vals.sc[i, 2]; atol=1e-6) + end + @test !isempty(res_dbrda.axis_tests) + @test res_dbrda.variance_partitioning.status == "computed" + + # Verify CCA, CAP, and scaling=1 on dune as well + res_cca = run_constrained_ordination( + dune_mat, + dune_meta, + OrdinationConfig(; + ordination_method = :cca, + ordination_distance = :chisq, + ordination_formula = "~ A1 + Condition(Management)", + ordination_scaling = 1, + permutations = 199, + seed = 777, + ); + beta_transform = "none", + ) + @test res_cca.method == "cca" + @test res_cca.scaling == 1 + @test res_cca.conditioned_inertia > 0.0 + @test res_cca.provenance["scaling_meaning"] == "distance" + + res_cap = run_constrained_ordination( + dune_mat, + dune_meta, + OrdinationConfig(; + ordination_method = :cap, + ordination_distance = :jaccard, + ordination_formula = "~ A1 + Management", + ordination_scaling = 2, + permutations = 199, + seed = 777, + ); + beta_transform = "none", + ) + @test res_cap.method == "cap" + @test res_cap.distance == "jaccard" + @test !isnothing(res_cap.omnibus_test.p_value) + + # 3. End-to-end HTTP POST with persistence, Aitchison CLR, and inline Newick UniFrac + _ord_with_study_root() do _ + runs_payload = [ + Dict("run" => "run1", "group" => "Control"), + Dict("run" => "run2", "group" => "Treatment"), + ] + resp_clr = _ord_http("POST", "/api/v1/studies/study1/analysis/ordination", Dict( + "method" => "rda", + "distance" => "aitchison", + "scaling" => 2, + "permutations" => 199, + "formula" => "~ group", + "persist" => true, + "show_features" => true, + "runs" => runs_payload, + )) + @test resp_clr.status == 200 + body_clr = ijson_read(resp_clr.body) + @test body_clr[:method] == "rda" + @test body_clr[:distance] == "aitchison" + @test body_clr[:n_samples] == 6 + @test haskey(body_clr, :biplot) + @test haskey(body_clr, :artifact) + aid = String(body_clr[:artifact][:analysis_id]) + resp_loaded = _ord_http("GET", "/api/v1/studies/study1/analysis/ordination/$aid") + @test resp_loaded.status == 200 + + # Weighted UniFrac via inline ASV Newick tree + t1 = CladeCumulus.asv_tip_name("ACGTACGTACGTACGT") + t2 = CladeCumulus.asv_tip_name("TGCAACGTACGTACGT") + t3 = CladeCumulus.asv_tip_name("GGCCACGTACGTACGT") + t4 = CladeCumulus.asv_tip_name("TTAAACGTACGTACGT") + nwk = "(($t1:0.15,$t2:0.15):0.25,($t3:0.12,$t4:0.12):0.28);" + resp_uf = _ord_http("POST", "/api/v1/studies/study1/analysis/ordination", Dict( + "method" => "dbrda", + "distance" => "weighted_unifrac", + "permutations" => 199, + "formula" => "~ group", + "tree_newick" => nwk, + "runs" => runs_payload, + )) + @test resp_uf.status == 200 + body_uf = ijson_read(resp_uf.body) + @test body_uf[:distance] == "weighted_unifrac" + @test body_uf[:provenance][:phylogeny_source] == "inline_newick" + end + end + else + @info "Skipping live R/vegan ordination execution tests (R or vegan unavailable in this environment)" + end +end