From 02ee915f7d73b6d10c6e7d57b633d28c7b71ca20 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sun, 18 Oct 2026 22:20:23 +0100 Subject: [PATCH 1/3] Add glmGamPoi to the pinned R dependencies The shrunken dispersion trend that analysis.differential.dispersion: glmGamPoi (issue #37) will use is fitted by glmGamPoi::glm_gp, so the package must be in the restored library and every list that enumerates it. renv.lock gains glmGamPoi 1.22.0 and its 11-package closure on R 4.5.0 (Bioc 3.22): assorthead 1.4.0, beachmat 2.26.0, DelayedMatrixStats 1.32.0, HDF5Array 1.38.0, h5mread 1.2.1, RcppArmadillo 15.6.2-1, rhdf5 2.54.1, rhdf5filters 1.22.0, Rhdf5lib 1.32.0, SingleCellExperiment 1.32.0 and sparseMatrixStats 1.22.0; DelayedArray, MatrixGenerics and SummarizedExperiment were already present and are unchanged. R/_renv_dependencies.R loads glmGamPoi, the runtime probe and install.jl's availability check require it, and both CI R-library installations install it, so the pinned closure and the required-package lists stay in step. Dependency plumbing only; no behavior change. Prepares for #37. Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- .github/workflows/ci.yml | 4 +- R/_renv_dependencies.R | 1 + install.jl | 4 +- renv.lock | 559 +++++++++++++++++++++++++++++++++++++++ src/core/r_runtime.jl | 1 + 5 files changed, 565 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 86fd896..573768b 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 1e91aca..589612f 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/install.jl b/install.jl index 964e00d..3be35d2 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 cfd72d4..9944f58 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/core/r_runtime.jl b/src/core/r_runtime.jl index 0fa982d..5fca35e 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)", ] From efe1df1a39cde30291b03acc20ba3d496f7e404a Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sun, 18 Oct 2026 22:20:24 +0100 Subject: [PATCH 2/3] Add glmGamPoi dispersion option to nb_glm differential abundance Implement #37. analysis.differential.dispersion chooses the fits and the test for nb_glm. The default glm_nb keeps MASS::glm.nb with its maximum-likelihood theta per taxon and a Wald test of the group coefficient, so #24's results are unchanged. The alternative glmGamPoi fits one glmGamPoi::glm_gp over the taxa passing min_prevalence with overdispersion shrinkage and reports glmGamPoi::test_de's quasi-likelihood F-test of the group contrast -- its own published inference -- with the returned log2 fold changes and their standard errors converted to the natural-log scale so every method reports one effect scale. The reported dispersion is the inverse of the shrunken quasi-likelihood dispersion the test used. A glm_gp or test_de error fails every taxon; a zero shrunken dispersion marks only its own taxon as a Poisson boundary and withholds inference there, and an unusable dispersion fails only its own taxon. Missing glmGamPoi raises GlmGamPoiUnavailable (HTTP 503 on the analysis route), never a silent fallback to another estimator or test. glmGamPoi is refused for clr_welch; the pair rule lives in DifferentialConfig and is reported at analysis time, while the config routes keep their ordinary per-key validation with the dispersion spellings registered alongside the other enums. The frontend offers the dispersion enum in method settings, and the tests compare the glmGamPoi path against direct glm_gp + test_de calls in R -- p-values, F statistics and log fold changes -- failing (not skipping) in CI when glmGamPoi is missing. Closes #37. Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- config/defaults/pipeline.yml | 9 + frontend/src/api/types.ts | 2 + .../src/components/PipelineStages.test.ts | 10 + frontend/src/components/PipelineStages.tsx | 2 + src/analysis/differential.jl | 316 ++++++++++++++---- src/server/routes/analysis.jl | 9 +- src/server/routes/config.jl | 1 + test/integration/test_server.jl | 15 + test/unit/test_differential.jl | 169 ++++++++++ 9 files changed, 459 insertions(+), 74 deletions(-) create mode 100644 frontend/src/components/PipelineStages.test.ts diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 8181ae6..2181a23 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -358,6 +358,15 @@ 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. A zero shrunken + # dispersion withholds inference for its taxon. 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 87457da..e6ee284 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -439,6 +439,8 @@ export interface DifferentialResult { min_prevalence: number /** nb_glm only. */ offset?: 'tss' | 'rle' + /** nb_glm only: where each taxon's theta comes from. */ + 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 0000000..dd950da --- /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 8512926..b938f77 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/src/analysis/differential.jl b/src/analysis/differential.jl index d2c78de..5e46c3f 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,43 @@ 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 each taxon's +# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa, with each taxon's +# overdispersion shrunk toward the cross-taxon trend, supplies glmGamPoi's own +# published inference -- glmGamPoi::test_de's quasi-likelihood F-test of the +# group coefficient. Its genes x samples convention, so +# the counts are transposed, and size_factors takes the linear size factor +# (the per-taxon fits below use its log as their offset). 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_gp <- tryCatch( + glmGamPoi::glm_gp(t(da_counts), design = model.matrix(~ da_g), + 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_term, 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 +335,98 @@ 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 + # per-taxon scale the test used is the shrunken quasi-likelihood + # dispersion `ql_disp_shrunken`, whose inverse is reported as the + # dispersion theta. Zero dispersion is the Poisson boundary, where the F + # ratio would report p = 0 for the wrong reason, so inference is withheld; + # an unusable dispersion fails only its own taxon. + da_ql <- da_gp[["overdispersion_shrinkage_list"]][["ql_disp_shrunken"]][da_j] + if (!is.finite(da_ql) || da_ql < 0) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the shrunken dispersion for this taxon is not usable" + next + } + if (da_ql == 0) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "zero dispersion: Poisson boundary; quasi-likelihood inference withheld" + 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_ql) + } 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 +439,70 @@ 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 `dispersion_theta` the inverse of the shrunken +quasi-likelihood dispersion the test used for the taxon. 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 shrunken dispersion +marks the Poisson boundary, where the F ratio would report p = 0 for the wrong +reason, so inference is withheld. 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_glp', 'da_glp_failed', 'da_gp', " * + "'da_j', 'da_y', 'da_absent', 'da_ql', '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 +743,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 @@ -612,8 +771,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}, @@ -651,7 +813,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, :] @@ -662,14 +825,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/server/routes/analysis.jl b/src/server/routes/analysis.jl index 0dd506a..3480bc4 100644 --- a/src/server/routes/analysis.jl +++ b/src/server/routes/analysis.jl @@ -1195,7 +1195,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) @@ -1203,7 +1206,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 """ @@ -1282,6 +1286,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 2f927ba..e434269 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 184e70f..ef63be8 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 3b03489..46a6e79 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,6 +267,92 @@ 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. + # 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$overdispersion_shrinkage_list$ql_disp_shrunken) + }) + """) + 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 dispersions 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 + MetaManifold.RRuntime.with_r_lock() do + reval(raw""" + test_original_gp <- get("glm_gp", asNamespace("glmGamPoi")) + test_ql <- 0 + assignInNamespace("glm_gp", function(...) { + fit <- test_original_gp(...) + fit$overdispersion_shrinkage_list$ql_disp_shrunken[1] <- test_ql + fit + }, ns="glmGamPoi") + """) + try + off = log.(tss_factors(effect, effect_samples)) + for (ql, status) in ((0.0, "boundary"), (-1.0, "failed"), (NaN, "failed")) + RCall.globalEnv[:test_ql] = ql + rows = fit_nb(effect, effect_groups, off; levels=("A", "B"), + dispersion="glmGamPoi") + @test rows.status[1] == status + @test isnothing(rows.pvalue[1]) + @test occursin("dispersion", 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_ql) + """) + end end end end @@ -353,6 +461,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 @@ -720,4 +838,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 From e28e832fd52466f5192406025fe6b16d9fd99ae1 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Tue, 20 Oct 2026 14:15:47 +0100 Subject: [PATCH 3/3] Prepare upstream PR 62 review fixes: theta, label-safe glmGamPoi contrast (#192) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ## Change and reason Prepares the next push to upstream [JoshuaJewell/MetaManifold-WebUI#62](https://github.com/JoshuaJewell/MetaManifold-WebUI/pull/62), the glmGamPoi dispersion option for JoshuaJewell/MetaManifold-WebUI#37. The branch is upstream PR 62's reviewed head (`efe1df1`) plus three new commits, so the update fast-forwards. Under upstream CONTRIBUTING item 8, nothing is rebased once review has started. 1. **Report glm_gp's negative-binomial overdispersion as theta.** Your `wip/pr62-theta` commit, cherry-picked unchanged. It covers both of Joshua's code requests from Oct 19 (theta is `1 / overdispersions`; the Poisson boundary is read from it). 2. **Pass the glmGamPoi contrast to test_de as numeric weights.** `glmGamPoi::test_de` parses and evaluates a text contrast as R code. The contrast label is built from run, group and sub-group names joined by `/`, and those names may contain `-`. So the glmGamPoi path failed every taxon for comparisons such as `run-A` vs `run-B` or `grpA/run_1` vs `grpB/run_2`, and label text that reads as R code was executed. The commit adds a regression test that covers both dispersions. 3. **State what glmGamPoi's theta, standard error and boundary mean.** These are precision fixes, checked against the glmGamPoi 1.22.0 source: - `overdispersions` is the per-taxon estimate before shrinkage. The coefficients and the F-test use the cross-taxon trend, and only `lfc_se` reads `overdispersions`. - The boundary note no longer says the fit is Poisson. - The README now lists `dispersion` with the other study-scoped settings. **Do not merge this PR on GitHub.** Merging would put unsigned commits, or a merge commit, onto upstream PR 62. Re-sign on your machine and push instead: ```sh git fetch origin claude/elegant-rubin-z92rx4 upstream/glmgampoi-r git switch -c pr62-next origin/claude/elegant-rubin-z92rx4 git rebase --exec 'git commit --amend --no-edit -S' efe1df1 # re-signs only the 3 new commits git log --format='%h %G? %s' efe1df1..HEAD # expect G or E, never N git push origin pr62-next:upstream/glmgampoi-r # fast-forward: no --force needed ``` When Joshua asks for a rebase (upstream `main` is now `ce2c645`), the five-commit stack rebases cleanly onto it. That was checked in a scratch worktree. ## Validation Environment: Julia 1.12.5 (the CI pin, from the official image), R 4.3.3, MASS 7.3-60.0.1, glmGamPoi 1.14.3, bun 1.4.2. glmGamPoi 1.14.3 lacks `compute_lfc_se`, so 1.22.0's pure-R `test_de`, `test_de_q` and `predict.glmGamPoi` were loaded over it. Diffing the 1.14.0 and 1.22.0 sources shows no other change on these paths. Upstream CI runs the real pins (R 4.5, glmGamPoi 1.22.0). - **`CI=true julia --project=. test/runtests.jl test_differential.jl`:** - **head:** 253 passed, 0 failed, 5 skipped. The skips are `clr_welch` tests that need zCompositions. `CI=true` makes a missing glmGamPoi fail rather than skip, so every glmGamPoi test ran. - **commits 1, 2 and 3 on their own:** 202, 253 and 253 passed, 0 failed. - **the new label test on the pre-fix code:** 5 failures, as expected. The four label sets fail every taxon, and the R-code label ran. - **Label experiments** (the `fit_nb` R block, run directly): - **Before:** `run-B`, `grpB/run_2` and `2024-02` failed every taxon, and a label of `[file.create(...)]` created the file. - **After:** every label set gives rows identical to `A`/`B`, the file is not created, and `glm_nb` is unchanged. - **Null calibration** (share of p < 0.05, binomial SE about 0.003): | samples per group | glm_nb (Wald) | glmGamPoi (QL F) | |---|---|---| | 4 | 0.132 | 0.046 | | 12 | 0.069 | 0.048 | - **Frontend:** `bun test src` 39 pass; typecheck clean; build succeeds. Not run: the integration and server suites, the DADA2 taxonomy and phylogeny paths, and the `clr_welch` tests (zCompositions unavailable; they skip loudly). ## Scope No lockfile change. The upstream PR description still has to be replaced; Joshua's third request asks for that. A draft follows, together with the summary comment for upstream issue 37.
Draft replacement description for upstream PR 62 ~~~~markdown ## Change and reason Adds an optional `analysis.differential.dispersion` setting for the `nb_glm` method (issue #37). `glm_nb` stays the default and is unchanged: `MASS::glm.nb` fits each taxon with its own maximum-likelihood theta, and the group coefficient gets a Wald test. `glmGamPoi` fits one `glmGamPoi::glm_gp` across the taxa that pass `min_prevalence`, with `overdispersion_shrinkage = TRUE`. It reports `glmGamPoi::test_de`'s quasi-likelihood F-test of the group contrast (Lund et al., 2012), which is glmGamPoi's own published inference. The trend-plus-Wald design in the first version of this PR has been replaced. In that path: - **estimate and se**: `lfc` and `lfc_se`, converted from log2 to the natural-log scale. `lfc_se` comes from `predict()` and does not feed the p-value. - **statistic and pvalue**: the F statistic and its p-value. The F-test refits the reduced model at the cross-taxon dispersion trend and divides by the shrunken quasi-likelihood dispersion. - **dispersion_theta**: `1 / overdispersions`. This is glm_gp's negative-binomial overdispersion, the taxon's own estimate before shrinkage, and the same parameter `MASS::glm.nb`'s theta inverts (requested in review). - **Poisson boundary**: based on that overdispersion (requested in review). A zero estimate marks the row `boundary` and keeps its statistics, since the test does not read that estimate. A negative or non-finite estimate fails only its own taxon. - **contrast**: passed to `test_de` as numeric weights on the design's columns. `test_de` parses and evaluates a text contrast as R code. Before this change, a contrast label containing `-` or `/` (e.g. `run-B`, or `grpB/run_2`, which the app builds when conditions differ by group and run) made every taxon fail. A label that reads as R code was executed. - **failures**: a missing glmGamPoi raises `GlmGamPoiUnavailable` (HTTP 503). A `glm_gp` or `test_de` error fails every taxon, with no fallback to `glm_nb`. `glmGamPoi` is refused for `clr_welch`. Config and UI: the config editor and `PATCH /api/v1/config` accept the new enum (`glm_nb`, `glmGamPoi`). The constructor refuses an invalid method/dispersion pair when the analysis runs (HTTP 400 `invalid_config`). There is no other config rewrite; the first version's description wrongly claimed one. The README lists the option with the other study-scoped settings. Commits, oldest first (the last three are new since your review): 1. `02ee915` Add glmGamPoi to the pinned R dependencies 2. `efe1df1` Add glmGamPoi dispersion option to nb_glm differential abundance 3. Report glm_gp's negative-binomial overdispersion as theta (review: theta and boundary) 4. Pass the glmGamPoi contrast to test_de as numeric weights (labels with `-` or `/`; label text evaluated as R) 5. State what glmGamPoi's theta, standard error and boundary mean (comments, note text, README, result type) Closes #37. ## Validation Run in a cloud container with Julia 1.12.5 (the CI pin), R 4.3.3, MASS 7.3-60.0.1, glmGamPoi 1.14.3, and bun 1.4.2. glmGamPoi 1.14.3 lacks `compute_lfc_se`, so 1.22.0's `test_de`, `test_de_q` and `predict.glmGamPoi` (pure R) were loaded over it. Diffing the 1.14.0 and 1.22.0 sources shows no other change on these code paths: `glm_gp_impl`, `parse_contrast`, the overdispersion estimate and the quasi-likelihood shrinkage are the same. CI on this PR runs the pinned R 4.5 with glmGamPoi 1.22.0. - `CI=true julia --project=. test/runtests.jl test_differential.jl`: - **branch head:** 253 passed, 0 failed, 5 skipped. The skips are the `clr_welch` tests, which need zCompositions. Under `CI=true` a missing glmGamPoi fails rather than skips, so the glmGamPoi tests all ran. - **each new commit on its own:** 202, 253 and 253 passed, 0 failed. - **the new label test against the pre-fix code:** it fails as intended. Each of the four label sets fails every taxon, and the R-code label's code runs. - The `fit_nb` R block, run directly in R on the `effect` fixture: - **labels**: `run-A`/`run-B`, `grpA/run_1`/`grpB/run_2` and `2024-01`/`2024-02` failed every taxon before the fix. With it, every label set, including one containing R code, gives rows identical to `A`/`B`. - **injection**: a label of `[file.create(...)]` created the file before the fix and does not after it. - **glm_nb**: rows are unchanged. - **field check**: setting `overdispersions[1] <- 0` changes `lfc_se[1]` only; `pval`, `f_statistic` and `lfc` stay the same. - Null calibration through the same R block (no taxon differs; TSS offsets; dispersion falling with the mean). Share of p < 0.05: | samples per group | taxa tested | glm_nb (Wald) | glmGamPoi (QL F) | |---|---|---|---| | 4 | ~5,900 | 0.132 | 0.046 | | 12 | 4,000 | 0.069 | 0.048 | The binomial SE is about 0.003. - Frontend (`frontend/`): `bun test src` 39 pass; `bun run typecheck` clean; `bun run build` succeeds. Not run here: the integration suite (`--integration`, `--server`), the DADA2 taxonomy and phylogeny paths, and the `clr_welch` tests (zCompositions was not available in the container; they skip loudly). ## Scope Left out: a native Julia gamma-Poisson fit (#37 settled on R), zero-inflated and occupancy models, any change to the `glm_nb` default or its results, and any UI beyond the one enum. Lockfile: `renv.lock` adds glmGamPoi 1.22.0 and the eleven packages it needs that were not already pinned. Existing pins are unchanged. The option cannot run without them, and CI and `install.jl` check that glmGamPoi loads. Noted for a separate issue, not this PR: in the calibration above, the default `glm_nb` Wald test rejects too often at small n. ## Contribution checks - [x] I followed [CONTRIBUTING.md](https://github.com/JoshuaJewell/MetaManifold-WebUI/blob/main/CONTRIBUTING.md). - [ ] This branch starts at current upstream main and contains no merges from fork main. It starts at `029fa61` and contains no merges. Since review has started, it gains commits rather than being rebased (item 8). It rebases cleanly onto `ce2c645` whenever you ask. - [ ] Every commit is signed and shows as Verified. *(tick after re-signing)* - [ ] CI has run and passes. *(tick after the push)* - [x] Existing licensing, citation and acknowledgement text is unchanged. 🤖 Generated with [Claude Code](https://claude.com/claude-code) https://claude.ai/code/session_01HtyPfZ1MWBB94RawgzQgz7 ~~~~
Draft summary comment for upstream issue 37 ~~~~markdown **Summary for the record** **Question (Oct 15):** take glmGamPoi's dispersion shrinkage for `nb_glm` (a) in R through RCall, or (b) as a native Julia gamma-Poisson fit with R kept as the oracle? **Answer (@JoshuaJewell, Oct 17):** (a). "glmGamPoi stays R, since there is no reputable Julia equivalent." **What follows from that, in #62:** - **Option:** `analysis.differential.dispersion: glmGamPoi` is optional. The default stays `glm_nb` (`MASS::glm.nb`), and its results are unchanged. - **Test:** inference is glmGamPoi's own quasi-likelihood F-test (`glm_gp` + `test_de`), not the trend-plus-Wald plan from Oct 15. - **Review fixes:** theta is `1 / overdispersions` (glm_gp's NB overdispersion), and the Poisson boundary is read from it, as requested. The contrast is now passed as numeric weights. Before that, any label containing `-` or `/` (e.g. `grpB/run_2`) made every taxon fail, and label text was evaluated as R code. - **Dependencies:** `renv.lock` adds glmGamPoi 1.22.0 and the eleven packages it needs; existing pins are unchanged. - **Evidence:** in a null simulation, the share of p < 0.05 was: | samples per group | `glm_nb` (Wald) | glmGamPoi (QL F) | |---|---|---| | 4 | 13.2% | 4.6% | | 12 | 6.9% | 4.8% | The `glm_nb` excess at small n is outside this issue's scope and is worth its own issue. **Good to go.** The R/RCall route is settled, and #62 is ready for re-review once the three new commits are signed and pushed. ~~~~
--- _Generated by [Claude Code](https://claude.ai/code/session_01HtyPfZ1MWBB94RawgzQgz7)_ --------- Co-authored-by: Claude --- README.md | 12 ++++- config/defaults/pipeline.yml | 7 +-- frontend/src/api/types.ts | 2 +- src/analysis/differential.jl | 80 ++++++++++++++++++++-------------- test/unit/test_differential.jl | 65 ++++++++++++++++++++++----- 5 files changed, 117 insertions(+), 49 deletions(-) diff --git a/README.md b/README.md index ecffd7d..d77f7c8 100644 --- a/README.md +++ b/README.md @@ -249,7 +249,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. @@ -259,7 +259,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 2181a23..debcff0 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -363,9 +363,10 @@ analysis: # 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. A zero shrunken - # dispersion withholds inference for its taxon. Requires the pinned - # glmGamPoi R package. Refused with clr_welch. + # 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 diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e6ee284..fff182b 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -439,7 +439,7 @@ export interface DifferentialResult { min_prevalence: number /** nb_glm only. */ offset?: 'tss' | 'rle' - /** nb_glm only: where each taxon's theta comes from. */ + /** 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/src/analysis/differential.jl b/src/analysis/differential.jl index 5e46c3f..f87d7bc 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -291,27 +291,32 @@ 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 each taxon's -# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa, with each taxon's -# overdispersion shrunk toward the cross-taxon trend, supplies glmGamPoi's own -# published inference -- glmGamPoi::test_de's quasi-likelihood F-test of the -# group coefficient. Its genes x samples convention, so -# the counts are transposed, and size_factors takes the linear size factor -# (the per-taxon fits below use its log as their offset). A fit- or test-level -# error fails every taxon. No fallback to another dispersion estimator or -# significance test is permitted. +# "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 = model.matrix(~ da_g), + 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_term, max_lfc = Inf, compute_lfc_se = TRUE), + 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)) @@ -349,20 +354,20 @@ for (da_j in seq_len(da_n)) { # 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 - # per-taxon scale the test used is the shrunken quasi-likelihood - # dispersion `ql_disp_shrunken`, whose inverse is reported as the - # dispersion theta. Zero dispersion is the Poisson boundary, where the F - # ratio would report p = 0 for the wrong reason, so inference is withheld; - # an unusable dispersion fails only its own taxon. - da_ql <- da_gp[["overdispersion_shrinkage_list"]][["ql_disp_shrunken"]][da_j] - if (!is.finite(da_ql) || da_ql < 0) { + # 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 shrunken dispersion for this taxon is not usable" - next - } - if (da_ql == 0) { - da_result[da_j, "status"] <- "boundary" - da_result[da_j, "note"] <- "zero dispersion: Poisson boundary; quasi-likelihood inference withheld" + 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) @@ -376,7 +381,11 @@ for (da_j in seq_len(da_n)) { next } da_result[da_j, c("estimate", "se", "statistic", "pvalue", "theta")] <- - c(da_est, da_se, da_stat, da_p, 1 / da_ql) + 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 @@ -453,17 +462,21 @@ 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 `dispersion_theta` the inverse of the shrunken -quasi-likelihood dispersion the test used for the taxon. Names are stripped -and matched without regard to case. +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. 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. 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 shrunken dispersion -marks the Poisson boundary, where the F ratio would report p = 0 for the wrong -reason, so inference is withheld. Constant counts, complete separation and +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"`. @@ -500,8 +513,9 @@ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:Abst finally RCall.reval("rm(list = intersect(ls(), c('da_counts', 'da_group', 'da_levels', " * "'da_offset', 'da_theta_upper', 'da_theta_lower', 'da_use_glp', 'da_g', " * - "'da_term', 'da_n', 'da_result', 'da_glp', 'da_glp_failed', 'da_gp', " * - "'da_j', 'da_y', 'da_absent', 'da_ql', 'da_est', 'da_se', 'da_stat', " * + "'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 diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 46a6e79..2c686a2 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -295,7 +295,9 @@ end # 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. - # One matrix column per quantity, one row per taxon. + # 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")) @@ -304,7 +306,7 @@ end 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$overdispersion_shrinkage_list$ql_disp_shrunken) + 1 / gp$overdispersions) }) """) for (j, t) in enumerate(effect_taxa) @@ -320,43 +322,86 @@ end end end - @testset "glmGamPoi unusable dispersions are isolated per taxon" begin + @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_ql <- 0 + test_od <- 0 assignInNamespace("glm_gp", function(...) { fit <- test_original_gp(...) - fit$overdispersion_shrinkage_list$ql_disp_shrunken[1] <- test_ql + fit$overdispersions[1] <- test_od fit }, ns="glmGamPoi") """) try off = log.(tss_factors(effect, effect_samples)) - for (ql, status) in ((0.0, "boundary"), (-1.0, "failed"), (NaN, "failed")) - RCall.globalEnv[:test_ql] = ql + # 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]) - @test occursin("dispersion", rows.note[1]) + @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_ql) + 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 + @testset "Degenerate and filtered taxa are reported, not dropped" begin if !HAVE_MASS @test_skip false