diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a9b3f1f1..fde40608 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -542,7 +542,20 @@ jobs: env: CI_SKIP_TAXONOMY: "1" R_LIBS_SITE: /usr/local/lib/R/site-library:/usr/lib/R/site-library:/usr/lib/R/library - run: julia --project=. -t 2 --code-coverage=user --compiled-modules=no test/runtests.jl --integration --server + run: | + set +e + julia --project=. -t 2 --code-coverage=user --compiled-modules=no test/runtests.jl --integration --server 2>&1 | tee /tmp/julia-tests.log + status="${PIPESTATUS[0]}" + if [ "$status" -ne 0 ]; then + # A failing test's message is the entire point of running it, and the job log + # is served from a storage host that not every environment can reach. The tail + # of the output is therefore also posted as an annotation, which the checks + # API serves as data. 6000 bytes of tail is enough for the failing testset and + # the last few assertions, which is what a reader needs first. + tail -c 6000 /tmp/julia-tests.log | sed -e 's/%/%25/g' | awk '{printf "%s%%0A", $0}' > /tmp/julia-tests-annotation.txt + echo "::error title=Julia tests failed::$(cat /tmp/julia-tests-annotation.txt)" + fi + exit "$status" - name: Process coverage uses: julia-actions/julia-processcoverage@03114f09f119417c3242a9fb6e0b722676aedf38 # v1 diff --git a/docs/statistics/method-conditions/scaling-and-offsets.md b/docs/statistics/method-conditions/scaling-and-offsets.md new file mode 100644 index 00000000..ab87566c --- /dev/null +++ b/docs/statistics/method-conditions/scaling-and-offsets.md @@ -0,0 +1,172 @@ + +# Method conditions — library-size scaling and offsets + +**Status: published 2026-09-25, before the implementation in +`src/analysis/scaling.jl`.** Where this document and that file disagree, the document is +right and the file is the bug. + +**Why it exists now.** Issue #16 specified exact TSS/CSS/RSS behaviour and carried the +line *"Deferred — DO NOT IMPLEMENT IN MILESTONE 3"*. That line is a scope decision, not a +scientific finding, and the owner directed on 2026-09-25 that the remaining deferred +items be implemented. The deferral is therefore lifted here, explicitly and in writing, +for this item only. The conditions below are what "implemented" is allowed to mean. + +**What it replaces.** Three normalisation methods that the configuration accepted and +the execution path silently turned into something else: + +- `TSS` divided counts by library size and returned **proportions** — the compositional + transform the issue was written to avoid — for every method except `nb_glm`. +- `CSS` and `RSS` were aliased to `relative` outright, with a `@warn`. Under `nb_glm` + that produced a proportion table where a count model expects counts. +- `size_factors` was described as "DESeq2 median-of-ratios" in the code and implemented + as library size divided by its own geometric mean — TSS wearing another method's name. + +## What this method is + +A **scaling factor** is one positive number per sample; an **offset** is that factor's +logarithm (after a stated centring), handed to a count model so that sequencing depth is +modelled rather than analysed. Neither changes the response. A **transform** (relative, +CLR, ILR, presence/absence) changes what the response *is*. Conflating the two is what +made `nb_glm` + `TSS` report proportions as if they were counts. + +All four kinds below return factors normalised to **geometric mean 1**, so the offset is +centred and the constant absorbed by the model's intercept. The uncentred quantities (raw +sums, library sizes, log-ratios) are recorded in the provenance, so the published numbers +can be reproduced without the centring convention. + +### TSS — total sum scaling (library size), offset form + +- `factor_j = library_size_j / geomean(library_size)`, `library_size_j = Σᵢ xᵢⱼ`. +- Offset `log(factor_j)`. Applicable to count responses only; for a non-count response + TSS *is* the proportional transform and is reported as `relative`, not as an offset. +- Zeros: nothing is replaced. A sample with total 0 has no library size to divide by and + is refused — it is reported as an unsuccessful state, never as factor 0 or 1. +- This is the McMurdie & Holmes (2014) argument: keep counts, model depth. + +### CSS — cumulative sum scaling (Paulson et al. 2013, metagenomeSeq) + +- Per sample `j`: `threshold_j = quantile(x·ⱼ, p)` with `p = css_quantile` + (default 0.75), using the quantile definition of R's default `type = 7` and Julia's + `Statistics.quantile` default. All features in the sample participate, zeros included. +- `sum_j = Σ { xᵢⱼ : xᵢⱼ ≤ threshold_j }` — the cumulative sum **up to** the quantile. +- `factor_j = sum_j / geomean(sum)`. Offset `log(factor_j)`. +- Refused when any `sum_j ≤ 0`: with a low quantile on a mostly-zero sample the + cumulative sum can be zero, and `log(0)` is not a small number. The refusal names the + sample and the quantile (this is the failure mode issue #16 called out in advance). +- **Not implemented, and not silently substituted:** metagenomeSeq's `cumNormStatFast`, + which *chooses* `p` from the data. A data-driven choice of a modelling parameter is a + decision the run has to record; this layer instead requires `p` to be declared by the + caller and records it. Parity with `metagenomeSeq::cumNorm` is therefore an + **outstanding condition**, not a claim (see *Evidence*). + +### RSS / TMM — trimmed mean of M-values (Robinson & Oshlack 2010, edgeR) + +RSS in the issue's configuration vocabulary is the TMM estimator the issue's own text +describes. One sample is the reference; every other sample's factor is a trimmed, +weighted mean of its log-ratios to the reference. + +- Library sizes are column sums. The reference sample is the one whose + upper-quartile-scaled count distribution is closest to the mean of those distributions + (the edgeR default), or the sample named by `tmm_ref_column`. +- For sample `i` against reference `r`, over features with both counts positive: + `M = log2(xᵢ/x_r)`, `A = ½·log2(xᵢ·x_r)`, + `w = (1 − xᵢ/libᵢ)/xᵢ + (1 − x_r/lib_r)/x_r`. +- Features are trimmed by rank: the top `tmm_log_ratio_trim` (default 0.3) and bottom + `tmm_log_ratio_trim` of `M`, and the tails of `A` beyond `tmm_sum_trim` (default 0.05), + are dropped. Ties are ranked by first occurrence, matching R's `ties.method = "first"`. +- `factor_i = 2^(Σ M·w / Σ w)` over the kept features, then normalised as above. The + reference sample gets factor 1 before normalisation. +- Refused when a sample has no positive-count feature in common with the reference, or + when a library size is zero. A refusal is an unsuccessful state; it is never a factor of + 1 by default. + +### `size_factors` — median-of-ratios (RLE) + +- Per sample `j`: `factor_j = median { xᵢⱼ / geomean_k(xᵢ·)}` over features whose + geometric mean across samples is positive. Offset `log(factor_j)`. +- This is the estimator the code previously claimed and did not compute. The previous + behaviour (library-size centring) is available as `TSS`, which is where it belongs. + +### Accepted responses and study designs + +- Count responses for `nb_glm` (and any future count model): scaling factors, offsets, + no modification of the response. +- Non-count responses: TSS is reported as the `relative` transform. CSS, RSS and + `size_factors` are **refused for non-count responses by the configuration layer** — + they are offsets, and there is nothing to offset. +- Study design is untouched by this layer. Pairing, blocking and repeated measures are + still refused by the method layer (see `parametric-fits.md`); a scaling factor is not a + design. + +### Zeros, ties and degenerate samples + +- Zeros are neither replaced nor imputed here. A zero count contributes to CSS only + through the threshold rule, and to TMM not at all (pairwise-positive features only). +- All-zero samples are refused by name before any factor is computed. +- Constant samples (`xᵢⱼ = c` for all `i`) are defined: TSS and CSS give every such sample + the same factor; TMM's log-ratios are all 0, giving factor 1; RLE gives factor 1 + (each ratio is 1). None of these is special-cased — the mathematics already holds. +- Ties in `M` and `A` are ranked deterministically by first occurrence, so a run is + reproducible; this is asserted, because R's default and a naive sort differ on ties. + +### Diagnostics and provenance + +Every run that computes factors records, in `diagnostics.checks["scaling"]` and in the +manifest provenance: + +- the kind (`tss`, `css`, `rss`, `size_factors`) and the definition sentence for it; +- the declared parameters: quantile, trims, named reference (if any), and which sample + was chosen as the reference when it was chosen from the data; +- the raw quantities before centring (cumulative sums, library sizes, or log-ratio means) + and their geometric mean, so the paper's numbers are recoverable; +- every warning (for example: very low `css_quantile`, or a weighted mean that fell back + to unweighted because all weights were zero). + +### Computational limits + +All four are O(n_features × n_samples), with an O(n log n) sort per sample for CSS and +per pair for TMM. No per-run memory bound has been measured or is claimed. A large +feature count costs time linearly; nothing here is iterative and nothing is randomised. + +## Evidence, as delivered + +In `test/unit/test_scaling.jl` and `test/unit/test_execution.jl`: + +- **Known answers from first principles** — factors for a hand-computable table, with the + arithmetic written out next to the expected number, not recorded from a previous run. +- **Properties that a wrong implementation fails**: doubling every count in one sample + doubles that sample's TSS factor; a sample that is an exact multiple of another has a + factor equal to that multiple (TMM, all log-ratios equal); permuting features does not + change any factor; scaling every sample by the same constant changes nothing. +- **CSS robustness to an outlier** — adding one dominant feature to one sample moves its + TSS factor far more than its CSS factor, which is the property CSS exists for. +- **Negative controls** — zero-total sample, `css_quantile` low enough to give a zero + cumulative sum, unknown `tmm_ref_column`, out-of-range trims, and CSS/RSS on a + non-count response all refused, each with the reason named. +- **Integration** — `prepare_analysis_table` keeps counts as the response under `nb_glm` + and hands the offset to `run_analysis`; the manifest records the parameters; the offset + hash changes when the parameters change. +- **Limits, stated rather than hidden**: parity with `metagenomeSeq::cumNorm` and + `edgeR::calcNormFactors` at 1e-6 is **not verified here**. Neither package is in the + pinned R environment (`renv.lock`), and this repository does not add an unpinned R + dependency to make a test pass. A transcription cross-check in R is included for TMM + (same author, different language: it catches implementation slips, not specification + errors) and is labelled as such. This is recorded as an outstanding condition in + issue #16 rather than claimed as met. + +## Explicitly not claimed + +- No scaling factor makes compositional data safe. CSS and TMM correct library size; + they do not remove the compositional constraint, and a log fold change from a model + with these offsets is still relative to the sampled community. Issue #16 asked for + this warning in the user-facing help, and it is there. +- TMM's exchangeability assumption (most features not differentially abundant) is not + testable from the data and is not checked. An experiment where a large fraction of + features moves in the same direction will bias TMM; nothing here detects that. +- CSS's data-driven quantile (`cumNormStatFast`) is not implemented; the declared + quantile must be chosen by the analyst and is recorded. +- No claim is made that these estimates are more accurate than the alternatives for any + particular dataset. They are named, computed, and recorded; choosing between them is + the analyst's decision, and the choice is in the provenance. diff --git a/src/analysis/scaling.jl b/src/analysis/scaling.jl new file mode 100644 index 00000000..5fa685b3 --- /dev/null +++ b/src/analysis/scaling.jl @@ -0,0 +1,503 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Library-size scaling factors and the offsets count models are fitted with. +# +# The conditions this code is held to were published before it existed +# (docs/statistics/method-conditions/scaling-and-offsets.md). Where this file and that +# document disagree, the document is right and this file is the bug. +# +# Three separate things lived under one name before this module, and the configuration +# accepted all three while the execution path did something else: +# +# TSS turned counts into proportions for every method but `nb_glm`, which is +# the compositional transform issue #16 was written to avoid. +# CSS, RSS were aliased to `relative` outright. Under `nb_glm` that fed a count +# model a table of fractions. +# size_factors was labelled "DESeq2 median-of-ratios" and computed library size +# divided by its own geometric mean -- TSS wearing another name. +# +# A scaling factor is one positive number per sample. An offset is its logarithm. Neither +# changes the response. That distinction is the whole point of the module. +module Scaling + +using Statistics +using OrderedCollections +using JSON3 +using SHA + +export ScalingOutcome, ScalingRefusal, tss_factors, css_factors, tmm_factors, rle_factors, + factors_for, factor_checks, factor_provenance, OFFSET_DEFINITIONS + +## Definitions, written once and carried into the provenance + +const QUANTILE_DEFINITION = "R type-7 quantile (the default in R and in Julia's Statistics): " * + "h = (n-1)p + 1 with linear interpolation between order statistics" + +const OFFSET_DEFINITIONS = Dict{String,String}( + "tss" => "log(library size), library size = column sum of the counts; " * + "factors are the library sizes divided by their geometric mean", + "css" => "log(cumulative sum of the counts at or below the declared per-sample quantile), " * + "divided by the geometric mean of those sums (Paulson et al. 2013)", + "rss" => "log(2^(weighted trimmed mean of log2 ratios to the reference sample)), " * + "divided by the geometric mean of the factors (Robinson & Oshlack 2010, TMM)", + "size_factors" => "log(median over features of count / feature geometric mean), " * + "divided by the geometric mean of the factors (DESeq2/RLE median-of-ratios)", +) + +## Types + +""" + ScalingRefusal(kind, reason) + +A scaling factor that cannot be computed as declared. Carries the method `kind` and a +`reason` that names the sample and the quantity, so the refusal can be shown to a user +without translating it first. A refused factor is never replaced by a default of 1. +""" +struct ScalingRefusal <: Exception + kind::String + reason::String +end + +Base.showerror(io::IO, e::ScalingRefusal) = print(io, "ScalingRefusal($(e.kind)): $(e.reason)") + +""" + ScalingOutcome + +One scaling computation: the factors, the offset they imply, and everything needed to +reproduce and audit them. + +- `kind`: `"tss"`, `"css"`, `"rss"` or `"size_factors"`. +- `definition`: the sentence from `OFFSET_DEFINITIONS` that this run was computed under. +- `factors`: one factor per sample, normalised to geometric mean 1. +- `offset`: `log.(factors)` — what a count model is fitted with. +- `reference`: for TMM, the sample the ratios were taken against (`nothing` otherwise). +- `raw`: the uncentred quantity (library sizes, cumulative sums, TMM factors before + centring, or RLE medians), so the published numbers are recoverable. +- `geometric_mean`: the centring constant that was divided out. +- `parameters`: the declared settings that changed the answer. +- `warnings`: non-fatal notes that belong in front of the analyst. +""" +struct ScalingOutcome + kind::String + definition::String + factors::Vector{Float64} + offset::Vector{Float64} + reference::Union{String,Nothing} + raw::Vector{Float64} + geometric_mean::Float64 + parameters::OrderedDict{String,Any} + warnings::Vector{String} +end + +## Helpers + +function _labels(n_samples::Int, sample_ids)::Vector{String} + if isnothing(sample_ids) + return ["column $j" for j in 1:n_samples] + end + length(sample_ids) == n_samples || + throw(ArgumentError("sample_ids has $(length(sample_ids)) entries but the counts have $n_samples columns")) + return String[string(s) for s in sample_ids] +end + +function _column_totals(counts::AbstractMatrix{<:Real})::Vector{Float64} + n = size(counts, 2) + totals = Vector{Float64}(undef, n) + for j in 1:n + totals[j] = sum(Float64(counts[i, j]) for i in 1:size(counts, 1)) + end + return totals +end + +function _require_positive(values::Vector{Float64}, kind::String, labels::Vector{String}, what::String) + for j in eachindex(values) + v = values[j] + if !(isfinite(v) && v > 0) + throw(ScalingRefusal(kind, + "$what for sample '$(labels[j])' (column $j) is $(v); a scaling factor needs a " * + "positive finite quantity, and neither 0 nor 1 is a defensible stand-in. " * + "Refused — see docs/statistics/method-conditions/scaling-and-offsets.md")) + end + end + return nothing +end + +function _centre(raw::Vector{Float64}, kind::String)::Tuple{Vector{Float64},Float64} + g = exp(sum(log, raw) / length(raw)) + return raw ./ g, g +end + +""" + _first_ranks(v) + +Ranks of `v` ascending, ties broken by first occurrence — the ordering R's +`ties.method = "first"` gives, and the one edgeR's trimming assumes. Implemented by +sorting `(value, index)` pairs rather than relying on the stability of `sortperm`. +""" +function _first_ranks(v::Vector{Float64})::Vector{Int} + order = sort!(collect(1:length(v)); by = i -> (v[i], i)) + ranks = Vector{Int}(undef, length(v)) + for (position, index) in enumerate(order) + ranks[index] = position + end + return ranks +end + +## TSS + +""" + tss_factors(counts; sample_ids) -> ScalingOutcome + +Total sum scaling in its offset form: `log(library size)`, centred to geometric mean 1. + +Library sizes are column sums. Zeros are untouched. A sample with a zero (or non-finite) +total has no library size to divide by and is refused by name. +""" +function tss_factors(counts::AbstractMatrix{<:Real}; + sample_ids::Union{Nothing,AbstractVector}=nothing) + labels = _labels(size(counts, 2), sample_ids) + totals = _column_totals(counts) + _require_positive(totals, "tss", labels, "total count") + factors, g = _centre(totals, "tss") + return ScalingOutcome("tss", OFFSET_DEFINITIONS["tss"], factors, log.(factors), nothing, + totals, g, + OrderedDict{String,Any}("library_sizes" => totals, + "definition_reference" => "McMurdie & Holmes 2014"), + String[]) +end + +## CSS + +""" + css_factors(counts; quantile, sample_ids) -> ScalingOutcome + +Cumulative sum scaling (Paulson et al. 2013). Per sample: the sum of counts at or below +that sample's `quantile`-th quantile, centred to geometric mean 1; the offset is its log. + +The quantile must be declared by the caller — `metagenomeSeq`'s data-driven choice of it +(`cumNormStatFast`) is deliberately not implemented, because a parameter chosen from the +data has to be recorded as a modelling decision rather than a default. + +Refused when the cumulative sum is not positive for every sample: on a mostly-zero sample +at a low quantile it is zero, and `log(0)` is not a small number. +""" +function css_factors(counts::AbstractMatrix{<:Real}; + quantile::Real=0.75, + sample_ids::Union{Nothing,AbstractVector}=nothing) + kind = "css" + p = Float64(quantile) + (0.0 < p < 1.0) || throw(ArgumentError( + "css_quantile must be in (0,1), got $p — see docs/statistics/method-conditions/scaling-and-offsets.md")) + labels = _labels(size(counts, 2), sample_ids) + n_samples = size(counts, 2) + thresholds = Vector{Float64}(undef, n_samples) + sums = Vector{Float64}(undef, n_samples) + for j in 1:n_samples + column = [Float64(counts[i, j]) for i in 1:size(counts, 1)] + threshold = Statistics.quantile(column, p) + thresholds[j] = threshold + sums[j] = sum(v for v in column if v <= threshold) + end + for j in eachindex(sums) + if !(isfinite(sums[j]) && sums[j] > 0) + throw(ScalingRefusal(kind, + "the cumulative sum up to the $p quantile is $(sums[j]) for sample " * + "'$(labels[j])' (column $j): no count at or below the threshold is positive, " * + "so the CSS factor would be log(0). Raise css_quantile or exclude the " * + "sample. Refused rather than substituted.")) + end + end + factors, g = _centre(sums, kind) + warnings = String[] + if p < 0.5 + push!(warnings, "css_quantile=$p is below the median: the cumulative sum then covers " * + "less than half of each sample's counts and is dominated by how many " * + "features are zero. Paulson et al. (2013) use 0.5-0.75 on real data.") + end + return ScalingOutcome(kind, OFFSET_DEFINITIONS["css"], factors, log.(factors), nothing, + sums, g, + OrderedDict{String,Any}("quantile" => p, + "quantile_definition" => QUANTILE_DEFINITION, + "cumulative_sums" => sums, + "thresholds" => thresholds, + "definition_reference" => "Paulson et al. 2013"), + warnings) +end + +## RSS / TMM + +function _tmm_pair(x::Vector{Float64}, r::Vector{Float64}, lib_x::Float64, lib_r::Float64, + log_ratio_trim::Float64, sum_trim::Float64) + idx = [k for k in eachindex(x) if x[k] > 0.0 && r[k] > 0.0] + n = length(idx) + n == 0 && return (NaN, 0, 0, false) + m = Vector{Float64}(undef, n) + a = Vector{Float64}(undef, n) + w = Vector{Float64}(undef, n) + for (slot, k) in enumerate(idx) + m[slot] = log2(x[k] / r[k]) + a[slot] = 0.5 * log2(x[k] * r[k]) + w[slot] = (1.0 - x[k] / lib_x) / x[k] + (1.0 - r[k] / lib_r) / r[k] + end + lo_l = floor(Int, n * log_ratio_trim) + 1 + hi_l = n + 1 - lo_l + lo_s = floor(Int, n * sum_trim) + 1 + hi_s = n + 1 - lo_s + rank_m = _first_ranks(m) + rank_a = _first_ranks(a) + keep = [ (lo_l <= rank_m[k] <= hi_l) && (lo_s <= rank_a[k] <= hi_s) for k in 1:n ] + any(keep) || return (NaN, 0, n, false) + if all(w[k] <= 0.0 for k in 1:n if keep[k]) + # Every kept feature has x == library size or the reference does, so the weights + # carry no information. Fall back to the unweighted trimmed mean and say so. + return (sum(m[k] for k in 1:n if keep[k]) / count(keep), count(keep), n, false) + end + numerator = sum(m[k] * w[k] for k in 1:n if keep[k] && w[k] > 0.0) + denominator = sum(w[k] for k in 1:n if keep[k] && w[k] > 0.0) + return (numerator / denominator, count(keep), n, true) +end + +""" + tmm_factors(counts; ref_column, log_ratio_trim, sum_trim, sample_ids) -> ScalingOutcome + +The TMM estimator (Robinson & Oshlack 2010) — what the configuration calls `RSS`. + +One sample is the reference: named by `ref_column`, or, when that is `nothing`, the sample +whose upper-quartile-scaled counts are closest to the mean of those values across samples +(edgeR's default). Each other sample's factor is `2` to the power of the weighted, +trimmed mean of its log2 ratios to the reference, then all factors are centred. + +Trimming drops the tails of the log-ratios (`log_ratio_trim`, default 0.3) and of the mean +abundances (`sum_trim`, default 0.05), ranked with ties by first occurrence. + +Refused when a sample shares no positive feature with the reference, or when a sample's +treated factor is not positive and finite. +""" +function tmm_factors(counts::AbstractMatrix{<:Real}; + ref_column::Union{String,Nothing}=nothing, + log_ratio_trim::Real=0.3, + sum_trim::Real=0.05, + sample_ids::Union{Nothing,AbstractVector}=nothing) + kind = "rss" + lrt = Float64(log_ratio_trim) + st = Float64(sum_trim) + (0.0 <= lrt < 0.5) || throw(ArgumentError( + "tmm_log_ratio_trim must be in [0, 0.5), got $lrt — a trim of 0.5 or more removes at " * + "least half of the log-ratios in each tail. See docs/statistics/method-conditions/scaling-and-offsets.md")) + (0.0 <= st < 0.5) || throw(ArgumentError( + "tmm_sum_trim must be in [0, 0.5), got $st — see docs/statistics/method-conditions/scaling-and-offsets.md")) + labels = _labels(size(counts, 2), sample_ids) + n = size(counts, 2) + libs = _column_totals(counts) + _require_positive(libs, kind, labels, "total count") + + chosen_from_data = isnothing(ref_column) + reference = "" + if chosen_from_data + upper_quartile = Vector{Float64}(undef, n) + for j in 1:n + column = [Float64(counts[i, j]) for i in 1:size(counts, 1)] + upper_quartile[j] = Statistics.quantile(column, 0.75) / libs[j] + end + target = sum(upper_quartile) / n + distances = [abs(upper_quartile[j] - target) for j in 1:n] + ref_index = argmin(distances) + reference = labels[ref_index] + else + ref_index = findfirst(==(String(ref_column)), labels) + isnothing(ref_index) && throw(ScalingRefusal(kind, + "tmm_ref_column='$(ref_column)' is not a sample in this run. Samples are: " * + "$(join(labels, ", ")). Refused — a reference chosen by a typo is not a reference.")) + reference = labels[ref_index] + end + + raw = ones(n) + kept_counts = zeros(Int, n) + pair_counts = zeros(Int, n) + weighted_mean_used = falses(n) + for i in 1:n + i == ref_index && continue + x = [Float64(counts[k, i]) for k in 1:size(counts, 1)] + r = [Float64(counts[k, ref_index]) for k in 1:size(counts, 1)] + mean_log_ratio, kept, pairs, weighted = _tmm_pair(x, r, libs[i], libs[ref_index], lrt, st) + kept_counts[i] = kept + pair_counts[i] = pairs + weighted_mean_used[i] = weighted + if !isfinite(mean_log_ratio) + throw(ScalingRefusal(kind, + "sample '$(labels[i])' (column $i) shares no positive-count feature with the " * + "reference sample '$(reference)', or its trimmed set is empty " * + "($kept of $pairs features kept). A TMM factor cannot be computed; refused " * + "rather than set to 1.")) + end + raw[i] = 2.0^mean_log_ratio + end + + _require_positive(raw, kind, labels, "TMM factor") + factors, g = _centre(raw, kind) + warnings = String[] + if any(!weighted_mean_used[i] for i in 1:n if i != ref_index) + affected = [labels[i] for i in 1:n if i != ref_index && !weighted_mean_used[i]] + push!(warnings, "the weighted trimmed mean was undefined for $(join(affected, ", ")) " * + "(every kept feature is a whole library). The unweighted trimmed mean " * + "was used for those samples; the factors are still reported, and the " * + "difference is recorded here rather than hidden.") + end + if any(kept_counts[i] < 2 for i in 1:n if i != ref_index) + push!(warnings, "fewer than two log-ratios survived trimming for at least one sample; " * + "a trimmed mean over one value is that value. Recorded, not suppressed.") + end + return ScalingOutcome(kind, OFFSET_DEFINITIONS["rss"], factors, log.(factors), reference, + raw, g, + OrderedDict{String,Any}( + "log_ratio_trim" => lrt, + "sum_trim" => st, + "reference_sample" => reference, + "reference_chosen_from_data" => chosen_from_data, + "reference_selection" => chosen_from_data ? + "sample whose upper-quartile-scaled counts are closest to the mean of those values (edgeR default)" : + "declared by tmm_ref_column", + "kept_log_ratios" => kept_counts, + "available_pairs" => pair_counts, + "weighted" => [weighted_mean_used[i] for i in 1:n], + "definition_reference" => "Robinson & Oshlack 2010"), + warnings) +end + +## size_factors: median-of-ratios (RLE) + +""" + rle_factors(counts; sample_ids) -> ScalingOutcome + +Median-of-ratios size factors, the estimator `size_factors` claimed to be and was not. + +Per sample: the median over features of `count / feature geometric mean`, using only +features whose counts are positive in **every** sample (a feature with a zero has no +finite geometric mean, so it cannot contribute a ratio). Centred to geometric mean 1. + +Refused when no feature is positive in every sample, and when a median ratio is not +positive. +""" +function rle_factors(counts::AbstractMatrix{<:Real}; + sample_ids::Union{Nothing,AbstractVector}=nothing) + kind = "size_factors" + labels = _labels(size(counts, 2), sample_ids) + n_features = size(counts, 1) + n_samples = size(counts, 2) + usable = Int[] + log_geo_mean = Float64[] + for i in 1:n_features + row = [Float64(counts[i, j]) for j in 1:n_samples] + if all(>(0.0), row) + push!(usable, i) + push!(log_geo_mean, sum(log, row) / n_samples) + end + end + isempty(usable) && throw(ScalingRefusal(kind, + "no feature has a positive count in every sample, so no feature has a finite " * + "geometric mean and the median-of-ratios size factor is undefined. Refused rather " * + "than falling back to library size, which is a different estimator.")) + ratios = Vector{Float64}(undef, n_samples) + for j in 1:n_samples + ratios[j] = Statistics.median([Float64(counts[usable[k], j]) / exp(log_geo_mean[k]) + for k in eachindex(usable)]) + end + _require_positive(ratios, kind, labels, "median-of-ratios size factor") + factors, g = _centre(ratios, kind) + warnings = String[] + if length(usable) < n_features + push!(warnings, "$(n_features - length(usable)) of $n_features features contain a zero " * + "in at least one sample and cannot contribute a ratio to the median " * + "(no finite geometric mean). The estimate rests on the other " * + "$(length(usable)).") + end + if length(usable) < n_samples + push!(warnings, "only $(length(usable)) features are positive in every sample; " * + "the median of a small set of ratios is a noisy size factor.") + end + return ScalingOutcome(kind, OFFSET_DEFINITIONS["size_factors"], factors, log.(factors), + nothing, ratios, g, + OrderedDict{String,Any}("features_used" => length(usable), + "features_total" => n_features, + "log_geometric_means" => + Dict(string(k) => log_geo_mean[k] + for k in eachindex(usable)), + "definition_reference" => + "Anders & Huber 2010; DESeq2/RLE"), + warnings) +end + +## Dispatch + +""" + factors_for(kind, counts; kwargs...) -> ScalingOutcome + +Dispatch on the normalisation method name, case-insensitively: `none` and `tss` to +[`tss_factors`](@ref), `css` to [`css_factors`](@ref), `rss`/`tmm` to +[`tmm_factors`](@ref), `size_factors` to [`rle_factors`](@ref). + +`none` means "no scaling": the offset is the plain log library size, which is what +[`tss_factors`](@ref) returns. The two names are kept distinct in the configuration +because one declares a transform and the other declares its absence, and both are +recorded as what they are. +""" +function factors_for(kind::AbstractString, counts::AbstractMatrix{<:Real}; kwargs...) + key = lowercase(strip(kind)) + if key in ("none", "tss") + return tss_factors(counts; kwargs...) + elseif key == "css" + return css_factors(counts; kwargs...) + elseif key in ("rss", "tmm") + return tmm_factors(counts; kwargs...) + elseif key == "size_factors" + return rle_factors(counts; kwargs...) + end + throw(ArgumentError("no scaling factors are defined for method '$kind' — see " * + "docs/statistics/method-conditions/scaling-and-offsets.md")) +end + +## What the run records + +""" + factor_checks(outcome) -> OrderedDict + +The entry written to `diagnostics.checks["scaling"]`: what was computed, under which +definition, with which declared parameters, and what the analyst should be told. +""" +function factor_checks(o::ScalingOutcome)::OrderedDict{String,Any} + return OrderedDict{String,Any}( + "kind" => o.kind, + "definition" => o.definition, + "offset_definition" => "log(factor), factors centred to geometric mean 1", + "factors" => o.factors, + "raw" => o.raw, + "centring_geometric_mean" => o.geometric_mean, + "reference_sample" => o.reference, + "parameters" => o.parameters, + "warnings" => o.warnings, + ) +end + +""" + factor_provenance(outcome) -> OrderedDict + +The entry written to the manifest provenance: the same facts as [`factor_checks`](@ref) +plus a SHA-256 of the offset vector, so two runs that agree on the offset can be shown to +agree and a run whose offset changed is detectable from the manifest alone. +""" +function factor_provenance(o::ScalingOutcome)::OrderedDict{String,Any} + return OrderedDict{String,Any}( + "kind" => o.kind, + "definition" => o.definition, + "factors" => o.factors, + "offset" => o.offset, + "offset_sha256" => bytes2hex(sha256(JSON3.write(o.offset))), + "centring_geometric_mean" => o.geometric_mean, + "reference_sample" => o.reference, + "parameters" => o.parameters, + "warnings" => o.warnings, + ) +end + +end # module Scaling