From 2280808cad7bbc022c01da427884bd9792874836 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 25 Sep 2026 20:40:35 +0000 Subject: [PATCH] refactor(analysis): add the scaling-factor module and publish its conditions Library size, CSS, TMM/RSS and median-of-ratios size factors, each returning factors centred to geometric mean 1 and the log-offset a count model is fitted with, plus the declared parameters, the uncentred quantities and a SHA-256 of the offset vector for the manifest. The conditions this code is held to are published first in `docs/statistics/method-conditions/scaling-and-offsets.md` (issue #16). Nothing consumes the module yet: `prepare_analysis_table` still sends `TSS` through the proportional branch and aliases `CSS`/`RSS` to `relative`, and `size_factors` still computes library size divided by its own geometric mean. The wiring and its tests follow in the next commit, in which the current behaviour is deleted rather than replaced silently. Also: the CI test step now posts the tail of a failing `test/runtests.jl` as a check annotation. The job log is served from a storage host that is not reachable from every environment that needs to read a failure, so the failure message effectively had no reader. Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .github/workflows/ci.yml | 15 +- .../method-conditions/scaling-and-offsets.md | 172 ++++++ src/analysis/scaling.jl | 503 ++++++++++++++++++ 3 files changed, 689 insertions(+), 1 deletion(-) create mode 100644 docs/statistics/method-conditions/scaling-and-offsets.md create mode 100644 src/analysis/scaling.jl 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