Root the indexed tree before computing Faith's PD (Codex P1 on #234)

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.

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

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,
@@ -292,6 +292,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}
 
@@ -347,12 +404,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}
@@ -376,16 +436,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},
@@ -394,6 +470,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, :])

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
 

test/unit/test_clade_cumulus.jl

@@ -184,6 +184,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