diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 48df92c7..a3c93d53 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -968,6 +968,15 @@ export interface CladeCumulusRequest { export interface CladeCumulusDiversity { samples: string[] mapped_tips: number + /** + * 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 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 unweighted_unifrac: number[][] weighted_unifrac: number[][] diff --git a/frontend/src/components/CladeCumulusPanel.tsx b/frontend/src/components/CladeCumulusPanel.tsx index 5a6a1613..dd59405b 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) @@ -249,6 +254,17 @@ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { return div.normalized_weighted_unifrac }, [result, unifracKind]) + // Only samples that actually lost mass are worth naming; a run with a + // complete placement file (or one produced before the field existed) reports + // nothing here. + const unplacedSamples = useMemo<[string, number][]>(() => { + const div = result?.diversity + if (!div?.unplaced_reads) return [] + return div.samples + .map(s => [s, div.unplaced_reads?.[s] ?? 0] as [string, number]) + .filter(([, n]) => n > 0) + }, [result]) + const exportMatrixCsv = () => { if (!result?.diversity || !activeMatrix) return const csv = formatDistanceMatrixCsv(result.diversity.samples, activeMatrix) @@ -700,6 +716,23 @@ 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 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 — 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(', ')} +

+ )} + {/* Pairwise UniFrac table */} {activeMatrix && (
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 14230418..4eb96fed 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -28,11 +28,31 @@ const DEFAULT_TAXONOMIC_RANKS = [ "Class", "Order", "Family", "Genus", "Species", ] -# Placement tools may truncate a query's candidate list, so retained -# like-weight ratios can legitimately sum to slightly less than 1 (a single -# retained placement at 0.98 is complete evidence). Sums far from 1 indicate -# a malformed or partial record and stay undecided. -const PLACEMENT_WEIGHT_SUM_TOL = 0.05 +# A multi-candidate placement partitions one query's probability mass across +# the edges it kept, so those like-weight ratios must sum to one: anything else +# is a malformed or truncated record. This is a floating-point tolerance on an +# exact identity, not slack for missing mass. +const PLACEMENT_WEIGHT_SUM_TOL = 1e-5 + +# A single-candidate placement is the other shape the same file can take: the +# tool reports one row whose `like_weight_ratio` is the share of the mass behind +# that site, and the remainder is mass it did not place anywhere. Such a record +# cannot be held to the sum-to-one identity above, so it is judged on its own +# share instead and is decisive while at least half the mass supports it. +# +# These two rules replace a single blanket tolerance: widening +# `PLACEMENT_WEIGHT_SUM_TOL` far enough to admit a lone 0.98 placement would +# also stop a genuinely malformed multi-candidate partition from being caught. +# +# Accepting a placement and trusting all of its mass are different questions, +# and the two consumers answer the second one differently on purpose. +# `_graft_jplace_queries` builds a count view of the observed ASV table and +# conserves reads: an unsplit query keeps its whole count at its one edge. +# The `/clade-cumulus` diversity block instead weights each tip by the raw +# like-weight ratio, so the mass the tool never placed stays off the tree and +# is reported as `unplaced_reads` rather than renormalised onto the edge that +# happened to be kept. +const PLACEMENT_SINGLE_MIN_LWR = 0.5 ## Phylogenetic tree structures & Newick / jplace parsers mutable struct PhyloNode @@ -208,9 +228,11 @@ function parse_newick(src::AbstractString; warnings::Vector{String} = String[]): 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} - # An NHX comment immediately after ')' precedes the (usually empty) - # label, so the capture refs must be threaded through this first skip - # or the support tag is consumed as an anonymous comment. + # 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`). That + # leading comment is consumed by this first skip, so the capture refs + # have to be threaded through it or the tag is read as an anonymous + # comment and the support is lost. skip_ws_and_leading_comments!(capture_nhx = capture_nhx, capture_edge = capture_edge) pos[] > n && return ("", false) if bytes[pos[]] == UInt8('\'') @@ -467,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) @@ -573,7 +602,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 @@ -737,11 +771,65 @@ 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. + +`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}, + 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`) +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. + +`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. + +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}, @@ -749,12 +837,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} @@ -791,17 +881,36 @@ 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}, 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} @@ -848,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}, @@ -855,11 +983,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) @@ -870,7 +1004,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 @@ -1783,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/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index dc331727..2fc52e4f 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 @@ -96,6 +96,13 @@ 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 = String(strip(string(get(body, :treeFile, "")))) inline_newick = String(strip(string(get(body, :newick, "")))) @@ -270,19 +277,29 @@ end tid = String(strip(string(r[id_col]))) !haskey(row_by_id, tid) && (row_by_id[tid] = ri) end - # Map split `.jplace` tip IDs back to `(base_query_name, share)` weights. + # Map split `.jplace` tip IDs back to the share of their query's reads + # that belongs on that edge. The share is the candidate's raw + # like-weight ratio, never a ratio renormalised over the retained + # candidates: a multi-candidate placement already partitions the mass + # (its ratios sum to 1 within `PLACEMENT_WEIGHT_SUM_TOL`), while a + # single-candidate placement accepted at `PLACEMENT_SINGLE_MIN_LWR` + # carries only the mass behind that one site. Dividing by the retained + # total would turn a lone 0.6 placement into a weight of 1.0 and push + # the 0.4 the tool left unplaced onto the reported edge, inflating that + # tip in weighted UniFrac and in Faith's PD. The shortfall is reported + # as `unplaced_reads` below instead. tip_weights = Dict{String,Float64}() if jplace_str !== nothing && !isempty(strip(jplace_str)) jd_div = MetaManifold.CladeCumulus.parse_jplace(jplace_str) for q in jd_div.queries q.status === :valid && !isempty(q.candidates) || continue - 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.like_weight / tw + # `parse_jplace` has already rejected any ratio outside + # [0, 1], so this needs no further clamping. + tip_weights[tid] = cand.like_weight end end end @@ -305,6 +322,36 @@ end mat[mi, ci] = ismissing(v) ? 0.0 : max(0.0, Float64(v)) * w end end + + # 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 (ci, sc) in enumerate(scols) + tot = 0.0 + for v in df[!, sc] + ismissing(v) && continue + tot += max(0.0, Float64(v)) + end + 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) for (ci, sc) in enumerate(scols) @@ -313,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 @@ -329,6 +379,7 @@ end diversity_payload = OrderedDict{String,Any}( "samples" => scols, "mapped_tips" => length(taxa_vec), + "unplaced_reads" => unplaced_reads, "faith_pd" => pd_by_sample, "unweighted_unifrac" => [uw[i, :] for i in 1:size(uw, 1)], "weighted_unifrac" => [w_raw[i, :] for i in 1:size(w_raw, 1)], 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/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 699f1c21..1a766344 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -16,6 +16,16 @@ const DM = MetaManifold.DiversityMetrics const PH = MetaManifold.Phylogeny const CLADE_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 = "") = CLADE_SV.Oxygen.internalrequest( CLADE_SV.HTTP.Request(method, path, ["Host" => "127.0.0.1:8080", "Content-Type" => "application/json"], body); @@ -110,6 +120,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[] @@ -149,6 +194,38 @@ _cc_req(method::String, path::String, body::String = "") = @test jp.queries[4].status === :undecided @test startswith(jp.queries[4].reason, "weight_sum_out_of_tolerance") + # Two acceptance rules, not one blanket tolerance. A single-candidate + # record is the tool naming one site and leaving the rest of the mass + # unplaced, so it is judged on its own share; a multi-candidate record + # partitions the mass and must still sum to one. Widening the sum + # tolerance far enough to admit a lone 0.98 would also have waved + # through a truncated 0.96 partition, which these cases pin down. + _lwr_doc(ps) = JSON3.write(OrderedDict( + "version" => 3, + "tree" => "((Ref_1:0.1{0},Ref_2:0.2{1})95:0.3{2},Ref_3:0.4{3}){4};", + "fields" => ["edge_num", "likelihood", "like_weight_ratio", "distal_length", "pendant_length"], + "placements" => ps, + )) + _single_at(w) = CC.parse_jplace(_lwr_doc([ + OrderedDict("n" => ["Q"], "p" => [[0, -100.0, w, 0.01, 0.02]]), + ])).queries[1] + _partition(a, b) = CC.parse_jplace(_lwr_doc([ + OrderedDict("n" => ["Q"], "p" => [[0, -100.0, a, 0.01, 0.02], [1, -100.4, b, 0.01, 0.03]]), + ])).queries[1] + + @test _single_at(0.60).status === :valid + @test _single_at(0.60).max_lwr ≈ 0.60 atol = 1e-12 + @test _single_at(CC.PLACEMENT_SINGLE_MIN_LWR).status === :valid + @test _single_at(0.49).status === :undecided + @test startswith(_single_at(0.49).reason, "weight_sum_out_of_tolerance") + + @test _partition(0.70, 0.30).status === :valid + @test _partition(0.80, 0.16).status === :undecided # 0.96: truncated partition + @test _partition(0.80, 0.24).status === :undecided # 1.04: over-complete partition + # A partition stays strict even where each part alone would clear the + # single-candidate bar, so the single rule cannot leak into the many. + @test _partition(0.60, 0.60).status === :undecided + # Placed query read counts project onto the reference tree edges without losing reads, # and `lwr_threshold = 0.95` marks `Q_borderline` (LWR = 0.85) as `present_in_some` # even when `support_threshold = 0.80`. @@ -252,6 +329,8 @@ 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) + # `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 @@ -275,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[] @@ -397,7 +507,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( @@ -414,6 +536,46 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 @test v_allow2.allowed == true @test String(v_allow2.proposal_file) != String(v_allow.proposal_file) + # 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 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};", + "fields" => ["edge_num", "likelihood", "like_weight_ratio", "distal_length", "pendant_length"], + "placements" => [ + OrderedDict("n" => ["ASV_1"], "p" => [[0, -100.0, 0.60, 0.05, 0.02]]), + OrderedDict("n" => ["ASV_2"], "p" => [[1, -99.0, 1.00, 0.05, 0.02]]), + ], + ))) + r_jp = _cc_req("POST", "/api/v1/studies/StudyCC/runs/run1/clade-cumulus", JSON3.write(Dict( + "table" => "merged", + "mode" => "combined", + "treeSource" => "study_tree", + "treeFile" => "placed.jplace", + "supportThreshold" => 0.85, + "lwrThreshold" => 0.95, + ))) + @test r_jp.status == 200 + 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 ≈ 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. + @test jp_body.root.cumulative_count ≈ 110.0 atol = 1e-9 + # Canonical DuckDB database is strictly untouched byte-for-byte! @test bytes2hex(open(sha256, db_path)) == before_db_sha diff --git a/test/unit/test_routes.jl b/test/unit/test_routes.jl index 62998944..f97b9bba 100644 --- a/test/unit/test_routes.jl +++ b/test/unit/test_routes.jl @@ -932,5 +932,30 @@ _saved_root = SV.ServerState._root[] end +@testset "Registered routes" begin + # Every route file is `include`d once into `Server`, and each routing macro + # appends to `Server.ROUTES`. A file included twice registers its routes + # twice, and no request test can reveal that: the second registration + # simply shadows the first. Assert the table holds no duplicate + # method+path pair. + registered = [(m, string(p)) for (m, p, _) in SV.ROUTES] + @test !isempty(registered) + + seen = Dict{Tuple{String,String},Int}() + for k in registered + seen[k] = get(seen, k, 0) + 1 + end + # Compared against an empty vector rather than asserted with `isempty`, so + # a failure names the route that doubled up instead of printing `false`. + duplicated = sort([k for (k, n) in seen if n > 1]) + @test duplicated == Tuple{String,String}[] + + # 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 "/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. SV.ServerState._root[] = _saved_root