From 5a92c065698bd8fa8f5fe8e587241354df98e379 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 15:24:25 +0000 Subject: [PATCH] Root the indexed tree before computing Faith's PD MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Codex (P1) and CodeRabbit (major) both flagged this on #234: IQ-TREE infers an *unrooted* tree for a reversible model, so the root of a `reference.treefile` is only whichever node happened to be serialized first. `analyse_phylogeny` used that root directly, while the legacy `faith_pd(counts, taxa, tree)` path midpoint-roots first via `_ensure_rooted_tree`. UniFrac is unaffected — it sums edge partitions and re-rooting preserves them — but Faith's PD is a path-to-root quantity, so for a sample with a single active tip the indexed path reported 1.0 where the legacy path reported 1.5 for the same tree, and an equivalent rerooted serialization could give a different number again. `ensure_rooted` now applies the legacy convention to `PhylogeneticTree`: a root with more than two children is midpoint-rooted by `midpoint_root` (the implementation already verified against scikit-bio reference outputs) and an already bifurcating tree is returned untouched. To reuse that audited code without reimplementing it, `_to_phylo`/`_from_phylo` convert to `PhyloNode` and back; the round trip preserves labels, tip names, lengths and the parents-precede-children ordering the roll-up relies on. `analyse_phylogeny`, `faith_pd(tree, counts, tip_names)` and `unifrac(tree, counts, tip_names)` root before rolling up; the roll-up overloads take a matrix, so they document that it must come from an already rooted tree. The cladistics `results.json` now records the convention in a `rooting` field, mirrored as an optional `CladisticsResult.rooting` in the frontend types. Tests pin the convention: an unrooted trifurcating tree roots to two children, PD for a lone tip matches the legacy `faith_pd` path exactly (1.5, versus the 1.0 the raw serialization gives), and a rooted and an unrooted serialization of the same inferred tree agree on `faith_pd`, both UniFrac matrices and `total_branch_length`. Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- docs/cladistics/clade-cumulus.md | 2 +- frontend/src/api/types.ts | 2 + src/analysis/clade_cumulus.jl | 89 +++++++++++++++++++++++++++++--- src/server/routes/cladistics.jl | 4 ++ test/unit/test_clade_cumulus.jl | 33 ++++++++++++ 5 files changed, 123 insertions(+), 7 deletions(-) diff --git a/docs/cladistics/clade-cumulus.md b/docs/cladistics/clade-cumulus.md index 7a386862..373e84a2 100644 --- a/docs/cladistics/clade-cumulus.md +++ b/docs/cladistics/clade-cumulus.md @@ -29,7 +29,7 @@ formal Agda proof suite, and React UI. - Supports `:phylogenetic`, `:taxonomic`, and `:combined` modes, surfacing non-monophyletic / dispersed taxa and discordant clades. - Strict raw `.jplace` evidence reader (`parse_jplace`) without defaulting missing `like_weight_ratio` values to `1.0`. - Exact `faith_pd`, `clade_faith_pd`, `unweighted_unifrac`, and `weighted_unifrac` (`normalized = false | true`), verified against `scikit-bio` reference outputs. - - Indexed study-tree API used by the de novo cladistics builds (`src/server/routes/cladistics.jl`): `asv_tip_name` (content-addressed tip names), `parse_phylogenetic_tree` (Newick → `PhylogeneticTree` with a `tips` index), `clade_counts` (nodes × samples roll-up), `unifrac` (either a raw samples × ASVs matrix or an existing roll-up), and `analyse_phylogeny`, the single call that produces the JSON payload the UI renders. + - Indexed study-tree API used by the de novo cladistics builds (`src/server/routes/cladistics.jl`): `asv_tip_name` (content-addressed tip names), `parse_phylogenetic_tree` (Newick → `PhylogeneticTree` with a `tips` index), `clade_counts` (nodes × samples roll-up), `unifrac` (either a raw samples × ASVs matrix or an existing roll-up), and `analyse_phylogeny`, the single call that produces the JSON payload the UI renders. Every entry point first passes the tree through `ensure_rooted`, which midpoint-roots an unrooted tree exactly as the legacy `PhyloNode` metrics do: IQ-TREE infers unrooted trees for reversible models, so the serialized root is arbitrary, and while UniFrac only sums edge partitions (and is root-invariant), Faith's PD is a path-to-root quantity. The stored `results.json` records the convention in its `rooting` field. 5. **Non-Destructive Server Routes & Evidence Mode UI (`src/server/routes/clade_cumulus.jl`, `frontend/src/components/EvidenceModeToggle.tsx`, `frontend/src/components/CladeCumulusPanel.tsx`)** - Gated behind an explicit `CladeCumulus Evidence Mode` switch on the Study `Trees` tab so standard workflows are never altered silently. diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 360d1121..158d8382 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -746,6 +746,8 @@ export interface CladisticsResult { unweighted_unifrac: number[][] weighted_unifrac: number[][] total_branch_length: number + /** How the inferred (unrooted) tree was rooted before computing PD. */ + rooting?: string normalisation: { method: CountNormalisation; depth: number; seed: number } excluded_samples: string[] weighted_unifrac_definition: string diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 34b80b43..e0d56652 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -18,7 +18,7 @@ export PhyloNode, IndexedPhyloNode, PhylogeneticTree, QueryCandidate, QueryPlacement, JplaceDocument, CladeNode, CladeTreeResult, parse_newick, parse_phylogenetic_tree, midpoint_root, parse_jplace, - asv_tip_name, clade_counts, unifrac, analyse_phylogeny, + asv_tip_name, clade_counts, unifrac, analyse_phylogeny, ensure_rooted, faith_pd, faith_pd_table, clade_faith_pd, unweighted_unifrac, weighted_unifrac, unifrac_matrix, branch_unifrac_contributions, @@ -278,6 +278,63 @@ function _tree_lengths(tree::PhylogeneticTree)::Vector{Float64} lengths end +## Rooting + +""" + ensure_rooted(tree::PhylogeneticTree) -> PhylogeneticTree + +Return `tree` with a biologically meaningful root, using the same convention as +the legacy `PhyloNode` metrics (`_ensure_rooted_tree`): an unrooted tree — more +than two children at the root — is midpoint-rooted by `midpoint_root`, and an +already bifurcating tree is returned unchanged. + +This matters because IQ-TREE infers an *unrooted* tree for a reversible +substitution model, so the root of a `reference.treefile` is only whichever +node happened to be serialized first. UniFrac is unaffected (it sums edge +partitions, and re-rooting preserves them), but Faith's PD is a path-to-root +quantity: for a sample with a single active tip the arbitrary root inflates PD +by the distance from that tip to the serialization point. Rooting here keeps +`analyse_phylogeny` consistent with `faith_pd(counts, taxa, tree)`. +""" +function ensure_rooted(tree::PhylogeneticTree)::PhylogeneticTree + length(tree.nodes[tree.root].children) <= 2 && return tree + rooted, changed = midpoint_root(_to_phylo(tree)) + changed || return tree # fewer than two tips: nothing to root against + _from_phylo(rooted) +end + +# The legacy rooting and metric code is audited against scikit-bio, so it is +# reused rather than reimplemented: convert to `PhyloNode`, root, convert back. +function _to_phylo(tree::PhylogeneticTree, id::Int = tree.root)::PhyloNode + node = tree.nodes[id] + is_root = id == tree.root + !is_root && isnothing(node.length) && throw(ArgumentError( + "Faith PD and UniFrac need a length on every non-root branch")) + children = PhyloNode[_to_phylo(tree, child) for child in node.children] + PhyloNode("n$id", node.name, is_root ? 0.0 : something(node.length, 0.0), + nothing, nothing, children) +end + +function _from_phylo(root::PhyloNode)::PhylogeneticTree + nodes = IndexedPhyloNode[] + tips = Dict{String,Int}() + function add!(node::PhyloNode, parent::Int)::Int + push!(nodes, IndexedPhyloNode(node.name, + parent == 0 ? nothing : node.branch_length, parent, Int[])) + id = length(nodes) + for child in node.children + push!(nodes[id].children, add!(child, id)) + end + if isempty(node.children) && !isempty(node.name) + haskey(tips, node.name) && throw(ArgumentError("duplicate tip '$(node.name)'")) + tips[node.name] = id + end + id + end + root_id = add!(root, 0) + PhylogeneticTree(nodes, root_id, tips) +end + """ unifrac(tree, cumulative; weighted=false) -> Matrix{Float64} @@ -333,12 +390,15 @@ end UniFrac over the natural samples × ASVs count matrix. The counts are rolled up through `clade_counts` first, so callers do not have to construct the internal nodes × samples matrix; `faith_pd(tree, counts, tip_names)` is the matching -raw-count entry point. +raw-count entry point. The tree is passed through `ensure_rooted` first, so an +unrooted IQ-TREE tree and a rerooted serialization of it agree; UniFrac itself +is root-invariant because it only sums edge partitions. """ function unifrac(tree::PhylogeneticTree, counts::AbstractMatrix{<:Real}, tip_names::AbstractVector{<:AbstractString}; weighted::Bool=false)::Matrix{Float64} - unifrac(tree, clade_counts(tree, counts, tip_names); weighted = weighted) + rooted = ensure_rooted(tree) + unifrac(rooted, clade_counts(rooted, counts, tip_names); weighted = weighted) end function _faith_pd_indexed(tree::PhylogeneticTree, cumulative::AbstractMatrix{<:Real})::Vector{Float64} @@ -362,16 +422,32 @@ function _faith_pd_indexed(tree::PhylogeneticTree, cumulative::AbstractMatrix{<: out end +""" + faith_pd(tree, cumulative) -> Vector{Float64} + faith_pd(tree, counts, tip_names) -> Vector{Float64} + +Faith's Phylogenetic Diversity per sample column. The roll-up overload sums the +length of every non-root edge whose clade holds reads, so `cumulative` must come +from `clade_counts` applied to a tree already passed through `ensure_rooted`; +the raw-count overload does that rooting for you and is the one to prefer. +""" faith_pd(tree::PhylogeneticTree, cumulative::AbstractMatrix{<:Real}) = _faith_pd_indexed(tree, cumulative) -faith_pd(tree::PhylogeneticTree, counts::AbstractMatrix{<:Real}, - tip_names::AbstractVector{<:AbstractString}) = - faith_pd(tree, clade_counts(tree, counts, tip_names)) +function faith_pd(tree::PhylogeneticTree, counts::AbstractMatrix{<:Real}, + tip_names::AbstractVector{<:AbstractString})::Vector{Float64} + rooted = ensure_rooted(tree) + faith_pd(rooted, clade_counts(rooted, counts, tip_names)) +end """ analyse_phylogeny(tree, counts, tip_names, sample_names) -> NamedTuple Compute the complete JSON-ready CladeCumulus result: cumulative clade counts and frequencies, Faith PD, and both UniFrac matrices. + +The tree is passed through `ensure_rooted` first, so an unrooted IQ-TREE +`reference.treefile` and any rerooted serialization of the same inferred tree +produce the same `faith_pd`, `unweighted_unifrac` and `weighted_unifrac` values. +The returned `nodes` describe the rooted tree. """ function analyse_phylogeny(tree::PhylogeneticTree, counts::AbstractMatrix{<:Real}, tip_names::AbstractVector{<:AbstractString}, @@ -380,6 +456,7 @@ function analyse_phylogeny(tree::PhylogeneticTree, counts::AbstractMatrix{<:Real "counts has $(size(counts, 1)) rows but $(length(sample_names)) sample names were supplied")) length(unique(String.(sample_names))) == length(sample_names) || throw(ArgumentError("sample names must be unique")) + tree = ensure_rooted(tree) cumulative = clade_counts(tree, counts, tip_names) samples = String.(sample_names) totals = vec(cumulative[tree.root, :]) diff --git a/src/server/routes/cladistics.jl b/src/server/routes/cladistics.jl index 14979026..5d13ee44 100644 --- a/src/server/routes/cladistics.jl +++ b/src/server/routes/cladistics.jl @@ -224,6 +224,10 @@ function _cc_compute(study::String, id::String) result["excluded_samples"] = counts_doc.excluded_samples result["weighted_unifrac_definition"] = "sum(branch * abs(a-b)) / sum(branch * (a+b)); clade relative abundances" + # IQ-TREE infers an unrooted tree for a reversible model, so PD is only + # meaningful relative to a stated rooting convention; record it next to the + # other metric definitions rather than leaving it implicit. + result["rooting"] = "midpoint (unrooted trees; Faith PD is root-dependent, UniFrac is not)" result end diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index 10f098ba..567c6897 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -135,6 +135,39 @@ _cc_req(method::String, path::String, body::String = "") = @test decoded.nodes[1].parent_id === nothing @test decoded.faith_pd[2] ≈ 5.0 @test decoded.unweighted_unifrac[1][2] ≈ 0.7 + + # IQ-TREE infers an unrooted tree for a reversible model, so the root of + # a `reference.treefile` is only whichever node was serialized first. + # `analyse_phylogeny` must therefore root the way the legacy metrics do + # (midpoint), or Faith's PD depends on that arbitrary choice. + unrooted_text = "(A:1.0,B:2.0,C:3.0);" + unrooted = CC.parse_phylogenetic_tree(unrooted_text) + @test length(unrooted.nodes[unrooted.root].children) == 3 + rooted = CC.ensure_rooted(unrooted) + @test length(rooted.nodes[rooted.root].children) == 2 + @test sort(collect(keys(rooted.tips))) == ["A", "B", "C"] + @test CC.ensure_rooted(rooted) === rooted + + tips3 = ["A", "B", "C"] + one_tip = [1.0 0.0 0.0] + counts3 = [1.0 0.0 0.0; + 0.0 2.0 3.0] + # The midpoint falls inside the root-C edge, so a lone A tip now has to + # reach it: 1.0 (A) + 0.5 (half of the B-root-C path remainder) = 1.5. + @test CC._faith_pd_indexed(unrooted, CC.clade_counts(unrooted, one_tip, tips3))[1] ≈ 1.0 + @test CC.faith_pd(unrooted, one_tip, tips3)[1] ≈ 1.5 + @test CC.faith_pd(unrooted, one_tip, tips3)[1] ≈ + CC.faith_pd([1.0, 0.0, 0.0], tips3, unrooted_text) + + # Rooted and unrooted serializations must agree on every payload metric: + # PD through the shared rooting convention, UniFrac through invariance. + unrooted_result = CC.analyse_phylogeny(unrooted, counts3, tips3, samples) + rooted_result = CC.analyse_phylogeny(rooted, counts3, tips3, samples) + @test unrooted_result.faith_pd ≈ rooted_result.faith_pd + @test unrooted_result.unweighted_unifrac ≈ rooted_result.unweighted_unifrac + @test unrooted_result.weighted_unifrac ≈ rooted_result.weighted_unifrac + @test unrooted_result.total_branch_length ≈ rooted_result.total_branch_length + @test unrooted_result.total_branch_length ≈ 6.0 end