From 6f7d5f5f37420c471ea92120ed6025c73b946b67 Mon Sep 17 00:00:00 2001 From: Arena Agent Date: Fri, 23 Oct 2026 07:30:37 +0000 Subject: [PATCH 1/7] fix(cladistics): load the DOI module, drop a duplicate route include, and repair the CladeCumulus readers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `main` still cannot load the package, and the CladeCumulus suite that never ran before now reports twelve defects. All of them are narrowed to: Package loading - src/core/doi.jl: #210 removed `module DOI` but left `end # module DOI`, a lowering error at include time ("end" at top level) — every test file fails before it starts. The opener is restored. - src/server/server.jl: routes/cladistics.jl was included twice, so every /cladistics route is registered twice. One include remains. CladeCumulus readers (test_clade_cumulus.jl, 3 failed / 14 errored) - routes/clade_cumulus.jl opened results.duckdb with `DuckDB.DB` and `DBInterface.connect` without importing either module, so the sources route answered 500 and every HTTP assertion downstream of it failed. Both are imported now, as the other route files already do. - `_graft_jplace_queries` and the diversity block in the route still used the pre-#207 field names (`weight`, `distal_length`, `pendant_length`) on `QueryCandidate` (whose fields are `like_weight`, `distal_len`, `pendant_len`), and `_graft_jplace_queries` used `.edge_num` on `PhyloNode` (whose field is `edge_number`). Five FieldErrors, now aligned with the structs. - parse_newick lost NHX support: the label reader consumed `[&&NHX:B=92]` before the call that records it, so `) [&&NHX:B=92] :0.1` parsed with support `nothing`. The reader now takes the capture refs. - parse_jplace rejected a single-candidate placement whose ratio is below 1.0, which is how placement tools report one site: the remainder is mass they did not place. A single candidate is now decisive when its ratio is at least the new `PLACEMENT_SINGLE_MIN_LWR = 0.5`; a multi-candidate placement must still partition the mass to within `PLACEMENT_WEIGHT_SUM_TOL`. A single candidate below half stays `undecided` with the same reason string. - test_clade_cumulus.jl reached the tool registry as `PH.Tools`, but `Phylogeny` imports `tool_bin` selectively and does not bind `Tools`; `MetaManifold.Tools` is the module that owns the registry. Guards against both merge-harm classes - scripts/check-julia-syntax.jl now rejects a bare `end` at the top level. Such a file parses cleanly and only fails at lowering, which is exactly how a missing `module` opener — the failure that stopped the whole suite — passed the parse gate. - test_routes.jl asserts that no method+path pair is registered twice and that both cladistics route files are present, which no request test could reveal (a duplicate registration silently shadows the first). Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- scripts/check-julia-syntax.jl | 13 +++++++++- src/analysis/clade_cumulus.jl | 40 ++++++++++++++++++++++-------- src/core/doi.jl | 2 ++ src/server/routes/clade_cumulus.jl | 6 ++--- src/server/server.jl | 1 - test/unit/test_clade_cumulus.jl | 10 +++++--- test/unit/test_routes.jl | 19 ++++++++++++++ 7 files changed, 71 insertions(+), 20 deletions(-) diff --git a/scripts/check-julia-syntax.jl b/scripts/check-julia-syntax.jl index bc8a9d60..f4a93eaf 100644 --- a/scripts/check-julia-syntax.jl +++ b/scripts/check-julia-syntax.jl @@ -20,6 +20,17 @@ for file in files issue = ex break end + # A bare `end` at the top level parses without error and then fails + # at lowering (`syntax: "end" at top level`), which is how a file + # that lost its `module Foo` opener reaches CI looking syntactically + # fine. An `end` balances exactly one opener, so top level cannot + # hold one: this is a missing or duplicated opener, and the file + # would not load. + if ex === :end + issue = "stray `end` at top level (a `module`/`function`/... opener " * + "is missing or was duplicated in this file)" + break + end next_position > position || error("Julia parser made no progress at code unit $position") position = next_position @@ -30,7 +41,7 @@ for file in files if !isnothing(issue) failed = true println(stderr, "Julia parse failed in $file:") - println(stderr, sprint(showerror, issue)) + println(stderr, issue isa AbstractString ? issue : sprint(showerror, issue)) end end failed && exit(1) diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index a702fde2..ca2850b4 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -30,6 +30,13 @@ const DEFAULT_TAXONOMIC_RANKS = [ const PLACEMENT_WEIGHT_SUM_TOL = 1e-5 +# Placement tools report a query placed at a single site as one candidate row +# whose `like_weight_ratio` is the share of the probability mass behind it; the +# rest is mass the tool did not place anywhere. Such a query is decisive while +# that share is at least a half. A multi-candidate placement instead partitions +# the mass, so its ratios must sum to one within `PLACEMENT_WEIGHT_SUM_TOL`. +const PLACEMENT_SINGLE_MIN_LWR = 0.5 + ## Phylogenetic tree structures & Newick / jplace parsers mutable struct PhyloNode id :: String @@ -202,8 +209,14 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]): end end - function read_label_and_quoted()::Tuple{String,Bool} - skip_ws_and_leading_comments!() + function read_label_and_quoted( + capture_nhx::Union{Ref{Union{Float64,Nothing}},Nothing} = nothing, + capture_edge::Union{Ref{Union{Int,Nothing}},Nothing} = nothing, + )::Tuple{String,Bool} + # A support tag may sit between the label and the branch length, or + # directly after ')' with no label at all (`)[&&NHX:B=92]:0.1`), so this + # skip has to record what it consumes. + skip_ws_and_leading_comments!(; capture_nhx = capture_nhx, capture_edge = capture_edge) pos[] > n && return ("", false) if bytes[pos[]] == UInt8('\'') pos[] += 1 @@ -260,7 +273,7 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]): nhx_sup = Ref{Union{Float64,Nothing}}(nothing) edge_num = Ref{Union{Int,Nothing}}(nothing) - raw_lbl, was_quoted = read_label_and_quoted() + raw_lbl, was_quoted = read_label_and_quoted(nhx_sup, edge_num) skip_ws_and_leading_comments!(capture_nhx = nhx_sup, capture_edge = edge_num) branch_len = 0.0 @@ -565,7 +578,12 @@ function parse_jplace(raw_json::AbstractString)::JplaceDocument if !isempty(bad_reason) push!(queries, QueryPlacement(q_ord, q_names, QueryCandidate[], :undecided, bad_reason, nothing)) - elseif abs(weight_sum - 1.0) > PLACEMENT_WEIGHT_SUM_TOL + elseif length(cands) == 1 && cands[1].like_weight < PLACEMENT_SINGLE_MIN_LWR + # One candidate but less than half the mass supports it: the rest of + # the mass is unexplained, so the placement is not decisive. + push!(queries, QueryPlacement(q_ord, q_names, cands, :undecided, + "weight_sum_out_of_tolerance ($(round(weight_sum; digits=6)))", nothing)) + elseif length(cands) > 1 && abs(weight_sum - 1.0) > PLACEMENT_WEIGHT_SUM_TOL push!(queries, QueryPlacement(q_ord, q_names, cands, :undecided, "weight_sum_out_of_tolerance ($(round(weight_sum; digits=6)))", nothing)) else @@ -991,7 +1009,7 @@ function _graft_jplace_queries( for q in jd.queries q.status === :valid && !isempty(q.candidates) || continue - total_w = sum((c.weight for c in q.candidates); init = 0.0) + total_w = sum((c.like_weight for c in q.candidates); init = 0.0) total_w > 0.0 || continue multi = length(q.candidates) > 1 for (qname, mult) in q.names @@ -1000,7 +1018,7 @@ function _graft_jplace_queries( tx = get(asv_deepest_taxon, qname, nothing) tconf = get(asv_tax_conf, qname, 1.0) for (ci, cand) in enumerate(q.candidates) - share = cand.weight / total_w + share = cand.like_weight / total_w tip_id = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname if multi asv_reads[tip_id] = get(asv_reads, tip_id, 0.0) + base_cnt * share @@ -1009,9 +1027,9 @@ function _graft_jplace_queries( elseif !haskey(asv_reads, tip_id) asv_reads[tip_id] = base_cnt end - query_lwr[tip_id] = cand.weight + query_lwr[tip_id] = cand.like_weight push!(get!(by_edge, cand.edge_num, _GraftedQuerySpec[]), - _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_length), max(0.0, cand.pendant_length), cand.weight)) + _GraftedQuerySpec(tip_id, qname, share, max(0.0, cand.distal_len), max(0.0, cand.pendant_len), cand.like_weight)) end end end @@ -1021,8 +1039,8 @@ function _graft_jplace_queries( function walk(u::PhyloNode)::PhyloNode rebuilt_kids = PhyloNode[walk(ch) for ch in u.children] - cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_num, rebuilt_kids) - specs = u.edge_num === nothing ? nothing : get(by_edge, u.edge_num, nothing) + cur = PhyloNode(u.id, u.name, u.branch_length, u.support, u.edge_number, rebuilt_kids) + specs = u.edge_number === nothing ? nothing : get(by_edge, u.edge_number, nothing) (specs === nothing || isempty(specs)) && return cur # Order grafts from most distal (closest to `cur`, smallest `distal`) to most proximal. @@ -1040,7 +1058,7 @@ function _graft_jplace_queries( sub = PhyloNode(string("jplace_clade#", seq_auto[]), "", 0.0, 1.0, nothing, PhyloNode[sub, qtip]) end rem_len = max(0.0, total_len - prev_distal) - PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_num, sub.children) + PhyloNode(sub.id, sub.name, rem_len, cur.support, cur.edge_number, sub.children) end walk(base_root) diff --git a/src/core/doi.jl b/src/core/doi.jl index 6a2f926a..f5b0bd95 100644 --- a/src/core/doi.jl +++ b/src/core/doi.jl @@ -16,6 +16,8 @@ # * Durable write-ahead journal; immutable DOI/provenance receipts # * DOI insertion never rewrites scientific objects +module DOI + using Dates, SHA, HTTP, JSON3, Logging export DOIConfig, Deposition, PublicationReceipt, diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index 45825b61..ec009353 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/src/server/routes/clade_cumulus.jl @@ -10,7 +10,7 @@ # Canonical run DuckDB tables are read-only here; staged reassignment proposals # write only study-local sidecars under `projects/{study}/trees/.clade_cumulus/`. -using DataFrames, OrderedCollections, SHA, Dates, JSON3 +using DataFrames, OrderedCollections, SHA, Dates, JSON3, DuckDB, DBInterface const _CC_SIDECAR_MAX_BYTES = 5 * 1024 * 1024 @@ -268,13 +268,13 @@ end jd_div = MetaManifold.CladeCumulus.parse_jplace(jplace_str) for q in jd_div.queries q.status === :valid && !isempty(q.candidates) || continue - tw = sum((c.weight for c in q.candidates); init = 0.0) + tw = sum((c.like_weight for c in q.candidates); init = 0.0) tw > 0.0 || continue multi = length(q.candidates) > 1 for (qname, _) in q.names for (ci, cand) in enumerate(q.candidates) tid = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname - tip_weights[tid] = cand.weight / tw + tip_weights[tid] = cand.like_weight / tw end end end diff --git a/src/server/server.jl b/src/server/server.jl index ebf5abb6..eca368a0 100644 --- a/src/server/server.jl +++ b/src/server/server.jl @@ -118,7 +118,6 @@ 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", "cladistics.jl")) ## Bad request middleware # Helpers deep in a route (e.g. `_req_group`) throw this to reject a request. diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index da09bc18..c05c8071 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -252,8 +252,10 @@ done echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_000004:0.08);" > "$OUT" """) chmod(mock_ft, 0o755) - old_ft = get(PH.Tools._tools, "fasttree", nothing) - PH.Tools._tools["fasttree"] = mock_ft + # `Phylogeny` imports `tool_bin` from `Tools` but does not bind the + # module itself, so the registry is reached through `MetaManifold`. + old_ft = get(MetaManifold.Tools._tools, "fasttree", nothing) + MetaManifold.Tools._tools["fasttree"] = mock_ft try s4 = [ "ASV_1" => "ACGTACGTACGTACGTACGT", @@ -267,9 +269,9 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 r["sha256"] == bytes2hex(open(sha256, mock_ft)), ft_res.tool_records) finally if old_ft === nothing - delete!(PH.Tools._tools, "fasttree") + delete!(MetaManifold.Tools._tools, "fasttree") else - PH.Tools._tools["fasttree"] = old_ft + MetaManifold.Tools._tools["fasttree"] = old_ft end end end diff --git a/test/unit/test_routes.jl b/test/unit/test_routes.jl index 62998944..2cb871ef 100644 --- a/test/unit/test_routes.jl +++ b/test/unit/test_routes.jl @@ -932,5 +932,24 @@ _saved_root = SV.ServerState._root[] end +@testset "Registered routes" begin + # Every route file is `include`d once into `Server`, and each macro call + # appends to `Server.ROUTES`. A file included twice would register its + # routes twice, which no request would reveal: the second registration just + # shadows the first. Assert the table has no duplicate method+path pair. + pairs = [(m, p) for (m, p, _) in SV.ROUTES] + @test !isempty(pairs) + @test allunique(pairs) + + dupes = unique([k for k in pairs if count(==(k), pairs) > 1]) + @test isempty(dupes) + + # The Cladistics tab and the run-level CladeCumulus evidence view are + # served by two different route files; both must be registered. + paths = Set(string(p) for (_, p) in pairs) + @test any(occursin("/cladistics", p) for p in paths) + @test any(occursin("/clade-cumulus", p) for p in paths) +end + # Testsets above point the server at temporary roots. SV.ServerState._root[] = _saved_root From d9a444c92415b5747b3eb8f796c09bbdd409f939 Mon Sep 17 00:00:00 2001 From: hyperpolymath Date: Fri, 23 Oct 2026 07:58:04 +0000 Subject: [PATCH 2/7] fix(cladistics): normalise treeFile at the request boundary CI on `main` (8df1041) shows `POST /clade-cumulus` answering 500 for every `study_tree` request: MethodError: no method matching _valid_tree_file(::SubString{String}) `JSON3` hands request fields back as `SubString`s and `strip` keeps them that way, but `_valid_tree_file` and `_tree_format` in `routes/trees.jl` are `String`-typed. The handler has no catch around `_resolve_tree_inputs_from_request`, so the MethodError surfaced as a 500 and took seven assertions down with it: `r_comb.status == 200` failed, and every `comb.*` access after it raised `KeyError` until the testset aborted. #222 fixed this class at `_load_run_table_df` -- "JSON3 hands request fields back as SubStrings; the DuckDB helpers are String-typed, so normalise once at this boundary" -- but `_resolve_tree_inputs_from_request` two functions below kept the unnormalised values. `cladistics.jl` already wraps its own call as `_valid_tree_file(String(published))`, so the convention is to normalise at the call site; this follows it. The existing `r_comb` assertion is the regression test: it sends `"treeFile" => "amplicon.nwk"` over JSON, which is exactly the SubString path that broke. With the `u.edge_num` -> `u.edge_number` repair already carried in from #223, this accounts for 8 of the 9 errors that run reported. --- src/server/routes/clade_cumulus.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index 03ced23f..8ca8d2f5 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/src/server/routes/clade_cumulus.jl @@ -96,9 +96,16 @@ function _resolve_tree_inputs_from_request( body, source_digests::OrderedDict{String,String}, ) + # JSON3 hands request fields back as SubStrings and `strip` keeps them that + # way, but `_valid_tree_file` and `_tree_format` are String-typed -- an + # unnormalised `treeFile` reached them as a SubString and raised a + # MethodError, which the handler had no catch for, so every `study_tree` + # request answered 500. Normalise once at this boundary, as + # `_load_run_table_df` does above and as `cladistics.jl` does at its own + # call into `_valid_tree_file`. (`lowercase` already returns a String.) tree_source = lowercase(strip(string(get(body, :treeSource, "auto")))) - tree_file = strip(string(get(body, :treeFile, ""))) - inline_newick = strip(string(get(body, :newick, ""))) + tree_file = String(strip(string(get(body, :treeFile, "")))) + inline_newick = String(strip(string(get(body, :newick, "")))) if tree_source == "study_tree" || (!isempty(tree_file) && tree_source == "auto") _valid_tree_file(tree_file) || return (nothing, nothing, :auto, From f5830a17416fd7e348f481934e5c61cb60db0326 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 08:23:57 +0000 Subject: [PATCH 3/7] fix(cladistics): divide weighted UniFrac by the library size, not the placed mass Codex raised two findings on the diversity block, and they are the same defect seen from two sides. Retaining the raw like-weight ratio stopped a lone 0.6 placement from being inflated to 1.0 in `mat`, but `weighted_unifrac` then divided by `sum(values(ab))` -- the placed-only mass -- which restored exactly the weight that had just been withheld. The same denominator also erased ASVs the tree does not cover: a sample could be 90% unmapped and still be compared as though every read had been placed. `_sample_total` now lets a caller declare the real library size, and `weighted_unifrac`, `branch_unifrac_contributions` and `unifrac_matrix` take it (`total_a`/`total_b`/`totals`). All default to `nothing` and the mass on the tree, so every existing call is unchanged. A declared total below the on-tree mass is rejected. The route computes one denominator per sample -- the sum of the table column -- and reports `library size - on-tree mass` as `unplaced_reads`, which now covers the like-weight shortfall and uncovered ASVs together. Faith's PD and unweighted UniFrac read presence only and are untouched. Two unrelated defects the suite could not previously reach: `_reaggregate_after_move` built `subtree_pd` by summing each child's branch length plus its subtree, which drops the branch subtending the node itself -- the opposite of `clade_faith_pd`, which the builder uses. Every proposed clade therefore came back short by its own branch and was compared against a current tree measured the other way: moving ASV_4 into a clade of PD 0.60 reported 0.55. It now mirrors `_pd_accum`, and the test pins the exact values plus total-PD conservation. `test_clade_cumulus.jl` bound `const SV`; the thirteen later unit files that bind `SV = MetaManifold.Server` plainly then threw `invalid assignment to constant Main.SV`. The throw escaped the include loop in `runtests.jl`, so roughly thirty unit files and the occupancy probe never ran at all. The binding is now plain, and each file is included inside its own testset so a top-level throw is attributed to its file and does not take the rest of the suite down with it. --- frontend/src/api/types.ts | 10 ++- frontend/src/components/CladeCumulusPanel.tsx | 14 ++-- src/analysis/clade_cumulus.jl | 82 +++++++++++++++++-- src/server/routes/clade_cumulus.jl | 53 ++++++------ test/runtests.jl | 9 +- test/unit/test_clade_cumulus.jl | 82 +++++++++++++++++-- 6 files changed, 200 insertions(+), 50 deletions(-) diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index be55f16a..249f9a12 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -969,10 +969,12 @@ export interface CladeCumulusDiversity { samples: string[] mapped_tips: number /** - * Reads per sample the `.jplace` file left on no edge: the like-weight - * ratios kept for a query need not reach 1, and the shortfall is counted - * here rather than redistributed onto the edges that survived. Absent for - * responses produced before this field existed. + * Reads per sample that no tip on the tree carries — a query whose kept + * `.jplace` like-weight ratios do not reach 1 contributes only that share, + * and an ASV the tree does not cover contributes nothing. This mass is held + * out of the weighted UniFrac denominators rather than renormalised onto the + * tips that were placed. Absent for responses produced before this field + * existed. */ unplaced_reads?: Record faith_pd: Record diff --git a/frontend/src/components/CladeCumulusPanel.tsx b/frontend/src/components/CladeCumulusPanel.tsx index 71ea0cb9..503d6ccc 100644 --- a/frontend/src/components/CladeCumulusPanel.tsx +++ b/frontend/src/components/CladeCumulusPanel.tsx @@ -712,15 +712,17 @@ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { {/* - Reads the placement file left on no edge. They are excluded - from the abundances above rather than renormalised onto the - edges a query did keep, so the shortfall is stated instead of - being silently absorbed. + Reads no tip carries: like-weight mass a placement left on no + edge, plus ASVs the tree does not cover. They are held out of + the weighted UniFrac denominators rather than renormalised onto + the tips that were placed, so the shortfall is stated here + instead of being silently absorbed. */} {unplacedSamples.length > 0 && (

- Unplaced reads (like-weight mass the placement file left on no edge, excluded - from the abundances above):{' '} + Unplaced reads (carried by no tip — like-weight mass left on no edge, plus + ASVs the tree does not cover — and held out of the weighted UniFrac + denominators):{' '} {unplacedSamples .map(([s, n]) => `${s} ${n.toFixed(2)}`) .join(', ')} diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 34e68f48..02efd2c6 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -764,11 +764,48 @@ function unweighted_unifrac( end """ - weighted_unifrac(counts_a, counts_b, taxa, tree; normalized=false) -> Float64 + _sample_total(abund, declared, name) -> Float64 + +The denominator for one sample's relative abundances. It defaults to the mass +that actually sits on the tree, which is the ordinary case. A caller that knows +part of the sample could not be placed on the tree -- `.jplace` like-weight mass +the tool left on no edge, or table ASVs the tree does not cover -- passes the +true library size instead, so that share stays off the tree rather than being +renormalised onto the tips that were placed. Without this, down-weighting a tip +to express placement uncertainty achieves nothing: dividing by the placed-only +total scales it straight back up. +""" +function _sample_total( + abund::Dict{String,Float64}, + declared::Union{Real,Nothing}, + name::AbstractString, +)::Float64 + on_tree = sum(values(abund); init = 0.0) + declared === nothing && return on_tree + t = Float64(declared) + (isfinite(t) && t >= 0.0) || + throw(ArgumentError("$name must be finite and non-negative, got $declared")) + # Below what the tree already carries would put a relative abundance above + # one. That is a caller error, not a discrepancy to round away -- but the + # equality case must survive floating-point summation, hence the slack. + t >= on_tree * (1 - 1e-9) - 1e-9 || + throw(ArgumentError( + "$name ($t) is below the $on_tree reads this sample already maps onto the tree")) + max(t, on_tree) +end + +""" + weighted_unifrac(counts_a, counts_b, taxa, tree; normalized=false, + total_a=nothing, total_b=nothing) -> Float64 Lozupone et al. (2007) raw (`normalized=false`) or normalized (`normalized=true`) weighted UniFrac distance: branch-length-weighted L1 distance between the two samples' CladeCumulus cumulative relative abundance trees. + +`total_a` / `total_b` override the per-sample denominator, which otherwise is +the mass on the tree. Pass the sample's true library size when part of it could +not be placed, so the unplaced share is absent from the tree instead of being +renormalised onto the tips that were. """ function weighted_unifrac( counts_a::AbstractVector{<:Real}, @@ -776,12 +813,14 @@ function weighted_unifrac( taxa::AbstractVector{<:AbstractString}, tree::Union{PhyloNode,AbstractString}; normalized::Bool = false, + total_a::Union{Real,Nothing} = nothing, + total_b::Union{Real,Nothing} = nothing, )::Float64 root = _ensure_rooted_tree(tree) ab_a = _validate_counts_and_taxa(counts_a, taxa, root) ab_b = _validate_counts_and_taxa(counts_b, taxa, root) - tot_a = sum(values(ab_a); init = 0.0) - tot_b = sum(values(ab_b); init = 0.0) + tot_a = _sample_total(ab_a, total_a, "total_a") + tot_b = _sample_total(ab_b, total_b, "total_b") (tot_a <= 0.0 && tot_b <= 0.0) && return 0.0 function walk(u::PhyloNode, depth_from_root::Float64, is_root::Bool)::Tuple{Float64,Float64,Float64,Float64} @@ -822,13 +861,17 @@ function branch_unifrac_contributions( counts_a::AbstractVector{<:Real}, counts_b::AbstractVector{<:Real}, taxa::AbstractVector{<:AbstractString}, - tree::Union{PhyloNode,AbstractString}, + tree::Union{PhyloNode,AbstractString}; + total_a::Union{Real,Nothing} = nothing, + total_b::Union{Real,Nothing} = nothing, )::Vector{OrderedDict{String,Any}} root = _ensure_rooted_tree(tree) ab_a = _validate_counts_and_taxa(counts_a, taxa, root) ab_b = _validate_counts_and_taxa(counts_b, taxa, root) - tot_a = sum(values(ab_a); init = 0.0) - tot_b = sum(values(ab_b); init = 0.0) + # Same denominators as `weighted_unifrac`, so the per-branch table adds up + # to the distance it is meant to explain. + tot_a = _sample_total(ab_a, total_a, "total_a") + tot_b = _sample_total(ab_b, total_b, "total_b") contribs = OrderedDict{String,Any}[] function walk(u::PhyloNode, is_root::Bool)::Tuple{Float64,Float64} @@ -882,11 +925,17 @@ function unifrac_matrix( tree::Union{PhyloNode,AbstractString}; metric::Symbol = :weighted, normalized::Bool = false, + totals::Union{AbstractVector{<:Real},Nothing} = nothing, )::Matrix{Float64} size(sample_matrix, 1) == length(taxa) || throw(ArgumentError("sample_matrix rows ($(size(sample_matrix, 1))) != length(taxa) ($(length(taxa)))")) size(sample_matrix, 2) == length(samples) || throw(ArgumentError("sample_matrix cols ($(size(sample_matrix, 2))) != length(samples) ($(length(samples)))")) + # Per-sample library sizes for the weighted metric, in `samples` order. See + # `_sample_total`: they keep mass that could not be placed off the tree + # rather than renormalised onto the tips that were. + isnothing(totals) || length(totals) == length(samples) || + throw(ArgumentError("totals ($(length(totals))) != length(samples) ($(length(samples)))")) root = _ensure_rooted_tree(tree) k = length(samples) M = zeros(Float64, k, k) @@ -897,7 +946,9 @@ function unifrac_matrix( d = if metric === :unweighted unweighted_unifrac(ci, cj, taxa, root) elseif metric === :weighted - weighted_unifrac(ci, cj, taxa, root; normalized = normalized) + weighted_unifrac(ci, cj, taxa, root; normalized = normalized, + total_a = isnothing(totals) ? nothing : totals[i], + total_b = isnothing(totals) ? nothing : totals[j]) else throw(ArgumentError("Unsupported UniFrac metric $(repr(metric))")) end @@ -1819,10 +1870,16 @@ function _reaggregate_after_move( end active_here = u.cumulative_count > 0.0 st = Epistemic.epistemic_status(case_obj, u.id) + # `subtree_pd` has to mean here exactly what it means on the tree + # this proposal is compared against. `clade_faith_pd`/`_pd_accum`, + # which the builder uses, counts a node's OWN branch when the + # subtree below it holds an observed tip -- so an observed tip + # contributes its own branch length, not zero. + leaf_pd = (active_here && !is_root) ? max(0.0, u.branch_length) : 0.0 leaf_node = CladeNode( u.id, u.name, u.rank, u.direct_count, u.cumulative_count, freq, u.residual_count, Epistemic.cloud_size_by_residual(u.residual_count), - u.branch_length, u.support, 0.0, u.dominant_taxon, u.taxon_purity, mono, + u.branch_length, u.support, leaf_pd, u.dominant_taxon, u.taxon_purity, mono, st, Epistemic.epistemic_colour(st), CladeNode[], ) return (leaf_node, tmap, 1, active_here) @@ -1844,13 +1901,20 @@ function _reaggregate_after_move( n_sub_tips += nt if ch_active any_active = true - sub_pd += max(0.0, rk.branch_length) + rk.subtree_pd + # `rk.subtree_pd` already carries the child's own branch length. + sub_pd += rk.subtree_pd end for (k, v) in ct merged_taxa[k] = get(merged_taxa, k, 0) + v end end + # Mirror `_pd_accum`: a non-root clade whose subtree is observed also + # contributes the branch that subtends it. + if any_active && !is_root + sub_pd += max(0.0, u.branch_length) + end + dom_tx = nothing dom_n = 0 for (tx, c) in merged_taxa diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index 8ca8d2f5..2dca43c3 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/src/server/routes/clade_cumulus.jl @@ -323,30 +323,34 @@ end end end - # Reads whose placement the file left on no edge at all: the - # like-weight ratios kept for a query need not reach 1, and the - # shortfall belongs to edges the placement tool discarded. Those - # reads are counted here rather than redistributed onto the edges - # that survived, so an abundance-weighted result can be read - # against the evidence that was actually placed. A row whose tips - # carry the full mass (every non-jplace ASV, and any partition - # complete within `PLACEMENT_WEIGHT_SUM_TOL`) contributes nothing. - row_weight = Dict{Int,Float64}() - for (ri, w) in taxa_rows - row_weight[ri] = get(row_weight, ri, 0.0) + w - end - short_rows = [(ri, 1.0 - wsum) for (ri, wsum) in row_weight - if wsum < 1.0 - MetaManifold.CladeCumulus.PLACEMENT_WEIGHT_SUM_TOL] + # Each sample's library size over the whole table, and the part of + # it that no tip carries. Reads go missing from `mat` two ways: a + # query whose kept like-weight ratios do not reach 1 contributes + # only that share, and an ASV the tree does not cover contributes + # nothing at all. Both are the same thing to a diversity metric -- + # mass that is not on the tree. + # + # `sample_totals` is then handed to the weighted metrics as their + # denominator. Without it they divide by the placed-only total, so + # a sample whose single mapped query was placed at LWR 0.6 comes + # back out at relative abundance 1.0 and the raw like-weight + # weighting above is undone at the last step. Faith's PD and + # unweighted UniFrac read presence only and are unaffected. + sample_totals = Float64[] unplaced_reads = OrderedDict{String,Float64}() - for sc in scols - col_vals = df[!, sc] + for (ci, sc) in enumerate(scols) tot = 0.0 - for (ri, lost_share) in short_rows - v = col_vals[ri] + for v in df[!, sc] ismissing(v) && continue - tot += max(0.0, Float64(v)) * lost_share + tot += max(0.0, Float64(v)) end - unplaced_reads[sc] = tot + placed = sum(@view mat[:, ci]) + # A multi-candidate partition only has to sum to one within + # `PLACEMENT_WEIGHT_SUM_TOL`, so the placed mass can sit a hair + # above the table total. The denominator must never fall below + # what the tree actually carries. + push!(sample_totals, max(tot, placed)) + unplaced_reads[sc] = max(0.0, tot - placed) end pd_by_sample = OrderedDict{String,Float64}( sc => MetaManifold.CladeCumulus.faith_pd(@view(mat[:, ci]), taxa_vec, cc_res.phylo_root) @@ -356,14 +360,17 @@ end mat, taxa_vec, scols, cc_res.phylo_root; metric = :unweighted ) w_raw = MetaManifold.CladeCumulus.unifrac_matrix( - mat, taxa_vec, scols, cc_res.phylo_root; metric = :weighted, normalized = false + mat, taxa_vec, scols, cc_res.phylo_root; metric = :weighted, normalized = false, + totals = sample_totals ) w_norm = MetaManifold.CladeCumulus.unifrac_matrix( - mat, taxa_vec, scols, cc_res.phylo_root; metric = :weighted, normalized = true + mat, taxa_vec, scols, cc_res.phylo_root; metric = :weighted, normalized = true, + totals = sample_totals ) top_branches = if length(scols) >= 2 first_n = MetaManifold.CladeCumulus.branch_unifrac_contributions( - @view(mat[:, 1]), @view(mat[:, 2]), taxa_vec, cc_res.phylo_root + @view(mat[:, 1]), @view(mat[:, 2]), taxa_vec, cc_res.phylo_root; + total_a = sample_totals[1], total_b = sample_totals[2] ) first_n[1:min(length(first_n), 15)] else diff --git a/test/runtests.jl b/test/runtests.jl index 88cdda1b..f199c368 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,8 +93,15 @@ using MetaManifold.Categories, MetaManifold.CompositionLibrary ## Unit tests (always run) @testset "MetaManifold" begin + # Each file gets its own testset so that a throw at a file's top level -- outside + # any of that file's own testsets -- is reported against the file that caused it + # and does not abort the loop. Bare `include` here meant one such throw silently + # skipped every remaining unit file (and the occupancy probe below) while the + # summary still looked like a single stray error. for f in _selected_files() - include(joinpath(@__DIR__, "unit", f)) + @testset "$f" begin + include(joinpath(@__DIR__, "unit", f)) + end end ## Occupancy probe and the two-thread disagreement test (MMW-01), with the diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index 14f34421..46708281 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -14,7 +14,22 @@ const EP = MetaManifold.Epistemic const CC = MetaManifold.CladeCumulus const DM = MetaManifold.DiversityMetrics const PH = MetaManifold.Phylogeny -const SV = MetaManifold.Server +# NOT `const`: every other unit file binds `SV = MetaManifold.Server` with a plain +# assignment, and all of them are `include`d into the same `Main`. Declaring it +# constant here -- this file runs sixth -- made the next plain assignment throw +# `invalid assignment to constant Main.SV`, which aborted the include loop and +# silently skipped the remaining unit files. +SV = MetaManifold.Server + +"""Depth-first lookup of a clade node by `id` in a decoded clade-tree payload.""" +function _cc_find(node, id::AbstractString) + String(node.id) == id && return node + for ch in node.children + hit = _cc_find(ch, id) + hit === nothing || return hit + end + nothing +end _cc_req(method::String, path::String, body::String = "") = SV.Oxygen.internalrequest( @@ -110,6 +125,41 @@ _cc_req(method::String, path::String, body::String = "") = @test CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; normalized = false) ≈ 2.3276502663333334 atol = 1e-11 @test CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; normalized = true) ≈ 0.6655215451489284 atol = 1e-11 + # An explicit per-sample total keeps mass that is not on the tree off + # it, instead of renormalising it onto the tips that are. Without this + # the weighted metric divides by the placed-only sum, so a sample whose + # placements carry a fraction of its reads is scaled straight back up + # to a full community and the down-weighting achieves nothing. + # + # Against an all-zero partner every `pb` is 0, so the raw distance is + # exactly `sum(branch_length * pa)` and is linear in `1 / total_a`: + # declaring twice the library size must halve it. + f_zero = zeros(Int, length(ft_taxa)) + d_placed_only = CC.weighted_unifrac(fa, f_zero, ft_taxa, ft_nwk; normalized = false) + d_true_total = CC.weighted_unifrac(fa, f_zero, ft_taxa, ft_nwk; + normalized = false, total_a = 2 * sum(fa)) + @test d_placed_only > 0.0 + @test d_true_total ≈ d_placed_only / 2 atol = 1e-12 + # Declaring the mass already on the tree is the default, to the bit. + @test CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; normalized = false, + total_a = sum(fa), total_b = sum(fb)) ≈ + CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; normalized = false) + # A total below what the tree already carries would push a relative + # abundance above one, so it is rejected rather than clamped. + @test_throws ArgumentError CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; total_a = 1.0) + @test_throws ArgumentError CC.weighted_unifrac(fa, fb, ft_taxa, ft_nwk; total_b = -1.0) + @test_throws ArgumentError CC.unifrac_matrix( + hcat(Float64.(fa), Float64.(fb)), ft_taxa, ["A", "B"], ft_nwk; + metric = :weighted, totals = [1.0]) + + # The per-branch breakdown uses the same denominators, so the table + # still explains the distance it sits next to. + contribs_default = CC.branch_unifrac_contributions(fa, f_zero, ft_taxa, ft_nwk) + contribs_total = CC.branch_unifrac_contributions(fa, f_zero, ft_taxa, ft_nwk; + total_a = 2 * sum(fa)) + @test sum(Float64(c["contribution"]) for c in contribs_total) ≈ + sum(Float64(c["contribution"]) for c in contribs_default) / 2 atol = 1e-12 + # Quoted labels, NHX support, negative branch clamping, non-finite branch rejection, # and midpoint split-support symmetry across both root children. warns = String[] @@ -428,7 +478,19 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 @test v_allow.proposal_file !== nothing @test isfile(joinpath(trees_dir, ".clade_cumulus", String(v_allow.proposal_file))) @test v_allow.proposed_root.cumulative_count ≈ 110.0 atol = 1e-9 - @test Float64(v_allow.proposed_root.children[1].subtree_pd) > before_target_pd + + # The proposal reports `subtree_pd` on exactly the same footing as the tree + # it is compared against (`clade_faith_pd`: a clade's own subtending branch + # counts whenever an observed tip sits below it). Moving ASV_4 hands its + # 0.25 branch to the target, so the target reads 0.30 + 0.10 + 0.20 + 0.25 + # and the clade it left falls to 0.40 + 0.15. + moved_target = _cc_find(v_allow.proposed_root, String(high_clade.id)) + @test moved_target !== nothing + @test Float64(moved_target.subtree_pd) ≈ 0.85 atol = 1e-9 + @test Float64(moved_target.subtree_pd) > before_target_pd + @test Float64(_cc_find(v_allow.proposed_root, String(low_clade.id)).subtree_pd) ≈ 0.55 atol = 1e-9 + # Relocating an observed tip cannot create or destroy phylogenetic diversity. + @test Float64(v_allow.proposed_root.subtree_pd) ≈ Float64(comb.root.subtree_pd) atol = 1e-9 # Re-saving under a different threshold writes a distinct immutable sidecar file r_allow2 = _cc_req("POST", "/api/v1/studies/StudyCC/runs/run1/clade-cumulus/validate-move", JSON3.write(Dict( @@ -448,9 +510,15 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 # 5. A `.jplace` study tree: the diversity block weights each tip by # the raw like-weight ratio, so mass the placement file left on no # edge is reported rather than renormalised onto the edge that was - # kept. ASV_1 is placed once at LWR 0.60, ASV_2 once at LWR 1.00, so - # ASV_1 contributes only 60% of its reads and the missing 40% - # (SampleA 40 x 0.4 = 16, SampleB 0) is counted as unplaced. + # kept. ASV_1 is placed once at LWR 0.60 and ASV_2 once at LWR 1.00; + # ASV_3 and ASV_4 are not placed at all. + # + # SampleA holds 40 + 10 + 0 + 5 = 55 reads, of which the tree + # carries 40 x 0.6 + 10 = 34, so 21 are unplaced: 16 of ASV_1's + # like-weight shortfall plus ASV_4's 5 unmapped reads. + # SampleB holds 0 + 20 + 25 + 10 = 55 and the tree carries 20, so + # 35 are unplaced (ASV_3's 25 and ASV_4's 10, ASV_1 having none + # here). write(joinpath(trees_dir, "placed.jplace"), JSON3.write(OrderedDict( "version" => 3, "tree" => "((Ref_1:0.10{0},Ref_2:0.20{1})0.96:0.30{2},Ref_3:0.40{3}){4};", @@ -472,8 +540,8 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 jp_body = JSON3.read(String(r_jp.body)) @test jp_body.diversity !== nothing @test jp_body.diversity.mapped_tips == 2 - @test jp_body.diversity.unplaced_reads.SampleA ≈ 16.0 atol = 1e-9 - @test jp_body.diversity.unplaced_reads.SampleB ≈ 0.0 atol = 1e-9 + @test jp_body.diversity.unplaced_reads.SampleA ≈ 21.0 atol = 1e-9 + @test jp_body.diversity.unplaced_reads.SampleB ≈ 35.0 atol = 1e-9 # The clade tree is a count view of the table and still conserves # every observed read; only the abundance-weighted diversity matrix # discounts the unplaced mass. From 7a6ead99535f75178948d624892fb1bd1603ebca Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 09:09:02 +0000 Subject: [PATCH 4/7] fix(stats): implement the missing erfc, and repair two never-run tests With the include loop no longer aborting, `test_differential.jl` and the server smoke tests ran for the first time and found three defects. `_wald` called `erfc`, which exists in no module this project loads -- `SpecialFunctions` is deliberately not a dependency, and the header of `compositional.jl` says so. Every multinomial and Dirichlet-multinomial fit therefore died with `UndefVarError: erfc` the moment it tried to report a p-value, which is to say the models could not be used at all. `compositional_erfc` now sits beside `compositional_loggamma` and `compositional_digamma` and is built the same way, from the regularised incomplete gamma function at a = 1/2: the power series for P below x^2 = 1.5 and the Legendre continued fraction for Q above it, needing only `exp` and `log`. Measured against the C library's erfc the relative error stays under 1e-14 wherever erfc(x) >= 1e-50, i.e. out to a two-sided deviate of |z| = 15, which matches the accuracy the other two already promise. The unit tests pin both sides of the crossover, the reflection, and the two-sided p-value that is its only caller. `min_prevalence` is the second positional argument of `DifferentialConfig`, not a keyword, so the "every taxon filtered out" case constructed a `MethodError` instead of reaching the engine. The same call also passed four samples that were all in group A, which the "each group needs at least one sample" guard would have answered long before the prevalence filter. Both are fixed, so the case now asserts what it says it does. --- src/analysis/compositional.jl | 75 +++++++++++++++++++++++++++++++--- test/unit/test_differential.jl | 32 +++++++++++++-- 2 files changed, 99 insertions(+), 8 deletions(-) diff --git a/src/analysis/compositional.jl b/src/analysis/compositional.jl index c05ce892..ab0ffac1 100644 --- a/src/analysis/compositional.jl +++ b/src/analysis/compositional.jl @@ -40,9 +40,10 @@ # # The special functions are implemented here rather than pulled from # SpecialFunctions, so the models add no dependency: `_loggamma` and `_digamma` -# are the Stirling series after the recurrence is applied up to x >= 10, -# accurate to about 1e-14 (checked against Julia's `lgamma` and known constants -# in the unit tests). +# are the Stirling series after the recurrence is applied up to x >= 10, and +# `compositional_erfc` is the regularised incomplete gamma function at a = 1/2. +# All three are accurate to about 1e-14 (checked against Julia's `lgamma` and +# known constants in the unit tests). module CompositionalModels export CompositionalFit, CompositionalRefusal, RuntimeRefusal, MAX_TAXA, @@ -50,7 +51,7 @@ export CompositionalFit, CompositionalRefusal, RuntimeRefusal, MAX_TAXA, CONCENTRATION_MAX, DM_METHODS, fit_multinomial, fit_dirichlet_multinomial, multinomial_loglik, dirichlet_multinomial_loglik, dm_concentration_mle, dm_concentration_mom, observed_information, composition_probabilities, - compositional_loggamma, compositional_digamma + compositional_loggamma, compositional_digamma, compositional_erfc # At most this many taxa are fitted. The observed information is # ((K-1) q) x ((K-1) q) doubles: 500 taxa with two covariates is about 8 MB, @@ -193,6 +194,70 @@ function compositional_digamma(x::Real) inv2 * (1 / 252 - inv2 * (1 / 240 - inv2 / 132)))) - shift end +# erfc through the regularised incomplete gamma function, which needs only +# `exp` and `log`: erfc(x) = Q(1/2, x^2) for x >= 0. The series for P and the +# continued fraction for Q each converge quickly on the side where the other +# one does not, and they cross over at x^2 = 1.5. Measured against the C +# library's erfc, the relative error is below 1e-14 wherever erfc(x) >= 1e-50 +# (|z| <= 15 as a two-sided normal deviate), which is the same order as +# `compositional_loggamma` and `compositional_digamma` above. +const _ERFC_LOG_SQRT_PI = 0.57236494292470008707 # log Gamma(1/2) +const _ERFC_EPS = 3e-16 +const _ERFC_FPMIN = 1e-300 +const _ERFC_MAX_ITER = 1000 + +# P(1/2, t) by its power series: sum t^n / ((1/2)(3/2)...(1/2 + n)). +function _lower_gamma_half(t::Float64) + ap = 0.5 + term = 2.0 # 1 / (1/2) + total = term + for _ in 1:_ERFC_MAX_ITER + ap += 1.0 + term *= t / ap + total += term + abs(term) < abs(total) * _ERFC_EPS && break + end + total * exp(-t + 0.5 * log(t) - _ERFC_LOG_SQRT_PI) +end + +# Q(1/2, t) by the Legendre continued fraction, evaluated with Lentz's method. +function _upper_gamma_half(t::Float64) + a = 0.5 + b = t + 1.0 - a + c = 1.0 / _ERFC_FPMIN + d = 1.0 / b + h = d + for i in 1:_ERFC_MAX_ITER + an = -i * (i - a) + b += 2.0 + d = an * d + b + abs(d) < _ERFC_FPMIN && (d = _ERFC_FPMIN) + c = b + an / c + abs(c) < _ERFC_FPMIN && (c = _ERFC_FPMIN) + d = 1.0 / d + delta = d * c + h *= delta + abs(delta - 1.0) <= _ERFC_EPS && break + end + h * exp(-t + 0.5 * log(t) - _ERFC_LOG_SQRT_PI) +end + +""" + compositional_erfc(x) -> Float64 + +The complementary error function. `compositional_erfc(abs(z) / sqrt(2))` is the +two-sided p-value of a standard normal deviate `z`, which is the only thing this +module asks of it. Throws `ArgumentError` for a non-finite argument. +""" +function compositional_erfc(x::Real) + xf = Float64(x) + isfinite(xf) || throw(ArgumentError("erfc needs a finite argument, not $x")) + xf < 0.0 && return 2.0 - compositional_erfc(-xf) + xf == 0.0 && return 1.0 + t = xf * xf + t < 1.5 ? 1.0 - _lower_gamma_half(t) : _upper_gamma_half(t) +end + ## Geometry """ @@ -860,7 +925,7 @@ function _wald(coefficients::Matrix{Float64}, variances::Vector{Float64}) se[l, k] = s if s > 0.0 && isfinite(s) z[l, k] = coefficients[l, k] / s - p[l, k] = erfc(abs(z[l, k]) / sqrt(2.0)) + p[l, k] = compositional_erfc(abs(z[l, k]) / sqrt(2.0)) else z[l, k] = NaN p[l, k] = NaN diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index d0f44e47..cbcc5831 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -1535,6 +1535,26 @@ end @test CM.compositional_digamma(1.0) ≈ -0.5772156649015329 rtol=1e-12 @test CM.compositional_digamma(10.0) ≈ 2.251752589066721 rtol=1e-12 @test CM.compositional_digamma(0.5) ≈ -1.9635100260214235 rtol=1e-12 + + # erfc is the regularised incomplete gamma at a = 1/2, with the series + # and the continued fraction swapping over at x^2 = 1.5. Both sides of + # that crossover are pinned here, against the C library's erfc. + @test CM.compositional_erfc(0.0) == 1.0 + @test CM.compositional_erfc(0.1) ≈ 0.8875370839817152 rtol=1e-13 + @test CM.compositional_erfc(0.5) ≈ 0.47950012218695348 rtol=1e-13 + @test CM.compositional_erfc(1.2) ≈ 0.08968602177036465 rtol=1e-13 # x^2 < 1.5 + @test CM.compositional_erfc(1.3) ≈ 0.06599205505934755 rtol=1e-13 # x^2 > 1.5 + @test CM.compositional_erfc(2.0) ≈ 0.0046777349810472654 rtol=1e-13 + @test CM.compositional_erfc(5.0) ≈ 1.5374597944280351e-12 rtol=1e-13 + # erfc(-x) = 2 - erfc(x), so the reflection has to be exact. + @test CM.compositional_erfc(-1.0) ≈ 1.842700792949715 rtol=1e-13 + @test CM.compositional_erfc(-3.0) + CM.compositional_erfc(3.0) ≈ 2.0 rtol=1e-15 + @test_throws ArgumentError CM.compositional_erfc(NaN) + @test_throws ArgumentError CM.compositional_erfc(Inf) + + # The only call site: the two-sided p-value of a normal deviate. + @test CM.compositional_erfc(1.959963984540054 / sqrt(2.0)) ≈ 0.05 rtol=1e-12 + @test CM.compositional_erfc(6.0 / sqrt(2.0)) ≈ 1.9731752900754024e-09 rtol=1e-13 end @testset "Multinomial: the DA engine reports a composition, not per-taxon fits" begin @@ -1644,11 +1664,17 @@ end @test e isa ArgumentError @test occursin("is not fitted", sprint(showerror, e)) # Every taxon filtered out is an explicit error, not an empty table. + # `min_prevalence` is the second positional argument of + # `DifferentialConfig`, not a keyword; and both groups have to be + # populated, or the earlier "each group needs at least one sample" + # guard answers first and the prevalence filter is never reached. + # Taxon 'a' is in 3 of 4 samples and 'b' in 2, so a threshold of 1.0 + # leaves nothing to fit. e = try differential_abundance([1 0; 2 3; 0 4; 5 0], ["s1", "s2", "s3", "s4"], - ["a", "b"], groups8[1:4]; - kw..., config=DifferentialConfig(; method="multinomial", - min_prevalence=1.0)) + ["a", "b"], ["A", "A", "B", "B"]; + kw..., config=DifferentialConfig("tss", 1.0; + method="multinomial")) nothing catch err err From 7c3661d3ebd5781588800846be29e37c2f65d89e Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 16:11:05 +0100 Subject: [PATCH 5/7] Update test/unit/test_routes.jl Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- test/unit/test_routes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/test_routes.jl b/test/unit/test_routes.jl index 2a9fa59d..f97b9bba 100644 --- a/test/unit/test_routes.jl +++ b/test/unit/test_routes.jl @@ -953,8 +953,8 @@ end # The Cladistics tab and the run-level CladeCumulus evidence view are # served by two different route files; both must be registered. paths = Set(p for (_, p) in registered) - @test any(occursin("/cladistics", p) for p in paths) - @test any(occursin("/clade-cumulus", p) for p in paths) + @test "/api/v1/studies/{study}/cladistics" in paths + @test "/api/v1/studies/{study}/runs/{run}/clade-cumulus/sources" in paths end # Testsets above point the server at temporary roots. From d2d63426d7ea2e2da34ed2a0f2dabdd0334d14de Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Fri, 23 Oct 2026 15:17:54 +0000 Subject: [PATCH 6/7] Document CladeCumulus parsing, UniFrac APIs, and panel behavior --- frontend/src/components/CladeCumulusPanel.tsx | 5 ++ src/analysis/clade_cumulus.jl | 68 +++++++++++++++++-- 2 files changed, 68 insertions(+), 5 deletions(-) diff --git a/frontend/src/components/CladeCumulusPanel.tsx b/frontend/src/components/CladeCumulusPanel.tsx index 503d6ccc..2a73c13c 100644 --- a/frontend/src/components/CladeCumulusPanel.tsx +++ b/frontend/src/components/CladeCumulusPanel.tsx @@ -33,6 +33,11 @@ export interface CladeCumulusPanelProps { runs: CladeCumulusRunOption[] } +/** + * Explore a study's run trees, diversity and unplaced reads, with CSV export + * and move previews that can be saved as sidecar proposals. Move checks use + * the last computed inputs until the user recomputes the result. + */ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { const toast = useToast() const reqCounter = useRef(0) diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 4f83a138..09ce8e54 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -489,13 +489,20 @@ end """ parse_jplace(raw_json::AbstractString) -> JplaceDocument -Strictly parse a `.jplace` document before any display normalization: +Strictly parse a `.jplace` document before any display normalisation: - Validates unique integer `{edge_num}` annotations across the tree - Validates explicit unique `"like_weight_ratio"` column in `"fields"` -- Validates each query's `p` candidates, finite weights in `[0, 1]`, edge - references present on non-root edges, and total weight within `1e-5` of `1.0` +- Validates each query's `p` candidates, finite weights in `[0, 1]`, and edge + references present on non-root edges. A single candidate needs weight at + least `0.5`; multiple candidates need total weight within `1e-5` of `1.0`. - Preserves raw query/candidate ordinals and `n` / `nm` multiplicities without - re-normalizing weights or collapsing tied maxima. + renormalising weights or collapsing tied maxima. + +Return the parsed tree, edge index, query placements and tree warnings. Missing +or invalid candidate data and insufficient weights produce `:undecided` queries +with a reason rather than an exception. Invalid document structure, edge +annotations, required columns or query names raise `ArgumentError`; Newick +parsing, JSON decoding and data conversion errors propagate. """ function parse_jplace(raw_json::AbstractString)::JplaceDocument doc = JSON3.read(raw_json) @@ -774,6 +781,11 @@ true library size instead, so that share stays off the tree rather than being renormalised onto the tips that were placed. Without this, down-weighting a tip to express placement uncertainty achieves nothing: dividing by the placed-only total scales it straight back up. + +`declared=nothing` returns the sum of `abund`; otherwise return at least that +sum, allowing a shortfall of `1e-9 * on_tree + 1e-9` for rounding. A declared +total that is non-finite, negative or further below the sum raises +`ArgumentError`. `name` identifies the offending parameter in these errors. """ function _sample_total( abund::Dict{String,Float64}, @@ -798,7 +810,7 @@ end weighted_unifrac(counts_a, counts_b, taxa, tree; normalized=false, total_a=nothing, total_b=nothing) -> Float64 -Lozupone et al. (2007) raw (`normalized=false`) or normalized (`normalized=true`) +Lozupone et al. (2007) raw (`normalized=false`) or normalised (`normalized=true`) weighted UniFrac distance: branch-length-weighted L1 distance between the two samples' CladeCumulus cumulative relative abundance trees. @@ -806,6 +818,18 @@ samples' CladeCumulus cumulative relative abundance trees. the mass on the tree. Pass the sample's true library size when part of it could not be placed, so the unplaced share is absent from the tree instead of being renormalised onto the tips that were. + +Counts correspond to `taxa`; repeated taxon names have their counts summed. +`tree` may be Newick text or a parsed tree and is midpoint-rooted when its root +has more than two children. The raw result uses the tree's branch-length units; +the normalised result divides by the sum of each tip's root distance weighted +by its relative abundance in both samples, returning `0.0` if that sum is zero. +Two samples with zero totals also return `0.0`. + +Raise `ArgumentError` for mismatched count and taxon lengths, non-finite or +negative counts, taxa absent from the tree, or explicit totals that are +non-finite, negative or below the mapped count sum beyond rounding tolerance. +Newick parsing errors propagate. """ function weighted_unifrac( counts_a::AbstractVector{<:Real}, @@ -857,6 +881,21 @@ function weighted_unifrac( den <= 0.0 ? 0.0 : num / den end +""" + branch_unifrac_contributions(counts_a, counts_b, taxa, tree; + total_a=nothing, total_b=nothing) + +Return positive contributions from non-root branches to raw weighted UniFrac, +sorted by decreasing contribution. Each record contains `clade_id`, `label`, +`is_tip`, `branch_length`, cumulative relative abundances `freq_a` and `freq_b`, +and `contribution` (branch length times the absolute frequency difference). +The contributions sum to [`weighted_unifrac`](@ref) with `normalized=false` +and the same totals; no positive contributions yields an empty vector. + +Counts, taxa, tree rooting and optional sample totals have the same meaning +and validation as in `weighted_unifrac`. Its `ArgumentError` conditions and +Newick parsing errors propagate here too. +""" function branch_unifrac_contributions( counts_a::AbstractVector{<:Real}, counts_b::AbstractVector{<:Real}, @@ -918,6 +957,25 @@ function branch_unifrac_contributions( contribs end +""" + unifrac_matrix(sample_matrix, taxa, samples, tree; + metric=:weighted, normalized=false, totals=nothing) -> Matrix{Float64} + +Return a symmetric pairwise UniFrac distance matrix with a zero diagonal, in +`samples` order. Input rows correspond to `taxa` and columns to `samples`. +`metric=:weighted` uses [`weighted_unifrac`](@ref), including its `normalized` +option; `metric=:unweighted` uses [`unweighted_unifrac`](@ref). + +For weighted distances, `totals` supplies library sizes in `samples` order, +including unplaced reads; by default each sample uses its mapped count sum. +For unweighted distances, `normalized` and the values of `totals` are ignored. + +Raise `ArgumentError` for mismatched matrix dimensions or `totals` length. +Newick parsing errors propagate. Pairwise evaluation also propagates the +chosen metric's input validation errors and rejects unsupported metrics with +`ArgumentError`. With fewer than two samples, return the zero matrix without +pairwise validation. +""" function unifrac_matrix( sample_matrix::AbstractMatrix{<:Real}, taxa::AbstractVector{<:AbstractString}, From cece7c714cddba98a4ab5cdd08f8b69d6bd91438 Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Fri, 23 Oct 2026 15:20:45 +0000 Subject: [PATCH 7/7] fix(cladistics): clear stale abundance and PD after moving a clade's last child and clarify unplaced-read normalization --- frontend/src/api/types.ts | 8 ++--- frontend/src/components/CladeCumulusPanel.tsx | 11 +++---- src/analysis/clade_cumulus.jl | 8 +++-- test/unit/test_clade_cumulus.jl | 31 +++++++++++++++++++ 4 files changed, 45 insertions(+), 13 deletions(-) diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 249f9a12..a3c93d53 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -971,10 +971,10 @@ export interface CladeCumulusDiversity { /** * Reads per sample that no tip on the tree carries — a query whose kept * `.jplace` like-weight ratios do not reach 1 contributes only that share, - * and an ASV the tree does not cover contributes nothing. This mass is held - * out of the weighted UniFrac denominators rather than renormalised onto the - * tips that were placed. Absent for responses produced before this field - * existed. + * and an ASV the tree does not cover contributes nothing. This mass remains + * in each sample's library-size denominator for weighted UniFrac but is + * excluded from tree-tip abundances. Absent for responses produced before + * this field existed. */ unplaced_reads?: Record faith_pd: Record diff --git a/frontend/src/components/CladeCumulusPanel.tsx b/frontend/src/components/CladeCumulusPanel.tsx index 2a73c13c..dd59405b 100644 --- a/frontend/src/components/CladeCumulusPanel.tsx +++ b/frontend/src/components/CladeCumulusPanel.tsx @@ -718,16 +718,15 @@ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { {/* Reads no tip carries: like-weight mass a placement left on no - edge, plus ASVs the tree does not cover. They are held out of - the weighted UniFrac denominators rather than renormalised onto - the tips that were placed, so the shortfall is stated here - instead of being silently absorbed. + edge, plus ASVs the tree does not cover. They remain in each + sample's library-size denominator for weighted UniFrac but are + excluded from tree-tip abundances. */} {unplacedSamples.length > 0 && (

Unplaced reads (carried by no tip — like-weight mass left on no edge, plus - ASVs the tree does not cover — and held out of the weighted UniFrac - denominators):{' '} + ASVs the tree does not cover — included in each sample’s library-size + denominator for weighted UniFrac but excluded from tree-tip abundances):{' '} {unplacedSamples .map(([s, n]) => `${s} ${n.toFixed(2)}`) .join(', ')} diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 09ce8e54..4eb96fed 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -1919,21 +1919,23 @@ function _reaggregate_after_move( function recompute(u::CladeNode, is_root::Bool)::Tuple{CladeNode,Dict{String,Int},Int,Bool} if isempty(u.children) - freq = total_reads > 0.0 ? u.cumulative_count / total_reads : 0.0 + # A clade emptied by the graft retains only its own direct reads. + cum = u.direct_count + freq = total_reads > 0.0 ? cum / total_reads : 0.0 tmap = Dict{String,Int}() mono = u.monophyletic if u.dominant_taxon !== nothing tmap[u.dominant_taxon] = 1 mono = get(total_tips_per_taxon, u.dominant_taxon, 0) == 1 end - active_here = u.cumulative_count > 0.0 + active_here = cum > 0.0 # `subtree_pd` is node-local, as in `clade_faith_pd`: a node's own # (non-root) branch counts once, here, and its parent adds only the # value it returns. own_pd = (active_here && !is_root) ? max(0.0, u.branch_length) : 0.0 st = Epistemic.epistemic_status(case_obj, u.id) leaf_node = CladeNode( - u.id, u.name, u.rank, u.direct_count, u.cumulative_count, freq, + u.id, u.name, u.rank, u.direct_count, cum, freq, u.residual_count, Epistemic.cloud_size_by_residual(u.residual_count), u.branch_length, u.support, own_pd, u.dominant_taxon, u.taxon_purity, mono, st, Epistemic.epistemic_colour(st), CladeNode[], diff --git a/test/unit/test_clade_cumulus.jl b/test/unit/test_clade_cumulus.jl index e8f6c762..1a766344 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -354,6 +354,37 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 end end + @testset "Moving the last child clears clade abundance and PD activity" begin + df = DataFrame(SeqName = ["ASV_1", "ASV_2", "ASV_3"], Sample1 = [10, 20, 0]) + tree = CC.build_clade_tree(df; mode = :phylogenetic, + newick = "(((ASV_1:0.10)0.99:0.20)0.99:0.30,(ASV_2:0.40,ASV_3:0.50)0.99:0.60)Root;") + ancestor = tree.root.children[1] + source = only(ancestor.children) + target = tree.root.children[2] + proposal = CC.validate_drag_drop(tree, "ASV_1", target.id) + @test proposal["allowed"] + proposed = JSON3.read(JSON3.write(proposal["proposed_root"])) + emptied = _cc_find(proposed, source.id) + @test isempty(emptied.children) + @test emptied.cumulative_count == 0.0 + @test emptied.cumulative_frequency == 0.0 + @test emptied.subtree_pd == 0.0 + inactive_ancestor = _cc_find(proposed, ancestor.id) + @test inactive_ancestor.cumulative_count == 0.0 + @test inactive_ancestor.cumulative_frequency == 0.0 + @test inactive_ancestor.subtree_pd == 0.0 + @test proposed.cumulative_count == tree.root.cumulative_count == 30.0 + @test proposed.cumulative_frequency == 1.0 + @test _cc_find(proposed, target.id).cumulative_count == 30.0 + @test _cc_find(proposed, "ASV_1").subtree_pd ≈ 0.10 atol = 1e-12 + @test _cc_find(proposed, "ASV_3").subtree_pd == 0.0 + @test proposed.subtree_pd ≈ 1.10 atol = 1e-12 + @test _cc_find(proposed, target.id).subtree_pd ≈ 1.10 atol = 1e-12 + @test source.cumulative_count == 10.0 + @test only(source.children).id == "ASV_1" + @test tree.root.subtree_pd ≈ 1.60 atol = 1e-12 + end + @testset "Dual-hierarchy CladeCumulus, polyphyly detection, drag-and-drop, and HTTP routes" begin tmp = mktempdir() old_root = CLADE_SV.ServerState._root[]