From c49527a15535881a1d3034522bd8f06cfe894495 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Tue, 22 Sep 2026 13:38:54 +0000 Subject: [PATCH 1/3] fix(analysis): heal zero-depth samples before the transform A sample with no reads at all reached the transform stage, where every transform divides by a sample total, so one empty sample poisoned whichever transform ran: * relative wrote 0.0 into every feature. That states that each feature's relative abundance is exactly 0 -- a value -- for a sample whose relative abundances do not exist. 0/0 is undefined and 0 is a value, and no downstream check can tell them apart, because 0.0 is an ordinary number. * rarefy took min_lib = minimum(lib_sizes), which is 0 when any sample is empty, and scaled EVERY sample by 0/lib. Measured before this change: four correctly kept samples and a 2x4 prepared table of zeros, reported as success. * clr took 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 columns were discarded only after the damage was computed; under 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 code said as much in a comment: "we re-apply transform after healing counts?" followed, for the impute case, by "for stub, just keep". All-zero samples are now healed before the transform, which removes the class rather than the instance: past that point no zero-depth column can reach a transform, under any policy. The visible result changes only where it was wrong. drop still drops, refuse still refuses (earlier), and impute now lets the transform see the imputed counts, so an imputed sample's relative abundances are uniform rather than a column of 0.0. Runs whose input has no zero-depth sample are untouched: the healing step returns immediately when there is nothing to heal. All-zero TAXA are deliberately still healed after the transform. 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 -- whereas moving the drop earlier would change the geometric mean CLR centres on and alter results for analyses that were already correct. Fixing a defect is not licence to change the numbers around it. The two branches that could only be reached by the old ordering now refuse by name instead of assigning a plausible-looking number. They are unreachable through the public entry point by construction, so they are defensive: their job is to turn a silent wrong number into a loud failure if a future refactor breaks the ordering, not to be exercised by a test that cannot reach them. This changes results that saved analyses were computed from, so it is recorded as a behaviour change rather than slipped in: docs/statistics/behaviour-change-zero-depth-samples.md states what changed, who is affected (only impute runs over inputs containing a zero-depth sample, under relative, rarefy or clr), how to tell from the diagnostics, and what to do. Owner-authorised 2026-09-22 as a narrow, explicit exception to the standing instruction that existing algorithms not change silently. Tests: 15 new assertions pinning the class -- impute proportions uniform and summing to one, drop still dropping with no all-zero column surviving, refuse still refusing, and rarefy not emptying the matrix. Held against the pre-fix file as a negative control, these fail in exactly 5 places (3 impute, 2 rarefy) while the file's 25 pre-existing prepare_analysis_table assertions still pass, so the change is pinned without disturbing what already worked. The whole file passes 87/87 on the corrected code. Refs #1 --- CHANGELOG.md | 15 +++ .../behaviour-change-zero-depth-samples.md | 83 +++++++++++++ src/analysis/Execution.jl | 106 +++++++++++----- test/unit/test_execution.jl | 115 ++++++++++++++++++ 4 files changed, 291 insertions(+), 28 deletions(-) create mode 100644 docs/statistics/behaviour-change-zero-depth-samples.md 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..45231523 --- /dev/null +++ b/docs/statistics/behaviour-change-zero-depth-samples.md @@ -0,0 +1,83 @@ + +# 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 + +**Only runs whose input contains a zero-depth sample, and only under +`drop_policy=impute`.** The matrix below is the whole of it: + +| `drop_policy` | Behaviour before | Behaviour now | +| --- | --- | --- | +| `drop` | sample dropped (after the damage was computed) | sample dropped (before) — **same visible result** | +| `refuse` | refused | refused, earlier — **same refusal** | +| `impute` | counts imputed to `epsilon`, but `prepared` kept a column of `0.0` | counts imputed, transform sees the imputed counts — **relative abundances are now uniform** | + +So a saved analysis is affected if, and only if, **all** of the following hold: it used +`normalization.method` of `relative`, `rarefy`, or `clr`, its input contained at least +one sample whose counts were all zero, and its `drop_policy` was `impute`. 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 an +`Imputed ... all-zero samples` healing entry. A run that contains one, under a relative, +rarefy or clr transform, is affected. + +## 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..2ea48283 100644 --- a/src/analysis/Execution.jl +++ b/src/analysis/Execution.jl @@ -903,6 +903,55 @@ 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 + 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 +1052,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 +1159,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 @@ -1115,8 +1181,9 @@ function prepare_analysis_table( (checks, warnings, errors) = self_diagnostics(filtered_counts, prepared, config; sample_metadata=sample_metadata, raw_counts=counts) - 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,31 +1219,14 @@ 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) diff --git a/test/unit/test_execution.jl b/test/unit/test_execution.jl index d694cb15..14cfe3f6 100644 --- a/test/unit/test_execution.jl +++ b/test/unit/test_execution.jl @@ -233,6 +233,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) From 21d12ee7201b5e4032b446c1b65193fd2aa1395c Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:47:18 +0100 Subject: [PATCH 2/3] Update docs/statistics/behaviour-change-zero-depth-samples.md Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- .../behaviour-change-zero-depth-samples.md | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/statistics/behaviour-change-zero-depth-samples.md b/docs/statistics/behaviour-change-zero-depth-samples.md index 45231523..a6371c51 100644 --- a/docs/statistics/behaviour-change-zero-depth-samples.md +++ b/docs/statistics/behaviour-change-zero-depth-samples.md @@ -35,19 +35,21 @@ face, and no downstream check can catch it, because `0.0` is an ordinary number. ## Who is affected -**Only runs whose input contains a zero-depth sample, and only under -`drop_policy=impute`.** The matrix below is the whole of it: +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` | sample dropped (after the damage was computed) | sample dropped (before) — **same visible result** | +| `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 imputed to `epsilon`, but `prepared` kept a column of `0.0` | counts imputed, transform sees the imputed counts — **relative abundances are now uniform** | - -So a saved analysis is affected if, and only if, **all** of the following hold: it used -`normalization.method` of `relative`, `rarefy`, or `clr`, its input contained at least -one sample whose counts were all zero, and its `drop_policy` was `impute`. Analyses whose -inputs had no zero-depth sample are unchanged — the healing step returns immediately when +| `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 From efd23755305f493ef670bc3cb6f6edd6fb78618b Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:05:23 +0000 Subject: [PATCH 3/3] fix: Preserve zero-taxon detection after imputation Capture all-zero taxa before healing zero-depth samples so diagnostics and taxon healing retain them. Add a regression test and clarify which analyses require reruns. --- .../behaviour-change-zero-depth-samples.md | 11 +++++-- src/analysis/Execution.jl | 30 +++++++++++++++---- test/unit/test_execution.jl | 17 +++++++++++ 3 files changed, 50 insertions(+), 8 deletions(-) diff --git a/docs/statistics/behaviour-change-zero-depth-samples.md b/docs/statistics/behaviour-change-zero-depth-samples.md index a6371c51..78e651b2 100644 --- a/docs/statistics/behaviour-change-zero-depth-samples.md +++ b/docs/statistics/behaviour-change-zero-depth-samples.md @@ -58,9 +58,14 @@ Re-run the affected analysis. The previous numbers for a zero-depth sample were 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 an -`Imputed ... all-zero samples` healing entry. A run that contains one, under a relative, -rarefy or clr transform, is affected. +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 diff --git a/src/analysis/Execution.jl b/src/analysis/Execution.jl index 2ea48283..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 @@ -944,6 +949,10 @@ function prepare_analysis_table( # 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, _) = @@ -1180,6 +1189,11 @@ 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` and `is_dangerous_diag` are initialised earlier, before the # all-zero sample healing that happens before the transform; re-initialising @@ -1229,7 +1243,13 @@ function prepare_analysis_table( # `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 14cfe3f6..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