diff --git a/frontend/src/api/events.ts b/frontend/src/api/events.ts index 3422641b..44992d87 100644 --- a/frontend/src/api/events.ts +++ b/frontend/src/api/events.ts @@ -13,12 +13,25 @@ export type SSEHandler = { export function openEventStream(handlers: SSEHandler & { onConnectedChange?: (connected: boolean) => void }): () => void { const es = new EventSource(apiUrl('/api/v1/events')) + // A malformed frame must not throw inside the listener: EventSource would + // keep the stream open but the exception surfaces as an uncaught error and + // the update is lost silently either way. Drop the bad frame instead. + const parseFrame = (data: string): unknown => { + try { + return JSON.parse(data) + } catch { + return undefined + } + } + es.addEventListener('job_update', (e: MessageEvent) => { - handlers.onJobUpdate?.(JSON.parse(e.data) as Job) + const job = parseFrame(e.data) + if (job !== undefined) handlers.onJobUpdate?.(job as Job) }) es.addEventListener('stage_update', (e: MessageEvent) => { - handlers.onStageUpdate?.(JSON.parse(e.data)) + const ev = parseFrame(e.data) + if (ev !== undefined) handlers.onStageUpdate?.(ev as Parameters>[0]) }) es.onopen = () => handlers.onConnectedChange?.(true) diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 360d1121..2039f8bc 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -968,6 +968,17 @@ 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. + * + * These reads are excluded from the tree-tip abundances, but they are still + * counted in each sample's library size, which is the denominator the + * weighted UniFrac metrics divide by. That is what keeps the shortfall off + * the tree instead of renormalising a partially placed sample back up to + * abundance 1. Absent for responses produced before this field existed. + */ unplaced_reads?: Record faith_pd: Record unweighted_unifrac: number[][] diff --git a/frontend/src/components/CladeCumulusPanel.tsx b/frontend/src/components/CladeCumulusPanel.tsx index c728daf4..5a666385 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 run's clade hierarchy, diversity and unplaced reads. Loads + * available sources and computes results on request; move checks can preview + * counterfactual trees or save sidecar proposals without changing canonical tables. + */ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { const toast = useToast() const reqCounter = useRef(0) @@ -249,6 +254,9 @@ 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 [] @@ -708,10 +716,25 @@ 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 absent from + the tree-tip abundances but they are *not* dropped from the + sample: each one still counts in its sample's library size, + which is the denominator the weighted UniFrac metrics divide + by. That is what keeps them off the tree instead of letting a + partially placed sample renormalise back up to abundance 1, so + the shortfall is stated here as well. + */} {unplacedSamples.length > 0 && (

- Unplaced reads (like-weight mass left off reference edges, excluded from weighted abundances above):{' '} - {unplacedSamples.map(([s, n]) => `${s} ${n.toFixed(2)}`).join(', ')} + Unplaced reads (carried by no tip — like-weight mass left on no edge, plus + ASVs the tree does not cover — excluded from the tree-tip abundances but + counted in each sample’s library size, so they lower its relative + abundances rather than being renormalised onto the placed tips):{' '} + {unplacedSamples + .map(([s, n]) => `${s} ${n.toFixed(2)}`) + .join(', ')}

)} diff --git a/frontend/src/tree/model.ts b/frontend/src/tree/model.ts index de557a87..646d6af5 100644 --- a/frontend/src/tree/model.ts +++ b/frontend/src/tree/model.ts @@ -112,14 +112,14 @@ export function parseNewick(text: string): SourceTree { skipSpace() if (s[i] === "'" || s[i] === '"') { const q = s[i++] - let out = '' + const parts: string[] = [] while (i < s.length) { if (s[i] === q) { - if (s[i + 1] === q) { out += q; i += 2; continue } + if (s[i + 1] === q) { parts.push(q); i += 2; continue } i++ - return out + return parts.join('') } - out += s[i++] + parts.push(s[i++]) } fail('unclosed quote') } @@ -209,7 +209,18 @@ export function parseNewick(text: string): SourceTree { } export function parseJplace(text: string): SourceTree { - const doc = JSON.parse(text) + let doc: { + tree?: unknown + fields?: string[] + placements?: Record[] + metadata?: Record + version?: unknown + } + try { + doc = JSON.parse(text) + } catch { + throw new Error('jplace file is not valid JSON') + } const tree = parseNewick(String(doc.tree)) const fields: string[] = doc.fields ?? [] const col = (f: string) => fields.indexOf(f) diff --git a/scripts/check-julia-syntax.jl b/scripts/check-julia-syntax.jl index 1ea4737a..f4a93eaf 100644 --- a/scripts/check-julia-syntax.jl +++ b/scripts/check-julia-syntax.jl @@ -20,8 +20,15 @@ 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 (missing or duplicated opener in this file)" + issue = "stray `end` at top level (a `module`/`function`/... opener " * + "is missing or was duplicated in this file)" break end next_position > position || diff --git a/src/analysis/clade_cumulus.jl b/src/analysis/clade_cumulus.jl index 34b80b43..0b00d549 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -37,6 +37,20 @@ const DEFAULT_TAXONOMIC_RANKS = [ # A sum over 1 is never legitimate: weights are mutually exclusive placement # ratios, so overweight evidence is malformed and is rejected by a tight upper # bound while the lower shortfall is tolerated. +# +# A single-candidate record is the shape that made the lower tolerance +# necessary, and it is judged on its own share instead +# (`PLACEMENT_SINGLE_MIN_LWR`): the tool named one site and left the rest of +# the mass unplaced. +# +# 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_WEIGHT_SUM_TOL = 0.05 const PLACEMENT_WEIGHT_OVER_TOL = 1e-9 const PLACEMENT_SINGLE_MIN_LWR = 0.5 @@ -576,9 +590,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('\'') @@ -835,13 +851,21 @@ 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 requires a weight + of at least `0.5`; multiple candidates must sum to 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. + +Candidate validation failures produce `:undecided` queries with a reason and +`max_lwr = nothing`; a lone weight below `0.5` gives the reason +`single_candidate_below_min_lwr ()`. Missing required fields, invalid +query names or multiplicities, duplicate edge numbers, and missing or duplicate +required field columns throw `ArgumentError`. JSON decoding, Newick parsing, +and value conversion errors propagate. """ function parse_jplace(raw_json::AbstractString)::JplaceDocument doc = JSON3.read(raw_json) @@ -942,8 +966,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 length(cands) == 1 && cands[1].like_weight < PLACEMENT_SINGLE_MIN_LWR + # One candidate and less than half the mass supports it: the rest of + # the mass is unexplained, so the placement is not decisive. The + # single-candidate sum is in tolerance by construction, so this + # branch reports the share it actually has, not a weight-sum error. push!(queries, QueryPlacement(q_ord, q_names, cands, :undecided, - "weight_sum_out_of_tolerance ($(round(weight_sum; digits=6)))", nothing)) + "single_candidate_below_min_lwr ($(round(cands[1].like_weight; digits=6)))", nothing)) elseif length(cands) > 1 && (weight_sum > 1.0 + PLACEMENT_WEIGHT_OVER_TOL || 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)) @@ -1108,11 +1136,63 @@ 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, throw `ArgumentError` +if the declared total is non-finite, negative, or below +`sum(values(abund)) * (1 - 1e-9) - 1e-9`. A smaller total within that tolerance +is raised to the on-tree sum. `name` identifies the total in error messages. +""" +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. + +Raw distances have the tree's branch-length units. Normalised distances divide +by the abundance-weighted sum of root-to-tip lengths, returning `0.0` when that +denominator is zero. A sample with zero total contributes zero abundance. + +Throw `ArgumentError` for counts that are non-finite, negative, or differ in +length from `taxa`, or for taxa absent from the tree. Duplicate taxa have their +counts added. Invalid totals throw `ArgumentError` as described by +`_sample_total`; errors parsing a Newick tree also propagate. Trees with more +than two root children are midpoint-rooted before calculating the distance. """ function weighted_unifrac( counts_a::AbstractVector{<:Real}, @@ -1120,12 +1200,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} @@ -1162,17 +1244,35 @@ 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 the positive contributions of 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 abundance difference). +Their contributions sum to `weighted_unifrac(...; normalized=false)` with the +same totals; return an empty vector when no branch contributes. + +Counts, taxa, tree rooting, and `total_a` / `total_b` follow `weighted_unifrac`, +including its `ArgumentError` validation and propagated Newick parsing errors. +""" 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} @@ -1219,6 +1319,22 @@ 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 matrix in `samples` order with a zero +diagonal. Rows of `sample_matrix` correspond to `taxa` and columns to `samples`. +Use `metric=:weighted` or `:unweighted`; `normalized` affects only the weighted +metric. Optional `totals` supplies library sizes in `samples` order using +`weighted_unifrac`'s denominator rules. The unweighted metric ignores their +values, but their length must still match `samples`. + +Throw `ArgumentError` for inconsistent matrix dimensions or totals length. +Pairwise metric validation errors (including unsupported metrics, invalid +counts, missing taxa, or invalid weighted totals) and tree parsing errors +propagate. With fewer than two samples, no pairwise validation is performed. +""" function unifrac_matrix( sample_matrix::AbstractMatrix{<:Real}, taxa::AbstractVector{<:AbstractString}, @@ -1226,11 +1342,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) @@ -1241,7 +1363,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 @@ -2110,6 +2234,20 @@ function _is_descendant(candidate_desc::AbstractString, ancestor::CladeNode)::Bo false end +""" + _reaggregate_after_move(root, moving, target_id, case_obj, total_reads) -> CladeNode + +Return a rebuilt tree with `moving` detached from its old parent and attached +beneath `target_id`, without mutating the input nodes. The caller must validate +the move first, as `validate_drag_drop` does. Recompute counts, residuals, +subtree PD and taxonomic summaries for the proposed topology; epistemic statuses +come from the supplied `case_obj`. Frequencies use `total_reads` as their +denominator, or zero when it is non-positive. A clade left childless carries +only its direct reads and has zero residual and PD when those reads are zero. + +Throw `ErrorException` if `moving` is the root. Invalid residuals propagate +`ArgumentError` from `Epistemic.cloud_size_by_residual`. +""" function _reaggregate_after_move( root::CladeNode, moving::CladeNode, @@ -2141,7 +2279,11 @@ function _reaggregate_after_move( # and purity are exact everywhere on the proposed tree. total_tips_per_taxon = Dict{String,Int}() function count_leaf_taxa!(u::CladeNode) - if isempty(u.children) + # A clade emptied by the graft is childless but is not a tree tip: + # counting it would inflate tip counts, purity and monophyly with a + # taxon whose reads have all moved elsewhere. Terminal `tip`/`ASV` + # nodes are real tips even when they carry no reads in this table. + if isempty(u.children) && (u.direct_count > 0.0 || u.rank in ("tip", "ASV")) u.dominant_taxon !== nothing && (total_tips_per_taxon[u.dominant_taxon] = get(total_tips_per_taxon, u.dominant_taxon, 0) + 1) else @@ -2154,23 +2296,41 @@ 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 + # The graft can leave a clade with no children: its only child was + # the clade that moved out. Its `cumulative_count` still counts the + # subtree that has just left, so trusting it would keep the emptied + # clade "active" and give it a nonzero `subtree_pd` and frequency + # for reads it no longer holds. A node with no children carries + # exactly its own reads, so recompute the count from + # `direct_count`, as the internal-node branch below recomputes its + # own from the rebuilt children. + own_count = u.direct_count + # `residual_count` comes from the same departed subtree: placement + # and support uncertainty belong to the reads that left with it, + # and they are recomputed at its new location. A clade left holding + # no reads has no residual uncertainty left to report, so the + # residual is derived from the same count as everything else here. + own_residual = own_count > 0.0 ? u.residual_count : 0.0 + freq = total_reads > 0.0 ? own_count / total_reads : 0.0 + # Same tip eligibility as `count_leaf_taxa!` above: an emptied + # internal clade contributes no taxon tip to its parent. + counts_as_tip = own_count > 0.0 || u.rank in ("tip", "ASV") tmap = Dict{String,Int}() mono = u.monophyletic - if u.dominant_taxon !== nothing + if counts_as_tip && 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 = own_count > 0.0 sub_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.residual_count, Epistemic.cloud_size_by_residual(u.residual_count), + u.id, u.name, u.rank, u.direct_count, own_count, freq, + own_residual, Epistemic.cloud_size_by_residual(own_residual), u.branch_length, u.support, sub_pd, u.dominant_taxon, u.taxon_purity, mono, st, Epistemic.epistemic_colour(st), CladeNode[], ) - return (leaf_node, tmap, 1, active_here) + return (leaf_node, tmap, counts_as_tip ? 1 : 0, active_here) end rebuilt_kids = CladeNode[] diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index d864b16f..c65a5a91 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 @@ -91,6 +91,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, "")))) @@ -198,6 +205,19 @@ end end end +""" + _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) + +Build a clade explorer JSON response from a run's results table and the requested +tree source, including a provenance receipt, without modifying canonical tables. +Diversity is `nothing` unless a phylogeny, sample columns and matching tips are +available. Placed `.jplace` abundances use raw like-weight ratios; unplaced reads +include unmapped ASVs and remain in the weighted UniFrac library-size denominators. + +Return HTTP error responses for request validation, missing inputs and failures +building the clade tree (the latter become HTTP 400). Other exceptions, including +diversity calculation failures, propagate to the route wrapper's HTTP 500 handler. +""" function _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) group = _req_group(req) err = _validate_run_request(study, run) @@ -271,9 +291,17 @@ function _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) tid = String(strip(string(r[id_col]))) !haskey(row_by_id, tid) && (row_by_id[tid] = ri) end - # Weight each placed `.jplace` tip by its raw like-weight ratio so any - # unplaced mass remains off the tree rather than renormalised onto the - # kept edge(s). + # 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) @@ -283,6 +311,8 @@ function _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) for (qname, _) in q.names for (ci, cand) in enumerate(q.candidates) tid = multi ? string(qname, "@edge", cand.edge_num, "#", ci) : qname + # `parse_jplace` has already rejected any ratio outside + # [0, 1], so this needs no further clamping. tip_weights[tid] = cand.like_weight end end @@ -306,29 +336,35 @@ function _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) mat[mi, ci] = ismissing(v) ? 0.0 : max(0.0, Float64(v)) * w end end - row_weight = Dict{Int,Float64}() - for (ri, w) in taxa_rows - row_weight[ri] = get(row_weight, ri, 0.0) + w - end - # `unplaced_reads` reports how much of each row's abundance never - # reached the tree, so its floor is floating-point noise, not the - # *validity* tolerance. `PLACEMENT_WEIGHT_SUM_TOL` (0.05) exists so - # that a candidate list truncated by the placement tool stays - # `:valid`; reusing it here reported 0 unplaced reads for a - # placement at LWR 0.99 whose remaining 1% really is off the tree. - noise_floor = MetaManifold.CladeCumulus.PLACEMENT_WEIGHT_OVER_TOL - short_rows = [(ri, 1.0 - wsum) for (ri, wsum) in row_weight - if 1.0 - wsum > noise_floor] + + # 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) @@ -338,14 +374,17 @@ function _clade_cumulus_post_impl(req::HTTP.Request, study::String, run::String) 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 10f098ba..d66fa476 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -16,6 +16,55 @@ 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 + +"""Depth-first lookup by `id` in a `clade_node_to_dict` payload (an `OrderedDict`).""" +function _cc_find_dict(node, id::AbstractString) + String(node["id"]) == id && return node + for ch in node["children"] + hit = _cc_find_dict(ch, id) + hit === nothing || return hit + end + nothing +end + +"""Every clade node below `root`, in depth-first order.""" +function _cc_nodes(root) + out = Any[root] + i = 1 + while i <= length(out) + append!(out, out[i].children) + i += 1 + end + out +end + +"""The first non-root clade holding exactly one child, or `nothing`.""" +function _cc_single_child_parent(root) + for u in _cc_nodes(root) + (u !== root && length(u.children) == 1) && return u + end + nothing +end + +"""The internal clade that directly contains the tip labelled `tipname`.""" +function _cc_parent_of_tip(root, tipname::AbstractString) + for u in _cc_nodes(root) + for ch in u.children + (isempty(ch.children) && startswith(String(ch.name), tipname)) && return u + end + 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); @@ -224,6 +273,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[] @@ -261,8 +345,13 @@ _cc_req(method::String, path::String, body::String = "") = @test jp.queries[3].status === :valid @test jp.queries[3].max_lwr ≈ 0.98 atol = 1e-12 @test jp.queries[4].status === :undecided - @test startswith(jp.queries[4].reason, "weight_sum_out_of_tolerance") + @test startswith(jp.queries[4].reason, "single_candidate_below_min_lwr") + # 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, whose shortfall may be the placement tool's + # own truncation but whose excess never is. _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};", @@ -280,9 +369,11 @@ _cc_req(method::String, path::String, body::String = "") = @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, "single_candidate_below_min_lwr") + @test _partition(0.70, 0.30).status === :valid - @test _partition(0.80, 0.16).status === :valid - @test _partition(0.80, 0.24).status === :undecided + @test _partition(0.80, 0.16).status === :valid # 0.96: a truncated candidate list + @test _partition(0.80, 0.24).status === :undecided # 1.04: overweight evidence @test _partition(0.60, 0.60).status === :undecided # Placed query read counts project onto the reference tree edges without losing reads, @@ -388,6 +479,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 @@ -566,6 +659,48 @@ 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 + @test jp_body.diversity.weighted_unifrac[1][2] ≈ + 0.15636363636363636 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 @@ -603,4 +738,81 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 end end + @testset "A move that empties a clade clears its activity and PD" begin + # `_reaggregate_after_move` detaches the moving clade from its old + # parent. When that parent's only child was the one that moved, the + # parent is left with no children but still carries the + # `cumulative_count` (and the `residual_count` behind its uncertainty + # cloud) of the subtree that has just left. Read as a leaf, that stale + # count would keep the emptied clade active and give it a nonzero + # `subtree_pd` and frequency for reads it no longer holds. + # Two domains, two ASVs each, so the proposed root's taxon purity + # distinguishes "four real tips" (0.5) from "four tips plus the + # emptied wrapper counted as a fifth Bacteria tip" (0.6). + df = DataFrame( + SeqName = ["ASV_1", "ASV_2", "ASV_3", "ASV_4"], + Domain = ["Bacteria", "Bacteria", "Archaea", "Archaea"], + SampleA = [10.0, 20.0, 30.0, 40.0], + ) + # `wrapper` has exactly one child (the ASV_1/ASV_2 clade); the other + # side of the root is a two-tip clade we can move that child into. + nwk = "(((ASV_1:0.10,ASV_2:0.20)0.96:0.30)0.94:0.05,(ASV_3:0.15,ASV_4:0.25)0.90:0.40)Root;" + cc_move = CC.build_clade_tree( + df; mode = :combined, newick = nwk, support_threshold = 0.5, lwr_threshold = 0.95, + ) + + wrapper = _cc_single_child_parent(cc_move.root) + @test wrapper !== nothing + only_child = wrapper === nothing ? nothing : wrapper.children[1] + @test only_child !== nothing + + # The target is the clade directly holding ASV_3: internal, well + # supported, and in the other half of the tree from the moving clade. + target = _cc_parent_of_tip(cc_move.root, "ASV_3") + @test target !== nothing + + if wrapper !== nothing && only_child !== nothing && target !== nothing + before_total = Float64(cc_move.root.cumulative_count) + vacated_reads = Float64(wrapper.cumulative_count) + vacated_residual = Float64(wrapper.residual_count) + moved_reads = Float64(only_child.cumulative_count) + moved_residual = Float64(only_child.residual_count) + @test vacated_reads > 0.0 + @test vacated_residual > 0.0 # so clearing it is observable + @test moved_reads > 0.0 + + verdict = CC.validate_drag_drop(cc_move, String(only_child.id), String(target.id)) + @test verdict["allowed"] == true + proposed_root = verdict["proposed_root"] + + vacated = _cc_find_dict(proposed_root, String(wrapper.id)) + @test vacated !== nothing + if vacated !== nothing + @test Float64(vacated["cumulative_count"]) == 0.0 + @test Float64(vacated["cumulative_frequency"]) == 0.0 + @test Float64(vacated["subtree_pd"]) == 0.0 + # The residual uncertainty belonged to the reads that left, so + # it leaves with them instead of lingering on an empty clade. + @test Float64(vacated["residual_count"]) == 0.0 + @test isempty(vacated["children"]) + end + + # The reads are relocated, not lost, and the clade that moved keeps + # every one of them. + relocated = _cc_find_dict(proposed_root, String(only_child.id)) + @test relocated !== nothing + if relocated !== nothing + @test Float64(relocated["cumulative_count"]) == moved_reads + @test Float64(relocated["residual_count"]) == moved_residual + end + @test Float64(proposed_root["cumulative_count"]) ≈ before_total atol = 1e-9 + # The emptied wrapper is childless but is not a tip: it must not + # add its departed subtree's taxon to the tip counts behind + # purity/monophyly. Four real tips, two per domain, give 0.5; + # counting the wrapper as a fifth Bacteria tip would give 0.6. + @test Float64(proposed_root["taxon_purity"]) ≈ 0.5 atol = 1e-9 + end + end + + end diff --git a/test/unit/test_routes.jl b/test/unit/test_routes.jl index aa61117d..e6ecf069 100644 --- a/test/unit/test_routes.jl +++ b/test/unit/test_routes.jl @@ -933,6 +933,11 @@ _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) @@ -940,12 +945,20 @@ end 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}[] - 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) + # The Cladistics tab and the run-level CladeCumulus evidence view are + # served by two different route files; both must be registered. These are + # whole method+path pairs, not substrings: `occursin` would also accept a + # misregistered `/api/v1/studies/{study}/cladistics-legacy` or a longer + # nested path, which is exactly the breakage this guards against. + @test ("GET", "/api/v1/studies/{study}/cladistics") in registered + @test ("GET", "/api/v1/studies/{study}/runs/{run}/clade-cumulus/sources") in registered + @test ("POST", "/api/v1/studies/{study}/runs/{run}/clade-cumulus") in registered + @test ("POST", "/api/v1/studies/{study}/runs/{run}/clade-cumulus/validate-move") in registered end # Testsets above point the server at temporary roots.