diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index 48df92c7..360d1121 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -968,6 +968,7 @@ export interface CladeCumulusRequest { export interface CladeCumulusDiversity { samples: string[] mapped_tips: number + 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..c728daf4 100644 --- a/frontend/src/components/CladeCumulusPanel.tsx +++ b/frontend/src/components/CladeCumulusPanel.tsx @@ -249,6 +249,14 @@ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { return div.normalized_weighted_unifrac }, [result, unifracKind]) + 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 +708,13 @@ export function CladeCumulusPanel({ study, runs }: CladeCumulusPanelProps) { ))} + {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(', ')} +

+ )} + {/* Pairwise UniFrac table */} {activeMatrix && (
diff --git a/scripts/check-julia-syntax.jl b/scripts/check-julia-syntax.jl index bc8a9d60..1ea4737a 100644 --- a/scripts/check-julia-syntax.jl +++ b/scripts/check-julia-syntax.jl @@ -20,6 +20,10 @@ for file in files issue = ex break end + if ex === :end + issue = "stray `end` at top level (missing or duplicated opener in this file)" + break + end next_position > position || error("Julia parser made no progress at code unit $position") position = next_position @@ -30,7 +34,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..ec92c0ae 100644 --- a/src/analysis/clade_cumulus.jl +++ b/src/analysis/clade_cumulus.jl @@ -28,11 +28,12 @@ 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 +# Multi-candidate placements partition probability mass across kept edges and +# must sum to 1 within floating-point tolerance; single-candidate placements +# leave unplaced mass off the tree and require at least half the mass at the +# reported site (`PLACEMENT_SINGLE_MIN_LWR`). +const PLACEMENT_WEIGHT_SUM_TOL = 1e-5 +const PLACEMENT_SINGLE_MIN_LWR = 0.5 ## Phylogenetic tree structures & Newick / jplace parsers mutable struct PhyloNode @@ -573,7 +574,10 @@ 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 + 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 @@ -1791,15 +1795,12 @@ function _reaggregate_after_move( mono = get(total_tips_per_taxon, u.dominant_taxon, 0) == 1 end active_here = u.cumulative_count > 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 + 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.branch_length, u.support, own_pd, u.dominant_taxon, u.taxon_purity, mono, + 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) @@ -1821,18 +1822,14 @@ function _reaggregate_after_move( n_sub_tips += nt if ch_active any_active = true - # The child's own branch is already inside `rk.subtree_pd`; - # adding `rk.branch_length` again would count every active - # internal edge twice. sub_pd += rk.subtree_pd end for (k, v) in ct merged_taxa[k] = get(merged_taxa, k, 0) + v end end - # Match clade_faith_pd: an active non-root clade adds the branch - # connecting it to its parent, once, on top of its descendants' values. - any_active && !is_root && (sub_pd += max(0.0, u.branch_length)) + + !is_root && any_active && (sub_pd += max(0.0, u.branch_length)) dom_tx = nothing dom_n = 0 diff --git a/src/pipeline/phylogeny.jl b/src/pipeline/phylogeny.jl index 7b0fa2af..0798491c 100644 --- a/src/pipeline/phylogeny.jl +++ b/src/pipeline/phylogeny.jl @@ -35,6 +35,7 @@ export REFERENCE_STEPS, PLACEMENT_STEPS, Step, placement_settings, run_workflow, using ..Validation import ..Tools using ..Tools: tool_bin, _sq, _run_logged, _safe_optional_args, _num + import ..Tools using ..RemoteExec: run_remote, checked_target import ..Provenance diff --git a/src/server/routes/clade_cumulus.jl b/src/server/routes/clade_cumulus.jl index dc331727..f7f81706 100644 --- a/src/server/routes/clade_cumulus.jl +++ b/src/server/routes/clade_cumulus.jl @@ -75,9 +75,7 @@ function _load_run_table_df(study::AbstractString, run::AbstractString, table::A db_path = joinpath(dir, "results.duckdb") isfile(db_path) || return (nothing, nothing, json_error(404, "no_results", "No results database for run '$run'")) - db = DuckDB.DB(db_path) - con = DBInterface.connect(db) - try + with_results_db(dir) do con cols = _duckdb_columns(con, table) isempty(cols) && return (nothing, nothing, json_error(404, "table_not_found", "Table '$table' not found in run '$run'")) @@ -85,9 +83,6 @@ function _load_run_table_df(study::AbstractString, run::AbstractString, table::A order = _stable_order_clause(cols) df = DataFrame(DBInterface.execute(con, "SELECT * FROM \"$table\" $order")) return (df, (db_path = db_path, sample_cols = scols, columns = cols), nothing) - finally - DBInterface.close!(con) - close(db) end end @@ -142,10 +137,9 @@ end dir = _require_duckdb(study, run; group) db_path = isnothing(dir) ? nothing : joinpath(dir, "results.duckdb") tables_info = OrderedDict{String,Any}[] - if !isnothing(db_path) && isfile(db_path) - db = DuckDB.DB(db_path) - con = DBInterface.connect(db) - try + if !isnothing(dir) && !isnothing(db_path) && isfile(db_path) + tables_info = with_results_db(dir) do con + available = OrderedDict{String,Any}[] # Real tables only: SHOW TABLES also lists views, and a view backed # by a session-registered DataFrame cannot be bound (let alone # counted) from this fresh connection. @@ -159,7 +153,7 @@ end nrows = Int(only(DataFrame(DBInterface.execute(con, "SELECT COUNT(*) AS n FROM \"$t\""))).n) has_seq = any(c -> c in ("Sequence", "sequence", "seq", "DNA"), cols) ranks = [r for r in MetaManifold.CladeCumulus.DEFAULT_TAXONOMIC_RANKS if r in cols] - push!(tables_info, OrderedDict{String,Any}( + push!(available, OrderedDict{String,Any}( "table" => t, "row_count" => nrows, "sample_count" => length(scols), @@ -167,9 +161,7 @@ end "ranks" => ranks, )) end - finally - DBInterface.close!(con) - close(db) + available end end @@ -270,19 +262,19 @@ 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. + # 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). 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 + tip_weights[tid] = cand.like_weight end end end @@ -305,6 +297,23 @@ end 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 + short_rows = [(ri, 1.0 - wsum) for (ri, wsum) in row_weight + if wsum < 1.0 - MetaManifold.CladeCumulus.PLACEMENT_WEIGHT_SUM_TOL] + unplaced_reads = OrderedDict{String,Float64}() + for sc in scols + col_vals = df[!, sc] + tot = 0.0 + for (ri, lost_share) in short_rows + v = col_vals[ri] + ismissing(v) && continue + tot += max(0.0, Float64(v)) * lost_share + end + unplaced_reads[sc] = tot + 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) @@ -329,6 +338,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/routes/trees.jl b/src/server/routes/trees.jl index f3cc8e07..beec713a 100644 --- a/src/server/routes/trees.jl +++ b/src/server/routes/trees.jl @@ -17,9 +17,9 @@ const _tree_view_lock = ReentrantLock() _trees_dir(study::String) = joinpath(ServerState.projects_dir(), study, "trees") _tree_view_path(study::String, file::String) = joinpath(_trees_dir(study), ".views", file * ".json") -_tree_format(file::String) = endswith(lowercase(file), ".jplace") ? "jplace" : "newick" +_tree_format(file::AbstractString) = endswith(lowercase(file), ".jplace") ? "jplace" : "newick" -function _valid_tree_file(file::String)::Bool +function _valid_tree_file(file::AbstractString)::Bool _valid_name(file) && any(ext -> endswith(lowercase(file), ext), _TREE_EXTENSIONS) 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 699f1c21..3b887de9 100644 --- a/test/unit/test_clade_cumulus.jl +++ b/test/unit/test_clade_cumulus.jl @@ -14,12 +14,12 @@ const EP = MetaManifold.Epistemic const CC = MetaManifold.CladeCumulus const DM = MetaManifold.DiversityMetrics const PH = MetaManifold.Phylogeny -const CLADE_SV = MetaManifold.Server +const _CLADE_SV = MetaManifold.Server _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); - middleware = CLADE_SV._MIDDLEWARE, + _CLADE_SV.Oxygen.internalrequest( + _CLADE_SV.HTTP.Request(method, path, ["Host" => "127.0.0.1:8080", "Content-Type" => "application/json"], body); + middleware = _CLADE_SV._MIDDLEWARE, ) @testset "CladeCumulus & Phylogenetic Integration" begin @@ -149,6 +149,28 @@ _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") + _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 _partition(0.70, 0.30).status === :valid + @test _partition(0.80, 0.16).status === :undecided + @test _partition(0.80, 0.24).status === :undecided + @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`. @@ -277,8 +299,8 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 @testset "Dual-hierarchy CladeCumulus, polyphyly detection, drag-and-drop, and HTTP routes" begin tmp = mktempdir() - old_root = CLADE_SV.ServerState._root[] - CLADE_SV.ServerState.set_root!(tmp) + old_root = _CLADE_SV.ServerState._root[] + _CLADE_SV.ServerState.set_root!(tmp) try data_dir = joinpath(tmp, "data", "StudyCC", "run1") mkpath(data_dir) @@ -310,8 +332,6 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 con = DBInterface.connect(db) DuckDB.register_data_frame(con, df, "src_df") DBInterface.execute(con, "CREATE TABLE merged AS SELECT * FROM src_df") - # Registration creates a persistent view backed by this Julia object; - # only the materialized table should remain in the on-disk fixture. DuckDB.unregister_data_frame(con, "src_df") DBInterface.close!(con) close(db) @@ -417,6 +437,28 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 # Canonical DuckDB database is strictly untouched byte-for-byte! @test bytes2hex(open(sha256, db_path)) == before_db_sha + # 6. POST combined mode against a study `.jplace` artifact + jplace_payload = JSON3.write(OrderedDict( + "version" => 3, + "tree" => "((ASV_1:0.05{0},ASV_2:0.05{1})0.96:0.20{2},ASV_3:0.08{3}){4};", + "fields" => ["edge_num", "likelihood", "like_weight_ratio", "distal_length", "pendant_length"], + "placements" => [ + OrderedDict("n" => ["ASV_4"], "p" => [[3, -10.0, 0.99, 0.02, 0.03]]), + ], + )) + write(joinpath(trees_dir, "placed.jplace"), jplace_payload) + 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", + ))) + @test r_jp.status == 200 + jp_body = JSON3.read(String(r_jp.body)) + @test jp_body.diversity !== nothing + @test jp_body.diversity.unplaced_reads.SampleA ≈ 0.05 atol = 1e-9 + @test jp_body.diversity.unplaced_reads.SampleB ≈ 0.10 atol = 1e-9 + # Path traversal attempts are rejected with 400! @test _cc_req("GET", "/api/v1/studies/..%2Fetc/runs/run1/clade-cumulus/sources").status == 400 @test _cc_req("POST", "/api/v1/studies/StudyCC/runs/run1/clade-cumulus", JSON3.write(Dict( @@ -424,7 +466,7 @@ echo "(MMTIP_000001:0.05,(MMTIP_000002:0.03,MMTIP_000003:0.04)0.910:0.12,MMTIP_0 "treeFile" => "../secret.nwk", ))).status == 400 finally - CLADE_SV.ServerState._root[] = old_root + _CLADE_SV.ServerState._root[] = old_root rm(tmp; recursive = true, force = true) end end diff --git a/test/unit/test_routes.jl b/test/unit/test_routes.jl index 62998944..aa61117d 100644 --- a/test/unit/test_routes.jl +++ b/test/unit/test_routes.jl @@ -932,5 +932,21 @@ _saved_root = SV.ServerState._root[] end +@testset "Registered routes" begin + 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 + 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) +end + # Testsets above point the server at temporary roots. SV.ServerState._root[] = _saved_root