diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 86fd896e..573768be 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -156,7 +156,7 @@ jobs: sudo chmod 777 "$RENV_PATHS_CACHE" sudo env R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'install.packages(c("renv", "BiocManager"), repos="https://cloud.r-project.org", lib=.Library)' sudo env RENV_PATHS_CACHE="$RENV_PATHS_CACHE" R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'renv::restore(project=".", library=.Library, prompt=FALSE)' - Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","MASS","zCompositions")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' + Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","MASS","zCompositions","glmGamPoi")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' - name: Check DADA2 CSV assembly run: Rscript --vanilla test/unit/test_dada2_tables.R @@ -250,7 +250,7 @@ jobs: julia --project=. -e ' using RCall R"cat(.libPaths(), sep=\"\n\")" - for pkg in ["dada2", "vegan", "MASS", "zCompositions"] + for pkg in ["dada2", "vegan", "MASS", "zCompositions", "glmGamPoi"] R"if (!require($pkg, character.only=TRUE, quietly=TRUE)) stop($pkg, \" not found\")" end ' diff --git a/R/_renv_dependencies.R b/R/_renv_dependencies.R index 1e91aca9..589612fe 100644 --- a/R/_renv_dependencies.R +++ b/R/_renv_dependencies.R @@ -3,6 +3,7 @@ if (FALSE) { library(dada2) library(vegan) + library(glmGamPoi) library(MASS) library(zCompositions) } diff --git a/README.md b/README.md index a79119f3..5c31ca30 100644 --- a/README.md +++ b/README.md @@ -250,7 +250,7 @@ is the reference. Counts are aggregated by taxon after applicable read-count bounds to the differential-abundance request, so those bounds do not exclude samples from this analysis. -All `analysis.differential` options below (`method`, `offset`, +All `analysis.differential` options below (`method`, `offset`, `dispersion`, `min_prevalence` and `replacement_delta`) are resolved at study scope. Set them in the study configuration; group and run overrides are not used by this endpoint, including when called from a group Analysis workspace. @@ -260,7 +260,15 @@ by this endpoint, including when called from a group Analysis workspace. - `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`, with a library-size offset. The effect is log2 fold change; positive values indicate greater abundance in the second condition. `offset` selects `tss` - (default) or `rle` scaling. + (default) or `rle` scaling. `dispersion` selects the fits and the test: + `glm_nb` (default) fits `MASS::glm.nb` per taxon with its own + maximum-likelihood theta and Wald-tests the group coefficient; `glmGamPoi` + fits one `glmGamPoi::glm_gp` across the taxa passing `min_prevalence`, + sharing a dispersion trend across them, and reports `glmGamPoi::test_de`'s + quasi-likelihood F-test. Under `glmGamPoi` the p-value comes from that + F-test, the standard error is `test_de`'s `lfc_se`, and theta is each + taxon's own estimate before shrinkage. It requires the pinned glmGamPoi R + package and is refused for `clr_welch`. - `clr_welch`: Welch's t-test on centred log-ratios after multiplicative zero replacement with `zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change. `replacement_delta` controls the replacement fraction diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 8181ae6d..debcff02 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -358,6 +358,16 @@ analysis: # rle is refused when no taxon is present in every sample. # Used by nb_glm only; clr_welch ignores it, since log-ratios carry no library size. offset: tss + # nb_glm only: the fits and the test. glm_nb (the default) fits + # MASS::glm.nb per taxon with its maximum-likelihood theta and Wald-tests + # the group coefficient. glmGamPoi fits one glmGamPoi::glm_gp across the + # taxa passing min_prevalence, with overdispersion shrinkage toward the + # cross-taxon trend, and reports glmGamPoi::test_de's quasi-likelihood + # F-test of the group contrast with its log2 fold changes; the reported + # theta inverts glm_gp's negative-binomial overdispersion, and a taxon + # whose overdispersion is zero is marked at the Poisson boundary, as under + # glm_nb. Requires the pinned glmGamPoi R package. Refused with clr_welch. + dispersion: glm_nb # Fraction of samples, 0 to 1, in which a taxon must have reads to be # tested. Taxa below it are listed as filtered, without a p-value. For # clr_welch it is applied before zero replacement, so on a sparse table it also diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index f516650e..8a452db0 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -440,6 +440,8 @@ export interface DifferentialResult { min_prevalence: number /** nb_glm only. */ offset?: 'tss' | 'rle' + /** nb_glm only: which fits and test produced the rows (per-taxon MASS::glm.nb with a Wald test, or one glmGamPoi::glm_gp with test_de's quasi-likelihood F-test). */ + dispersion?: 'glm_nb' | 'glmGamPoi' /** clr_welch only. */ replacement_delta?: number } diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts new file mode 100644 index 00000000..dd950da3 --- /dev/null +++ b/frontend/src/components/PipelineStages.test.ts @@ -0,0 +1,10 @@ +// SPDX-License-Identifier: AGPL-3.0-only +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +import { expect, test } from 'bun:test' +import { CONFIG_TYPES } from './PipelineStages' + +test('dispersion uses the enum editor with both supported models', () => { + expect(CONFIG_TYPES['analysis.differential.dispersion']).toEqual({ + kind: 'enum', options: ['glm_nb', 'glmGamPoi'], + }) +}) diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index 85129267..b938f77a 100644 --- a/frontend/src/components/PipelineStages.tsx +++ b/frontend/src/components/PipelineStages.tsx @@ -152,6 +152,7 @@ export const CONFIG_DESCRIPTIONS: Record = { 'analysis.beta.depth': 'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.', 'analysis.differential.method': 'Model fitted to each taxon: nb_glm (negative-binomial GLM on read counts; effect is a log2 fold change) or clr_welch (Welch\'s t-test on centred log-ratios over the taxa passing min_prevalence, after zCompositions::cmultRepl zero replacement; effect is a CLR difference, not a fold change).', 'analysis.differential.offset': 'Library-size offset in each negative-binomial model: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none). Used by nb_glm only.', + 'analysis.differential.dispersion': 'nb_glm only: glm_nb fits MASS::glm.nb per taxon with its own maximum-likelihood theta; glmGamPoi fits one glmGamPoi::glm_gp across taxa with overdispersion shrinkage and reports glmGamPoi::test_de’s quasi-likelihood F-test with its log2 fold changes.', 'analysis.differential.min_prevalence': 'Fraction of samples, 0 to 1, in which a taxon must have reads to be tested; taxa below it are listed as filtered.', 'analysis.differential.replacement_delta': 'clr_welch only: cmultRepl\'s frac, between 0 and 1; an imputed proportion above its taxon\'s smallest observed proportion is lowered to this fraction of it. Replaced zeros are not measurements.', 'analysis.alpha.show_points': 'Show individual samples overlaid on alpha comparison boxplots.', @@ -252,6 +253,7 @@ export const CONFIG_TYPES: Record = { 'analysis.beta.depth': { kind: 'int' }, 'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch'] }, 'analysis.differential.offset': { kind: 'enum', options: ['tss', 'rle'] }, + 'analysis.differential.dispersion': { kind: 'enum', options: ['glm_nb', 'glmGamPoi'] }, 'analysis.differential.min_prevalence': { kind: 'float', min: 0, max: 1, step: 0.05 }, 'analysis.differential.replacement_delta': { kind: 'float', min: 0.01, max: 0.99, step: 0.05 }, 'analysis.alpha.show_points': { kind: 'boolean' }, diff --git a/install.jl b/install.jl index 964e00d9..3be35d24 100644 --- a/install.jl +++ b/install.jl @@ -879,8 +879,8 @@ function setup_r_packages(; rebuild::Bool=false)::Tuple{Symbol,String} # Verdict is whether the packages the pipeline loads are actually usable, not # renv's exit code: restore can fail on one leaf package and still leave a # working library, or "succeed" against a stale cache. - # MASS and zCompositions back the differential abundance methods. - pkgs = ["dada2", "vegan", "MASS", "zCompositions"] + # MASS, zCompositions and glmGamPoi back the differential abundance methods. + pkgs = ["dada2", "vegan", "MASS", "zCompositions", "glmGamPoi"] probe = """ pkgs <- c($(join(repr.(pkgs), ", "))) cat(paste(pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)], diff --git a/renv.lock b/renv.lock index cfd72d4b..9944f58c 100644 --- a/renv.lock +++ b/renv.lock @@ -3011,6 +3011,565 @@ "BugReports": "https://github.com/Japal/zCompositions/issues", "NeedsCompilation": "no", "Author": "Javier Palarea-Albaladejo [cre, aut] (ORCID: ), Josep Antoni Martin-Fernandez [aut] (ORCID: )" + }, + "DelayedMatrixStats": { + "Package": "DelayedMatrixStats", + "Version": "1.32.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Functions that Apply to Rows and Columns of 'DelayedMatrix' Objects", + "Date": "2025-01-09", + "Authors@R": "c(person(\"Peter\", \"Hickey\", role = c(\"aut\", \"cre\"), email = \"peter.hickey@gmail.com\", comment = c(ORCID = \"0000-0002-8153-6258\")), person(\"Hervé\", \"Pagès\", role = \"ctb\", email = \"hpages.on.github@gmail.com\"), person(\"Aaron\", \"Lun\", role = \"ctb\", email = \"infinite.monkeys.with.keyboards@gmail.com\"))", + "Description": "A port of the 'matrixStats' API for use with DelayedMatrix objects from the 'DelayedArray' package. High-performing functions operating on rows and columns of DelayedMatrix objects, e.g. col / rowMedians(), col / rowRanks(), and col / rowSds(). Functions optimized per data type and for subsetted calculations such that both memory usage and processing time is minimized.", + "License": "MIT + file LICENSE", + "Encoding": "UTF-8", + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.2", + "Depends": [ + "MatrixGenerics (>= 1.15.1)", + "DelayedArray (>= 0.31.7)" + ], + "Imports": [ + "methods", + "sparseMatrixStats (>= 1.13.2)", + "Matrix (>= 1.5-0)", + "S4Vectors (>= 0.17.5)", + "IRanges (>= 2.25.10)", + "SparseArray (>= 1.5.19)" + ], + "Suggests": [ + "testthat", + "knitr", + "rmarkdown", + "BiocStyle", + "microbenchmark", + "profmem", + "HDF5Array", + "matrixStats (>= 1.0.0)" + ], + "VignetteBuilder": "knitr", + "URL": "https://github.com/PeteHaitch/DelayedMatrixStats", + "BugReports": "https://github.com/PeteHaitch/DelayedMatrixStats/issues", + "biocViews": "Infrastructure, DataRepresentation, Software", + "git_url": "https://git.bioconductor.org/packages/DelayedMatrixStats", + "git_branch": "RELEASE_3_22", + "git_last_commit": "cbf7d75", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Peter Hickey [aut, cre] (ORCID: ), Hervé Pagès [ctb], Aaron Lun [ctb]", + "Maintainer": "Peter Hickey " + }, + "HDF5Array": { + "Package": "HDF5Array", + "Version": "1.38.0", + "Source": "Bioconductor", + "Title": "HDF5 datasets as array-like objects in R", + "Description": "The HDF5Array package is an HDF5 backend for DelayedArray objects. It implements the HDF5Array, H5SparseMatrix, H5ADMatrix, and TENxMatrix classes, 4 convenient and memory-efficient array-like containers for representing and manipulating either: (1) a conventional (a.k.a. dense) HDF5 dataset, (2) an HDF5 sparse matrix (stored in CSR/CSC/Yale format), (3) the central matrix of an h5ad file (or any matrix in the /layers group), or (4) a 10x Genomics sparse matrix. All these containers are DelayedArray extensions and thus support all operations (delayed or block-processed) supported by DelayedArray objects.", + "biocViews": "Infrastructure, DataRepresentation, DataImport, Sequencing, RNASeq, Coverage, Annotation, GenomeAnnotation, SingleCell, ImmunoOncology", + "URL": "https://bioconductor.org/packages/HDF5Array", + "BugReports": "https://github.com/Bioconductor/HDF5Array/issues", + "License": "Artistic-2.0", + "Encoding": "UTF-8", + "Authors@R": "person(\"Hervé\", \"Pagès\", role=c(\"aut\", \"cre\"), email=\"hpages.on.github@gmail.com\", comment=c(ORCID=\"0009-0002-8272-4522\"))", + "Depends": [ + "R (>= 3.4)", + "methods", + "SparseArray (>= 1.7.5)", + "DelayedArray (>= 0.33.5)", + "h5mread (>= 0.99.4)" + ], + "Imports": [ + "utils", + "stats", + "tools", + "Matrix", + "BiocGenerics (>= 0.51.2)", + "S4Vectors", + "IRanges", + "S4Arrays (>= 1.1.1)", + "rhdf5" + ], + "Suggests": [ + "BiocParallel", + "GenomicRanges", + "SummarizedExperiment (>= 1.15.1)", + "h5vcData", + "ExperimentHub", + "TENxBrainData", + "zellkonverter", + "GenomicFeatures", + "SingleCellExperiment", + "DelayedMatrixStats", + "genefilter", + "RSpectra", + "RUnit", + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "Collate": "utils.R h5utils.R HDF5ArraySeed-class.R HDF5Array-class.R ReshapedHDF5ArraySeed-class.R ReshapedHDF5Array-class.R dump-management.R writeHDF5Array.R saveHDF5SummarizedExperiment.R H5SparseMatrixSeed-class.R H5SparseMatrix-class.R H5ADMatrixSeed-class.R H5ADMatrix-class.R TENxMatrixSeed-class.R TENxMatrix-class.R writeTENxMatrix.R zzz.R", + "git_url": "https://git.bioconductor.org/packages/HDF5Array", + "git_branch": "RELEASE_3_22", + "git_last_commit": "9bca08f", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Hervé Pagès [aut, cre] (ORCID: )", + "Maintainer": "Hervé Pagès " + }, + "RcppArmadillo": { + "Package": "RcppArmadillo", + "Version": "15.6.2-1", + "Source": "Repository", + "Type": "Package", + "Title": "'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library", + "Date": "2026-10-16", + "Authors@R": "c(person(\"Dirk\", \"Eddelbuettel\", role = c(\"aut\", \"cre\"), email = \"edd@debian.org\", comment = c(ORCID = \"0000-0001-6419-907X\")), person(\"Romain\", \"Francois\", role = \"aut\", comment = c(ORCID = \"0000-0002-2444-4226\")), person(\"Doug\", \"Bates\", role = \"aut\", comment = c(ORCID = \"0000-0001-8316-9503\")), person(\"Binxiang\", \"Ni\", role = \"aut\"), person(\"Conrad\", \"Sanderson\", role = \"aut\", comment = c(ORCID = \"0000-0002-0049-4501\")))", + "Description": "'Armadillo' is a templated C++ linear algebra library aiming towards a good balance between speed and ease of use. It provides high-level syntax and functionality deliberately similar to Matlab. It is useful for algorithm development directly in C++, or quick conversion of research code into production environments. It provides efficient classes for vectors, matrices and cubes where dense and sparse matrices are supported. Integer, floating point and complex numbers are supported. A sophisticated expression evaluator (based on template meta-programming) automatically combines several operations to increase speed and efficiency. Dynamic evaluation automatically chooses optimal code paths based on detected matrix structures. Matrix decompositions are provided through integration with LAPACK, or one of its high performance drop-in replacements (such as 'MKL' or 'OpenBLAS'). It can automatically use 'OpenMP' multi-threading (parallelisation) to speed up computationally expensive operations. . The 'RcppArmadillo' package includes the header files from the 'Armadillo' library; users do not need to install 'Armadillo' itself in order to use 'RcppArmadillo'. Starting from release 15.0.0, the minimum compilation standard is C++14. . Since release 7.800.0, 'Armadillo' is licensed under Apache License 2; previous releases were under licensed as MPL 2.0 from version 3.800.0 onwards and LGPL-3 prior to that; 'RcppArmadillo' (the 'Rcpp' bindings/bridge to Armadillo) is licensed under the GNU GPL version 2 or later, as is the rest of 'Rcpp'.", + "License": "GPL (>= 2)", + "LazyLoad": "yes", + "Depends": [ + "R (>= 3.3.0)" + ], + "LinkingTo": [ + "Rcpp" + ], + "Imports": [ + "Rcpp (>= 1.1.1)", + "stats", + "utils", + "methods" + ], + "Suggests": [ + "tinytest", + "Matrix (>= 1.3.0)", + "pkgKitten", + "reticulate", + "slam" + ], + "URL": "https://github.com/RcppCore/RcppArmadillo, https://dirk.eddelbuettel.com/code/rcpp.armadillo.html", + "BugReports": "https://github.com/RcppCore/RcppArmadillo/issues", + "RoxygenNote": "6.0.1", + "Encoding": "UTF-8", + "VignetteBuilder": "Rcpp", + "NeedsCompilation": "yes", + "Author": "Dirk Eddelbuettel [aut, cre] (ORCID: ), Romain Francois [aut] (ORCID: ), Doug Bates [aut] (ORCID: ), Binxiang Ni [aut], Conrad Sanderson [aut] (ORCID: )", + "Maintainer": "Dirk Eddelbuettel ", + "Repository": "CRAN" + }, + "Rhdf5lib": { + "Package": "Rhdf5lib", + "Version": "1.32.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "hdf5 library as an R package", + "Authors@R": "c( person( \"Mike\", \"Smith\", role = \"ctb\", comment = c(ORCID = \"0000-0002-7800-3848\") ), person( \"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\")) , person( given = \"The HDF Group\", role = \"cph\" ) )", + "Description": "Provides C and C++ hdf5 libraries.", + "License": "Artistic-2.0", + "Copyright": "src/hdf5/COPYING", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Depends": [ + "R (>= 4.2.0)" + ], + "Suggests": [ + "BiocStyle", + "knitr", + "rmarkdown", + "tinytest", + "mockery" + ], + "URL": "https://github.com/Huber-group-EMBL/Rhdf5lib", + "BugReports": "https://github.com/Huber-group-EMBL/Rhdf5lib/issues", + "SystemRequirements": "GNU make", + "Encoding": "UTF-8", + "biocViews": "Infrastructure", + "RoxygenNote": "7.1.2", + "git_url": "https://git.bioconductor.org/packages/Rhdf5lib", + "git_branch": "RELEASE_3_22", + "git_last_commit": "f62ae28", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Mike Smith [ctb] (ORCID: ), Hugo Gruson [cre] (ORCID: ), The HDF Group [cph]", + "Maintainer": "Hugo Gruson " + }, + "SingleCellExperiment": { + "Package": "SingleCellExperiment", + "Version": "1.32.0", + "Source": "Bioconductor", + "Date": "2025-06-27", + "Title": "S4 Classes for Single Cell Data", + "Authors@R": "c( person(\"Aaron\", \"Lun\", role=c(\"aut\", \"cph\"), email=\"infinite.monkeys.with.keyboards@gmail.com\"), person(\"Davide\",\"Risso\", role=c(\"aut\",\"cre\", \"cph\"), email=\"risso.davide@gmail.com\"), person(\"Keegan\", \"Korthauer\", role=\"ctb\"), person(\"Kevin\", \"Rue-Albrecht\", role=\"ctb\"), person(\"Luke\", \"Zappia\", role = \"ctb\", comment = c(ORCID = \"0000-0001-7744-8565\", github = \"lazappi\")))", + "Depends": [ + "SummarizedExperiment" + ], + "Imports": [ + "methods", + "utils", + "stats", + "S4Vectors", + "BiocGenerics", + "GenomicRanges", + "DelayedArray" + ], + "Suggests": [ + "testthat", + "BiocStyle", + "knitr", + "rmarkdown", + "Matrix", + "scRNAseq (>= 2.9.1)", + "Rtsne" + ], + "biocViews": "ImmunoOncology, DataRepresentation, DataImport, Infrastructure, SingleCell", + "Description": "Defines a S4 class for storing data from single-cell experiments. This includes specialized methods to store and retrieve spike-in information, dimensionality reduction coordinates and size factors for each cell, along with the usual metadata for genes and libraries.", + "License": "GPL-3", + "VignetteBuilder": "knitr", + "RoxygenNote": "7.3.2", + "git_url": "https://git.bioconductor.org/packages/SingleCellExperiment", + "git_branch": "RELEASE_3_22", + "git_last_commit": "db7133e", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Aaron Lun [aut, cph], Davide Risso [aut, cre, cph], Keegan Korthauer [ctb], Kevin Rue-Albrecht [ctb], Luke Zappia [ctb] (ORCID: , github: lazappi)", + "Maintainer": "Davide Risso " + }, + "assorthead": { + "Package": "assorthead", + "Version": "1.4.0", + "Source": "Bioconductor", + "Date": "2025-10-27", + "Title": "Assorted Header-Only C++ Libraries", + "Authors@R": "person(\"Aaron\", \"Lun\", role=c(\"cre\", \"aut\"), email=\"infinite.monkeys.with.keyboards@gmail.com\")", + "Description": "Vendors an assortment of useful header-only C++ libraries. Bioconductor packages can use these libraries in their own C++ code by LinkingTo this package without introducing any additional dependencies. The use of a central repository avoids duplicate vendoring of libraries across multiple R packages, and enables better coordination of version updates across cohorts of interdependent C++ libraries.", + "License": "MIT + file LICENSE", + "Suggests": [ + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "URL": "https://github.com/LTLA/assorthead", + "BugReports": "https://github.com/LTLA/assorthead/issues", + "Encoding": "UTF-8", + "biocViews": "SingleCell, QualityControl, Normalization, DataRepresentation, DataImport, DifferentialExpression, Alignment", + "git_url": "https://git.bioconductor.org/packages/assorthead", + "git_branch": "RELEASE_3_22", + "git_last_commit": "9255f2f", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Aaron Lun [cre, aut]", + "Maintainer": "Aaron Lun " + }, + "beachmat": { + "Package": "beachmat", + "Version": "2.26.0", + "Source": "Bioconductor", + "Date": "2025-08-19", + "Title": "Compiling Bioconductor to Handle Each Matrix Type", + "Authors@R": "c(person(\"Aaron\", \"Lun\", role = c(\"aut\", \"cre\"), email = \"infinite.monkeys.with.keyboards@gmail.com\"), person(\"Hervé\", \"Pagès\", role=\"aut\"), person(\"Mike\", \"Smith\", role=\"aut\"))", + "Imports": [ + "methods", + "DelayedArray (>= 0.27.2)", + "SparseArray", + "BiocGenerics", + "Matrix", + "Rcpp" + ], + "Suggests": [ + "testthat", + "BiocStyle", + "knitr", + "rmarkdown", + "rcmdcheck", + "BiocParallel", + "HDF5Array", + "beachmat.hdf5" + ], + "LinkingTo": [ + "Rcpp", + "assorthead (>= 1.3.3)" + ], + "biocViews": "DataRepresentation, DataImport, Infrastructure", + "Description": "Provides a consistent C++ class interface for reading from a variety of commonly used matrix types. Ordinary matrices and several sparse/dense Matrix classes are directly supported, along with a subset of the delayed operations implemented in the DelayedArray package. All other matrix-like objects are supported by calling back into R.", + "License": "GPL-3", + "NeedsCompilation": "yes", + "VignetteBuilder": "knitr", + "SystemRequirements": "C++17", + "URL": "https://github.com/tatami-inc/beachmat", + "BugReports": "https://github.com/tatami-inc/beachmat/issues", + "RoxygenNote": "7.3.2", + "Encoding": "UTF-8", + "git_url": "https://git.bioconductor.org/packages/beachmat", + "git_branch": "RELEASE_3_22", + "git_last_commit": "016c55e", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "Author": "Aaron Lun [aut, cre], Hervé Pagès [aut], Mike Smith [aut]", + "Maintainer": "Aaron Lun " + }, + "glmGamPoi": { + "Package": "glmGamPoi", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Fit a Gamma-Poisson Generalized Linear Model", + "Authors@R": "c(person(\"Constantin\", \"Ahlmann-Eltze\", email = \"artjom31415@googlemail.com\", role = c(\"aut\", \"cre\"), comment = c(ORCID = \"0000-0002-3762-068X\")), person(\"Nathan\", \"Lubock\", email=\"lubock.nathan@gmail.com\", role = \"ctb\", comment = c(ORCID = \"0000-0001-8064-2465\")), person(\"Michael\", \"Love\", email=\"michaelisaiahlove@gmail.com\", role = \"ctb\"))", + "Description": "Fit linear models to overdispersed count data. The package can estimate the overdispersion and fit repeated models for matrix input. It is designed to handle large input datasets as they typically occur in single cell RNA-seq experiments.", + "License": "GPL-3", + "Encoding": "UTF-8", + "SystemRequirements": "C++17", + "Suggests": [ + "testthat (>= 2.1.0)", + "zoo", + "DESeq2", + "edgeR", + "limma", + "MASS", + "statmod", + "ggplot2", + "bench", + "BiocParallel", + "knitr", + "rmarkdown", + "BiocStyle", + "TENxPBMCData", + "muscData", + "scran", + "dplyr" + ], + "LinkingTo": [ + "Rcpp", + "RcppArmadillo", + "beachmat", + "assorthead" + ], + "Imports": [ + "Rcpp", + "beachmat", + "DelayedMatrixStats", + "matrixStats", + "MatrixGenerics", + "SparseArray (>= 1.5.21)", + "S4Vectors", + "DelayedArray", + "HDF5Array", + "Matrix", + "SummarizedExperiment", + "SingleCellExperiment", + "BiocGenerics", + "methods", + "stats", + "utils", + "splines", + "rlang", + "vctrs" + ], + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.2", + "URL": "https://github.com/const-ae/glmGamPoi", + "BugReports": "https://github.com/const-ae/glmGamPoi/issues", + "biocViews": "Regression, RNASeq, Software, SingleCell", + "VignetteBuilder": "knitr", + "git_url": "https://git.bioconductor.org/packages/glmGamPoi", + "git_branch": "RELEASE_3_22", + "git_last_commit": "908cde9", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Constantin Ahlmann-Eltze [aut, cre] (ORCID: ), Nathan Lubock [ctb] (ORCID: ), Michael Love [ctb]", + "Maintainer": "Constantin Ahlmann-Eltze ", + "Depends": [ + "R (>= 4.1.0)" + ] + }, + "h5mread": { + "Package": "h5mread", + "Version": "1.2.1", + "Source": "Bioconductor", + "Title": "A fast HDF5 reader", + "Description": "The main function in the h5mread package is h5mread(), which allows reading arbitrary data from an HDF5 dataset into R, similarly to what the h5read() function from the rhdf5 package does. In the case of h5mread(), the implementation has been optimized to make it as fast and memory-efficient as possible.", + "biocViews": "Infrastructure, DataRepresentation, DataImport", + "URL": "https://bioconductor.org/packages/h5mread", + "BugReports": "https://github.com/Bioconductor/h5mread/issues", + "License": "Artistic-2.0", + "Encoding": "UTF-8", + "Authors@R": "person(\"Hervé\", \"Pagès\", role=c(\"aut\", \"cre\"), email=\"hpages.on.github@gmail.com\", comment=c(ORCID=\"0009-0002-8272-4522\"))", + "Depends": [ + "R (>= 4.5)", + "methods", + "rhdf5", + "BiocGenerics", + "SparseArray" + ], + "Imports": [ + "stats", + "tools", + "rhdf5filters", + "S4Vectors", + "IRanges", + "S4Arrays" + ], + "LinkingTo": [ + "Rhdf5lib", + "S4Vectors" + ], + "SystemRequirements": "GNU make", + "Suggests": [ + "BiocParallel", + "ExperimentHub", + "TENxBrainData", + "HDF5Array", + "testthat", + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "Collate": "utils.R h5dim.R H5File-class.R h5ls.R H5DSetDescriptor-class.R uaselection.R h5mread.R h5summarize.R h5mread_from_reshaped.R h5dimscales.R h5writeDimnames.R zzz.R", + "git_url": "https://git.bioconductor.org/packages/h5mread", + "git_branch": "RELEASE_3_22", + "git_last_commit": "b377076", + "git_last_commit_date": "2025-11-21", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Hervé Pagès [aut, cre] (ORCID: )", + "Maintainer": "Hervé Pagès " + }, + "rhdf5": { + "Package": "rhdf5", + "Version": "2.54.1", + "Source": "Bioconductor", + "Type": "Package", + "Title": "R Interface to HDF5", + "Authors@R": "c( person(\"Bernd\", \"Fischer\", role = \"aut\"), person( \"Mike\", \"Smith\", role= \"aut\", comment = c(ORCID = \"0000-0002-7800-3848\", \"Maintainer from 2017 to 2025\") ), person(\"Gregoire\", \"Pau\", role=\"aut\"), person(\"Martin\", \"Morgan\", role = \"ctb\"), person(\"Daniel\", \"van Twisk\", role = \"ctb\"), person( \"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\") ) )", + "Description": "This package provides an interface between HDF5 and R. HDF5's main features are the ability to store and access very large and/or complex datasets and a wide variety of metadata on mass storage (disk) through a completely portable file format. The rhdf5 package is thus suited for the exchange of large and/or complex datasets between R and other software package, and for letting R applications work on datasets that are larger than the available RAM.", + "License": "Artistic-2.0", + "URL": "https://github.com/Huber-group-EMBL/rhdf5", + "BugReports": "https://github.com/Huber-group-EMBL/rhdf5/issues", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Imports": [ + "Rhdf5lib (>= 1.13.4)", + "rhdf5filters (>= 1.15.5)" + ], + "Depends": [ + "R (>= 4.0.0)", + "methods" + ], + "Suggests": [ + "bit64", + "BiocStyle", + "knitr", + "rmarkdown", + "testthat", + "bench", + "dplyr", + "ggplot2", + "mockery", + "BiocParallel" + ], + "LinkingTo": [ + "Rhdf5lib" + ], + "SystemRequirements": "GNU make", + "biocViews": "Infrastructure, DataImport", + "Encoding": "UTF-8", + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.3", + "git_url": "https://git.bioconductor.org/packages/rhdf5", + "git_branch": "RELEASE_3_22", + "git_last_commit": "7f691e4", + "git_last_commit_date": "2025-12-02", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Bernd Fischer [aut], Mike Smith [aut] (ORCID: , Maintainer from 2017 to 2025), Gregoire Pau [aut], Martin Morgan [ctb], Daniel van Twisk [ctb], Hugo Gruson [cre] (ORCID: )", + "Maintainer": "Hugo Gruson " + }, + "rhdf5filters": { + "Package": "rhdf5filters", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "HDF5 Compression Filters", + "Authors@R": "c( person(\"Mike\", \"Smith\", role = c(\"aut\", \"ccp\"), comment = c(ORCID = \"0000-0002-7800-3848\") ), person(\"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\") ) )", + "Description": "Provides a collection of additional compression filters for HDF5 datasets. The package is intended to provide seemless integration with rhdf5, however the compiled filters can also be used with external applications.", + "License": "BSD_2_clause + file LICENSE", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Suggests": [ + "BiocStyle", + "knitr", + "rmarkdown", + "tinytest", + "rhdf5 (>= 2.47.7)" + ], + "SystemRequirements": "GNU make", + "URL": "https://github.com/Huber-group-EMBL/rhdf5filters", + "BugReports": "https://github.com/Huber-group-EMBL/rhdf5filters/issues", + "LinkingTo": [ + "Rhdf5lib" + ], + "RoxygenNote": "7.3.2", + "Encoding": "UTF-8", + "biocViews": "Infrastructure, DataImport", + "git_url": "https://git.bioconductor.org/packages/rhdf5filters", + "git_branch": "RELEASE_3_22", + "git_last_commit": "3465c24", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Mike Smith [aut, ccp] (ORCID: ), Hugo Gruson [cre] (ORCID: )", + "Maintainer": "Hugo Gruson " + }, + "sparseMatrixStats": { + "Package": "sparseMatrixStats", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Summary Statistics for Rows and Columns of Sparse Matrices", + "Authors@R": "person(\"Constantin\", \"Ahlmann-Eltze\", email = \"artjom31415@googlemail.com\", role = c(\"aut\", \"cre\"), comment = c(ORCID = \"0000-0002-3762-068X\"))", + "Description": "High performance functions for row and column operations on sparse matrices. For example: col / rowMeans2, col / rowMedians, col / rowVars etc. Currently, the optimizations are limited to data in the column sparse format. This package is inspired by the matrixStats package by Henrik Bengtsson.", + "License": "MIT + file LICENSE", + "Encoding": "UTF-8", + "RoxygenNote": "7.3.1", + "LinkingTo": [ + "Rcpp" + ], + "Imports": [ + "Rcpp", + "Matrix", + "matrixStats (>= 0.60.0)", + "methods" + ], + "Depends": [ + "MatrixGenerics (>= 1.5.3)" + ], + "Suggests": [ + "testthat (>= 2.1.0)", + "knitr", + "bench", + "rmarkdown", + "BiocStyle" + ], + "SystemRequirements": "C++11", + "VignetteBuilder": "knitr", + "URL": "https://github.com/const-ae/sparseMatrixStats", + "BugReports": "https://github.com/const-ae/sparseMatrixStats/issues", + "biocViews": "Infrastructure, Software, DataRepresentation", + "git_url": "https://git.bioconductor.org/packages/sparseMatrixStats", + "git_branch": "RELEASE_3_22", + "git_last_commit": "3acf348", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Constantin Ahlmann-Eltze [aut, cre] (ORCID: )", + "Maintainer": "Constantin Ahlmann-Eltze " } } } diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index a4722ba6..b3fe3938 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -22,13 +22,23 @@ using DataFrames, RCall using ..RRuntime: with_r_lock using ..Analysis: _palette_hex, R_WAIT_SECONDS export DifferentialConfig, ScalingRefusal, MASSUnavailable, ZCompositionsUnavailable, - OFFSET_METHODS, METHODS, tss_factors, rle_factors, size_factors, bh_adjust, - validate_counts, replace_zeros, clr_transform, fit_nb, fit_welch, - differential_abundance, volcano_chart + GlmGamPoiUnavailable, OFFSET_METHODS, METHODS, DISPERSIONS, tss_factors, + rle_factors, size_factors, bh_adjust, validate_counts, replace_zeros, + clr_transform, fit_nb, fit_welch, differential_abundance, volcano_chart # Methods accepted under `analysis.differential.method`. const METHODS = ("nb_glm", "clr_welch") +# What `nb_glm` does and how it tests, accepted under +# `analysis.differential.dispersion`: `"glm_nb"` fits each taxon with +# MASS::glm.nb and its maximum-likelihood theta; `"glmGamPoi"` fits one +# glmGamPoi::glm_gp across the taxa and reports glmGamPoi::test_de's own +# quasi-likelihood F-test on that fit, with each taxon's dispersion shrunk +# toward the cross-taxon trend. The second shares information across taxa, +# which steadies low-count taxa and small studies. `"glmGamPoi"` is refused +# for `clr_welch`. +const DISPERSIONS = ("glm_nb", "glmGamPoi") + # Offset (size factor) methods accepted under `analysis.differential.offset`. const OFFSET_METHODS = ("tss", "rle") @@ -43,7 +53,7 @@ const DEFAULT_DELTA = 0.65 """ DifferentialConfig(offset="tss", min_prevalence=0.0; method="nb_glm", - replacement_delta=0.65) + replacement_delta=0.65, dispersion="glm_nb") Settings for one differential abundance analysis, read from `analysis.differential` in the pipeline config. `method` is `"nb_glm"` or @@ -54,26 +64,49 @@ log-ratios, which are already free of library size. have at least one read to be tested. `replacement_delta`, in (0, 1), is the `frac` passed to `zCompositions::cmultRepl` by `clr_welch`: an imputed proportion above its taxon's smallest observed proportion is lowered to this fraction of it. +`dispersion` names what `nb_glm` does and how it tests (`"glm_nb"` or +`"glmGamPoi"`): `"glm_nb"` (the default) fits `MASS::glm.nb` per taxon with +its maximum-likelihood theta and Wald-tests the group coefficient; +`"glmGamPoi"` fits one `glmGamPoi::glm_gp` across the taxa passing +`min_prevalence`, with overdispersion shrinkage toward the cross-taxon trend, +and reports `glmGamPoi::test_de`'s quasi-likelihood F-test of the group +contrast — glmGamPoi's own published inference (Lund et al., 2012) — with its +log2 fold changes. `"glmGamPoi"` is refused for `clr_welch`. +Method and offset names are stripped and lowercased; dispersion names are +stripped and matched without regard to case, then stored as `"glm_nb"` or +`"glmGamPoi"`. Throws `ArgumentError` for unsupported names, incompatible +method and dispersion choices, or non-finite or out-of-range fractions, +including settings unused by the selected method. """ struct DifferentialConfig method :: String offset :: String min_prevalence :: Float64 replacement_delta :: Float64 + dispersion :: String function DifferentialConfig(offset::AbstractString="tss", min_prevalence::Real=0.0; method::AbstractString="nb_glm", - replacement_delta::Real=DEFAULT_DELTA) + replacement_delta::Real=DEFAULT_DELTA, + dispersion::AbstractString="glm_nb") m = lowercase(strip(String(method))) m in METHODS || throw(ArgumentError( "analysis.differential.method must be one of $(join(METHODS, ", ")), not '$method'")) o = lowercase(strip(String(offset))) o in OFFSET_METHODS || throw(ArgumentError( "analysis.differential.offset must be one of $(join(OFFSET_METHODS, ", ")), not '$offset'")) + d = strip(String(dispersion)) + di = findfirst(s -> lowercase(s) == lowercase(d), DISPERSIONS) + di === nothing && throw(ArgumentError( + "analysis.differential.dispersion must be one of $(join(DISPERSIONS, ", ")), not '$dispersion'")) + d = DISPERSIONS[di] + (d == "glm_nb" || m == "nb_glm") || throw(ArgumentError( + "analysis.differential.dispersion is a negative-binomial setting, but method is " * + "'$m'; only nb_glm has a dispersion to choose")) (isfinite(min_prevalence) && 0 <= min_prevalence <= 1) || throw(ArgumentError( "analysis.differential.min_prevalence must be a fraction in [0, 1], not $min_prevalence")) (isfinite(replacement_delta) && 0 < replacement_delta < 1) || throw(ArgumentError( "analysis.differential.replacement_delta must be in (0, 1), not $replacement_delta")) - new(m, o, Float64(min_prevalence), Float64(replacement_delta)) + new(m, o, Float64(min_prevalence), Float64(replacement_delta), d) end end @@ -115,6 +148,21 @@ Base.showerror(io::IO, ::ZCompositionsUnavailable) = print(io, "the R package zCompositions is not available, so zeros cannot be replaced for " * "clr_welch; restore the R library from renv.lock (Rscript -e 'renv::restore()')") +""" + GlmGamPoiUnavailable() + +Raised when the R package glmGamPoi, which provides the shrunken dispersion +trend that `analysis.differential.dispersion: glmGamPoi` asks for, cannot be +loaded. There is no fallback to another dispersion estimator: the method is +refused outright. +""" +struct GlmGamPoiUnavailable <: Exception end + +Base.showerror(io::IO, ::GlmGamPoiUnavailable) = print(io, + "the R package glmGamPoi is not available, so the shrunken dispersion cannot be " * + "fitted; restore the pinned R library (Rscript -e 'renv::restore()') " * + "or use analysis.differential.dispersion: glm_nb") + """ _geomean(x) -> Float64 @@ -243,7 +291,48 @@ da_result <- data.frame(status = rep("ok", da_n), note = rep("", da_n), estimate = rep(NA_real_, da_n), se = rep(NA_real_, da_n), statistic = rep(NA_real_, da_n), pvalue = rep(NA_real_, da_n), theta = rep(NA_real_, da_n), stringsAsFactors = FALSE) +# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa supplies glmGamPoi's +# own published inference -- glmGamPoi::test_de's quasi-likelihood F-test of +# the group coefficient. glm_gp refits the coefficients at its cross-taxon +# dispersion trend and shrinks each taxon's quasi-likelihood dispersion by +# empirical Bayes. glm_gp works genes x samples, so the counts are transposed, +# and size_factors takes the linear size factor (the per-taxon fits below use +# its log as their offset). The contrast reaches test_de as numeric weights on +# the design's columns, never as text: test_de parses and evaluates a text +# contrast as R code, so a condition label such as "run-B" or "grpB/run_2" +# would stop it, and a label that reads as R code would run. A fit- or +# test-level error fails every taxon. No fallback to another dispersion +# estimator or significance test is permitted. +da_glp <- NULL +da_glp_failed <- "" +if (da_use_glp) { + da_design <- model.matrix(~ da_g) + da_contrast <- as.numeric(colnames(da_design) == da_term) + da_gp <- tryCatch( + glmGamPoi::glm_gp(t(da_counts), design = da_design, + size_factors = exp(da_offset), overdispersion_shrinkage = TRUE), + error = function(e) e) + if (inherits(da_gp, "error")) { + da_glp_failed <- paste("glm_gp stopped with an error:", conditionMessage(da_gp)) + } else { + da_glp <- tryCatch( + glmGamPoi::test_de(da_gp, contrast = da_contrast, max_lfc = Inf, compute_lfc_se = TRUE), + error = function(e) e) + if (inherits(da_glp, "error")) { + da_glp_failed <- paste("test_de stopped with an error:", conditionMessage(da_glp)) + da_glp <- NULL + } else if (nrow(da_glp) != da_n) { + da_glp_failed <- "test_de returned no row for each taxon" + da_glp <- NULL + } + } +} for (da_j in seq_len(da_n)) { + if (nzchar(da_glp_failed)) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- da_glp_failed + next + } da_y <- da_counts[, da_j] if (length(unique(da_y)) < 2L) { da_result[da_j, "status"] <- "failed" @@ -251,62 +340,102 @@ for (da_j in seq_len(da_n)) { next } # Complete separation: with no reads in every sample of one group the group - # coefficient runs off to infinity, its standard error with it, and the Wald - # p-value collapses towards 1 although the fit reports convergence. That p is + # coefficient runs off to infinity and its test loses its footing. That p is # not evidence of no difference, so the taxon is refused with its reason. da_absent <- levels(da_g)[vapply(levels(da_g), function(l) all(da_y[da_g == l] == 0), logical(1))] if (length(da_absent) > 0L) { da_result[da_j, "status"] <- "failed" da_result[da_j, "note"] <- paste0("no reads in any '", da_absent[[1]], "' sample (complete ", - "separation): the fold change is unbounded and the Wald test ", + "separation): the fold change is unbounded and its test ", "is undefined, so no p-value is reported") next } - da_warns <- character(0) - da_fit <- tryCatch( - withCallingHandlers( - MASS::glm.nb(da_y ~ da_g + offset(da_offset), control = glm.control(maxit = 100)), - warning = function(w) { - da_warns <<- c(da_warns, conditionMessage(w)) - invokeRestart("muffleWarning") - }), - error = function(e) e) - if (inherits(da_fit, "error")) { - da_result[da_j, "status"] <- "failed" - da_result[da_j, "note"] <- paste("glm.nb stopped with an error:", conditionMessage(da_fit)) - next - } - da_warns <- unique(da_warns) - if (!isTRUE(da_fit[["converged"]]) || !is.null(da_fit[["th.warn"]]) || - any(grepl("iteration limit|alternation limit|did not converge|NaNs produced", da_warns))) { - da_result[da_j, "status"] <- "failed" - da_result[da_j, "note"] <- paste("the fit did not converge:", - paste(c(da_fit[["th.warn"]], da_warns), collapse = "; ")) - next - } - da_co <- summary(da_fit)[["coefficients"]] - if (!(da_term %in% rownames(da_co))) { - da_result[da_j, "status"] <- "failed" - da_result[da_j, "note"] <- "the group coefficient is not estimable (aliased)" - next - } - da_row <- da_co[da_term, ] - if (!all(is.finite(da_row)) || da_row[[4]] < 0 || da_row[[4]] > 1) { - da_result[da_j, "status"] <- "failed" - da_result[da_j, "note"] <- "the fit returned a non-finite estimate, standard error or p-value" - next - } - da_theta <- da_fit[["theta"]] - da_result[da_j, c("estimate", "se", "statistic", "pvalue")] <- unname(da_row[1:4]) - da_result[da_j, "theta"] <- da_theta - if (!is.finite(da_theta) || da_theta >= da_theta_upper) { - da_result[da_j, "status"] <- "boundary" - da_result[da_j, "note"] <- "the dispersion parameter theta reached its upper bound: these counts show no overdispersion, so the fit is effectively Poisson" - } else if (da_theta <= da_theta_lower) { - da_result[da_j, "status"] <- "boundary" - da_result[da_j, "note"] <- "the dispersion parameter theta reached its lower bound: the variance is extreme relative to the mean" - } else if (length(da_warns) > 0L) { - da_result[da_j, "note"] <- paste("R warned:", paste(da_warns, collapse = "; ")) + if (da_use_glp) { + # glmGamPoi's own test on its own fit. `lfc` is the log2 fold change of + # the group contrast and `lfc_se` its standard error; both are converted + # to the natural-log scale so every method reports one effect scale. The + # reported theta is the inverse of glm_gp's negative-binomial + # overdispersion for the taxon (`overdispersions`: variance = mu + mu^2 * + # overdispersion), the same parameter glm.nb's theta inverts. It is the + # taxon's own estimate before shrinkage: the coefficients and the F-test's + # reduced fit use the cross-taxon trend, and the F ratio divides by the + # shrunken quasi-likelihood dispersion; only lfc_se, from predict(), is + # weighted by this estimate. Zero overdispersion is the Poisson boundary: + # the row is marked, as under glm_nb, and keeps its statistics, since the + # test does not read it. An unusable overdispersion fails only its own + # taxon. + da_od <- da_gp[["overdispersions"]][da_j] + if (!is.finite(da_od) || da_od < 0) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the negative-binomial overdispersion for this taxon is not usable" + next + } + da_est <- da_glp[["lfc"]][da_j] * log(2) + da_se <- da_glp[["lfc_se"]][da_j] * log(2) + da_stat <- da_glp[["f_statistic"]][da_j] + da_p <- da_glp[["pval"]][da_j] + if (!is.finite(da_est) || !is.finite(da_se) || !is.finite(da_stat) || + !is.finite(da_p) || da_p < 0 || da_p > 1) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the fit returned a non-finite estimate, standard error or p-value" + next + } + da_result[da_j, c("estimate", "se", "statistic", "pvalue", "theta")] <- + c(da_est, da_se, da_stat, da_p, 1 / da_od) + if (da_od == 0) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "the taxon's negative-binomial overdispersion estimate is zero: these counts show no overdispersion beyond Poisson, so theta is unbounded; the test used the cross-taxon trend, so its statistics stand" + } + } else { + # "glm_nb" estimates theta per taxon inside MASS::glm.nb's alternating + # fit, then Wald-tests the group coefficient at fixed scale 1, which is + # summary.negbin's own default. + da_warns <- character(0) + da_fit <- tryCatch( + withCallingHandlers( + MASS::glm.nb(da_y ~ da_g + offset(da_offset), control = glm.control(maxit = 100)), + warning = function(w) { + da_warns <<- c(da_warns, conditionMessage(w)) + invokeRestart("muffleWarning") + }), + error = function(e) e) + if (inherits(da_fit, "error")) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- paste("glm.nb stopped with an error:", conditionMessage(da_fit)) + next + } + da_warns <- unique(da_warns) + if (!isTRUE(da_fit[["converged"]]) || !is.null(da_fit[["th.warn"]]) || + any(grepl("iteration limit|alternation limit|did not converge|NaNs produced", da_warns))) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- paste("the fit did not converge:", + paste(c(da_fit[["th.warn"]], da_warns), collapse = "; ")) + next + } + da_co <- summary(da_fit, dispersion = 1)[["coefficients"]] + if (!(da_term %in% rownames(da_co))) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the group coefficient is not estimable (aliased)" + next + } + da_row <- da_co[da_term, ] + if (!all(is.finite(da_row)) || da_row[[4]] < 0 || da_row[[4]] > 1) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the fit returned a non-finite estimate, standard error or p-value" + next + } + da_theta <- da_fit[["theta"]] + da_result[da_j, c("estimate", "se", "statistic", "pvalue")] <- unname(da_row[1:4]) + da_result[da_j, "theta"] <- da_theta + if (!is.finite(da_theta) || da_theta >= da_theta_upper) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "the dispersion parameter theta reached its upper bound: these counts show no overdispersion, so the fit is effectively Poisson" + } else if (da_theta <= da_theta_lower) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "the dispersion parameter theta reached its lower bound: the variance is extreme relative to the mean" + } else if (length(da_warns) > 0L) { + da_result[da_j, "note"] <- paste("R warned:", paste(da_warns, collapse = "; ")) + } } } """ @@ -319,36 +448,75 @@ A finite number as `Float64`, or `nothing` for R's NA, NaN or an infinity. _num(x) = (ismissing(x) || isnothing(x) || !isfinite(x)) ? nothing : Float64(x) """ - fit_nb(counts, groups, offset; levels) -> DataFrame + fit_nb(counts, groups, offset; levels, dispersion="glm_nb") -> DataFrame + +Fit a negative-binomial GLM with a log link to each column of the samples x +taxa integer matrix `counts`. `offset` contains the natural logarithm of each +sample's size factor. `levels` is `(reference, contrast)`, so the estimate is +the natural-log fold change of `contrast` over `reference`. +`dispersion` names what the fits do and how they test: `"glm_nb"` (the +default) takes `MASS::glm.nb`'s maximum-likelihood theta per taxon and +Wald-tests the group coefficient at fixed scale 1; `"glmGamPoi"` fits one +`glmGamPoi::glm_gp` across all supplied taxa with overdispersion shrinkage +toward the cross-taxon trend and reports `glmGamPoi::test_de`'s +quasi-likelihood F-test of the group contrast — its `lfc` (log2 fold change) +and `lfc_se`, both converted to the natural-log scale so the columns carry one +effect scale. In that path `statistic` is the quasi-likelihood F, `pvalue` +its p-value, and `theta` the inverse of `glm_gp`'s negative-binomial +overdispersion for the taxon (its `overdispersions` element, the taxon's own +estimate before shrinkage), the same parameter `MASS::glm.nb`'s theta inverts; +it is `nothing` where that overdispersion is zero. The contrast is passed to +`test_de` as numeric weights, so condition labels are never parsed as R code. +Names are stripped and matched without regard to case. -Fit `MASS::glm.nb(y ~ group + offset(offset))` to each column of the samples x -taxa integer matrix `counts`. `levels` is `(reference, contrast)`, so the -estimate is the natural-log fold change of `contrast` over `reference`. Returns one row per taxon with `status` (`ok`, `boundary` or `failed`), `note`, `estimate`, `se`, `statistic`, `pvalue` and `theta`; a failed taxon has -`nothing` for every statistic. Throws `MASSUnavailable` when MASS cannot be -loaded and `RBusyError` when the R runtime stays busy. +`nothing` for every statistic. Boundary fits retain their statistics: under +`"glm_nb"` theta at or above 1e7 marks the Poisson end and theta at or below +1e-8 marks extreme variance; under `"glmGamPoi"` a zero negative-binomial +overdispersion marks the Poisson boundary, and the row keeps its statistics, +since the F-test divides by the shrunken quasi-likelihood dispersion, which +stays defined there. Constant counts, complete separation and +unusable coefficients produce failed rows. Errors caught from the model-fitting +calls also produce failed rows; a caught `glm_gp` or `test_de` error fails every +taxon, with no fallback to `"glm_nb"`. + +Throws `ArgumentError` for an unsupported dispersion name or mismatched +sample, group and offset lengths, `MASSUnavailable` when MASS cannot be loaded, +`GlmGamPoiUnavailable` when `"glmGamPoi"` is selected and its package cannot +be loaded, and `RBusyError` when the R runtime stays busy. R errors outside the +caught model-fitting calls propagate. """ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:AbstractString}, - offset::AbstractVector{<:Real}; levels::NTuple{2,String}) + offset::AbstractVector{<:Real}; levels::NTuple{2,String}, + dispersion::AbstractString="glm_nb") + di = findfirst(s -> lowercase(s) == lowercase(strip(String(dispersion))), DISPERSIONS) + di === nothing && throw(ArgumentError( + "dispersion must be one of $(join(DISPERSIONS, ", ")), not '$dispersion'")) + use_glp = DISPERSIONS[di] == "glmGamPoi" size(counts, 1) == length(groups) == length(offset) || throw(ArgumentError( "counts have $(size(counts, 1)) samples, groups $(length(groups)), offset $(length(offset))")) raw = with_r_lock(; timeout=R_WAIT_SECONDS[]) do RCall.rcopy(RCall.reval("requireNamespace('MASS', quietly = TRUE)")) || throw(MASSUnavailable()) + use_glp && !RCall.rcopy(RCall.reval("requireNamespace('glmGamPoi', quietly = TRUE)")) && + throw(GlmGamPoiUnavailable()) RCall.globalEnv[:da_counts] = Matrix{Int}(counts) RCall.globalEnv[:da_group] = String.(groups) RCall.globalEnv[:da_levels] = collect(levels) RCall.globalEnv[:da_offset] = Float64.(offset) RCall.globalEnv[:da_theta_upper] = THETA_UPPER RCall.globalEnv[:da_theta_lower] = THETA_LOWER + RCall.globalEnv[:da_use_glp] = use_glp try RCall.reval(_FIT_R) DataFrame(RCall.rcopy(RCall.reval("da_result"))) finally RCall.reval("rm(list = intersect(ls(), c('da_counts', 'da_group', 'da_levels', " * - "'da_offset', 'da_theta_upper', 'da_theta_lower', 'da_g', 'da_term', " * - "'da_n', 'da_result', 'da_j', 'da_y', 'da_warns', 'da_fit', 'da_co', " * - "'da_row', 'da_theta', 'da_absent'))); invisible(gc())") + "'da_offset', 'da_theta_upper', 'da_theta_lower', 'da_use_glp', 'da_g', " * + "'da_term', 'da_n', 'da_result', 'da_design', 'da_contrast', 'da_glp', " * + "'da_glp_failed', 'da_gp', " * + "'da_j', 'da_y', 'da_absent', 'da_od', 'da_est', 'da_se', 'da_stat', " * + "'da_p', 'da_warns', 'da_fit', 'da_co', 'da_row', 'da_theta'))); invisible(gc())") end end DataFrame(status = String.(raw.status), @@ -589,7 +757,12 @@ labelled `reference` in `groups`, by `config.method`. `nb_glm`: size factors are computed on the full matrix by `config.offset` and enter each negative-binomial model as `log(factor)`; the estimate is a -natural-log fold change. +natural-log fold change. `config.dispersion` selects the fits and the test: +`"glm_nb"` (per-taxon `MASS::glm.nb` with its maximum-likelihood theta) or +`"glmGamPoi"` (one `glmGamPoi::glm_gp` fit across the taxa passing the +prevalence filter, with overdispersion shrinkage, and `glmGamPoi::test_de`'s +quasi-likelihood F-test of the group contrast). The result's `config` +includes the selected dispersion name. `clr_welch`, in this order: (1) `min_prevalence` chooses the taxa, and taxa with no reads in any sample are always left out; (2) the zeros of those taxa only @@ -614,8 +787,11 @@ effect size and its label. Throws `ArgumentError` for unusable input (including a zero replacement that cannot be made), `ScalingRefusal` when the offset cannot be computed, -`MASSUnavailable` or `ZCompositionsUnavailable` when the R package a method -needs is missing, and an `ErrorException` when no taxon could be fitted at all. +`MASSUnavailable`, `GlmGamPoiUnavailable` or `ZCompositionsUnavailable` when a +required R package is missing, `RBusyError` when the R runtime stays busy, and +an `ErrorException` when no taxon could be fitted at all, including when the +`glm_gp` fit fails. Uncaught R errors from the fitting and zero-replacement +helpers propagate. """ function differential_abundance(counts::AbstractMatrix{<:Real}, samples::AbstractVector{<:AbstractString}, @@ -664,7 +840,8 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, factors = size_factors(counts, samples, config.offset) tested = findall(>=(config.min_prevalence), prevalence) fits = isempty(tested) ? nothing : - fit_nb(Int.(counts[:, tested]), groups, log.(factors); levels) + fit_nb(Int.(counts[:, tested]), groups, log.(factors); levels, + dispersion=config.dispersion) if !isnothing(fits) for (k, j) in enumerate(tested) f = fits[k, :] @@ -675,14 +852,23 @@ function differential_abundance(counts::AbstractMatrix{<:Real}, "pvalue" => f.pvalue, "dispersion_theta" => f.theta)) end end - method_text = "Negative-binomial GLM per taxon (MASS::glm.nb, log link) with a " * - "$(uppercase(config.offset)) size-factor offset; Wald test of the group " * - "coefficient; Benjamini-Hochberg adjustment over the fitted taxa" + nb_fit = config.dispersion == "glm_nb" ? "MASS::glm.nb" : "glmGamPoi::glm_gp" + if config.dispersion == "glm_nb" + method_text = "Negative-binomial GLM per taxon ($nb_fit, log link) with a " * + "$(uppercase(config.offset)) size-factor offset; Wald test of the group " * + "coefficient; Benjamini-Hochberg adjustment over the fitted taxa" + else + method_text = "Negative-binomial GLM per taxon ($nb_fit, log link) with a " * + "$(uppercase(config.offset)) size-factor offset; quasi-likelihood F-test " * + "of the group contrast (glmGamPoi::test_de) with overdispersion shrinkage; " * + "Benjamini-Hochberg adjustment over the fitted taxa" + end effect = Dict("key" => "log2_fold_change", "label" => "log2 fold change") size_factor_rows = [Dict("sample" => String(samples[i]), "group" => String(groups[i]), "factor" => factors[i]) for i in eachindex(samples)] config_echo = Dict{String,Any}("method" => config.method, "offset" => config.offset, - "min_prevalence" => config.min_prevalence) + "min_prevalence" => config.min_prevalence, + "dispersion" => config.dispersion) else (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( "Welch's t-test needs at least 2 samples in each group to estimate its variance; " * diff --git a/src/core/r_runtime.jl b/src/core/r_runtime.jl index 0fa982da..5fca35e3 100644 --- a/src/core/r_runtime.jl +++ b/src/core/r_runtime.jl @@ -25,6 +25,7 @@ const _R_LOCK = ReentrantLock() const REQUIRED_R_PACKAGES = [ "dada2" => "the DADA2 stages that run on this machine", "vegan" => "NMDS, PERMANOVA and PERMDISP", + "glmGamPoi" => "shared dispersion-trend differential abundance", "MASS" => "negative-binomial differential abundance", "zCompositions" => "CLR differential abundance (clr_welch)", ] diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 9c1f0129..29875bff 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1196,7 +1196,10 @@ end _differential_config(study) -> DifferentialConfig The study's `analysis.differential` settings, resolved through the config -cascade. +cascade. Missing or `nothing` values use the defaults in `DifferentialConfig`, +including `dispersion="glm_nb"`. Invalid settings propagate the constructor's +`ArgumentError`; values that cannot be converted to strings or numbers also +raise their conversion errors. """ function _differential_config(study::String) resolved = _resolved_run_config(study, nothing, nothing, nothing) @@ -1204,7 +1207,8 @@ function _differential_config(study::String) DifferentialConfig(String(val("analysis.differential.offset", "tss")), Float64(val("analysis.differential.min_prevalence", 0.0)); method=String(val("analysis.differential.method", "nb_glm")), - replacement_delta=Float64(val("analysis.differential.replacement_delta", 0.65))) + replacement_delta=Float64(val("analysis.differential.replacement_delta", 0.65)), + dispersion=String(val("analysis.differential.dispersion", "glm_nb"))) end """ @@ -1283,6 +1287,7 @@ function _differential(study::String, body) e isa ScalingRefusal && return json_error(422, "scaling_refused", sprint(showerror, e)) e isa MASSUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) e isa ZCompositionsUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) + e isa GlmGamPoiUnavailable && return json_error(503, "r_unavailable", sprint(showerror, e)) e isa ErrorException && return json_error(422, "no_taxon_fitted", e.msg) rethrow() end diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl index 2f927ba1..e434269f 100644 --- a/src/server/routes/config.jl +++ b/src/server/routes/config.jl @@ -108,6 +108,7 @@ const _ENUM_VALUES = Dict( "analysis.beta.normalisation" => ("hellinger", "none", "rarefy", "srs"), "analysis.differential.method" => Differential.METHODS, "analysis.differential.offset" => Differential.OFFSET_METHODS, + "analysis.differential.dispersion" => Differential.DISPERSIONS, "phylogeny.reference.align.strategy" => Validation.PHYLOGENY_ALIGN_STRATEGIES, "phylogeny.placement.align.strategy" => Validation.PHYLOGENY_ALIGN_STRATEGIES, "phylogeny.reference.trim.method" => Validation.PHYLOGENY_TRIM_METHODS, diff --git a/test/integration/test_server.jl b/test/integration/test_server.jl index 184e70fa..ef63be81 100644 --- a/test/integration/test_server.jl +++ b/test/integration/test_server.jl @@ -70,6 +70,21 @@ status_exception=false) @test r5.status in (403, 404) + ## PATCH /api/v1/config stores the dispersion spelling over the + ## real listener and refuses the pre-release "ml". + r6 = HTTP.patch("http://localhost:$port/api/v1/config", + ["Content-Type" => "application/json"], + JSON3.write(Dict("analysis.differential.dispersion" => "glmGamPoi")); + status_exception=false) + @test r6.status == 200 + saved = YAML.load_file(joinpath(tmp_root, "config", "pipeline.yml")) + @test saved["analysis"]["differential"]["dispersion"] == "glmGamPoi" + r7 = HTTP.patch("http://localhost:$port/api/v1/config", + ["Content-Type" => "application/json"], + JSON3.write(Dict("analysis.differential.dispersion" => "ml")); + status_exception=false) + @test r7.status == 400 + @info "Server smoke tests passed on port $port" end finally diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 0c7968f9..868b96b3 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -31,6 +31,15 @@ const HAVE_MASS = HAVE_R && _r_has("MASS") HAVE_MASS || @info "R or MASS unavailable - the negative-binomial fit tests are SKIPPED" const HAVE_ZCOMP = HAVE_R && _r_has("zCompositions") HAVE_ZCOMP || @info "R or zCompositions unavailable - the clr_welch tests are SKIPPED" +const HAVE_GLMGAMPOI = HAVE_MASS && _r_has("glmGamPoi") + +# Issue #37's test contract: the direct glmGamPoi comparisons fail, not skip, +# in CI when glmGamPoi is missing. Locally they skip loudly. This guard is the +# first statement of every glmGamPoi-dependent testset. +const _GLMGAMPOI_IN_CI = get(ENV, "CI", "") == "true" +HAVE_GLMGAMPOI || @info (_GLMGAMPOI_IN_CI ? + "glmGamPoi unavailable IN CI - the shrunken-dispersion tests FAIL" : + "R, MASS or glmGamPoi unavailable - the shrunken-dispersion fit tests are SKIPPED") # Rows per run: Up is split over two ASVs so the route must pool them by genus. const _DA_ROWS = Dict( @@ -172,6 +181,16 @@ end @test_throws ArgumentError DifferentialConfig(; method="ancombc") # The pre-review name, never released, is not accepted. @test_throws ArgumentError DifferentialConfig(; method="clr_lm") + # The dispersion option: glm_nb is the default and keeps nb_glm's + # per-taxon maximum-likelihood theta; glmGamPoi is the alternative and + # is a negative-binomial setting only. + @test DifferentialConfig().dispersion == "glm_nb" + @test DifferentialConfig(dispersion="glmGamPoi").dispersion == "glmGamPoi" + @test DifferentialConfig(dispersion=" GLMGAMPOI ").dispersion == "glmGamPoi" + @test DifferentialConfig(dispersion=" Glm_Nb ").dispersion == "glm_nb" + @test_throws ArgumentError DifferentialConfig(dispersion="shrinkage") + @test_throws ArgumentError DifferentialConfig(method="clr_welch", dispersion="glmGamPoi") + @test DifferentialConfig(method="clr_welch").dispersion == "glm_nb" for d in (0, 1, -0.5, NaN) @test_throws ArgumentError DifferentialConfig(; method="clr_welch", replacement_delta=d) end @@ -212,6 +231,9 @@ end byt = Dict(r["taxon"] => r for r in res["rows"]) @test res["status"] == "ok" @test res["diagnostics"]["n_tested"] == 3 + @test res["method"] == "Negative-binomial GLM per taxon (MASS::glm.nb, log link) " * + "with a TSS size-factor offset; Wald test of the group " * + "coefficient; Benjamini-Hochberg adjustment over the fitted taxa" @test all(r -> r["status"] == "ok", res["rows"]) @test isapprox(byt["t_up"]["estimate"], log(4); atol=0.5) @test isapprox(byt["t_up"]["log2_fold_change"], 2; atol=0.8) @@ -245,7 +267,138 @@ end @test r["pvalue"] ≈ direct[3] rtol=1e-6 @test r["dispersion_theta"] ≈ direct[4] rtol=1e-6 end + @test res["config"]["dispersion"] == "glm_nb" + end + end + end + + @testset "glmGamPoi path matches a direct glmGamPoi::test_de call" begin + if !HAVE_GLMGAMPOI + # Issue #37's contract: these comparisons fail, not skip, in CI. + @test !_GLMGAMPOI_IN_CI + @test_skip false + else + res = differential_abundance(effect, effect_samples, effect_taxa, effect_groups; + reference="A", contrast="B", + config=DifferentialConfig(dispersion="glmGamPoi")) + @test res["status"] == "ok" + @test res["method"] == "Negative-binomial GLM per taxon (glmGamPoi::glm_gp, log link) " * + "with a TSS size-factor offset; quasi-likelihood F-test of the " * + "group contrast (glmGamPoi::test_de) with overdispersion shrinkage; " * + "Benjamini-Hochberg adjustment over the fitted taxa" + @test res["config"]["dispersion"] == "glmGamPoi" + byt = Dict(r["taxon"] => r for r in res["rows"]) + off = log.(tss_factors(effect, effect_samples)) + g = effect_groups + # The reference is a direct glmGamPoi::test_de call here in R, not + # the fitted path under test: the same model and the same published + # quasi-likelihood F-test, computed independently. lfc is log2 and + # lfc_se its standard error, so both are scaled by log(2) to meet + # the natural-log convention of every method's estimate and se. + # theta inverts glm_gp's negative-binomial overdispersion, not the + # quasi-likelihood scale the F-test divides by. One matrix column + # per quantity, one row per taxon. + direct = rcopy(R""" + local({ + g <- factor($g, levels = c("A", "B")) + gp <- glmGamPoi::glm_gp(t($effect), design = model.matrix(~ g), + size_factors = exp($off), overdispersion_shrinkage = TRUE) + td <- glmGamPoi::test_de(gp, contrast = "gB", max_lfc = Inf, + compute_lfc_se = TRUE) + cbind(td$pval, td$f_statistic, td$lfc * log(2), td$lfc_se * log(2), + 1 / gp$overdispersions) + }) + """) + for (j, t) in enumerate(effect_taxa) + r = byt[t] + @test r["status"] == "ok" + @test r["pvalue"] ≈ direct[j, 1] rtol=1e-6 + @test r["statistic"] ≈ direct[j, 2] rtol=1e-6 + @test r["estimate"] ≈ direct[j, 3] rtol=1e-6 + @test r["standard_error"] ≈ direct[j, 4] rtol=1e-6 + @test r["log2_fold_change"] ≈ direct[j, 3] / log(2) rtol=1e-6 + @test r["dispersion_theta"] ≈ direct[j, 5] rtol=1e-6 + end + end + end + + @testset "glmGamPoi unusable overdispersions are isolated per taxon" begin + if !HAVE_GLMGAMPOI + # Issue #37's contract: these comparisons fail, not skip, in CI. + @test !_GLMGAMPOI_IN_CI + @test_skip false + else + # test_de refits the reduced model from the shrinkage list's + # dispersion trend and divides by ql_disp_shrunken, so altering + # `overdispersions` after the fit changes theta, the boundary check + # and lfc_se (predict() weights by it), never the F statistic or + # its p-value. + MetaManifold.RRuntime.with_r_lock() do + reval(raw""" + test_original_gp <- get("glm_gp", asNamespace("glmGamPoi")) + test_od <- 0 + assignInNamespace("glm_gp", function(...) { + fit <- test_original_gp(...) + fit$overdispersions[1] <- test_od + fit + }, ns="glmGamPoi") + """) + try + off = log.(tss_factors(effect, effect_samples)) + # A zero overdispersion is the Poisson boundary: marked, with + # the statistics kept and no finite theta, as under glm_nb. + # A negative or non-finite one fails its taxon alone. + for (od, status, tested) in ((0.0, "boundary", true), (-1.0, "failed", false), + (NaN, "failed", false)) + RCall.globalEnv[:test_od] = od + rows = fit_nb(effect, effect_groups, off; levels=("A", "B"), + dispersion="glmGamPoi") + @test rows.status[1] == status + @test isnothing(rows.pvalue[1]) == !tested + @test isnothing(rows.theta[1]) + @test occursin("overdispersion", rows.note[1]) + @test all(rows.status[2:end] .== "ok") + end + finally + reval(raw""" + assignInNamespace("glm_gp", test_original_gp, ns="glmGamPoi") + rm(test_original_gp, test_od) + """) + end + end + end + end + + @testset "Condition labels never change the fits or reach R as code" begin + if !HAVE_GLMGAMPOI + # Issue #37's contract: these comparisons fail, not skip, in CI. + @test !_GLMGAMPOI_IN_CI + @test_skip false + else + # Condition labels join group, run and sub-group names with '/', and + # those names may carry '-' and '.'. test_de parses and evaluates a + # text contrast as R code, so the contrast goes to it as numeric + # weights. These labels, and one that reads as R code, must give the + # A/B results under both dispersions, and that code must never run. + # The coded label opens with '[' so that the text the old path built, + # "da_g[...]", indexed the live factor da_g and so ran the code. + off = log.(tss_factors(effect, effect_samples)) + evaluated = "mm_contrast_label_evaluated" + labels = (("run-A", "run-B"), ("grpA/run_1", "grpB/run_2"), ("2024-01", "2024-02"), + ("A", "[assign(\"$evaluated\", TRUE, envir = globalenv())]")) + for dispersion in ("glm_nb", "glmGamPoi") + base = fit_nb(effect, effect_groups, off; levels=("A", "B"), dispersion) + @test all(base.status .== "ok") + for (ref, con) in labels + groups = [g == "A" ? ref : con for g in effect_groups] + rows = fit_nb(effect, groups, off; levels=(ref, con), dispersion) + @test rows.status == base.status + for col in (:estimate, :se, :statistic, :pvalue, :theta) + @test all(isapprox.(rows[!, col], base[!, col]; rtol=1e-10)) + end + end end + @test !rcopy(R"exists($evaluated, envir = globalenv())") end end @@ -384,6 +537,16 @@ end @test isnothing(SV._config_value_error("analysis.differential.min_prevalence", 0.1)) @test isnothing(SV._config_value_error("analysis.differential.method", "clr_welch")) @test !isnothing(SV._config_value_error("analysis.differential.method", "ancombc")) + # Like every other enum key, the config route takes exact spellings. + for dispersion in ("glm_nb", "glmGamPoi") + @test isnothing(SV._config_value_error("analysis.differential.dispersion", dispersion)) + end + # "ml" was the pre-release spelling of glm_nb and is not accepted; + # case and spaces are DifferentialConfig's lenience, not the route's. + for dispersion in ("shrinkage", "ml", "", nothing, 1, true, ["glm_nb"], + " GLM_NB ", " GLMGAMPOI ", "GLMGAMPOI") + @test !isnothing(SV._config_value_error("analysis.differential.dispersion", dispersion)) + end end @testset "CLR transform" begin @@ -751,4 +914,55 @@ end @test rcopy(Bool, reval("requireNamespace('zCompositions', quietly = TRUE)")) == HAVE_ZCOMP end end + + ## Without glmGamPoi, dispersion: glmGamPoi refuses with + # GlmGamPoiUnavailable and the route answers 503. As above, a + # `requireNamespace` bound in R's global environment reports glmGamPoi + # missing and defers to base for every other package, so the production + # probe sees the package as absent on any machine. MASS is probed first + # and must be real, so this testset needs it. + @testset "nb_glm without glmGamPoi refuses with GlmGamPoiUnavailable and a 503" begin + if !HAVE_MASS + @test_skip false # without MASS the glm_nb path's own guard fires first + else + reval("""requireNamespace <- function(package, ...) { + if (identical(package, "glmGamPoi")) FALSE + else base::requireNamespace(package, ...) + }""") + try + @test !rcopy(Bool, reval("requireNamespace('glmGamPoi', quietly = TRUE)")) + e = try + fit_nb([5 6; 7 3; 9 4; 4 8], ["A", "A", "B", "B"], [0.0, 0.0, 0.0, 0.0]; + levels=("A", "B"), dispersion="glmGamPoi") + nothing + catch err + err + end + @test e isa GlmGamPoiUnavailable + @test occursin("renv::restore()", sprint(showerror, e)) + @test occursin("dispersion: glm_nb", sprint(showerror, e)) + _da_fixture() do + root = MetaManifold.Server.ServerState._root[] + write(joinpath(root, "data", "studyD", "pipeline.yml"), + "analysis:\n differential:\n dispersion: glmGamPoi\n") + response = _da_post("studyD", Dict("runs" => [Dict("run" => "runA"), Dict("run" => "runB")], + "table" => "merged", "rank" => "Genus")) + body = JSON3.read(String(response.body)) + @test response.status == 503 + @test body.error == "r_unavailable" + @test occursin("glmGamPoi", body.message) + end + # The override defers to base for everything else, so the + # default dispersion still fits normally with it in place: + # the injected-effect dataset, measured with MASS::glm.nb, + # reaches status ok on its own. + r = fit_nb(effect, effect_groups, log.(tss_factors(effect, effect_samples)); + levels=("A", "B")) + @test all(r.status .== "ok") + finally + reval("rm(requireNamespace)") + end + @test !rcopy(Bool, reval("exists('requireNamespace', envir = globalenv(), inherits = FALSE)")) + end + end end