From 8671521f64a683c7cfcebebcb12a4fa509195ceb Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 25 Sep 2026 20:52:03 +0000 Subject: [PATCH] fix(analysis): method names compare in lower case; offsets wired in MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two things, one cause: the layer that claimed to compute TSS/CSS/RSS offsets did not, and the layer that accepted them refused them. The refusal. NormalizationConfig canonicalises its method to lower case while VALID_NORMALIZATION_FOR_METHOD listed "TSS", "CSS", "RSS", so the membership test failed for those three spellings and AnalysisConfig raised "normalization.method 'tss' incompatible with method 'nb_glm'". That is what reddened the CI run for the merged TSS-offsets commit — the existing test_execution.jl testset constructs method="TSS" exactly as documented. The table and both comparisons are lower case now, and the new tests assert that every admissible spelling of every method name is accepted. The wiring. prepare_analysis_table now computes scaling factors through src/analysis/scaling.jl and hands over the offset, rather than aliasing TSS/CSS/RSS to `relative` (which fed a count model a table of fractions) and calling library-size centring "DESeq2 median-of-ratios": tss log library size (counts stay counts) css log cumulative sum at the declared css_quantile rss (TMM) log trimmed weighted mean of log-ratios to a reference sample size_factors median-of-ratios (RLE), the estimator the name claimed none the plain log library size, as it always meant for a count model css and rss are refused for a response with no counts to offset, in the configuration layer and again in the execution path. css_quantile, tmm_ref_column, tmm_log_ratio_trim and tmm_sum_trim are validated in the constructor, part of the config hash and canonical JSON, round-tripped through to_json/from_json/Nickel/DEED, mirrored in the JSON schema and Nickel contracts, exposed through the server route, and explained in context_help and the frontend help. What was computed is recorded in diagnostics.checks["scaling"] and the manifest provenance, including the parameters, the raw quantities before centring and a SHA-256 of the offset vector. Evidence: test/unit/test_scaling.jl — known answers with the arithmetic written beside them, properties a wrong implementation fails (feature permutation, sample-multiple invariance, shared-constant invariance), the CSS outlier robustness property, negative controls for every refusal, the integration path through prepare_analysis_table, and a base-R transcription of the same conditions. Parity with metagenomeSeq::cumNorm and edgeR::calcNormFactors is not claimed: neither is in renv.lock, and the outstanding condition stays recorded in issue #16 and in docs/statistics/method-conditions/scaling-and-offsets.md. Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- CHANGELOG.md | 49 +++ config/schemas/analysis_config.ncl | 44 +- config/schemas/analysis_config.schema.json | 32 +- config/templates/analysis_config_chora.deed | 11 +- frontend/src/types/analysis_config.ts | 60 ++- src/MetaManifold.jl | 5 + src/analysis/AnalysisConfig.jl | 171 ++++++-- src/analysis/Execution.jl | 126 ++++-- src/server/routes/analysis_config.jl | 7 + test/runtests.jl | 1 + test/unit/test_scaling.jl | 423 ++++++++++++++++++++ 11 files changed, 852 insertions(+), 77 deletions(-) create mode 100644 test/unit/test_scaling.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 8b9b6971..a70b2156 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,55 @@ types, tests, infrastructure, and alignment. ## [Unreleased] +### Fixed — `method = "TSS"` was refused by the layer that claimed to have shipped it (2026-09-25) + +- **The allowed-normalisation table held the three names in a different case from the + one the configuration stored.** `NormalizationConfig` canonicalises its method to + lower case (`"TSS"` becomes `"tss"`), while `VALID_NORMALIZATION_FOR_METHOD` listed + `"TSS"`, `"CSS"`, `"RSS"`. The membership test therefore failed for every one of those + spellings, and `AnalysisConfig` raised + `normalization.method 'tss' incompatible with method 'nb_glm'` — including for + `test/unit/test_execution.jl`'s `TSS offset for NB_GLM` testset, which constructed + `method="TSS"` exactly as the documentation instructs. That failure is what reddened + the CI run for the merged TSS-offsets commit; it was not a flaky test. +- Both the table and the two comparisons (`AnalysisConfig` constructor, `validate_config`) + are lower case now, `test/unit/test_scaling.jl` asserts that every admissible spelling + of every method name is accepted, and the JSON schema enum still accepts the upper-case + spellings for existing documents. + +### Added — TSS/CSS/RSS are computed as offsets and recorded as such (2026-09-25) + +- **`src/analysis/scaling.jl` is wired into the execution path** (issue #16). Under + `nb_glm` the response stays the counts and the offset is the declared scaling factor's + logarithm: `tss` = log library size, `css` = log cumulative sum at the declared + `css_quantile` (Paulson et al. 2013), `rss` = log of the trimmed weighted mean of + log-ratios to a reference sample (TMM, Robinson & Oshlack 2010), `size_factors` = + median-of-ratios (DESeq2/RLE) — the estimator the name claimed all along and the code + did not compute. `none` keeps the plain log library size it has always meant. +- **The declared parameters travel with the configuration**: `css_quantile` (default + 0.75), `tmm_ref_column` (default: chosen the way edgeR chooses it, and recorded), + `tmm_log_ratio_trim` (0.3) and `tmm_sum_trim` (0.05) are validated in the constructor, + included in the config hash and canonical JSON, round-tripped through `to_json`/ + `from_json`/Nickel/DEED, mirrored in the JSON schema and the Nickel contracts, exposed + through the server's configuration route, and documented in `context_help` and the + frontend help. A run that asked for a different quantile is a different run. +- **What was computed is recorded**: `diagnostics.checks["scaling"]` and the manifest + provenance carry the kind, the definition sentence, the parameters, the raw quantities + before centring and a SHA-256 of the offset vector, so two runs can be shown to agree + and a run whose offset changed is detectable from the manifest alone. +- **Refusals instead of substitutions**: `css` and `rss` are refused for a response with + no counts to offset (a non-count response is what the old alias silently produced), + a zero or non-finite sample total is refused by name, a CSS cumulative sum of zero is + refused as `log(0)` rather than replaced, and an unknown `tmm_ref_column` is refused + with the real sample names instead of falling back to the data-driven choice. +- Evidence: `test/unit/test_scaling.jl` (known answers with the arithmetic written + beside them, properties a wrong implementation fails, the CSS robustness property, + negative controls, the integration path through `prepare_analysis_table`, and a + base-R transcription of the same conditions), plus the conditions document + `docs/statistics/method-conditions/scaling-and-offsets.md`. Parity with + `metagenomeSeq::cumNorm` and `edgeR::calcNormFactors` is **not** claimed: neither + package is in `renv.lock`, and the outstanding condition is recorded in issue #16. + ### Fixed — the analysis path no longer returns placeholder statistics (2026-09-25) - **`run_analysis` computed nothing and returned numbers anyway.** It derived diff --git a/config/schemas/analysis_config.ncl b/config/schemas/analysis_config.ncl index 1e0b7685..b6fffa29 100644 --- a/config/schemas/analysis_config.ncl +++ b/config/schemas/analysis_config.ncl @@ -3,7 +3,7 @@ # AnalysisConfig Nickel contract — hyperpolymath/standards style Milestone 3 # From 1-formats/k9/*.ncl and .machine_readable/contractiles/_base.ncl # Implements BH mandatory, DANGER banner, advanced validation for pseudocount/epsilon/zero_policy/etc. -# TSS/CSS/RSS deferred as alias to relative with warning, see GitHub issues +# tss/css/rss are exact offsets since 2026-09-25 (issue #16), see docs/statistics/method-conditions/scaling-and-offsets.md let AnalysisMethod = std.enum.TagOrString & [| 'nb_glm, 'clr_lm, 'ilr_lm, 'logistic |] in let NormalizationMethod = std.enum.TagOrString & [| 'none, 'rarefy, 'relative, 'size_factors, 'clr, 'ilr, 'presence_absence, 'TSS, 'CSS, 'RSS, 'tss, 'css, 'rss |] in @@ -40,6 +40,20 @@ let PseudocountContract = fun label value => 'Ok value in +let CssQuantileContract = fun label value => + if value <= 0 || value >= 1 then + 'Error { message = "css_quantile must be in (0,1), got %{std.to_string value}. It is the per-sample quantile of the count distribution whose cumulative sum becomes the CSS scaling factor. See context_help('normalization.css_quantile')" } + else + 'Ok value +in + +let TmmTrimContract = fun label value => + if value < 0 || value >= 0.5 then + 'Error { message = "a TMM trim must be in [0,0.5), got %{std.to_string value}. See context_help('normalization.tmm_log_ratio_trim')" } + else + 'Ok value +in + let EpsilonContract = fun label value => if value <= 0 || value >= 1 then 'Error { message = "epsilon must be in (0,1) for numerical stability, got %{std.to_string value}. Typical 1e-6. See context_help('advanced.epsilon')" } @@ -94,7 +108,7 @@ let MethodNormalizationCompatibility = fun label value => let norm = value.normalization.method in if method == 'nb_glm then if norm == 'clr || norm == 'ilr then - 'Error { message = "NB_GLM expects count data, not CLR/ILR transforms. Use clr_lm/ilr_lm for compositional, or change normalization to none/size_factors/relative/TSS/CSS/RSS. Refusing meaningless combination. See context_help('normalization.method')" } + 'Error { message = "NB_GLM expects count data, not CLR/ILR transforms. Use clr_lm/ilr_lm for compositional, or change normalization to none/size_factors/relative/tss/css/rss. Refusing meaningless combination. See context_help('normalization.method')" } else 'Ok value else if method == 'clr_lm then @@ -182,7 +196,7 @@ in | { method | NormalizationMethod - | doc "Normalization / transform, must be compatible with method. TSS/CSS/RSS deferred alias to relative with warning, see GitHub issue 01.", + | doc "Normalization / transform, must be compatible with method. tss/css/rss are exact offsets (issue #16), not aliases to relative; see docs/statistics/method-conditions/scaling-and-offsets.md.", pseudocount | Number @@ -211,7 +225,29 @@ in tss_css_rss_note | std.option.String - | doc "Note for TSS/CSS/RSS deferred features — currently aliased to relative with warning. See GitHub issue 01.", + | doc "Free-form note recorded with the configuration. Historical name: TSS/CSS/RSS are not deferred any more (issue #16), and this field is kept so existing documents keep parsing.", + + css_quantile + | Number + | CssQuantileContract + | default = 0.75 + | doc "CSS: per-sample quantile whose cumulative sum is the scaling factor. Used only by method='css'. Paulson et al. 2013 use 0.5-0.75. See docs/statistics/method-conditions/scaling-and-offsets.md.", + + tmm_ref_column + | std.option.String + | doc "RSS/TMM: reference sample, or omitted to choose it the way edgeR chooses it (recorded either way). Used only by method='rss'.", + + tmm_log_ratio_trim + | Number + | TmmTrimContract + | default = 0.3 + | doc "RSS/TMM: fraction trimmed from each log-ratio tail before the weighted mean; edgeR default 0.3.", + + tmm_sum_trim + | Number + | TmmTrimContract + | default = 0.05 + | doc "RSS/TMM: fraction trimmed from each mean-abundance tail; edgeR default 0.05.", }, correction diff --git a/config/schemas/analysis_config.schema.json b/config/schemas/analysis_config.schema.json index b31b9689..45a25522 100644 --- a/config/schemas/analysis_config.schema.json +++ b/config/schemas/analysis_config.schema.json @@ -2,7 +2,7 @@ "$schema": "https://json-schema.org/draft/2020-12/schema", "$id": "https://hyperpolymath.github.io/MetaManifold-WebUI/schemas/analysis_config.schema.json", "title": "AnalysisConfig — versioned, explicit, provenance-rich analysis configuration — Milestone 3", - "description": "Safe, explicit, versioned AnalysisConfig layer for parametric and nonparametric analyses (NB GLM, CLR/ILR+Gaussian LM, logistic in v1; BH mandatory; DANGER banner on overrides; Advanced Analysis section heavy validation/help/warnings for custom pseudocount/epsilon/zero_policy/etc.; DOI-ready bundles). From hyperpolymath/standards JSON + Nickel + DEED schemes. TSS/CSS/RSS deferred as alias to relative with warning, see GitHub issues.", + "description": "Safe, explicit, versioned AnalysisConfig layer for parametric and nonparametric analyses (NB GLM, CLR/ILR+Gaussian LM, logistic in v1; BH mandatory; DANGER banner on overrides; Advanced Analysis section heavy validation/help/warnings for custom pseudocount/epsilon/zero_policy/etc.; DOI-ready bundles). From hyperpolymath/standards JSON + Nickel + DEED schemes. tss/css/rss are exact offsets (issue #16, 2026-09-25), not aliases to relative; see docs/statistics/method-conditions/scaling-and-offsets.md.", "type": "object", "required": ["schema_version", "id", "method", "formula", "metadata_columns", "normalization", "correction"], "properties": { @@ -60,7 +60,7 @@ "method": { "type": "string", "enum": ["none", "rarefy", "relative", "size_factors", "clr", "ilr", "presence_absence", "TSS", "CSS", "RSS", "tss", "css", "rss"], - "description": "Normalization / transform. Must be compatible with method: nb_glm allows none/rarefy/size_factors/relative/TSS/CSS/RSS (TSS/CSS/RSS deferred alias to relative with warning), clr_lm requires clr, ilr_lm requires ilr, logistic allows none/relative/rarefy/presence_absence/TSS. See MethodNormalizationCompatibility Nickel contract and GitHub issues for TSS/CSS/RSS exact offsets." + "description": "Normalization / transform. Must be compatible with method: nb_glm allows none/rarefy/size_factors/relative/tss/css/rss, clr_lm requires clr, ilr_lm requires ilr, logistic allows none/relative/rarefy/presence_absence/tss. tss/css/rss are offsets computed exactly (issue #16, docs/statistics/method-conditions/scaling-and-offsets.md), not aliases to relative; css and rss apply to count responses only, because an offset needs something to offset. See MethodNormalizationCompatibility Nickel contract." }, "pseudocount": { "type": "number", @@ -93,7 +93,33 @@ }, "tss_css_rss_note": { "type": ["string", "null"], - "description": "Note for TSS/CSS/RSS deferred features — currently aliased to relative with warning. See GitHub issue 01-tss-css-rss-offsets." + "description": "Free-form note recorded with the configuration. The field name is historical: TSS/CSS/RSS are not deferred any more (issue #16, 2026-09-25) and are exact offsets. Retained so existing documents keep parsing." + }, + "css_quantile": { + "type": "number", + "exclusiveMinimum": 0, + "exclusiveMaximum": 1, + "default": 0.75, + "description": "CSS: the per-sample quantile of the count distribution whose cumulative sum becomes the scaling factor (Paulson et al. 2013 use 0.5-0.75). Used only by normalization.method='css'. Refused at run time when that cumulative sum is zero for any sample, because log(0) is not a small number. metagenomeSeq's data-driven cumNormStatFast choice of the quantile is deliberately not implemented: a parameter chosen from the data must be recorded, not defaulted. See docs/statistics/method-conditions/scaling-and-offsets.md." + }, + "tmm_ref_column": { + "type": ["string", "null"], + "default": null, + "description": "RSS/TMM: the sample every other sample's log-ratios are taken against. Null chooses it the way edgeR chooses it (the sample whose upper-quartile-scaled counts are closest to the mean of those values) and records which one. An unknown name is refused at run time, never silently replaced by the data-driven choice. Used only by normalization.method='rss'." + }, + "tmm_log_ratio_trim": { + "type": "number", + "minimum": 0, + "exclusiveMaximum": 0.5, + "default": 0.3, + "description": "RSS/TMM: fraction trimmed from each tail of the log-ratios before the weighted mean; edgeR's default is 0.3. Zero disables the trimming, which removes the robustness TMM is used for — both values are recorded, so an untrimmed run is visible rather than assumed. Used only by normalization.method='rss'." + }, + "tmm_sum_trim": { + "type": "number", + "minimum": 0, + "exclusiveMaximum": 0.5, + "default": 0.05, + "description": "RSS/TMM: fraction trimmed from each tail of the mean abundances before the weighted mean; edgeR's default is 0.05. Used only by normalization.method='rss'." } }, "allOf": [ diff --git a/config/templates/analysis_config_chora.deed b/config/templates/analysis_config_chora.deed index 945ceafe..961735e3 100644 --- a/config/templates/analysis_config_chora.deed +++ b/config/templates/analysis_config_chora.deed @@ -24,7 +24,11 @@ :zero-policy "pseudocount" :ilr-basis "" :multiplicative-replacement-delta 0 - :tss-css-rss-note "TSS/CSS/RSS deferred alias to relative, see GitHub issue 01") + :tss-css-rss-note "" + :css-quantile 0.75 + :tmm-ref-column "" + :tmm-log-ratio-trim 0.3 + :tmm-sum-trim 0.05) (correction :method "BH" @@ -69,7 +73,7 @@ :method "Analysis method — explicit, no auto-selection v1 nb_glm clr_lm ilr_lm logistic. See JSON schema for scientific context and GitHub issues for deferred multinomial/occupancy/ordination." :formula "R-style formula, e.g. '~ group' or 'disease ~ group + batch'. Must reference only metadata_columns. Forbidden ; backtick dollar. See ValidFormula Nickel contract." :correction "BH mandatory in v1. Any override triggers DANGER banner and requires acknowledgment token I_UNDERSTAND_THE_RISK_AND_WANT_TO_OVERRIDE_BH. Scary banner for paper writers." - :normalization "Normalization must be compatible with method: nb_glm allows none/rarefy/size_factors/relative/TSS/CSS/RSS (TSS/CSS/RSS deferred alias), clr_lm requires clr, ilr_lm requires ilr. Pseudocount/epsilon/zero_policy heavy validation." + :normalization "Normalization must be compatible with method: nb_glm allows none/rarefy/size_factors/relative/tss/css/rss, clr_lm requires clr, ilr_lm requires ilr. tss/css/rss are exact offsets (tss = log library size, css = cumulative sum at css-quantile, rss = trimmed mean of log-ratios to a reference sample) and are refused for non-count responses; see docs/statistics/method-conditions/scaling-and-offsets.md. None of them removes the compositional constraint. Pseudocount/epsilon/zero_policy heavy validation." :advanced "All advanced options behind Advanced Analysis expander, hidden unless Evidence Mode enabled. Heavy validation, refusal of meaningless inputs, warnings for custom pseudocount/epsilon/zero_policy/min_prevalence/max_features/min_samples_per_group.")) ;; SPDX-License-Identifier: AGPL-3.0-only @@ -81,4 +85,5 @@ ;; ;; ║ Config ID: 00000000-0000-0000-0000-000000000000 ║ ;; ;; ║ If you are writing a paper, you MUST disclose these overrides in Methods║ ;; ;; ╚════════════════════════════════════════════════════════════════════════════╝ -;; Deferred features: TSS/CSS/RSS offsets, multinomial/DM, occupancy, constrained ordinations, ILR basis phylogenetic/SBP, glmGamPoi/Bayesian multiplicative — see docs/issues/milestone3/ +;; Deferred features: multinomial/DM, occupancy, constrained ordinations, ILR basis phylogenetic/SBP, glmGamPoi/Bayesian multiplicative — see docs/issues/milestone3/ +;; Implemented since this template was written: TSS/CSS/RSS offsets (issue #16), src/analysis/scaling.jl. diff --git a/frontend/src/types/analysis_config.ts b/frontend/src/types/analysis_config.ts index 8574b767..8a752dcb 100644 --- a/frontend/src/types/analysis_config.ts +++ b/frontend/src/types/analysis_config.ts @@ -19,6 +19,16 @@ export interface NormalizationConfig { ilr_basis?: string | null multiplicative_replacement_delta?: number | null tss_css_rss_note?: string | null + // CSS: quantile of each sample's count distribution whose cumulative sum is the + // scaling factor. Only used by method='css'. Default 0.75 (Paulson et al. 2013). + css_quantile?: number + // RSS/TMM: reference sample every other sample's log-ratios are taken against. + // Omitted means "chosen the way edgeR chooses it", and the choice is recorded. + tmm_ref_column?: string | null + // RSS/TMM: fraction trimmed from each log-ratio tail before the weighted mean (0.3). + tmm_log_ratio_trim?: number + // RSS/TMM: fraction trimmed from each mean-abundance tail (0.05). + tmm_sum_trim?: number } export interface CorrectionConfig { @@ -163,15 +173,55 @@ See JSON schema pattern ^[^;\`$]+$ and Nickel ValidFormula. `, 'normalization.method': `Normalization / Transform (method-dependent) — must be compatible with method -- For NB_GLM: none, size_factors (DESeq2 default, preferred), relative, rarefy (discouraged, use with caution), TSS (alias for relative, deferred exact TSS), CSS (deferred), RSS (deferred) +- For NB_GLM: none, size_factors (median-of-ratios), relative, rarefy (discouraged), tss, css, rss - For CLR_LM: must be clr — Centered Log-Ratio, requires pseudocount >0 - For ILR_LM: must be ilr — Isometric Log-Ratio, requires pseudocount >0 and ilr_basis -- For LOGISTIC: presence_absence, none, relative, rarefy, TSS +- For LOGISTIC: presence_absence, none, relative, rarefy, tss -TSS/CSS/RSS offsets are deferred features (see GitHub issues) — currently aliased to relative. For exact TSS/CSS/RSS offsets, see deferred issue with value/difficulty/risk. +How the count-model choices differ (each is an offset, not a transform; counts stay counts): -Refuses meaningless: NB_GLM + clr/ilr (counts vs compositional), CLR_LM + none, etc. See Nickel MethodNormalizationCompatibility contract. -See JSON schema enum and DEED (normalization :method). +- tss (total sum scaling): offset = log(library size). The McMurdie & Holmes (2014) answer to rarefaction: model depth, do not divide by it. +- css (cumulative sum scaling, Paulson et al. 2013): offset = log of the sum of counts at or below each sample's own quantile, set by css_quantile (default 0.75). Robust to a few dominant taxa. Refused when that sum is zero for any sample. +- rss (TMM, Robinson & Oshlack 2010): offset = log of a trimmed weighted mean of log-ratios against a reference sample; set by tmm_log_ratio_trim (0.3), tmm_sum_trim (0.05) and tmm_ref_column. Assumes most features are not differentially abundant. +- size_factors: median-of-ratios (DESeq2/RLE), in the offset form. +- relative: proportions. A transform, not an offset: it discards the count nature of the data. +- none: no scaling; for a count model the offset is the plain log library size. + +None of these removes the compositional constraint. They correct for sequencing depth; a log fold change from a model with these offsets is still relative to the sampled community. css and rss are refused for a response with no counts to offset. + +Refuses meaningless: NB_GLM + clr/ilr (counts vs compositional), CLR_LM + none, css/rss on a non-count response, etc. See Nickel MethodNormalizationCompatibility contract. +See docs/statistics/method-conditions/scaling-and-offsets.md, JSON schema enum and DEED (normalization :method). +`, + 'normalization.css_quantile': `CSS quantile (normalization.method = css only) + +- The per-sample quantile of the count distribution whose cumulative sum becomes the scaling factor. Paulson et al. (2013) use 0.5-0.75; the default here is 0.75. +- Must be in (0,1). Below 0.5 the cumulative sum covers less than half of a sample's counts and is dominated by how many features are zero — a warning is issued. +- Refused at run time when that cumulative sum is zero for any sample (log(0) is not a small number). Raise the quantile or exclude the sample. +- metagenomeSeq's data-driven choice of this quantile (cumNormStatFast) is deliberately not implemented: a parameter chosen from the data is a decision the run has to record. + +Recorded in the config hash and in the manifest provenance. +`, + 'normalization.tmm_ref_column': `TMM reference sample (normalization.method = rss only) + +- Names the sample every other sample's log-ratios are taken against. +- Omitted (the default): the reference is chosen the way edgeR chooses it — the sample whose upper-quartile-scaled counts are closest to the mean of those values. Which sample was chosen is recorded in the provenance. +- The name is not checked until the run has the sample list; an unknown name is refused there, by name, rather than silently falling back to the data-driven choice. + +Recorded in the config hash. +`, + 'normalization.tmm_log_ratio_trim': `TMM log-ratio trimming (normalization.method = rss only) + +- Fraction trimmed from each tail of the log-ratios before the weighted mean: 0.3 by default, as in edgeR. +- Must be in [0,0.5). Zero means no trimming, which removes the robustness TMM is used for; both values are recorded, so an untrimmed run is visible rather than assumed. + +Recorded in the config hash. +`, + 'normalization.tmm_sum_trim': `TMM abundance trimming (normalization.method = rss only) + +- Fraction trimmed from each tail of the mean abundances before the weighted mean: 0.05 by default, as in edgeR. +- Must be in [0,0.5). Zero means no trimming of the abundance tails. + +Recorded in the config hash. `, 'normalization.pseudocount': `Pseudocount for zero replacement (CLR/ILR mandatory, NB_GLM optional but warned) diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index 4f609440..e8c4b006 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -37,6 +37,11 @@ include("analysis/AnalysisConfig.jl") include("analysis/clade_cumulus.jl") # Estimation fits what AnalysisConfig declares and Execution runs it, so it sits between them. include("analysis/estimation.jl") +# Library-size scaling factors and the offsets count models are fitted with. Included +# before Execution, which is the only consumer: prepare_analysis_table decides what the +# response is and which offset it gets, and Scaling computes both under the conditions +# published in docs/statistics/method-conditions/scaling-and-offsets.md. +include("analysis/scaling.jl") include("analysis/Execution.jl") end diff --git a/src/analysis/AnalysisConfig.jl b/src/analysis/AnalysisConfig.jl index cca4176f..5f812c4f 100644 --- a/src/analysis/AnalysisConfig.jl +++ b/src/analysis/AnalysisConfig.jl @@ -98,11 +98,18 @@ const ZERO_POLICY_STRINGS = Dict{String,ZeroPolicy}( const VALID_DISPERSION_METHODS = ("parametric", "local", "mean", "pooled", "glmGamPoi") const VALID_ZERO_HANDLING = ("pseudocount", "multiplicative_replacement", "bayesian_multiplicative", "refuse") const VALID_ILR_BASIS = ("default", "phylogenetic", "sequential_binary_partition", "balance_dendrogram") +# Method names are canonicalised to lower case by `NormalizationConfig`, so the allowed +# names are held in the same case and compared in it. They were not always: the entries +# for TSS/CSS/RSS arrived upper case while the constructor stored `"tss"`, so every one of +# those combinations raised "normalization.method ... incompatible" from `AnalysisConfig` +# -- including the exact TSS offset the CHANGELOG claimed had shipped. The comparison is +# now case-insensitive in both directions so neither side can drift again, and the tests +# assert that each admissible spelling of every method name is accepted. const VALID_NORMALIZATION_FOR_METHOD = Dict{AnalysisMethod, Vector{String}}( - NB_GLM => ["none", "rarefy", "size_factors", "relative", "TSS", "CSS", "RSS"], # TSS/CSS/RSS deferred but allowed as alias for relative + NB_GLM => ["none", "rarefy", "size_factors", "relative", "tss", "css", "rss"], CLR_LM => ["clr"], ILR_LM => ["ilr"], - LOGISTIC => ["none", "relative", "rarefy", "presence_absence", "TSS"], + LOGISTIC => ["none", "relative", "rarefy", "presence_absence", "tss"], ) # -------------------------------------------------------------------------- @@ -125,7 +132,11 @@ struct NormalizationConfig zero_policy::ZeroPolicy ilr_basis::Union{String,Nothing} multiplicative_replacement_delta::Union{Float64,Nothing} - tss_css_rss_note::Union{String,Nothing} # deferred feature note + tss_css_rss_note::Union{String,Nothing} # free-form note, recorded with the config + css_quantile::Float64 # CSS: quantile of each sample's counts + tmm_ref_column::Union{String,Nothing} # RSS/TMM: reference sample, or nothing + tmm_log_ratio_trim::Float64 # RSS/TMM: trim of the log-ratio tail + tmm_sum_trim::Float64 # RSS/TMM: trim of the abundance tail function NormalizationConfig(; method::String, @@ -134,7 +145,11 @@ struct NormalizationConfig zero_policy::String="pseudocount", ilr_basis::Union{String,Nothing}=nothing, multiplicative_replacement_delta::Union{Float64,Nothing}=nothing, - tss_css_rss_note::Union{String,Nothing}=nothing + tss_css_rss_note::Union{String,Nothing}=nothing, + css_quantile::Real=0.75, + tmm_ref_column::Union{String,Nothing}=nothing, + tmm_log_ratio_trim::Real=0.3, + tmm_sum_trim::Real=0.05 ) method_clean = lowercase(strip(method)) isempty(method_clean) && throw(ArgumentError("normalization.method must be non-empty (e.g. 'clr', 'size_factors', 'TSS') — see context_help('normalization.method')")) @@ -189,14 +204,43 @@ struct NormalizationConfig (delta <= 0 || delta >= 1) && throw(ArgumentError("multiplicative_replacement_delta must be in (0,1), got $delta — see context_help('advanced.zero_policy')")) end - # TSS/CSS/RSS note — deferred feature - if method_clean in ("tss", "css", "rss") - if isnothing(tss_css_rss_note) - @warn "TSS/CSS/RSS are deferred features (see GitHub issues). Currently aliased to relative/TSS. For exact TSS/CSS/RSS offsets, see deferred issue with value/difficulty/risk." + # TSS/CSS/RSS are exact offsets as of 2026-09-25 (issue #16), implemented in + # src/analysis/scaling.jl. Their parameters are validated here and recorded in the + # config hash, because each one changes the numbers a count model reports. + css_q = Float64(css_quantile) + if !(0 < css_q < 1) + throw(ArgumentError("css_quantile must be in (0,1), got $css_q — it is the quantile of each sample's counts that the cumulative sum is taken up to. See context_help('normalization.css_quantile') and docs/statistics/method-conditions/scaling-and-offsets.md")) + end + if css_q < 0.5 + @warn "css_quantile=$css_q 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." css_quantile + end + lrt = Float64(tmm_log_ratio_trim) + if !(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 context_help('normalization.tmm_log_ratio_trim')")) + end + st = Float64(tmm_sum_trim) + if !(0 <= st < 0.5) + throw(ArgumentError("tmm_sum_trim must be in [0,0.5), got $st — see context_help('normalization.tmm_sum_trim')")) + end + if lrt == 0.0 && st == 0.0 + @warn "tmm_log_ratio_trim=0 and tmm_sum_trim=0: the trimmed mean becomes an untrimmed mean over every feature, which is the robustness the TMM estimator is used for. Proceeding, and recording both values." tmm_log_ratio_trim tmm_sum_trim + end + ref_col = isnothing(tmm_ref_column) ? nothing : strip(tmm_ref_column) + if !isnothing(ref_col) + isempty(ref_col) && throw(ArgumentError("tmm_ref_column must name a sample, or be omitted — got an empty string. See context_help('normalization.tmm_ref_column')")) + occursin(r"[;`$]", ref_col) && throw(ArgumentError("tmm_ref_column contains a forbidden ; ` \$ — got '$ref_col'")) + end + if !(method_clean in ("css", "rss")) + if css_q != 0.75 + @warn "css_quantile is only used by normalization.method='css'; it is recorded but has no effect for '$method_clean'." + end + if !isnothing(ref_col) || lrt != 0.3 || st != 0.05 + @warn "tmm_ref_column/tmm_log_ratio_trim/tmm_sum_trim are only used by normalization.method='rss'; they are recorded but have no effect for '$method_clean'." end end - new(method_clean, pseudocount, epsilon, zp, ilr_basis, multiplicative_replacement_delta, tss_css_rss_note) + new(method_clean, pseudocount, epsilon, zp, ilr_basis, multiplicative_replacement_delta, + tss_css_rss_note, css_q, ref_col, lrt, st) end end @@ -479,8 +523,11 @@ struct AnalysisConfig end end - # Normalization compatibility with method - norm_method = normalization.method + # Normalization compatibility with method. Both sides are lowercased here: the + # table is the canonical list and the config has already canonicalised its own + # name, and a comparison that depends on which side was canonicalised is how the + # "tss" vs "TSS" mismatch shipped. + norm_method = lowercase(strip(normalization.method)) allowed_norms = get(VALID_NORMALIZATION_FOR_METHOD, method_enum, String[]) if !(norm_method in allowed_norms) throw(ArgumentError("normalization.method '$norm_method' incompatible with method '$(METHOD_TO_STRING[method_enum])'. Allowed for $(METHOD_TO_STRING[method_enum]): $(join(allowed_norms, ", ")). See context_help('normalization.method') and MethodNormalizationCompatibility contract in Nickel. Refusing meaningless combination.")) @@ -537,7 +584,11 @@ struct AnalysisConfig "epsilon" => normalization.epsilon, "zero_policy" => string(normalization.zero_policy), "ilr_basis" => normalization.ilr_basis, - "multiplicative_replacement_delta" => normalization.multiplicative_replacement_delta + "multiplicative_replacement_delta" => normalization.multiplicative_replacement_delta, + "css_quantile" => normalization.css_quantile, + "tmm_ref_column" => normalization.tmm_ref_column, + "tmm_log_ratio_trim" => normalization.tmm_log_ratio_trim, + "tmm_sum_trim" => normalization.tmm_sum_trim ), "correction" => OrderedDict( "method" => correction.method, @@ -684,7 +735,7 @@ function validate_config(config::AnalysisConfig, available_columns::Vector{Strin # Check normalization compatibility already done in constructor, but re-check for strict if strict - norm_method = config.normalization.method + norm_method = lowercase(strip(config.normalization.method)) allowed = get(VALID_NORMALIZATION_FOR_METHOD, config.method, String[]) if !(norm_method in allowed) push!(errors, "Incompatible normalization.method '$norm_method' for method '$(METHOD_TO_STRING[config.method])' — allowed $(join(allowed, ", "))") @@ -726,7 +777,11 @@ function canonical_json(config::AnalysisConfig) "pseudocount" => config.normalization.pseudocount, "epsilon" => config.normalization.epsilon, "zero_policy" => string(config.normalization.zero_policy), - "ilr_basis" => config.normalization.ilr_basis + "ilr_basis" => config.normalization.ilr_basis, + "css_quantile" => config.normalization.css_quantile, + "tmm_ref_column" => config.normalization.tmm_ref_column, + "tmm_log_ratio_trim" => config.normalization.tmm_log_ratio_trim, + "tmm_sum_trim" => config.normalization.tmm_sum_trim ), "correction" => OrderedDict( "method" => config.correction.method, @@ -805,15 +860,69 @@ function context_help(field_path::String) "normalization.method" => """ Normalization / Transform (method-dependent) — must be compatible with method - - For NB_GLM: none, size_factors (DESeq2 default, preferred), relative, rarefy (discouraged, use with caution), TSS (alias for relative, deferred exact TSS), CSS (deferred), RSS (deferred) + - For NB_GLM: none, size_factors (DESeq2/RLE median-of-ratios), relative, rarefy (discouraged), tss, css, rss - For CLR_LM: must be clr — Centered Log-Ratio, requires pseudocount >0 - For ILR_LM: must be ilr — Isometric Log-Ratio, requires pseudocount >0 and ilr_basis - - For LOGISTIC: presence_absence, none, relative, rarefy, TSS + - For LOGISTIC: presence_absence, none, relative, rarefy, tss + + How the count-model choices differ (each is an offset, not a transform; counts stay counts): + + - tss (total sum scaling): offset = log(library size). The McMurdie & Holmes (2014) answer to rarefaction: model depth, do not divide by it. + - css (cumulative sum scaling, Paulson et al. 2013): offset = log of the sum of counts at or below each sample's own quantile, set by css_quantile (default 0.75). Robust to a few dominant taxa. Refused when that sum is zero for any sample. + - rss (TMM, Robinson & Oshlack 2010): offset = log of a trimmed weighted mean of log-ratios against a reference sample; set by tmm_log_ratio_trim (0.3), tmm_sum_trim (0.05) and tmm_ref_column. Assumes most features are not differentially abundant. + - size_factors: median-of-ratios (DESeq2/RLE) as a size factor, in the offset form. + - relative: proportions. A transform, not an offset: it discards the count nature of the data. + - none: no scaling; the offset is the plain log library size. + + None of these removes the compositional constraint. They correct for sequencing depth; a log fold change from a model with these offsets is still relative to the sampled community. + + Refuses meaningless: NB_GLM + clr/ilr (counts vs compositional), CLR_LM + none, css/rss on a non-count response (an offset needs something to offset), etc. See Nickel MethodNormalizationCompatibility contract. + See docs/statistics/method-conditions/scaling-and-offsets.md and JSON schema enum. + """, + "normalization.css_quantile" => """ + CSS quantile (normalization.method = css only) + + - The per-sample quantile of the count distribution whose cumulative sum becomes the + scaling factor. Paulson et al. (2013) use 0.5-0.75; the default here is 0.75. + - Must be in (0,1). Below 0.5 the cumulative sum covers less than half of a sample's + counts and is dominated by how many features are zero — a warning is issued. + - Refused at run time when the cumulative sum is zero for any sample (log(0) is not a + small number). Raise the quantile or exclude the sample. + - metagenomeSeq's data-driven choice of this quantile (cumNormStatFast) is deliberately + not implemented: a parameter chosen from the data is a decision the run must record. + + Recorded in the config hash and in the manifest provenance. + """, + "normalization.tmm_ref_column" => """ + TMM reference sample (normalization.method = rss only) + + - Names the sample every other sample's log-ratios are taken against. + - Omitted (the default): the reference is chosen the way edgeR chooses it — the sample + whose upper-quartile-scaled counts are closest to the mean of those values. Which + sample was chosen is recorded in the provenance. + - The name is not checked until the run has the sample list; an unknown name is refused + there, by name, rather than silently falling back to the data-driven choice. + + Recorded in the config hash. + """, + "normalization.tmm_log_ratio_trim" => """ + TMM log-ratio trimming (normalization.method = rss only) + + - Fraction trimmed from each tail of the log-ratios before the weighted mean: 0.3 by + default, as in edgeR. + - Must be in [0,0.5). Zero means no trimming, which removes the robustness TMM is used + for; both values are recorded, so an untrimmed run is visible rather than assumed. + + Recorded in the config hash. + """, + "normalization.tmm_sum_trim" => """ + TMM abundance trimming (normalization.method = rss only) - TSS/CSS/RSS offsets are deferred features (see GitHub issues) — currently aliased to relative. For exact TSS/CSS/RSS offsets, see deferred issue with value/difficulty/risk. + - Fraction trimmed from each tail of the mean abundances before the weighted mean: 0.05 + by default, as in edgeR. + - Must be in [0,0.5). Zero means no trimming of the abundance tails. - Refuses meaningless: NB_GLM + clr/ilr (counts vs compositional), CLR_LM + none, etc. See Nickel MethodNormalizationCompatibility contract. - See JSON schema enum and DEED (normalization :method). + Recorded in the config hash. """, "normalization.pseudocount" => """ Pseudocount for zero replacement (CLR/ILR mandatory, NB_GLM optional but warned) @@ -1077,7 +1186,11 @@ function to_json(config::AnalysisConfig) "epsilon" => config.normalization.epsilon, "zero_policy" => string(config.normalization.zero_policy), "ilr_basis" => config.normalization.ilr_basis, - "multiplicative_replacement_delta" => config.normalization.multiplicative_replacement_delta + "multiplicative_replacement_delta" => config.normalization.multiplicative_replacement_delta, + "css_quantile" => config.normalization.css_quantile, + "tmm_ref_column" => config.normalization.tmm_ref_column, + "tmm_log_ratio_trim" => config.normalization.tmm_log_ratio_trim, + "tmm_sum_trim" => config.normalization.tmm_sum_trim ), "correction" => OrderedDict( "method" => config.correction.method, @@ -1121,7 +1234,11 @@ function from_json(json_str::String) epsilon=get(norm_data, :epsilon, 1e-6), zero_policy=get(norm_data, :zero_policy, "pseudocount"), ilr_basis=get(norm_data, :ilr_basis, nothing), - multiplicative_replacement_delta=get(norm_data, :multiplicative_replacement_delta, nothing) + multiplicative_replacement_delta=get(norm_data, :multiplicative_replacement_delta, nothing), + css_quantile=Float64(get(norm_data, :css_quantile, 0.75)), + tmm_ref_column=get(norm_data, :tmm_ref_column, nothing), + tmm_log_ratio_trim=Float64(get(norm_data, :tmm_log_ratio_trim, 0.3)), + tmm_sum_trim=Float64(get(norm_data, :tmm_sum_trim, 0.05)) ) corr_data = data.correction @@ -1206,6 +1323,10 @@ function to_nickel(config::AnalysisConfig) epsilon = $(config.normalization.epsilon), zero_policy = '$(string(config.normalization.zero_policy))', ilr_basis = $(isnothing(config.normalization.ilr_basis) ? "null" : "'$(config.normalization.ilr_basis)'"), + css_quantile = $(config.normalization.css_quantile) | CssQuantileContract, + tmm_ref_column = $(isnothing(config.normalization.tmm_ref_column) ? "null" : "\"$(config.normalization.tmm_ref_column)\""), + tmm_log_ratio_trim = $(config.normalization.tmm_log_ratio_trim) | TmmTrimContract, + tmm_sum_trim = $(config.normalization.tmm_sum_trim) | TmmTrimContract, } | MethodNormalizationCompatibility, correction = { @@ -1327,7 +1448,11 @@ function to_deed(config::AnalysisConfig) :epsilon $(config.normalization.epsilon) :zero-policy "$(string(config.normalization.zero_policy))" :ilr-basis "$(isnothing(config.normalization.ilr_basis) ? "" : config.normalization.ilr_basis)" - :multiplicative-replacement-delta $(isnothing(config.normalization.multiplicative_replacement_delta) ? "0" : string(config.normalization.multiplicative_replacement_delta))) + :multiplicative-replacement-delta $(isnothing(config.normalization.multiplicative_replacement_delta) ? "0" : string(config.normalization.multiplicative_replacement_delta)) + :css-quantile $(config.normalization.css_quantile) + :tmm-ref-column "$(isnothing(config.normalization.tmm_ref_column) ? "" : config.normalization.tmm_ref_column)" + :tmm-log-ratio-trim $(config.normalization.tmm_log_ratio_trim) + :tmm-sum-trim $(config.normalization.tmm_sum_trim)) (correction :method "$(config.correction.method)" @@ -1372,7 +1497,7 @@ function to_deed(config::AnalysisConfig) :method "Analysis method explicit no auto-selection v1 nb_glm clr_lm ilr_lm logistic" :formula "R-style formula e.g. ~ group must reference only metadata_columns forbids ; backtick dollar" :correction "BH mandatory in v1 any override triggers DANGER banner requires acknowledgment token $(DANGER_ACK_TOKEN)" - :normalization "Normalization must be compatible with method nb_glm allows none/rarefy/size_factors/relative/TSS/CLR/ILR" + :normalization "Normalization must be compatible with method: nb_glm allows none/rarefy/size_factors/relative/tss/css/rss (tss/css/rss are exact offsets, see docs/statistics/method-conditions/scaling-and-offsets.md), clr_lm requires clr, ilr_lm requires ilr" :advanced "All advanced options behind Advanced Analysis expander hidden unless Evidence Mode heavy validation refusal meaningless")) """ end diff --git a/src/analysis/Execution.jl b/src/analysis/Execution.jl index 0cf6336d..2376bcad 100644 --- a/src/analysis/Execution.jl +++ b/src/analysis/Execution.jl @@ -22,9 +22,14 @@ Design: - Counts: Matrix{Float64} where rows=taxa, cols=samples (or transposed? We use rows=taxa, cols=samples as per DESeq2 style) - Sample metadata: Dict or DataFrame-like (we use OrderedDict and DataFrames if available) - Taxa metadata: optional with avec_fibre, epistemic_status, residual_count -- Transforms: none, relative, size_factors, clr, ilr, presence_absence, rarefy (discouraged), TSS/CSS/RSS deferred alias to relative +- Transforms: none, relative, clr, ilr, presence_absence, rarefy (discouraged) +- Offsets (counts stay counts): tss = log library size, css = cumulative sum at the declared + quantile, rss (TMM) = trimmed mean of log-ratios to a reference sample, size_factors = + median-of-ratios. See src/analysis/scaling.jl and + docs/statistics/method-conditions/scaling-and-offsets.md - Zero policies: pseudocount (default safe), multiplicative_replacement, bayesian_multiplicative, refuse (DANGEROUS) -- Offsets: log library size or size_factors for NB_GLM +- Offsets: log library size or size factors for count responses; scaling factors and + their provenance are recorded in `diagnostics.checks["scaling"]` and the manifest - Epistemic filtering: avec_fibre true, epistemic_status present_in_every_admissible_world, min_prevalence, min_abundance, max_features - Self-diagnostics: NaN/Inf, zero variance, all-zero samples/taxa, library size outliers, batch confounding, prevalence/abundance, etc. - Safe self-healing: heal NaN/Inf with epsilon, drop all-zero samples/taxa with warning, record healing in diagnostics, never silent @@ -38,6 +43,7 @@ using UUIDs using JSON3 using OrderedCollections using Logging +import ..Scaling using Statistics # Use AnalysisConfig from parent module @@ -1045,52 +1051,83 @@ function prepare_analysis_table( # Transform and offset # ---------------------------------------------------------------------- - transform_method = config.normalization.method + # The scaling factors and offsets are computed in Scaling, held to the conditions + # published in docs/statistics/method-conditions/scaling-and-offsets.md. This block + # decides only what the response is and which offset it gets. + transform_method = lowercase(strip(config.normalization.method)) prepared = copy(counts_after_zero) offset = nothing + scaling = nothing + count_response = config.method == AnalysisConfig.NB_GLM + + # Proportions of each sample. Used by `relative`, and by `tss` on a non-count + # response, where the total-sum transform and the proportional transform are the same + # operation and the analyst's declared depth handling is recorded in the provenance. + function _proportions_of(counts_matrix, sample_ids_local, config_local) + lib = vec(sum(counts_matrix, dims=1)) + out = zeros(Float64, size(counts_matrix)) + for j in 1:size(counts_matrix, 2) + if lib[j] > 0 + out[:, j] = counts_matrix[:, j] ./ lib[j] + else + sample_name = j <= length(sample_ids_local) ? sample_ids_local[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_local.id).")) + end + end + return out + end if transform_method == "none" - # No transform, keep counts as is + # No transform: the response stays the counts. For a count model the offset is the + # log library size, which is what `none` has always meant in offset form; for + # anything else there is no offset at all. prepared = counts_after_zero - # For NB_GLM, offset is log library size or size_factors - if config.method == AnalysisConfig.NB_GLM - lib_sizes = vec(sum(counts_after_zero, dims=1)) - offset = log.(lib_sizes .+ effective_epsilon) + if count_response + scaling = Scaling.tss_factors(counts_after_zero; sample_ids = filtered_sample_ids) + offset = scaling.offset end - elseif transform_method in ("relative", "tss", "css", "rss", "TSS", "CSS", "RSS") - # Relative abundance: counts / sum per sample - # TSS/CSS/RSS currently alias to relative with warning (deferred exact) - if transform_method in ("tss", "TSS", "css", "CSS", "rss", "RSS") - @warn "TSS/CSS/RSS are deferred features (see GitHub issues 01-tss-css-rss-offsets). Currently aliased to relative with warning. For exact TSS/CSS/RSS offsets, see deferred issue." + elseif transform_method == "tss" + if count_response + # Total sum scaling as an offset: counts stay counts, depth is modelled + # (McMurdie & Holmes 2014) instead of divided out. + prepared = counts_after_zero + scaling = Scaling.tss_factors(counts_after_zero; sample_ids = filtered_sample_ids) + offset = scaling.offset + else + prepared = _proportions_of(counts_after_zero, filtered_sample_ids, config) end - lib_sizes = vec(sum(counts_after_zero, dims=1)) - # For NB_GLM with TSS, preserve raw counts as response and supply log(lib_sizes) as offset, - # preserving count distribution (McMurdie & Holmes 2014) instead of converting to proportions. - if config.method == AnalysisConfig.NB_GLM && transform_method in ("tss", "TSS") - prepared = copy(counts_after_zero) - offset = log.(lib_sizes .+ effective_epsilon) + elseif transform_method in ("css", "rss") + # Offsets for a count model. A non-count response has nothing to offset, and + # silently handing it proportions under the name `css` is exactly the substitution + # this work removed; the configuration layer refuses the pair as well, and this is + # the second door. + count_response || throw(ArgumentError( + "normalization.method='$transform_method' is an offset for a count model, and " * + "method '$(AnalysisConfig.METHOD_TO_STRING[config.method])' has no counts to offset. " * + "Use 'relative' for proportions, or 'clr'/'ilr' for a compositional transform. " * + "Nothing was computed.")) + prepared = counts_after_zero + scaling = if transform_method == "css" + Scaling.css_factors(counts_after_zero; + quantile = config.normalization.css_quantile, + sample_ids = filtered_sample_ids) else - for j in 1:size(counts_after_zero, 2) - if lib_sizes[j] > 0 - prepared[:, j] = counts_after_zero[:, j] ./ lib_sizes[j] - else - 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 + Scaling.tmm_factors(counts_after_zero; + ref_column = config.normalization.tmm_ref_column, + log_ratio_trim = config.normalization.tmm_log_ratio_trim, + sum_trim = config.normalization.tmm_sum_trim, + sample_ids = filtered_sample_ids) end + offset = scaling.offset + elseif transform_method == "relative" + prepared = _proportions_of(counts_after_zero, filtered_sample_ids, config) elseif transform_method == "size_factors" - # Size factors: DESeq2 median-of-ratios (stub: use median ratio or simple) - # For stub, compute size_factors as median of counts / geometric mean per taxon - # Simplified: size_factors = colSums / mean(colSums) or median ratio - lib_sizes = vec(sum(counts_after_zero, dims=1)) - # Geometric mean per taxon (row) across samples, ignoring zeros - # For stub, use lib_sizes / exp(mean(log(lib_sizes))) - geo_mean_lib = exp(mean(log.(lib_sizes .+ effective_epsilon))) - size_factors = lib_sizes ./ geo_mean_lib - # For NB_GLM, we keep counts but store size_factors as offset? Actually DESeq2 uses size_factors as normalization, not offset, but for GLM we can use log(size_factors) as offset - offset = log.(size_factors .+ effective_epsilon) - prepared = counts_after_zero # keep counts, offset stored separately + # Median-of-ratios (RLE), in the offset form. This used to compute library size + # divided by its own geometric mean -- that is `tss`, and calling it DESeq2's + # size factor was the kind of substitution this repository refuses. + scaling = Scaling.rle_factors(counts_after_zero; sample_ids = filtered_sample_ids) + offset = scaling.offset + prepared = counts_after_zero # counts stay counts, offset stored separately elseif transform_method == "clr" # Centered Log-Ratio: log(x) - mean(log(x)) per sample # Requires pseudocount>0 already applied @@ -1189,6 +1226,16 @@ function prepare_analysis_table( # ---------------------------------------------------------------------- (checks, warnings, errors) = self_diagnostics(filtered_counts, prepared, config; sample_metadata=sample_metadata, raw_counts=counts) + + # Record what was computed, and only what was computed: a run that declared `css` gets + # a `scaling` entry with the quantile, the thresholds and the cumulative sums; a run + # that declared `none` gets the plain log library size, labelled as such. This is + # written *after* self_diagnostics, which rebuilds `checks` and `warnings` from + # scratch — put before it, the entry would have been silently overwritten. + if !isnothing(scaling) + checks["scaling"] = Scaling.factor_checks(scaling) + append!(warnings, scaling.warnings) + end checks["all_zero_taxa"] = OrderedDict{String,Any}( "all_zero_taxa_indices" => all_zero_taxa_indices, "all_zero_taxa_count" => length(all_zero_taxa_indices), @@ -1352,7 +1399,8 @@ function prepare_analysis_table( "transform" => transform_method, "zero_policy" => string(zero_policy), "pseudocount" => effective_pseudocount, - "epsilon" => effective_epsilon + "epsilon" => effective_epsilon, + "scaling" => isnothing(scaling) ? nothing : Scaling.factor_provenance(scaling) ) ) diff --git a/src/server/routes/analysis_config.jl b/src/server/routes/analysis_config.jl index dd2d888e..7f0f1cf4 100644 --- a/src/server/routes/analysis_config.jl +++ b/src/server/routes/analysis_config.jl @@ -71,6 +71,13 @@ end pseudocount=Float64(get(norm_body, "pseudocount", 0.5)), ilr_basis=get(norm_body, "ilr_basis", nothing) isa Nothing ? nothing : String(get(norm_body, "ilr_basis", nothing)), multiplicative_replacement_delta=get(norm_body, "multiplicative_replacement_delta", nothing) isa Nothing ? nothing : Float64(get(norm_body, "multiplicative_replacement_delta", nothing)), + # The declared scaling parameters travel with the request, are validated by the + # same constructor the rest of the configuration uses, and end up in the config + # hash — a run that asked for a different CSS quantile is a different run. + css_quantile=Float64(get(norm_body, "css_quantile", 0.75)), + tmm_ref_column=get(norm_body, "tmm_ref_column", nothing) isa Nothing ? nothing : String(get(norm_body, "tmm_ref_column", nothing)), + tmm_log_ratio_trim=Float64(get(norm_body, "tmm_log_ratio_trim", 0.3)), + tmm_sum_trim=Float64(get(norm_body, "tmm_sum_trim", 0.05)), ) corr_body = get(body, "correction", Dict{String,Any}()) diff --git a/test/runtests.jl b/test/runtests.jl index a1e4f1ea..8eef0f69 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,7 @@ using MetaManifold: AnalysisConfig include("unit/test_migrate_composition.jl") include("unit/test_analysis_config.jl") include("unit/test_execution.jl") + include("unit/test_scaling.jl") include("unit/test_estimation.jl") ## Integration tests (opt-in) diff --git a/test/unit/test_scaling.jl b/test/unit/test_scaling.jl new file mode 100644 index 00000000..97c1e7b3 --- /dev/null +++ b/test/unit/test_scaling.jl @@ -0,0 +1,423 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# +# Evidence for issue #16 — TSS, CSS and RSS/TMM as exact offsets — held to the conditions +# published before the implementation in +# docs/statistics/method-conditions/scaling-and-offsets.md: +# +# * known answers whose arithmetic is written next to the expected number, so a reader +# can check the implementation with a pencil rather than against a previous run; +# * an independent-in-language transcription in R (base R only) for TMM, CSS and RLE, +# which catches implementation slips — it is NOT a reference implementation from the +# field, and it does not catch a shared specification error; +# * properties a wrong implementation fails: invariance to permuting features, a sample +# that is an exact multiple of another, and the outlier robustness CSS exists for; +# * negative controls: every refusal is asserted to fire, and to name what it refused; +# * the integration path, from AnalysisConfig through prepare_analysis_table to the +# offset the manifest hands the fit. +# +# What is NOT claimed here: parity with metagenomeSeq::cumNorm or edgeR::calcNormFactors. +# 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; the outstanding condition is recorded +# in issue #16 and in the conditions document instead of being asserted as met. + +@testset "Scaling — TSS, CSS, RSS/TMM and size factors (issue #16)" begin + + using OrderedCollections + using RCall + using MetaManifold.Scaling + using MetaManifold.RRuntime: with_r_lock + + # ------------------------------------------------------------------ + # The table the known answers are computed on + # + # columns: [10 20 40; [60 ] column sums + # 20 20 10; [100] — deliberately not sorted, so a + # 30 60 50] [100] column-sum implementation that assumed + # order would be caught + # ------------------------------------------------------------------ + table_a = [10.0 20 40; + 20.0 20 10; + 30.0 60 50] + ids_a = ["s1", "s2", "s3"] + + @testset "TSS — library sizes, centred" begin + outcome = tss_factors(table_a; sample_ids = ids_a) + # 60, 100, 100; geometric mean = (60*100*100)^(1/3) = 84.34326653017494 + @test outcome.raw ≈ [60.0, 100.0, 100.0] atol = 1e-12 + @test outcome.factors ≈ [0.711378660898012, 1.185631101496687, 1.185631101496687] atol = 1e-10 + @test outcome.offset ≈ log.(outcome.factors) atol = 1e-15 + @test exp(sum(log, outcome.factors) / 3) ≈ 1.0 atol = 1e-12 + @test outcome.kind == "tss" + @test occursin("library size", outcome.definition) + @test outcome.parameters["library_sizes"] == [60.0, 100.0, 100.0] + end + + @testset "CSS — the cumulative sum at the declared quantile" begin + outcome = css_factors(table_a; quantile = 0.75, sample_ids = ids_a) + # column 1: quantile([10,20,30], 0.75) = 25 (type 7: 20 + 0.5*(30-20)); sum of + # counts at or below 25 = 10 + 20 = 30 + # column 2: quantile([20,20,60], 0.75) = 40; sum = 20 + 20 = 40 + # column 3: quantile([10,40,50], 0.75) = 45; sum = 10 + 40 = 50 + @test outcome.raw ≈ [30.0, 40.0, 50.0] atol = 1e-12 + @test outcome.parameters["thresholds"] ≈ [25.0, 40.0, 45.0] atol = 1e-12 + @test outcome.parameters["quantile"] == 0.75 + # geometric mean of [30,40,50] = 39.14867641168864 + @test outcome.factors ≈ [0.766309432393553, 1.021745909858071, 1.277182387322588] atol = 1e-10 + @test outcome.kind == "css" + @test isempty(outcome.warnings) + end + + @testset "CSS is robust to a dominant feature, and TSS is not" begin + # sample 2 has one taxon carrying almost everything + table_c = [10.0 10; + 10.0 10; + 10.0 1000] + tss = tss_factors(table_c) + css = css_factors(table_c; quantile = 0.75) + # thresholds: column 1 = 10 (all counts equal), column 2 = 505 (type 7 between + # the 2nd and 3rd of [10,10,1000]), so the second cumulative sum is 10 + 10 = 20 + @test css.raw ≈ [30.0, 20.0] atol = 1e-12 + @test log(css.factors[2] / css.factors[1]) ≈ log(20 / 30) atol = 1e-12 + @test log(tss.factors[2] / tss.factors[1]) ≈ log(1020 / 30) atol = 1e-12 + # The property CSS exists for: the dominant feature moves it far less than TSS. + @test abs(log(css.factors[2] / css.factors[1])) < abs(log(tss.factors[2] / tss.factors[1])) / 5 + end + + @testset "RSS/TMM — a sample that is a multiple of another" begin + table_b = [100.0 200; + 200.0 400; + 300.0 600] # column 2 = 2 * column 1, so every log ratio is log2(2) + outcome = tmm_factors(table_b) + @test outcome.reference == "column 1" # chosen from the data + @test outcome.raw ≈ [1.0, 2.0] atol = 1e-12 + # centred: 1/sqrt(2), sqrt(2) + @test outcome.factors ≈ [0.7071067811865476, 1.4142135623730951] atol = 1e-10 + @test outcome.offset[2] - outcome.offset[1] ≈ log(2) atol = 1e-12 + @test outcome.parameters["reference_chosen_from_data"] == true + end + + @testset "RSS/TMM — known answer on the reference table" begin + outcome = tmm_factors(table_a; sample_ids = ids_a) + # reference selection: upper-quartile / library size is 25/60, 40/100, 45/100; + # the mean is 0.4222222222222222 and 25/60 is closest to it, so s1 is the reference + @test outcome.reference == "s1" + # raw factors before centring: sample 1 against itself is 1 by definition + @test outcome.raw ≈ [1.0, 1.5874010519681996, 1.296335601619136] atol = 1e-10 + @test outcome.factors ≈ [0.786198050576376, 1.248011612540287, 1.019176522885718] atol = 1e-10 + @test outcome.parameters["log_ratio_trim"] == 0.3 + @test outcome.parameters["sum_trim"] == 0.05 + @test length(outcome.parameters["kept_log_ratios"]) == 3 + end + + @testset "size_factors — median-of-ratios, which is what the name claimed" begin + outcome = rle_factors(table_a; sample_ids = ids_a) + # feature geometric means across samples: row 1 = 20.0, row 2 = 15.874010519681996, + # row 3 = 44.82219709646434 (all rows positive: every feature used) + @test outcome.parameters["features_used"] == 3 + # column 1: median(10/20, 20/15.874010519681996, 30/44.82219709646434) + # = median(0.5, 1.2599210498948732, 0.6694329500821695) = 0.6694329500821695 + @test outcome.raw ≈ [0.6694329500821695, 1.259921049894873, 1.1157215834702825] atol = 1e-10 + @test outcome.factors ≈ [0.683132584572384, 1.285704749170459, 1.13855430762064] atol = 1e-10 + # The old implementation returned library size / its geometric mean, which is the + # TSS factor. Assert the two are different on data where they differ, so a silent + # revert to the old behaviour fails here. + tss = tss_factors(table_a) + @test !isapprox(outcome.factors, tss.factors; atol = 1e-6) + end + + @testset "properties a wrong implementation fails" begin + base = css_factors(table_a; quantile = 0.75) + permuted = css_factors(table_a[[3, 1, 2], :]; quantile = 0.75) + @test permuted.factors ≈ base.factors atol = 1e-12 + + # Doubling one sample's counts doubles its own scaling factor relative to the others. + doubled = hcat(table_a[:, 1], 2 .* table_a[:, 2], table_a[:, 3]) + t0 = tss_factors(table_a) + t1 = tss_factors(doubled) + @test (t1.factors[2] / t1.factors[1]) ≈ 2 * (t0.factors[2] / t0.factors[1]) atol = 1e-12 + + # Scaling every sample by the same constant changes nothing: the factors are + # ratios, and the constant is absorbed by the centring. + scaled = 7.5 .* table_a + @test tss_factors(scaled).factors ≈ t0.factors atol = 1e-12 + @test tmm_factors(scaled).factors ≈ tmm_factors(table_a).factors atol = 1e-12 + @test rle_factors(scaled).factors ≈ rle_factors(table_a).factors atol = 1e-12 + end + + @testset "refusals name what they refused" begin + # A sample with no counts has no library size to divide by. + err = try + tss_factors([0.0 1.0; 0.0 2.0]) + catch e + e + end + @test err isa Scaling.ScalingRefusal + @test occursin("column 1", err.reason) + @test !occursin("log(0)", err.reason) # TSS says library size, not log(0) + + # CSS on a mostly-zero sample at a low quantile: the cumulative sum is zero, and + # log(0) is not a small number. This is the failure mode issue #16 predicted. + mostly_zero = [0.0 100; + 0.0 100; + 1.0 100] + err_css = try + css_factors(mostly_zero; quantile = 0.25) + catch e + e + end + @test err_css isa Scaling.ScalingRefusal + @test occursin("column 1", err_css.reason) + @test occursin("log(0)", err_css.reason) + @test occursin("0.25", err_css.reason) + + # A reference that does not exist is refused, and the refusal lists the real names + # rather than silently choosing one. + err_ref = try + tmm_factors(table_a; ref_column = "s99", sample_ids = ids_a) + catch e + e + end + @test err_ref isa Scaling.ScalingRefusal + @test occursin("s99", err_ref.reason) + @test occursin("s1", err_ref.reason) + + # Samples that share no positive feature cannot be compared. + err_disjoint = try + tmm_factors([1.0 0.0; 2.0 0.0; 0.0 3.0; 0.0 4.0]) + catch e + e + end + @test err_disjoint isa Scaling.ScalingRefusal + @test occursin("shares no positive-count feature", err_disjoint.reason) + + # No feature positive in every sample: the median-of-ratios is undefined, and the + # refusal says so instead of falling back to library size. + err_rle = try + rle_factors([1.0 0.0; 2.0 0.0]) + catch e + e + end + @test err_rle isa Scaling.ScalingRefusal + @test occursin("geometric mean", err_rle.reason) + + # Declared parameters outside their domain are refused before anything is computed. + @test_throws ArgumentError css_factors(table_a; quantile = 0.0) + @test_throws ArgumentError css_factors(table_a; quantile = 1.0) + @test_throws ArgumentError css_factors(table_a; quantile = 1.5) + @test_throws ArgumentError tmm_factors(table_a; log_ratio_trim = 0.5) + @test_throws ArgumentError tmm_factors(table_a; log_ratio_trim = -0.1) + @test_throws ArgumentError tmm_factors(table_a; sum_trim = 0.5) + end + + @testset "dispatch and provenance" begin + @test factors_for("TSS", table_a).kind == "tss" + @test factors_for("none", table_a).kind == "tss" + @test factors_for("css", table_a).kind == "css" + @test factors_for("tmm", table_a).kind == "rss" + @test factors_for("RSS", table_a).kind == "rss" + @test factors_for("size_factors", table_a).kind == "size_factors" + @test_throws ArgumentError factors_for("clr", table_a) + + checks = factor_checks(css_factors(table_a; quantile = 0.75)) + @test checks["kind"] == "css" + @test checks["parameters"]["quantile"] == 0.75 + @test haskey(checks, "raw") + + p1 = factor_provenance(css_factors(table_a; quantile = 0.75)) + p2 = factor_provenance(css_factors(table_a; quantile = 0.3)) + @test p1["offset_sha256"] != p2["offset_sha256"] + @test p1["offset_sha256"] == + factor_provenance(css_factors(table_a; quantile = 0.75))["offset_sha256"] + @test length(p1["offset_sha256"]) == 64 + end + + @testset "a quantile below the median is warned about, and still computed" begin + outcome = css_factors(table_a; quantile = 0.4) + @test outcome.kind == "css" + @test any(w -> occursin("below the median", w), outcome.warnings) + @test outcome.parameters["quantile"] == 0.4 + end + + # ------------------------------------------------------------------ + # Integration: the configuration layer, the execution path, and the manifest + # ------------------------------------------------------------------ + @testset "every admissible spelling of a method name is accepted" begin + # Regression: the allowed-normalisation table arrived upper case while + # NormalizationConfig stored lower case, so `method="TSS"` raised + # "incompatible with method" from AnalysisConfig — including for the exact TSS + # offset the CHANGELOG said had shipped, and the existing test execution testset + # that constructs method="TSS". Both sides are compared in lower case now. + adv = AnalysisConfig.AdvancedConfig(min_prevalence = 0.0, min_abundance = 0.0, + min_samples_per_group = 2) + for spelling in ["TSS", "tss", "CSS", "css", "RSS", "rss", "size_factors"] + norm = AnalysisConfig.NormalizationConfig(method = spelling) + @test norm.method == lowercase(spelling) + config = AnalysisConfig.AnalysisConfig( + method = "nb_glm", formula = "~ group", metadata_columns = ["group"], + normalization = norm, advanced = adv, created_by = "test_scaling") + @test config.normalization.method == lowercase(spelling) + end + + # CSS/RSS are offsets and stay refused for a response with no counts to offset. + @test_throws ArgumentError AnalysisConfig.AnalysisConfig( + method = "clr_lm", formula = "~ group", metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = "css"), + advanced = adv, created_by = "test_scaling") + @test_throws ArgumentError AnalysisConfig.AnalysisConfig( + method = "logistic", formula = "group ~ group", outcome_column = "group", + metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = "rss"), + advanced = adv, created_by = "test_scaling") + end + + @testset "declared scaling parameters are validated and hashed" begin + @test_throws ArgumentError AnalysisConfig.NormalizationConfig(method = "css", css_quantile = 0.0) + @test_throws ArgumentError AnalysisConfig.NormalizationConfig(method = "css", css_quantile = 1.2) + @test_throws ArgumentError AnalysisConfig.NormalizationConfig(method = "rss", tmm_log_ratio_trim = 0.5) + @test_throws ArgumentError AnalysisConfig.NormalizationConfig(method = "rss", tmm_sum_trim = -0.1) + @test_throws ArgumentError AnalysisConfig.NormalizationConfig(method = "rss", tmm_ref_column = " ") + + adv = AnalysisConfig.AdvancedConfig(min_prevalence = 0.0, min_abundance = 0.0, + min_samples_per_group = 2) + make_config = q -> AnalysisConfig.AnalysisConfig( + method = "nb_glm", formula = "~ group", metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = "css", css_quantile = q), + advanced = adv, created_by = "test_scaling") + # The parameter is part of the config hash: a run that declares a different + # quantile is a different run, and the manifest has to say so. + @test make_config(0.75).hash != make_config(0.5).hash + end + + @testset "prepare_analysis_table keeps counts and hands over the offset" begin + adv = AnalysisConfig.AdvancedConfig(min_prevalence = 0.0, min_abundance = 0.0, + min_samples_per_group = 2) + expectations = [ + ("TSS", "tss", tss_factors(table_a).factors), + ("CSS", "css", css_factors(table_a; quantile = 0.75).factors), + ("RSS", "rss", tmm_factors(table_a).factors), + ("size_factors", "size_factors", rle_factors(table_a).factors), + ] + for (method_name, kind, expected) in expectations + config = AnalysisConfig.AnalysisConfig( + method = "nb_glm", formula = "~ group", metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = method_name), + advanced = adv, created_by = "test_scaling") + (prepared, diagnostics, manifest, sample_ids, _) = + Execution.prepare_analysis_table(config, table_a; + sample_ids = ids_a, + taxa_ids = ["t1", "t2", "t3"], + drop_policy = "drop") + @test prepared == table_a # the count response is untouched + @test manifest.offset ≈ log.(expected) atol = 1e-10 + @test diagnostics.checks["scaling"]["kind"] == kind + @test manifest.provenance["scaling"]["kind"] == kind + @test manifest.provenance["scaling"]["offset_sha256"] isa String + @test length(manifest.provenance["scaling"]["offset_sha256"]) == 64 + @test sample_ids == ids_a + end + + # `none` is not "no offset": for a count model it is the plain log library size, + # and it is recorded as such rather than left to be inferred. + config_none = AnalysisConfig.AnalysisConfig( + method = "nb_glm", formula = "~ group", metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = "none"), + advanced = adv, created_by = "test_scaling") + (_, diag_none, manifest_none, _, _) = Execution.prepare_analysis_table( + config_none, table_a; sample_ids = ids_a, taxa_ids = ["t1", "t2", "t3"], + drop_policy = "drop") + @test diag_none.checks["scaling"]["kind"] == "tss" + @test manifest_none.offset ≈ log.(tss_factors(table_a).factors) atol = 1e-10 + + # A run that declares `relative` has no offset, and does not get one by accident. + config_rel = AnalysisConfig.AnalysisConfig( + method = "nb_glm", formula = "~ group", metadata_columns = ["group"], + normalization = AnalysisConfig.NormalizationConfig(method = "relative"), + advanced = adv, created_by = "test_scaling") + (_, diag_rel, manifest_rel, _, _) = Execution.prepare_analysis_table( + config_rel, table_a; sample_ids = ids_a, taxa_ids = ["t1", "t2", "t3"], + drop_policy = "drop") + @test isnothing(manifest_rel.offset) + @test !haskey(diag_rel.checks, "scaling") + end + + # ------------------------------------------------------------------ + # The R transcription + # + # Base R only, written from the same conditions document. It is the same author and a + # different language: it catches transcription slips (an off-by-one in the trim, a + # quantile type, the wrong centring), not errors in the specification itself. The + # real reference packages are not in renv.lock, and this does not pretend to be them. + # ------------------------------------------------------------------ + r_reference_script = """ + r_A <- matrix(c(10, 20, 30, 20, 20, 60, 40, 10, 50), nrow = 3, byrow = FALSE) + r_libs <- colSums(r_A) + r_tss <- r_libs / exp(mean(log(r_libs))) + r_q <- 0.75 + r_css_raw <- sapply(seq_len(ncol(r_A)), function(j) { + r_col <- r_A[, j] + r_thr <- quantile(r_col, r_q, type = 7) + sum(r_col[r_col <= r_thr]) + }) + r_css <- r_css_raw / exp(mean(log(r_css_raw))) + r_uq <- sapply(seq_len(ncol(r_A)), function(j) quantile(r_A[, j], 0.75, type = 7) / r_libs[j]) + r_ref <- which.min(abs(r_uq - mean(r_uq))) + r_raw <- rep(1, ncol(r_A)) + for (i in seq_len(ncol(r_A))) { + if (i == r_ref) next + r_x <- r_A[, i] + r_r <- r_A[, r_ref] + r_keep <- r_x > 0 & r_r > 0 + r_m <- log2(r_x[r_keep] / r_r[r_keep]) + r_a <- 0.5 * log2(r_x[r_keep] * r_r[r_keep]) + r_w <- (1 - r_x[r_keep] / r_libs[i]) / r_x[r_keep] + (1 - r_r[r_keep] / r_libs[r_ref]) / r_r[r_keep] + r_n <- length(r_m) + r_lo <- floor(r_n * 0.3) + 1 + r_hi <- r_n + 1 - r_lo + r_los <- floor(r_n * 0.05) + 1 + r_his <- r_n + 1 - r_los + r_rm <- rank(r_m, ties.method = "first") + r_ra <- rank(r_a, ties.method = "first") + r_sel <- r_rm >= r_lo & r_rm <= r_hi & r_ra >= r_los & r_ra <= r_his + r_raw[i] <- 2^(sum(r_m[r_sel] * r_w[r_sel]) / sum(r_w[r_sel])) + } + r_tmm <- r_raw / exp(mean(log(r_raw))) + r_gm <- function(v) exp(mean(log(v))) + r_rle_raw <- sapply(seq_len(ncol(r_A)), function(j) median(r_A[, j] / apply(r_A, 1, r_gm))) + r_rle <- r_rle_raw / exp(mean(log(r_rle_raw))) + """ + + r_reference = try + with_r_lock(; timeout = 60.0) do + RCall.globalEnv[:scaling_A] = table_a + RCall.reval(r_reference_script) + ( + tss = convert(Vector{Float64}, RCall.reval("r_tss")), + css = convert(Vector{Float64}, RCall.reval("r_css")), + tmm = convert(Vector{Float64}, RCall.reval("r_tmm")), + rle = convert(Vector{Float64}, RCall.reval("r_rle")), + reference = RCall.rcopy(RCall.reval("r_ref")), + ) + end + catch + nothing + end + + if r_reference === nothing + # Loud, by name, with the reason: a silent skip is how a validation claim goes + # missing. CI installs R and restores renv.lock, so a skip here is a regression. + @info "SKIPPED (by name): the R transcription cross-check for TSS/CSS/RSS/RLE — R is not reachable in this environment." + @test_skip css_factors(table_a; quantile = 0.75) + else + @testset "an R transcription of the same conditions agrees" begin + @test r_reference.tss ≈ tss_factors(table_a).factors atol = 1e-12 + @test r_reference.css ≈ css_factors(table_a; quantile = 0.75).factors atol = 1e-12 + @test r_reference.tmm ≈ tmm_factors(table_a).factors atol = 1e-12 + @test r_reference.rle ≈ rle_factors(table_a).factors atol = 1e-12 + # And R reproduces the published numbers, so the constants above are not + # Julia-version artefacts. + @test r_reference.tss ≈ [0.711378660898012, 1.185631101496687, 1.185631101496687] atol = 1e-12 + @test r_reference.tmm ≈ [0.786198050576376, 1.248011612540287, 1.019176522885718] atol = 1e-12 + @test r_reference.reference == 1 + end + end +end