diff --git a/CHANGELOG.md b/CHANGELOG.md index 843a433e..93f1aaa2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,21 @@ types, tests, infrastructure, and alignment. - CI: pipeline extended to licence-header, formatting, lint, and commit convention checks ahead of typecheck/test/bench/build. +### Fixed + +- **Zero-depth samples poisoned the transform stage.** A sample with no reads + reached the transform, where every transform divides by a sample total: + `relative` wrote `0.0` for every feature (a value where the answer is + *undefined*), `rarefy` took `min_lib = minimum(lib_sizes) = 0` and scaled + **every** sample by zero, and `clr` took `log(0) = -Inf` that the NaN/Inf + healing later replaced with `epsilon`. All-zero samples are now healed + *before* the transform rather than after, which also removes the + inconsistency where `prepared` and `filtered_counts` described different + data under `drop_policy=impute`. Visible results change only for `impute` + runs over inputs containing a zero-depth sample; `drop` and `refuse` are + unchanged. Owner-authorised behaviour change, recorded in + `docs/statistics/behaviour-change-zero-depth-samples.md`. + ## [0.1.0] — 2026-05-21 (upstream baseline) Initial public state of the application as inherited from upstream diff --git a/docs/statistics/behaviour-change-zero-depth-samples.md b/docs/statistics/behaviour-change-zero-depth-samples.md new file mode 100644 index 00000000..78e651b2 --- /dev/null +++ b/docs/statistics/behaviour-change-zero-depth-samples.md @@ -0,0 +1,90 @@ + +# Behaviour change — zero-depth samples in the relative-abundance transform + +**Status: decided and implemented on 2026-09-22.** This note exists because the change +alters results that saved analyses were computed from, which is not a thing to do quietly. + +## What changed + +A sample with **no reads at all** (a *zero-depth* sample: every count zero) used to reach +the transform stage inside `prepare_analysis_table`. Every transform divides by a sample +total, so such a sample poisoned whichever transform ran. All-zero samples are now healed +**before** the transform instead of after it. + +## Why it had to change + +The pipeline wrote a value where the truth is that no value exists: + +| Transform | What a zero-depth sample produced | +| --- | --- | +| `relative` | `0.0` for every feature — the claim *"this feature's relative abundance is exactly 0"* about a sample whose relative abundances do not exist. `0/0` is undefined; `0.0` is a value. | +| `rarefy` | `min_lib = minimum(lib_sizes)` is `0`, so **every** sample was scaled by `0/lib`. The result was a prepared table of zeros for all samples, reported as success. | +| `clr` | `log(0) = -Inf`, which centring turns into `NaN`, which the later NaN/Inf healing replaced with `epsilon` — a wrong number presented as a healed one. | + +Healing afterwards could not repair any of it. Under `drop_policy=drop` the affected +columns were discarded, but only after the damage was computed; under +`drop_policy=impute` the imputed counts were assigned while `prepared` kept the values +computed from the empty column, so **`prepared` and `filtered_counts` described different +data** — the one inconsistency a pipeline must never have. + +This is a case of the principle stated in `method-conditions/exact-descriptive-summaries.md`: +reporting `0` for a zero-depth sample's relative abundance is a lie told with a straight +face, and no downstream check can catch it, because `0.0` is an ordinary number. + +## Who is affected + +Runs whose input contains a zero-depth sample can have different prepared +results when healing moves before the transform: + +| `drop_policy` | Behaviour before | Behaviour now | +| --- | --- | --- | +| `drop` | `rarefy` included the zero-depth sample, so retained columns were scaled to zero before the sample was dropped | sample dropped before `rarefy` — **prepared result changes for `rarefy`** | +| `refuse` | refused | refused, earlier — **same refusal** | +| `impute` | counts were transformed before healing, so `prepared` could retain values for the unhealed zero-depth column | counts imputed, transform sees the imputed counts — **prepared output can change** | + +For a zero-depth input, `drop_policy=impute` can change prepared output for +`none`, `relative`, `tss`, `css`, `rss`, `size_factors`, `presence_absence`, or +`rarefy`. `drop_policy=drop` changes prepared output for `rarefy`. Successful +`clr` and `ilr` runs using pseudocount replacement retain equal log-ratio values +for the zero-depth column. Analyses whose inputs had no zero-depth sample are +unchanged — the healing step returns immediately when +there is nothing to heal. + +## What to do + +Re-run the affected analysis. The previous numbers for a zero-depth sample were not +measurements of a quantity that exists, so there is no correction to apply: the +computation itself was refusing to say "undefined" and saying "zero" instead. + +If you are unsure whether an analysis is affected, check its diagnostics for the +recorded all-zero-sample healing and its `drop_policy`. Re-run `rarefy` analyses +with a `Dropped ... all-zero samples` entry (`drop_policy=drop`), and re-run +`none`, `relative`, `tss`, `css`, `rss`, `size_factors`, `presence_absence`, or +`rarefy` analyses with an `Imputed ... all-zero samples` entry +(`drop_policy=impute`). Successful pseudocount-based `clr` and `ilr` runs retain +equal log-ratio values, as recorded in the affected method/policy combinations +above. + +## What deliberately did NOT change + +- **All-zero *taxa* are still healed after the transform.** The asymmetry is intentional. + A zero-count taxon in a sample that *has* reads has a relative abundance of exactly + `0.0` — that value is true, so there is no lie to remove — while moving the drop earlier + would change the geometric mean that CLR centres on, altering results for analyses that + were already correct. Fixing a defect is not a licence to change the numbers around it. +- **The legacy path is otherwise byte-identical.** The exact-precision layer + (`exact_summaries.jl`) is a separate, opt-in module that nothing in the pipeline calls; + this change does not touch it. +- **No method is switched silently.** The healing is recorded in `healings`, surfaced in + diagnostics and the manifest, and raises the DANGER self-healing banner as before. + +## Provenance + +The repository owner (`@hyperpolymath`) authorised this change on 2026-09-22, choosing +"fix it now, change the results" over leaving the defect documented or hiding it behind an +opt-in flag, with this note as the accompanying migration record. This is an explicit, +narrow exception to the standing instruction that existing pipeline algorithms and +previously saved analyses not be silently changed: the change is not silent, and the +results it changes were wrong. diff --git a/src/analysis/Execution.jl b/src/analysis/Execution.jl index 58624e70..4c835d40 100644 --- a/src/analysis/Execution.jl +++ b/src/analysis/Execution.jl @@ -636,13 +636,18 @@ function heal_all_zero_samples(counts::Matrix{Float64}, sample_ids::Vector{Strin end end -function heal_all_zero_taxa(counts::Matrix{Float64}, taxa_ids::Vector{String}, drop_policy::DropPolicy, epsilon::Float64) - check = check_all_zero_taxa(counts) - if !check["has_all_zero_taxa"] +function heal_all_zero_taxa( + counts::Matrix{Float64}, + taxa_ids::Vector{String}, + drop_policy::DropPolicy, + epsilon::Float64; + indices::Union{Nothing,Vector{Int}}=nothing +) + indices = isnothing(indices) ? check_all_zero_taxa(counts)["all_zero_taxa_indices"] : indices + if isempty(indices) return (counts, taxa_ids, String[], Int[]) end - indices = check["all_zero_taxa_indices"] healings = String[] if drop_policy == DROP @@ -903,6 +908,59 @@ function prepare_analysis_table( effective_pseudocount = config.normalization.method in ("clr", "ilr") ? pseudocount : adv_pseudocount effective_epsilon = epsilon # could also use adv_epsilon, but use normalization epsilon for now + # ---------------------------------------------------------------------- + # All-zero SAMPLES are healed here — before the transform, not after + # ---------------------------------------------------------------------- + # This was the defect. Healing used to run after the transform, so `prepared` + # was computed from counts that still contained zero-depth samples, and every + # transform below divides by a sample total: + # + # * relative wrote `0.0` into each feature — the statement "this feature's + # relative abundance is exactly 0" about a sample whose relative + # abundances do not exist at all. 0/0 is undefined; 0 is a value. + # * rarefy took `min_lib = minimum(lib_sizes)`, which is 0 when any sample + # is empty, and then scaled EVERY sample by 0/lib — emptying the + # whole prepared table while reporting success. + # * clr took log(0) = -Inf, which centring turns into NaN, which the + # NaN/Inf healing later replaced with epsilon: a wrong number + # presented as a healed one. + # + # Healing afterwards could not repair any of this. Under drop_policy=drop the + # affected columns were discarded, but only after the damage was computed; + # under drop_policy=impute the imputed counts were assigned while `prepared` + # kept the values computed from the empty column, so `prepared` and + # `filtered_counts` described different data — the exact inconsistency a + # pipeline must never have. + # + # Healing first removes the class rather than the instance: past this point no + # zero-depth column can reach a transform, under any policy. + # + # Visible results change only where they were wrong. Under drop_policy=drop the + # sample is still dropped; under refuse the same refusal is raised, earlier; + # under impute the sample is imputed and the transform now sees the imputed + # counts, so its relative abundances are a uniform distribution rather than a + # column of 0.0. + # + # All-zero TAXA are deliberately still healed after the transform, and the + # asymmetry is intentional. A zero-count taxon in a sample that has reads has a + # relative abundance of exactly 0.0 — that value is true — so there is no lie + # to remove, while moving the drop earlier would change the geometric mean CLR + # centres on and therefore change results for analyses that were already + # correct. Fixing a defect is not licence to change the numbers around it. + healings = String[] + is_dangerous_diag = false + # Preserve the original all-zero taxa before imputing zero-depth samples. + # Otherwise an epsilon inserted for an empty sample makes an all-zero taxon + # appear nonzero when diagnostics and taxon healing run after the transform. + all_zero_taxa_indices = check_all_zero_taxa(filtered_counts)["all_zero_taxa_indices"] + zero_depth = check_all_zero_samples(filtered_counts) + if zero_depth["has_all_zero_samples"] + (filtered_counts, filtered_sample_ids, sample_healings, _) = + heal_all_zero_samples(filtered_counts, filtered_sample_ids, dp, effective_epsilon) + append!(healings, sample_healings) + is_dangerous_diag = true + end + # Apply zero policy counts_after_zero = copy(filtered_counts) @@ -1003,7 +1061,18 @@ function prepare_analysis_table( if lib_sizes[j] > 0 prepared[:, j] = counts_after_zero[:, j] ./ lib_sizes[j] else - prepared[:, j] .= 0.0 + # Unreachable since all-zero samples are healed before the transform, + # and deliberately a refusal rather than an assignment. This branch + # used to write `prepared[:, j] .= 0.0`, which states that every + # feature has a relative abundance of exactly 0 -- a value -- for a + # sample whose relative abundances do not exist. 0/0 is undefined, + # and "0" and "undefined" are different claims; the old value was + # wrong in a way no downstream check could detect, because 0.0 is a + # perfectly ordinary number. If a future path reintroduces a + # zero-depth column here, that is a bug to surface, not to paper + # over with a plausible-looking number. + sample_name = j <= length(filtered_sample_ids) ? filtered_sample_ids[j] : string(j) + throw(ArgumentError("INTERNAL: sample '$sample_name' (column $j) has zero total counts at the relative-abundance transform. All-zero samples are healed before this point, so reaching here means the healing was bypassed — the relative abundances of an empty sample are undefined and must not be reported as 0. Please report this with the config id $(config.id).")) end end # Offset for NB_GLM if needed? For relative, no offset, but for TSS offset we would use log(lib_sizes) @@ -1099,6 +1168,12 @@ function prepare_analysis_table( # Rarefy to min library size lib_sizes = vec(sum(counts_after_zero, dims=1)) min_lib = minimum(lib_sizes) + # A single zero-depth sample would make min_lib 0, and every sample would then + # be scaled by 0/lib — emptying the entire prepared table while reporting + # success. All-zero samples are healed before this point; this refusal exists + # so that if that guarantee is ever broken the failure is loud rather than a + # matrix of zeros that still passes every shape check downstream. + min_lib <= 0 && throw(ArgumentError("INTERNAL: rarefaction target is $(min_lib) — at least one sample has zero total counts. All-zero samples are healed before this point, so reaching here means the healing was bypassed. Rarefying to an empty sample would set every value to 0. Please report this with the config id $(config.id).")) # For stub, rarefy by subsampling proportionally to min_lib (not exact, just scaling) for j in 1:size(counts_after_zero, 2) if lib_sizes[j] > 0 @@ -1114,9 +1189,15 @@ function prepare_analysis_table( # ---------------------------------------------------------------------- (checks, warnings, errors) = self_diagnostics(filtered_counts, prepared, config; sample_metadata=sample_metadata, raw_counts=counts) + checks["all_zero_taxa"] = OrderedDict{String,Any}( + "all_zero_taxa_indices" => all_zero_taxa_indices, + "all_zero_taxa_count" => length(all_zero_taxa_indices), + "has_all_zero_taxa" => !isempty(all_zero_taxa_indices) + ) - healings = String[] - is_dangerous_diag = false + # `healings` and `is_dangerous_diag` are initialised earlier, before the + # all-zero sample healing that happens before the transform; re-initialising + # them here would silently discard those entries. banner = nothing # If errors present, hard-stop with DANGER banner @@ -1152,34 +1233,23 @@ function prepare_analysis_table( is_dangerous_diag = true end - # Heal all-zero samples/taxa based on drop_policy - # For counts (not prepared), check all-zero - if checks["all_zero_samples"]["has_all_zero_samples"] - try - (healed_counts, healed_sample_ids, heal_s, _) = heal_all_zero_samples(filtered_counts, filtered_sample_ids, dp, effective_epsilon) - # Need to also heal prepared table accordingly — for simplicity, we re-apply transform after healing counts? - # For stub, we just record healing and drop corresponding columns from prepared if DROP - if dp == DROP - # Drop columns from prepared - keep_cols = [j for j in 1:size(prepared,2) if !(j in checks["all_zero_samples"]["all_zero_samples_indices"])] - prepared = prepared[:, keep_cols] - filtered_sample_ids = healed_sample_ids - filtered_counts = healed_counts - else - # Impute: keep prepared but healing already applied to counts, need to re-prepare? For stub, just keep - filtered_counts = healed_counts - end - append!(healings, heal_s) - is_dangerous_diag = true - catch e - # REFUSE policy throws - throw(e) - end - end - + # All-zero SAMPLES were healed before the transform, above. The block that used to + # stand here healed them afterwards, and in doing so had to guess how to reconcile + # `prepared` with the healed counts -- the code said so: "for simplicity, we + # re-apply transform after healing counts?" and, for the impute case, "for stub, + # just keep". That guess is why a zero-depth sample could keep a column of 0.0 + # while its counts were imputed. Nothing is removed by deleting it: the guard + # below cannot fire, because after healing there are no all-zero samples left for + # `check_all_zero_samples` to find. if checks["all_zero_taxa"]["has_all_zero_taxa"] try - (healed_counts, healed_taxa_ids, heal_t, _) = heal_all_zero_taxa(filtered_counts, filtered_taxa_ids, dp, effective_epsilon) + (healed_counts, healed_taxa_ids, heal_t, _) = heal_all_zero_taxa( + filtered_counts, + filtered_taxa_ids, + dp, + effective_epsilon; + indices=all_zero_taxa_indices + ) if dp == DROP keep_rows = [i for i in 1:size(prepared,1) if !(i in checks["all_zero_taxa"]["all_zero_taxa_indices"])] # For ILR, prepared has n-1 rows, so need to handle differently — for stub, skip if ILR diff --git a/test/unit/test_execution.jl b/test/unit/test_execution.jl index d694cb15..fc83c4e2 100644 --- a/test/unit/test_execution.jl +++ b/test/unit/test_execution.jl @@ -217,6 +217,23 @@ @test size(prepared_drop_azt, 1) == 2 # dropped t2 @test "t2" ∉ tids_drop_azt + # Snapshot all-zero taxa before zero-depth sample imputation so the + # imputed epsilon cannot hide a taxon that must still be healed. + counts_with_all_zero_taxon_and_sample = [ + 10.0 20.0 0.0 30.0 25.0; + 0.0 0.0 0.0 0.0 0.0; + 30.0 20.0 0.0 10.0 15.0 + ] + (_, diag_imp_azt, _, _, _) = Execution.prepare_analysis_table( + config, counts_with_all_zero_taxon_and_sample; + sample_ids=sample_ids_az, + taxa_ids=taxa_ids_az, + drop_policy="impute", + impute_policy="epsilon" + ) + @test diag_imp_azt.checks["all_zero_taxa"]["all_zero_taxa_indices"] == [2] + @test any(h -> occursin("Imputed 1 all-zero taxa", h), diag_imp_azt.healings) + # Test CLR with pseudocount=1 and epsilon=1e-6 specifically # pseudocount=1 should be added to all counts for CLR, so 0->1, 1->2, etc. # Then log(1)=0, log(2)=0.693, etc., CLR centered @@ -233,6 +250,121 @@ @test !any(isnan, prepared_simple) end + @testset "zero-depth samples are healed before the transform, not after" begin + # The defect these tests pin, in the owner's words: "Execution.jl still writes + # a zero-depth sample's relative abundance as 0.0 — the exact thing item 1's + # conditions forbid." A sample with no reads at all has NO relative + # abundances: 0/0 is undefined, while 0.0 is a value, and reporting the second + # as the first is a claim the data does not support. No downstream check could + # catch it, because 0.0 is a perfectly ordinary number. + # + # The cause was ordering. `prepared` was computed from counts that still + # contained such samples, and every transform divides by a sample total, so + # one empty sample poisoned whichever transform ran: relative wrote 0.0 into + # every feature; rarefy took min_lib = minimum(lib_sizes) = 0 and scaled EVERY + # sample by 0, emptying the whole prepared table while reporting success; clr + # took log(0) = -Inf, which centring turns into NaN, which the later NaN/Inf + # healing then replaced with epsilon — a wrong number presented as a healed + # one. Healing afterwards could not repair any of it: under drop_policy=impute + # the imputed counts were assigned while `prepared` kept the values computed + # from the empty column, so `prepared` and `filtered_counts` described + # different data. + # + # Only drop_policy=impute changes its visible result, and only where it was + # wrong: the sample is imputed and the transform now sees the imputed counts, + # so its proportions are uniform rather than a column of 0.0. The drop and + # refuse results are unchanged, which is why the assertion that carries the + # weight here is the impute one. + norm = AnalysisConfig.NormalizationConfig(method="relative", pseudocount=1.0, + epsilon=1e-6, zero_policy="pseudocount") + adv = AnalysisConfig.AdvancedConfig(pseudocount=1.0, epsilon=1e-6, + min_prevalence=0.0, min_abundance=0.0, + min_samples_per_group=2) + config = AnalysisConfig.AnalysisConfig( + method="nb_glm", + formula="~ group", + metadata_columns=["group"], + normalization=norm, + advanced=adv, + created_by="test_zero_depth" + ) + + # s3 has no reads at all; the other four samples do. + counts_zero_depth = [10.0 20.0 0.0 30.0 25.0; + 30.0 20.0 0.0 10.0 15.0] + sids_zd = ["s1", "s2", "s3", "s4", "s5"] + tids_zd = ["t1", "t2"] + + @testset "impute — the imputed sample's proportions are uniform, not 0.0" begin + (prepared_zd, diag_zd, _, kept_zd, _) = Execution.prepare_analysis_table( + config, counts_zero_depth; + sample_ids=sids_zd, taxa_ids=tids_zd, + drop_policy="impute", impute_policy="epsilon" + ) + @test size(prepared_zd, 2) == 5 + @test "s3" in kept_zd + @test any(h -> occursin("Imputed", h), diag_zd.healings) + + # The assertion the old code fails. An equal epsilon in every feature is + # an equal share, so this column is uniform and sums to one; the old code + # left it at 0.0 while claiming the counts had been imputed. + col = prepared_zd[:, 3] + @test all(col .≈ 1 / length(tids_zd)) + @test !all(col .== 0.0) + @test sum(col) ≈ 1.0 + end + + @testset "drop — the empty sample is still dropped" begin + (prepared_dz, diag_dz, _, kept_dz, _) = Execution.prepare_analysis_table( + config, counts_zero_depth; + sample_ids=sids_zd, taxa_ids=tids_zd, + drop_policy="drop", impute_policy="epsilon" + ) + @test size(prepared_dz, 2) == 4 + @test "s3" ∉ kept_dz + @test any(h -> occursin("Dropped", h), diag_dz.healings) + # The property, whatever the policy: no column of the prepared table + # stands for a sample with no reads. + @test !any(all(c .== 0.0) for c in eachcol(prepared_dz)) + end + + @testset "refuse — the same refusal, raised earlier" begin + @test_throws ArgumentError Execution.prepare_analysis_table( + config, counts_zero_depth; + sample_ids=sids_zd, taxa_ids=tids_zd, + drop_policy="refuse", impute_policy="epsilon" + ) + end + end + + @testset "rarefy does not empty the matrix when a sample is empty" begin + # min_lib = minimum(lib_sizes) is 0 whenever any sample is empty, and every + # sample was then scaled by 0/lib. Measured before the fix: four correctly + # kept samples and a 2x4 matrix of zeros, reported as success. + norm_rar = AnalysisConfig.NormalizationConfig(method="rarefy", epsilon=1e-6, + zero_policy="pseudocount") + adv_rar = AnalysisConfig.AdvancedConfig(min_samples_per_group=2) + config_rar = AnalysisConfig.AnalysisConfig( + method="nb_glm", + formula="~ group", + metadata_columns=["group"], + normalization=norm_rar, + advanced=adv_rar, + created_by="test_zero_depth_rarefy" + ) + counts_zero_depth = [10.0 20.0 0.0 30.0 25.0; + 30.0 20.0 0.0 10.0 15.0] + (prepared_rar, _, _, kept_rar, _) = Execution.prepare_analysis_table( + config_rar, counts_zero_depth; + sample_ids=["s1", "s2", "s3", "s4", "s5"], taxa_ids=["t1", "t2"], + drop_policy="drop", impute_policy="epsilon" + ) + @test size(prepared_rar, 2) == 4 # s3 dropped + @test "s3" ∉ kept_rar + @test !all(prepared_rar .== 0.0) # the old code produced all zeros + @test any(prepared_rar .> 0) + end + @testset "run_analysis — stub with manifest and diagnostics, hard-stop DANGER banner" begin norm = AnalysisConfig.NormalizationConfig(method="size_factors", epsilon=1e-6) adv = AnalysisConfig.AdvancedConfig(min_samples_per_group=2)