diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 86fd896e..573768be 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -156,7 +156,7 @@ jobs: sudo chmod 777 "$RENV_PATHS_CACHE" sudo env R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'install.packages(c("renv", "BiocManager"), repos="https://cloud.r-project.org", lib=.Library)' sudo env RENV_PATHS_CACHE="$RENV_PATHS_CACHE" R_CLI_DYNAMIC=false TERM=dumb Rscript --no-init-file -e 'renv::restore(project=".", library=.Library, prompt=FALSE)' - Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","MASS","zCompositions")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' + Rscript --no-init-file -e 'for (pkg in c("dada2","Biostrings","ShortRead","vegan","MASS","zCompositions","glmGamPoi")) if (!require(pkg,character.only=TRUE,quietly=TRUE)) stop(pkg, " failed to install")' - name: Check DADA2 CSV assembly run: Rscript --vanilla test/unit/test_dada2_tables.R @@ -250,7 +250,7 @@ jobs: julia --project=. -e ' using RCall R"cat(.libPaths(), sep=\"\n\")" - for pkg in ["dada2", "vegan", "MASS", "zCompositions"] + for pkg in ["dada2", "vegan", "MASS", "zCompositions", "glmGamPoi"] R"if (!require($pkg, character.only=TRUE, quietly=TRUE)) stop($pkg, \" not found\")" end ' diff --git a/R/_renv_dependencies.R b/R/_renv_dependencies.R index 1e91aca9..589612fe 100644 --- a/R/_renv_dependencies.R +++ b/R/_renv_dependencies.R @@ -3,6 +3,7 @@ if (FALSE) { library(dada2) library(vegan) + library(glmGamPoi) library(MASS) library(zCompositions) } diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 8181ae6d..db28ddcd 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: where each taxon's negative-binomial theta comes from. ml + # (the default) takes the maximum-likelihood estimate MASS::glm.nb fits per + # taxon. glmGamPoi uses glmGamPoi::glm_gp + # dispersion_trend across the taxa passing min_prevalence, then fits each + # taxon through MASS::negative.binomial(theta = 1/dispersion). This is a + # plug-in NB Wald test with scale 1, not glmGamPoi quasi-likelihood testing; + # trend-estimation uncertainty is not propagated. A zero trend withholds + # NB inference. Requires the pinned glmGamPoi R package. + dispersion: ml # 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 87457dad..4ea20add 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?: 'ml' | 'glmGamPoi' /** clr_welch only. */ replacement_delta?: number } diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts new file mode 100644 index 00000000..37fa1f24 --- /dev/null +++ b/frontend/src/components/PipelineStages.test.ts @@ -0,0 +1,9 @@ +// SPDX-License-Identifier: AGPL-3.0-only +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: ['ml', 'glmGamPoi'], + }) +}) diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx index 85129267..39fcecdf 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: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits. Trend-estimation uncertainty is not propagated; this is not glmGamPoi quasi-likelihood testing.', '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: ['ml', 'glmGamPoi'] }, 'analysis.differential.min_prevalence': { kind: 'float', min: 0, max: 1, step: 0.05 }, 'analysis.differential.replacement_delta': { kind: 'float', min: 0.01, max: 0.99, step: 0.05 }, 'analysis.alpha.show_points': { kind: 'boolean' }, diff --git a/install.jl b/install.jl index 964e00d9..3be35d24 100644 --- a/install.jl +++ b/install.jl @@ -879,8 +879,8 @@ function setup_r_packages(; rebuild::Bool=false)::Tuple{Symbol,String} # Verdict is whether the packages the pipeline loads are actually usable, not # renv's exit code: restore can fail on one leaf package and still leave a # working library, or "succeed" against a stale cache. - # MASS and zCompositions back the differential abundance methods. - pkgs = ["dada2", "vegan", "MASS", "zCompositions"] + # MASS, zCompositions and glmGamPoi back the differential abundance methods. + pkgs = ["dada2", "vegan", "MASS", "zCompositions", "glmGamPoi"] probe = """ pkgs <- c($(join(repr.(pkgs), ", "))) cat(paste(pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)], diff --git a/renv.lock b/renv.lock index cfd72d4b..9944f58c 100644 --- a/renv.lock +++ b/renv.lock @@ -3011,6 +3011,565 @@ "BugReports": "https://github.com/Japal/zCompositions/issues", "NeedsCompilation": "no", "Author": "Javier Palarea-Albaladejo [cre, aut] (ORCID: ), Josep Antoni Martin-Fernandez [aut] (ORCID: )" + }, + "DelayedMatrixStats": { + "Package": "DelayedMatrixStats", + "Version": "1.32.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Functions that Apply to Rows and Columns of 'DelayedMatrix' Objects", + "Date": "2025-01-09", + "Authors@R": "c(person(\"Peter\", \"Hickey\", role = c(\"aut\", \"cre\"), email = \"peter.hickey@gmail.com\", comment = c(ORCID = \"0000-0002-8153-6258\")), person(\"Hervé\", \"Pagès\", role = \"ctb\", email = \"hpages.on.github@gmail.com\"), person(\"Aaron\", \"Lun\", role = \"ctb\", email = \"infinite.monkeys.with.keyboards@gmail.com\"))", + "Description": "A port of the 'matrixStats' API for use with DelayedMatrix objects from the 'DelayedArray' package. High-performing functions operating on rows and columns of DelayedMatrix objects, e.g. col / rowMedians(), col / rowRanks(), and col / rowSds(). Functions optimized per data type and for subsetted calculations such that both memory usage and processing time is minimized.", + "License": "MIT + file LICENSE", + "Encoding": "UTF-8", + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.2", + "Depends": [ + "MatrixGenerics (>= 1.15.1)", + "DelayedArray (>= 0.31.7)" + ], + "Imports": [ + "methods", + "sparseMatrixStats (>= 1.13.2)", + "Matrix (>= 1.5-0)", + "S4Vectors (>= 0.17.5)", + "IRanges (>= 2.25.10)", + "SparseArray (>= 1.5.19)" + ], + "Suggests": [ + "testthat", + "knitr", + "rmarkdown", + "BiocStyle", + "microbenchmark", + "profmem", + "HDF5Array", + "matrixStats (>= 1.0.0)" + ], + "VignetteBuilder": "knitr", + "URL": "https://github.com/PeteHaitch/DelayedMatrixStats", + "BugReports": "https://github.com/PeteHaitch/DelayedMatrixStats/issues", + "biocViews": "Infrastructure, DataRepresentation, Software", + "git_url": "https://git.bioconductor.org/packages/DelayedMatrixStats", + "git_branch": "RELEASE_3_22", + "git_last_commit": "cbf7d75", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Peter Hickey [aut, cre] (ORCID: ), Hervé Pagès [ctb], Aaron Lun [ctb]", + "Maintainer": "Peter Hickey " + }, + "HDF5Array": { + "Package": "HDF5Array", + "Version": "1.38.0", + "Source": "Bioconductor", + "Title": "HDF5 datasets as array-like objects in R", + "Description": "The HDF5Array package is an HDF5 backend for DelayedArray objects. It implements the HDF5Array, H5SparseMatrix, H5ADMatrix, and TENxMatrix classes, 4 convenient and memory-efficient array-like containers for representing and manipulating either: (1) a conventional (a.k.a. dense) HDF5 dataset, (2) an HDF5 sparse matrix (stored in CSR/CSC/Yale format), (3) the central matrix of an h5ad file (or any matrix in the /layers group), or (4) a 10x Genomics sparse matrix. All these containers are DelayedArray extensions and thus support all operations (delayed or block-processed) supported by DelayedArray objects.", + "biocViews": "Infrastructure, DataRepresentation, DataImport, Sequencing, RNASeq, Coverage, Annotation, GenomeAnnotation, SingleCell, ImmunoOncology", + "URL": "https://bioconductor.org/packages/HDF5Array", + "BugReports": "https://github.com/Bioconductor/HDF5Array/issues", + "License": "Artistic-2.0", + "Encoding": "UTF-8", + "Authors@R": "person(\"Hervé\", \"Pagès\", role=c(\"aut\", \"cre\"), email=\"hpages.on.github@gmail.com\", comment=c(ORCID=\"0009-0002-8272-4522\"))", + "Depends": [ + "R (>= 3.4)", + "methods", + "SparseArray (>= 1.7.5)", + "DelayedArray (>= 0.33.5)", + "h5mread (>= 0.99.4)" + ], + "Imports": [ + "utils", + "stats", + "tools", + "Matrix", + "BiocGenerics (>= 0.51.2)", + "S4Vectors", + "IRanges", + "S4Arrays (>= 1.1.1)", + "rhdf5" + ], + "Suggests": [ + "BiocParallel", + "GenomicRanges", + "SummarizedExperiment (>= 1.15.1)", + "h5vcData", + "ExperimentHub", + "TENxBrainData", + "zellkonverter", + "GenomicFeatures", + "SingleCellExperiment", + "DelayedMatrixStats", + "genefilter", + "RSpectra", + "RUnit", + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "Collate": "utils.R h5utils.R HDF5ArraySeed-class.R HDF5Array-class.R ReshapedHDF5ArraySeed-class.R ReshapedHDF5Array-class.R dump-management.R writeHDF5Array.R saveHDF5SummarizedExperiment.R H5SparseMatrixSeed-class.R H5SparseMatrix-class.R H5ADMatrixSeed-class.R H5ADMatrix-class.R TENxMatrixSeed-class.R TENxMatrix-class.R writeTENxMatrix.R zzz.R", + "git_url": "https://git.bioconductor.org/packages/HDF5Array", + "git_branch": "RELEASE_3_22", + "git_last_commit": "9bca08f", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Hervé Pagès [aut, cre] (ORCID: )", + "Maintainer": "Hervé Pagès " + }, + "RcppArmadillo": { + "Package": "RcppArmadillo", + "Version": "15.6.2-1", + "Source": "Repository", + "Type": "Package", + "Title": "'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library", + "Date": "2026-10-16", + "Authors@R": "c(person(\"Dirk\", \"Eddelbuettel\", role = c(\"aut\", \"cre\"), email = \"edd@debian.org\", comment = c(ORCID = \"0000-0001-6419-907X\")), person(\"Romain\", \"Francois\", role = \"aut\", comment = c(ORCID = \"0000-0002-2444-4226\")), person(\"Doug\", \"Bates\", role = \"aut\", comment = c(ORCID = \"0000-0001-8316-9503\")), person(\"Binxiang\", \"Ni\", role = \"aut\"), person(\"Conrad\", \"Sanderson\", role = \"aut\", comment = c(ORCID = \"0000-0002-0049-4501\")))", + "Description": "'Armadillo' is a templated C++ linear algebra library aiming towards a good balance between speed and ease of use. It provides high-level syntax and functionality deliberately similar to Matlab. It is useful for algorithm development directly in C++, or quick conversion of research code into production environments. It provides efficient classes for vectors, matrices and cubes where dense and sparse matrices are supported. Integer, floating point and complex numbers are supported. A sophisticated expression evaluator (based on template meta-programming) automatically combines several operations to increase speed and efficiency. Dynamic evaluation automatically chooses optimal code paths based on detected matrix structures. Matrix decompositions are provided through integration with LAPACK, or one of its high performance drop-in replacements (such as 'MKL' or 'OpenBLAS'). It can automatically use 'OpenMP' multi-threading (parallelisation) to speed up computationally expensive operations. . The 'RcppArmadillo' package includes the header files from the 'Armadillo' library; users do not need to install 'Armadillo' itself in order to use 'RcppArmadillo'. Starting from release 15.0.0, the minimum compilation standard is C++14. . Since release 7.800.0, 'Armadillo' is licensed under Apache License 2; previous releases were under licensed as MPL 2.0 from version 3.800.0 onwards and LGPL-3 prior to that; 'RcppArmadillo' (the 'Rcpp' bindings/bridge to Armadillo) is licensed under the GNU GPL version 2 or later, as is the rest of 'Rcpp'.", + "License": "GPL (>= 2)", + "LazyLoad": "yes", + "Depends": [ + "R (>= 3.3.0)" + ], + "LinkingTo": [ + "Rcpp" + ], + "Imports": [ + "Rcpp (>= 1.1.1)", + "stats", + "utils", + "methods" + ], + "Suggests": [ + "tinytest", + "Matrix (>= 1.3.0)", + "pkgKitten", + "reticulate", + "slam" + ], + "URL": "https://github.com/RcppCore/RcppArmadillo, https://dirk.eddelbuettel.com/code/rcpp.armadillo.html", + "BugReports": "https://github.com/RcppCore/RcppArmadillo/issues", + "RoxygenNote": "6.0.1", + "Encoding": "UTF-8", + "VignetteBuilder": "Rcpp", + "NeedsCompilation": "yes", + "Author": "Dirk Eddelbuettel [aut, cre] (ORCID: ), Romain Francois [aut] (ORCID: ), Doug Bates [aut] (ORCID: ), Binxiang Ni [aut], Conrad Sanderson [aut] (ORCID: )", + "Maintainer": "Dirk Eddelbuettel ", + "Repository": "CRAN" + }, + "Rhdf5lib": { + "Package": "Rhdf5lib", + "Version": "1.32.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "hdf5 library as an R package", + "Authors@R": "c( person( \"Mike\", \"Smith\", role = \"ctb\", comment = c(ORCID = \"0000-0002-7800-3848\") ), person( \"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\")) , person( given = \"The HDF Group\", role = \"cph\" ) )", + "Description": "Provides C and C++ hdf5 libraries.", + "License": "Artistic-2.0", + "Copyright": "src/hdf5/COPYING", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Depends": [ + "R (>= 4.2.0)" + ], + "Suggests": [ + "BiocStyle", + "knitr", + "rmarkdown", + "tinytest", + "mockery" + ], + "URL": "https://github.com/Huber-group-EMBL/Rhdf5lib", + "BugReports": "https://github.com/Huber-group-EMBL/Rhdf5lib/issues", + "SystemRequirements": "GNU make", + "Encoding": "UTF-8", + "biocViews": "Infrastructure", + "RoxygenNote": "7.1.2", + "git_url": "https://git.bioconductor.org/packages/Rhdf5lib", + "git_branch": "RELEASE_3_22", + "git_last_commit": "f62ae28", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Mike Smith [ctb] (ORCID: ), Hugo Gruson [cre] (ORCID: ), The HDF Group [cph]", + "Maintainer": "Hugo Gruson " + }, + "SingleCellExperiment": { + "Package": "SingleCellExperiment", + "Version": "1.32.0", + "Source": "Bioconductor", + "Date": "2025-06-27", + "Title": "S4 Classes for Single Cell Data", + "Authors@R": "c( person(\"Aaron\", \"Lun\", role=c(\"aut\", \"cph\"), email=\"infinite.monkeys.with.keyboards@gmail.com\"), person(\"Davide\",\"Risso\", role=c(\"aut\",\"cre\", \"cph\"), email=\"risso.davide@gmail.com\"), person(\"Keegan\", \"Korthauer\", role=\"ctb\"), person(\"Kevin\", \"Rue-Albrecht\", role=\"ctb\"), person(\"Luke\", \"Zappia\", role = \"ctb\", comment = c(ORCID = \"0000-0001-7744-8565\", github = \"lazappi\")))", + "Depends": [ + "SummarizedExperiment" + ], + "Imports": [ + "methods", + "utils", + "stats", + "S4Vectors", + "BiocGenerics", + "GenomicRanges", + "DelayedArray" + ], + "Suggests": [ + "testthat", + "BiocStyle", + "knitr", + "rmarkdown", + "Matrix", + "scRNAseq (>= 2.9.1)", + "Rtsne" + ], + "biocViews": "ImmunoOncology, DataRepresentation, DataImport, Infrastructure, SingleCell", + "Description": "Defines a S4 class for storing data from single-cell experiments. This includes specialized methods to store and retrieve spike-in information, dimensionality reduction coordinates and size factors for each cell, along with the usual metadata for genes and libraries.", + "License": "GPL-3", + "VignetteBuilder": "knitr", + "RoxygenNote": "7.3.2", + "git_url": "https://git.bioconductor.org/packages/SingleCellExperiment", + "git_branch": "RELEASE_3_22", + "git_last_commit": "db7133e", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Aaron Lun [aut, cph], Davide Risso [aut, cre, cph], Keegan Korthauer [ctb], Kevin Rue-Albrecht [ctb], Luke Zappia [ctb] (ORCID: , github: lazappi)", + "Maintainer": "Davide Risso " + }, + "assorthead": { + "Package": "assorthead", + "Version": "1.4.0", + "Source": "Bioconductor", + "Date": "2025-10-27", + "Title": "Assorted Header-Only C++ Libraries", + "Authors@R": "person(\"Aaron\", \"Lun\", role=c(\"cre\", \"aut\"), email=\"infinite.monkeys.with.keyboards@gmail.com\")", + "Description": "Vendors an assortment of useful header-only C++ libraries. Bioconductor packages can use these libraries in their own C++ code by LinkingTo this package without introducing any additional dependencies. The use of a central repository avoids duplicate vendoring of libraries across multiple R packages, and enables better coordination of version updates across cohorts of interdependent C++ libraries.", + "License": "MIT + file LICENSE", + "Suggests": [ + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "URL": "https://github.com/LTLA/assorthead", + "BugReports": "https://github.com/LTLA/assorthead/issues", + "Encoding": "UTF-8", + "biocViews": "SingleCell, QualityControl, Normalization, DataRepresentation, DataImport, DifferentialExpression, Alignment", + "git_url": "https://git.bioconductor.org/packages/assorthead", + "git_branch": "RELEASE_3_22", + "git_last_commit": "9255f2f", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "no", + "Author": "Aaron Lun [cre, aut]", + "Maintainer": "Aaron Lun " + }, + "beachmat": { + "Package": "beachmat", + "Version": "2.26.0", + "Source": "Bioconductor", + "Date": "2025-08-19", + "Title": "Compiling Bioconductor to Handle Each Matrix Type", + "Authors@R": "c(person(\"Aaron\", \"Lun\", role = c(\"aut\", \"cre\"), email = \"infinite.monkeys.with.keyboards@gmail.com\"), person(\"Hervé\", \"Pagès\", role=\"aut\"), person(\"Mike\", \"Smith\", role=\"aut\"))", + "Imports": [ + "methods", + "DelayedArray (>= 0.27.2)", + "SparseArray", + "BiocGenerics", + "Matrix", + "Rcpp" + ], + "Suggests": [ + "testthat", + "BiocStyle", + "knitr", + "rmarkdown", + "rcmdcheck", + "BiocParallel", + "HDF5Array", + "beachmat.hdf5" + ], + "LinkingTo": [ + "Rcpp", + "assorthead (>= 1.3.3)" + ], + "biocViews": "DataRepresentation, DataImport, Infrastructure", + "Description": "Provides a consistent C++ class interface for reading from a variety of commonly used matrix types. Ordinary matrices and several sparse/dense Matrix classes are directly supported, along with a subset of the delayed operations implemented in the DelayedArray package. All other matrix-like objects are supported by calling back into R.", + "License": "GPL-3", + "NeedsCompilation": "yes", + "VignetteBuilder": "knitr", + "SystemRequirements": "C++17", + "URL": "https://github.com/tatami-inc/beachmat", + "BugReports": "https://github.com/tatami-inc/beachmat/issues", + "RoxygenNote": "7.3.2", + "Encoding": "UTF-8", + "git_url": "https://git.bioconductor.org/packages/beachmat", + "git_branch": "RELEASE_3_22", + "git_last_commit": "016c55e", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "Author": "Aaron Lun [aut, cre], Hervé Pagès [aut], Mike Smith [aut]", + "Maintainer": "Aaron Lun " + }, + "glmGamPoi": { + "Package": "glmGamPoi", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Fit a Gamma-Poisson Generalized Linear Model", + "Authors@R": "c(person(\"Constantin\", \"Ahlmann-Eltze\", email = \"artjom31415@googlemail.com\", role = c(\"aut\", \"cre\"), comment = c(ORCID = \"0000-0002-3762-068X\")), person(\"Nathan\", \"Lubock\", email=\"lubock.nathan@gmail.com\", role = \"ctb\", comment = c(ORCID = \"0000-0001-8064-2465\")), person(\"Michael\", \"Love\", email=\"michaelisaiahlove@gmail.com\", role = \"ctb\"))", + "Description": "Fit linear models to overdispersed count data. The package can estimate the overdispersion and fit repeated models for matrix input. It is designed to handle large input datasets as they typically occur in single cell RNA-seq experiments.", + "License": "GPL-3", + "Encoding": "UTF-8", + "SystemRequirements": "C++17", + "Suggests": [ + "testthat (>= 2.1.0)", + "zoo", + "DESeq2", + "edgeR", + "limma", + "MASS", + "statmod", + "ggplot2", + "bench", + "BiocParallel", + "knitr", + "rmarkdown", + "BiocStyle", + "TENxPBMCData", + "muscData", + "scran", + "dplyr" + ], + "LinkingTo": [ + "Rcpp", + "RcppArmadillo", + "beachmat", + "assorthead" + ], + "Imports": [ + "Rcpp", + "beachmat", + "DelayedMatrixStats", + "matrixStats", + "MatrixGenerics", + "SparseArray (>= 1.5.21)", + "S4Vectors", + "DelayedArray", + "HDF5Array", + "Matrix", + "SummarizedExperiment", + "SingleCellExperiment", + "BiocGenerics", + "methods", + "stats", + "utils", + "splines", + "rlang", + "vctrs" + ], + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.2", + "URL": "https://github.com/const-ae/glmGamPoi", + "BugReports": "https://github.com/const-ae/glmGamPoi/issues", + "biocViews": "Regression, RNASeq, Software, SingleCell", + "VignetteBuilder": "knitr", + "git_url": "https://git.bioconductor.org/packages/glmGamPoi", + "git_branch": "RELEASE_3_22", + "git_last_commit": "908cde9", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Constantin Ahlmann-Eltze [aut, cre] (ORCID: ), Nathan Lubock [ctb] (ORCID: ), Michael Love [ctb]", + "Maintainer": "Constantin Ahlmann-Eltze ", + "Depends": [ + "R (>= 4.1.0)" + ] + }, + "h5mread": { + "Package": "h5mread", + "Version": "1.2.1", + "Source": "Bioconductor", + "Title": "A fast HDF5 reader", + "Description": "The main function in the h5mread package is h5mread(), which allows reading arbitrary data from an HDF5 dataset into R, similarly to what the h5read() function from the rhdf5 package does. In the case of h5mread(), the implementation has been optimized to make it as fast and memory-efficient as possible.", + "biocViews": "Infrastructure, DataRepresentation, DataImport", + "URL": "https://bioconductor.org/packages/h5mread", + "BugReports": "https://github.com/Bioconductor/h5mread/issues", + "License": "Artistic-2.0", + "Encoding": "UTF-8", + "Authors@R": "person(\"Hervé\", \"Pagès\", role=c(\"aut\", \"cre\"), email=\"hpages.on.github@gmail.com\", comment=c(ORCID=\"0009-0002-8272-4522\"))", + "Depends": [ + "R (>= 4.5)", + "methods", + "rhdf5", + "BiocGenerics", + "SparseArray" + ], + "Imports": [ + "stats", + "tools", + "rhdf5filters", + "S4Vectors", + "IRanges", + "S4Arrays" + ], + "LinkingTo": [ + "Rhdf5lib", + "S4Vectors" + ], + "SystemRequirements": "GNU make", + "Suggests": [ + "BiocParallel", + "ExperimentHub", + "TENxBrainData", + "HDF5Array", + "testthat", + "knitr", + "rmarkdown", + "BiocStyle" + ], + "VignetteBuilder": "knitr", + "Collate": "utils.R h5dim.R H5File-class.R h5ls.R H5DSetDescriptor-class.R uaselection.R h5mread.R h5summarize.R h5mread_from_reshaped.R h5dimscales.R h5writeDimnames.R zzz.R", + "git_url": "https://git.bioconductor.org/packages/h5mread", + "git_branch": "RELEASE_3_22", + "git_last_commit": "b377076", + "git_last_commit_date": "2025-11-21", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Hervé Pagès [aut, cre] (ORCID: )", + "Maintainer": "Hervé Pagès " + }, + "rhdf5": { + "Package": "rhdf5", + "Version": "2.54.1", + "Source": "Bioconductor", + "Type": "Package", + "Title": "R Interface to HDF5", + "Authors@R": "c( person(\"Bernd\", \"Fischer\", role = \"aut\"), person( \"Mike\", \"Smith\", role= \"aut\", comment = c(ORCID = \"0000-0002-7800-3848\", \"Maintainer from 2017 to 2025\") ), person(\"Gregoire\", \"Pau\", role=\"aut\"), person(\"Martin\", \"Morgan\", role = \"ctb\"), person(\"Daniel\", \"van Twisk\", role = \"ctb\"), person( \"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\") ) )", + "Description": "This package provides an interface between HDF5 and R. HDF5's main features are the ability to store and access very large and/or complex datasets and a wide variety of metadata on mass storage (disk) through a completely portable file format. The rhdf5 package is thus suited for the exchange of large and/or complex datasets between R and other software package, and for letting R applications work on datasets that are larger than the available RAM.", + "License": "Artistic-2.0", + "URL": "https://github.com/Huber-group-EMBL/rhdf5", + "BugReports": "https://github.com/Huber-group-EMBL/rhdf5/issues", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Imports": [ + "Rhdf5lib (>= 1.13.4)", + "rhdf5filters (>= 1.15.5)" + ], + "Depends": [ + "R (>= 4.0.0)", + "methods" + ], + "Suggests": [ + "bit64", + "BiocStyle", + "knitr", + "rmarkdown", + "testthat", + "bench", + "dplyr", + "ggplot2", + "mockery", + "BiocParallel" + ], + "LinkingTo": [ + "Rhdf5lib" + ], + "SystemRequirements": "GNU make", + "biocViews": "Infrastructure, DataImport", + "Encoding": "UTF-8", + "Roxygen": "list(markdown = TRUE)", + "RoxygenNote": "7.3.3", + "git_url": "https://git.bioconductor.org/packages/rhdf5", + "git_branch": "RELEASE_3_22", + "git_last_commit": "7f691e4", + "git_last_commit_date": "2025-12-02", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Bernd Fischer [aut], Mike Smith [aut] (ORCID: , Maintainer from 2017 to 2025), Gregoire Pau [aut], Martin Morgan [ctb], Daniel van Twisk [ctb], Hugo Gruson [cre] (ORCID: )", + "Maintainer": "Hugo Gruson " + }, + "rhdf5filters": { + "Package": "rhdf5filters", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "HDF5 Compression Filters", + "Authors@R": "c( person(\"Mike\", \"Smith\", role = c(\"aut\", \"ccp\"), comment = c(ORCID = \"0000-0002-7800-3848\") ), person(\"Hugo\", \"Gruson\", role = \"cre\", email = \"hugo.gruson@embl.de\", comment = c(ORCID = \"0000-0002-4094-1476\") ) )", + "Description": "Provides a collection of additional compression filters for HDF5 datasets. The package is intended to provide seemless integration with rhdf5, however the compiled filters can also be used with external applications.", + "License": "BSD_2_clause + file LICENSE", + "LazyLoad": "true", + "VignetteBuilder": "knitr", + "Suggests": [ + "BiocStyle", + "knitr", + "rmarkdown", + "tinytest", + "rhdf5 (>= 2.47.7)" + ], + "SystemRequirements": "GNU make", + "URL": "https://github.com/Huber-group-EMBL/rhdf5filters", + "BugReports": "https://github.com/Huber-group-EMBL/rhdf5filters/issues", + "LinkingTo": [ + "Rhdf5lib" + ], + "RoxygenNote": "7.3.2", + "Encoding": "UTF-8", + "biocViews": "Infrastructure, DataImport", + "git_url": "https://git.bioconductor.org/packages/rhdf5filters", + "git_branch": "RELEASE_3_22", + "git_last_commit": "3465c24", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Mike Smith [aut, ccp] (ORCID: ), Hugo Gruson [cre] (ORCID: )", + "Maintainer": "Hugo Gruson " + }, + "sparseMatrixStats": { + "Package": "sparseMatrixStats", + "Version": "1.22.0", + "Source": "Bioconductor", + "Type": "Package", + "Title": "Summary Statistics for Rows and Columns of Sparse Matrices", + "Authors@R": "person(\"Constantin\", \"Ahlmann-Eltze\", email = \"artjom31415@googlemail.com\", role = c(\"aut\", \"cre\"), comment = c(ORCID = \"0000-0002-3762-068X\"))", + "Description": "High performance functions for row and column operations on sparse matrices. For example: col / rowMeans2, col / rowMedians, col / rowVars etc. Currently, the optimizations are limited to data in the column sparse format. This package is inspired by the matrixStats package by Henrik Bengtsson.", + "License": "MIT + file LICENSE", + "Encoding": "UTF-8", + "RoxygenNote": "7.3.1", + "LinkingTo": [ + "Rcpp" + ], + "Imports": [ + "Rcpp", + "Matrix", + "matrixStats (>= 0.60.0)", + "methods" + ], + "Depends": [ + "MatrixGenerics (>= 1.5.3)" + ], + "Suggests": [ + "testthat (>= 2.1.0)", + "knitr", + "bench", + "rmarkdown", + "BiocStyle" + ], + "SystemRequirements": "C++11", + "VignetteBuilder": "knitr", + "URL": "https://github.com/const-ae/sparseMatrixStats", + "BugReports": "https://github.com/const-ae/sparseMatrixStats/issues", + "biocViews": "Infrastructure, Software, DataRepresentation", + "git_url": "https://git.bioconductor.org/packages/sparseMatrixStats", + "git_branch": "RELEASE_3_22", + "git_last_commit": "3acf348", + "git_last_commit_date": "2025-10-29", + "Repository": "Bioconductor 3.22", + "NeedsCompilation": "yes", + "Author": "Constantin Ahlmann-Eltze [aut, cre] (ORCID: )", + "Maintainer": "Constantin Ahlmann-Eltze " } } } diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index d2c78de0..9cbd392f 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -22,13 +22,18 @@ 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, 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") +# The per-taxon dispersion `nb_glm` uses: maximum likelihood (MASS::glm.nb's +# own theta) or glmGamPoi's shrinkage estimator, which shares information +# across taxa. `"ml"` is the default and keeps #24's results unchanged. +const DISPERSIONS = ("ml", "glmGamPoi") + # Offset (size factor) methods accepted under `analysis.differential.offset`. const OFFSET_METHODS = ("tss", "rle") @@ -43,7 +48,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="ml") Settings for one differential abundance analysis, read from `analysis.differential` in the pipeline config. `method` is `"nb_glm"` or @@ -54,26 +59,50 @@ 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 the per-taxon dispersion `nb_glm` uses (`"ml"` or +`"glmGamPoi"`): `"ml"` takes the maximum-likelihood theta from `MASS::glm.nb`, +`"glmGamPoi"` uses `overdispersion_shrinkage_list.dispersion_trend` from +`glmGamPoi::glm_gp` for +the taxa passing the prevalence filter, with `theta = 1 / dispersion` fixed +in each negative-binomial GLM. Both options use a Wald test of the group +coefficient with fixed scale 1; `"glmGamPoi"` is refused for `clr_welch`. +The trend is plugged into an NB Wald model; its estimation uncertainty is not +propagated. This option does not implement glmGamPoi's quasi-likelihood test. +Method and offset names are stripped and lowercased; dispersion names are +stripped and matched without regard to case, then stored as `"ml"` 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="ml") 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 == "ml" || 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 +144,25 @@ 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 +`analysis.differential.dispersion: glmGamPoi` asks for, cannot be loaded. +""" +struct GlmGamPoiUnavailable <: Exception end + +""" + Base.showerror(io::IO, ::GlmGamPoiUnavailable) + +Write the missing-package error and guidance to restore the pinned R library or select +`dispersion: ml` to `io`, without a trailing newline. +""" +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: ml") + """ _geomean(x) -> Float64 @@ -243,7 +291,37 @@ 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's shrunken dispersion, fitted across all taxa in one call -- its +# genes x samples convention, so the counts are transposed, and size_factors +# is the linear size factor (the per-taxon fits below take its log as the +# offset). A global fit error or missing trend fails the batch. A zero trend +# is a per-taxon Poisson boundary; malformed values fail only their own taxa. +# No fallback to a different dispersion estimator is permitted. +da_phi <- rep(NA_real_, da_n) +da_gp_failed <- "" +if (da_use_gp) { + 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_gp_failed <- paste("glm_gp stopped with an error:", conditionMessage(da_gp)) + } else { + # glm_gp's overdispersions retain the unshrunk MLEs. Its NB refit uses + # dispersion_trend; ql_disp_shrunken is a distinct quasi-likelihood scale. + da_phi <- as.numeric(da_gp[["overdispersion_shrinkage_list"]][["dispersion_trend"]]) + if (length(da_phi) != da_n) { + da_gp_failed <- "glm_gp returned no dispersion trend for each taxon" + da_phi <- rep(NA_real_, da_n) + } + } +} for (da_j in seq_len(da_n)) { + if (nzchar(da_gp_failed)) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- da_gp_failed + next + } da_y <- da_counts[, da_j] if (length(unique(da_y)) < 2L) { da_result[da_j, "status"] <- "failed" @@ -263,18 +341,58 @@ for (da_j in seq_len(da_n)) { 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 + # "ml": theta is estimated per taxon by glm.nb's own alternating iterations. + # "glmGamPoi": theta is TAKEN from the shrunken dispersion fitted above + # (theta = 1 / phi) and the per-taxon fit goes through glm() with MASS's + # negative.binomial family, so the Wald test and every check below are the + # same machinery on a given dispersion instead of an estimated one. + if (da_use_gp) { + if (!is.finite(da_phi[da_j]) || da_phi[da_j] < 0) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the dispersion trend for this taxon is invalid" + next + } + if (da_phi[da_j] == 0) { + da_result[da_j, "status"] <- "boundary" + da_result[da_j, "note"] <- "zero dispersion: Poisson boundary; NB Wald inference withheld" + next + } + da_theta <- 1 / da_phi[da_j] + if (!is.finite(da_theta) || da_theta <= 0) { + da_result[da_j, "status"] <- "failed" + da_result[da_j, "note"] <- "the shrunken dispersion for this taxon is not usable" + next + } + da_fit <- tryCatch( + withCallingHandlers( + glm(da_y ~ da_g + offset(da_offset), + family = MASS::negative.binomial(theta = da_theta), + 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("the negative-binomial fit stopped with an error:", + conditionMessage(da_fit)) + next + } + } else { + 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"]]) || @@ -284,7 +402,7 @@ for (da_j in seq_len(da_n)) { paste(c(da_fit[["th.warn"]], da_warns), collapse = "; ")) next } - da_co <- summary(da_fit)[["coefficients"]] + 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)" @@ -296,15 +414,23 @@ for (da_j in seq_len(da_n)) { da_result[da_j, "note"] <- "the fit returned a non-finite estimate, standard error or p-value" next } - da_theta <- da_fit[["theta"]] + # "glmGamPoi" set da_theta before the fit; "ml" takes it from the fit. + if (!da_use_gp) 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" + da_result[da_j, "note"] <- if (da_use_gp) + paste("the shrunken dispersion puts these counts at the Poisson end:", + "no overdispersion is left, so the fit is effectively Poisson") + else + "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" + da_result[da_j, "note"] <- if (da_use_gp) + "the shrunken dispersion puts these counts at the extreme-variance end" + else + "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 +445,64 @@ 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="ml") -> 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 where each taxon's theta comes from: `"ml"` (the default) +takes `MASS::glm.nb`'s maximum-likelihood estimate; `"glmGamPoi"` uses the +`overdispersion_shrinkage_list.dispersion_trend` from `glmGamPoi::glm_gp` on +all supplied taxa in one call, then fits each taxon with `theta = 1 / dispersion` fixed. Both options +use a Wald test of the group coefficient. 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: theta +at or above 1e7 marks the Poisson end, and theta at or below 1e-8 marks extreme +variance. Constant counts, complete separation, non-convergence and unusable +coefficients produce failed rows. Errors caught from the model-fitting calls +also produce failed rows; a caught `glm_gp` error or an invalid dispersion +vector fails every taxon, with no fallback to `"ml"`. A zero trend marks only +its taxon as a Poisson boundary and withholds NB Wald statistics. Invalid +negative/non-finite trends fail only the affected taxa. + +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="ml") + d = strip(String(dispersion)) + lowercase(d) in lowercase.(DISPERSIONS) || throw(ArgumentError( + "dispersion must be one of $(join(DISPERSIONS, ", ")), not '$dispersion'")) + use_gp = lowercase(d) == "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_gp && !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_gp] = use_gp 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_gp', 'da_g', 'da_term', " * + "'da_n', 'da_result', 'da_phi', 'da_gp_failed', 'da_gp', 'da_j', 'da_y', " * + "'da_warns', 'da_fit', 'da_co', 'da_row', 'da_theta', 'da_absent'))); invisible(gc())") end end DataFrame(status = String.(raw.status), @@ -589,7 +743,10 @@ 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 per-taxon maximum-likelihood +theta (`"ml"`) or fixed theta from `glmGamPoi::glm_gp`'s cross-taxon dispersion trend +(`"glmGamPoi"`), fitted in one call on the taxa passing the prevalence filter. +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 +769,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 +811,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 +823,21 @@ 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 " * + nb_fit = config.dispersion == "ml" ? "MASS::glm.nb" : + "MASS::negative.binomial with dispersion trend from glmGamPoi::glm_gp" + 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" + if config.dispersion == "glmGamPoi" + method_text *= "; the trend is treated as fixed, without propagating its " * + "estimation uncertainty; this is not glmGamPoi's quasi-likelihood test" + end effect = Dict("key" => "log2_fold_change", "label" => "log2 fold change") size_factor_rows = [Dict("sample" => String(samples[i]), "group" => String(groups[i]), "factor" => factors[i]) for i in eachindex(samples)] config_echo = Dict{String,Any}("method" => config.method, "offset" => config.offset, - "min_prevalence" => config.min_prevalence) + "min_prevalence" => config.min_prevalence, + "dispersion" => config.dispersion) else (n_ref >= 2 && n_con >= 2) || throw(ArgumentError( "Welch's t-test needs at least 2 samples in each group to estimate its variance; " * diff --git a/src/core/r_runtime.jl b/src/core/r_runtime.jl index 0fa982da..5fca35e3 100644 --- a/src/core/r_runtime.jl +++ b/src/core/r_runtime.jl @@ -25,6 +25,7 @@ const _R_LOCK = ReentrantLock() const REQUIRED_R_PACKAGES = [ "dada2" => "the DADA2 stages that run on this machine", "vegan" => "NMDS, PERMANOVA and PERMDISP", + "glmGamPoi" => "shared dispersion-trend differential abundance", "MASS" => "negative-binomial differential abundance", "zCompositions" => "CLR differential abundance (clr_welch)", ] diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl index 0dd506a3..2d3af55f 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="ml"`. 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", "ml"))) 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 2f927ba1..f4f2ba0b 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, @@ -130,6 +131,11 @@ factory default; a key whose default is ~ accepts a string, number or boolean. """ function _config_value_error(key::String, value) v = value isa JSON3.Array ? collect(value) : value + if key == "analysis.differential.dispersion" && v isa AbstractString + v = lowercase(strip(v)) + return v in lowercase.(Differential.DISPERSIONS) ? nothing : + "$key must be one of $(join(Differential.DISPERSIONS, ", "))" + end if haskey(_ENUM_VALUES, key) allowed = _ENUM_VALUES[key] v in allowed || return "$key must be one of $(join(filter(!isnothing, collect(allowed)), ", "))" @@ -177,50 +183,172 @@ function _config_value_error(key::String, value) nothing end +# Validate the effective pair after merging a PATCH, before writing any key. +# Resolved entries carry source annotations; global defaults are raw values. +function _differential_pair_error(resolved::AbstractDict, changes::AbstractDict) + keys_ = ("analysis.differential.method", "analysis.differential.dispersion") + any(k -> haskey(changes, k), keys_) || return nothing + value(k, default) = if haskey(changes, k) + changes[k] + else + current = get(resolved, k, default) + current isa NamedTuple ? current.value : current + end + method = value(keys_[1], "nb_glm") + dispersion = value(keys_[2], "ml") + if method isa AbstractString && dispersion isa AbstractString && + lowercase(strip(method)) == "clr_welch" && lowercase(strip(dispersion)) == "glmgampoi" + return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; use dispersion: ml with clr_welch" + end + nothing +end + +# Resolve every affected study/group/run against the prospective document. +# Enumerate runnable descendants too: a run need not own a pipeline.yml file. +# This is validation evidence, so an unreadable cascade must refuse the save. +function _differential_cascade_error(path::String, document) + data = ServerState.data_dir() + root = dirname(data) + global_path = joinpath(root, "config", "pipeline.yml") + prefix = path == global_path ? String[] : splitpath(relpath(dirname(path), data)) + cache = Dict{String,Any}(path => _flatten(document)) + flat(file) = get!(cache, file) do + _flatten(_load_yml(file)) + end + pair = ("analysis.differential.method", "analysis.differential.dispersion") + function check(values, label) + effective = Dict(pair[1]=>get(values, pair[1], "nb_glm"), + pair[2]=>get(values, pair[2], "ml")) + bad = _differential_pair_error(effective, effective) + isnothing(bad) ? nothing : "$label: $bad" + end + try + base = merge(flat(joinpath(root, "config", "defaults", "pipeline.yml")), flat(global_path)) + if isempty(prefix) + bad = check(base, "global config") + isnothing(bad) || return bad + end + targets = Set{String}() + isempty(prefix) || push!(targets, joinpath(prefix...)) + if length(prefix) >= 2 + # A run edit has no descendants; a group edit can affect its runs. + # Inspect only the edited subtree, so malformed siblings cannot block it. + if length(prefix) == 2 && !get(document, "pool_children", false) + for run in _group_run_names(prefix[1], prefix[2]) + push!(targets, joinpath(prefix..., run)) + end + end + else + for study in _study_names() + isempty(prefix) || first(prefix) == study || continue + push!(targets, study) + for run in _run_names(study) + push!(targets, joinpath(study, run)) + end + for group in _group_names(study) + push!(targets, joinpath(study, group)) + for run in _group_run_names(study, group) + push!(targets, joinpath(study, group, run)) + end + end + end + end + for label in sort!(collect(targets)) + parts = splitpath(label) + length(parts) >= length(prefix) && parts[1:length(prefix)] == prefix || continue + effective = copy(base) + for i in eachindex(parts) + merge!(effective, flat(joinpath(data, parts[1:i]..., "pipeline.yml"))) + end + bad = check(effective, label) + isnothing(bad) || return bad + end + nothing + catch e + "Cannot validate the affected config cascade: $(sprint(showerror, e))" + end +end + +# Validate, prepare the complete document, and publish once under one lock. +# Readers do not take this lock, so per-key renames would expose partial PATCHes. +function _write_config_updates(path::String, changes::AbstractDict; remove=()) + lock(_config_file_lock) do + native = Dict{String,Any}(string(k) => _to_plain(v) for (k, v) in changes) + key = "analysis.differential.dispersion" + if haskey(native, key) && native[key] isa AbstractString + i = findfirst(v -> lowercase(v) == lowercase(strip(native[key])), Differential.DISPERSIONS) + isnothing(i) || (native[key] = Differential.DISPERSIONS[i]) + end + isempty(native) && isempty(remove) && return nothing + cfg = _load_yml(path) + for (k, v) in native + _set_config_value!(cfg, k, v) + end + removed = false + for k in remove + removed |= _remove_config_value!(cfg, k) + end + isempty(native) && !removed && return nothing + pair = ("analysis.differential.method", "analysis.differential.dispersion") + if any(k -> haskey(native, k) || k in remove, pair) + bad = _differential_cascade_error(path, cfg) + isnothing(bad) || return json_error(400, "invalid_config_value", bad) + end + _atomic_write_yaml(path, cfg) + nothing + end +end + function _validate_config_key(dotted_key::String) dotted_key in _ALLOWED_CONFIG_KEYS || error( "Config key '$dotted_key' is not a recognised pipeline option") end -function _write_override(path::String, dotted_key::String, value) +# Shared in-memory setter. A bad later key must not publish earlier changes. +function _set_config_value!(cfg::AbstractDict, dotted_key::String, value) _validate_config_key(dotted_key) - lock(_config_file_lock) do - cfg = _load_yml(path) - parts = split(dotted_key, ".") - d = cfg - for p in parts[1:end-1] - child = get(d, p, nothing) - if child isa Dict - d = child - elseif isnothing(child) - new_dict = Dict() - d[p] = new_dict - d = new_dict - else - error("Config key segment '$p' is not a section; cannot set sub-key") - end + parts = split(dotted_key, ".") + d = cfg + for p in parts[1:end-1] + child = get(d, p, nothing) + if child isa Dict + d = child + elseif isnothing(child) + new_dict = Dict() + d[p] = new_dict + d = new_dict + else + error("Config key segment '$p' is not a section; cannot set sub-key") end - d[parts[end]] = value - _atomic_write_yaml(path, cfg) end + d[parts[end]] = value + cfg end -function _delete_override(path::String, dotted_key::String) - isfile(path) || return +function _write_override(path::String, dotted_key::String, value) lock(_config_file_lock) do - cfg = _load_yml(path) - parts = split(dotted_key, ".") - d = cfg - for p in parts[1:end-1] - haskey(d, p) || return - d = d[p] - d isa Dict || return - end - delete!(d, parts[end]) + cfg = _load_yml(path) + _set_config_value!(cfg, dotted_key, value) _atomic_write_yaml(path, cfg) end end +function _remove_config_value!(cfg::AbstractDict, dotted_key::String) + parts = split(dotted_key, ".") + d = cfg + for p in parts[1:end-1] + haskey(d, p) || return false + d = d[p] + d isa AbstractDict || return false + end + haskey(d, parts[end]) || return false + delete!(d, parts[end]) + true +end + +_delete_override(path::String, dotted_key::String) = + _write_config_updates(path, Dict(); remove=(dotted_key,)) + ## Routes # Default (global) config endpoints @get "/api/v1/config" function(req) @@ -238,16 +366,16 @@ end bad = _config_value_error(string(k), v) isnothing(bad) || return json_error(400, "invalid_config_value", bad) end - for (k, v) in body - _write_override(user_path, string(k), v) - end + bad = _write_config_updates(user_path, body) + isnothing(bad) || return bad json(Dict(k => (; value=v, source="default") for (k, v) in _default_cfg(root))) end @delete "/api/v1/config/{key}" function(req, key::String) root = dirname(ServerState.data_dir()) user_path = joinpath(root, "config", "pipeline.yml") - _delete_override(user_path, key) + bad = _delete_override(user_path, key) + isnothing(bad) || return bad json(Dict(k => (; value=v, source="default") for (k, v) in _default_cfg(root))) end @@ -502,9 +630,8 @@ end isnothing(bad) || return json_error(400, "invalid_config_value", bad) end path = joinpath(ServerState.data_dir(), study, "pipeline.yml") - for (k, v) in body - _write_override(path, string(k), v) - end + bad = _write_config_updates(path, body) + isnothing(bad) || return bad json(_resolve_config(study)) end @@ -512,7 +639,8 @@ end study in _study_names() || return json_error(404, "study_not_found", "Study '$study' not found") path = joinpath(ServerState.data_dir(), study, "pipeline.yml") - _delete_override(path, key) + bad = _delete_override(path, key) + isnothing(bad) || return bad json(_resolve_config(study)) end @@ -586,9 +714,8 @@ end isnothing(bad) || return json_error(400, "invalid_config_value", bad) end path = joinpath(ServerState.data_dir(), study, group, "pipeline.yml") - for (k, v) in body - _write_override(path, string(k), v) - end + bad = _write_config_updates(path, body) + isnothing(bad) || return bad json(_resolve_config(study, nothing, group)) end @@ -600,7 +727,8 @@ end group in _group_names(study) || return json_error(404, "group_not_found", "Group '$group' not found") path = joinpath(ServerState.data_dir(), study, group, "pipeline.yml") - _delete_override(path, key) + bad = _delete_override(path, key) + isnothing(bad) || return bad json(_resolve_config(study, nothing, group)) end @@ -634,9 +762,8 @@ end isnothing(bad) || return json_error(400, "invalid_config_value", bad) end path = joinpath(ServerState.data_dir(), study, run_path, "pipeline.yml") - for (k, v) in body - _write_override(path, string(k), v) - end + bad = _write_config_updates(path, body) + isnothing(bad) || return bad json(_resolve_config(study, run, group)) end @@ -648,6 +775,7 @@ end group = let g = _req_group(req); isnothing(g) ? _run_group(study, run) : g end run_path = isnothing(group) ? run : joinpath(group, run) path = joinpath(ServerState.data_dir(), study, run_path, "pipeline.yml") - _delete_override(path, key) + bad = _delete_override(path, key) + isnothing(bad) || return bad json(_resolve_config(study, run, group)) end diff --git a/test/integration/test_server.jl b/test/integration/test_server.jl index 184e70fa..4cc02029 100644 --- a/test/integration/test_server.jl +++ b/test/integration/test_server.jl @@ -65,6 +65,22 @@ r4 = HTTP.get("http://localhost:$port/api/v1/databases"; status_exception=false) @test r4.status == 200 + ## Configuration persistence through the real listener/middleware. + config_url = "http://localhost:$port/api/v1/config" + headers = ["Content-Type"=>"application/json", "Origin"=>"http://localhost:$port"] + for (submitted, canonical) in ((" ML ", "ml"), (" GLMGAMPOI ", "glmGamPoi")) + response = HTTP.patch(config_url, headers, + JSON3.write(Dict("analysis.differential.method"=>"nb_glm", + "analysis.differential.dispersion"=>submitted)); status_exception=false) + @test response.status == 200 + @test JSON3.read(String(response.body))["analysis.differential.dispersion"].value == canonical + response = HTTP.get(config_url; status_exception=false) + @test response.status == 200 + @test JSON3.read(String(response.body))["analysis.differential.dispersion"].value == canonical + saved = YAML.load_file(joinpath(tmp_root, "config", "pipeline.yml")) + @test saved["analysis"]["differential"]["dispersion"] == canonical + end + ## Path traversal guard -> 403 r5 = HTTP.get("http://localhost:$port/files/smoke_study/runs/run_A/../../../etc/passwd"; status_exception=false) diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 3b034896..e5c12d20 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -10,7 +10,8 @@ using MetaManifold using MetaManifold.Differential using RCall using Test -using DuckDB, DBInterface, JSON3 +using DuckDB, DBInterface, JSON3, HTTP, YAML +using OrderedCollections: OrderedDict """ _r_has(pkg) -> Bool @@ -29,6 +30,8 @@ end const HAVE_R = _r_has("stats") const HAVE_MASS = HAVE_R && _r_has("MASS") HAVE_MASS || @info "R or MASS unavailable - the negative-binomial fit tests are SKIPPED" +const HAVE_GLMGAMPOI = HAVE_MASS && _r_has("glmGamPoi") +HAVE_GLMGAMPOI || @info "R, MASS or glmGamPoi unavailable - the shrunken-dispersion 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" @@ -172,6 +175,14 @@ 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: ml is the default, glmGamPoi is the alternative, + # and it is a negative-binomial setting only. + @test DifferentialConfig().dispersion == "ml" + @test DifferentialConfig(dispersion="glmGamPoi").dispersion == "glmGamPoi" + @test DifferentialConfig(dispersion=" GLMGAMPOI ").dispersion == "glmGamPoi" + @test_throws ArgumentError DifferentialConfig(dispersion="shrinkage") + @test_throws ArgumentError DifferentialConfig(method="clr_welch", dispersion="glmGamPoi") + @test DifferentialConfig(method="clr_welch").dispersion == "ml" for d in (0, 1, -0.5, NaN) @test_throws ArgumentError DifferentialConfig(; method="clr_welch", replacement_delta=d) end @@ -211,6 +222,9 @@ end reference="A", contrast="B") byt = Dict(r["taxon"] => r for r in res["rows"]) @test res["status"] == "ok" + @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 res["diagnostics"]["n_tested"] == 3 @test all(r -> r["status"] == "ok", res["rows"]) @test isapprox(byt["t_up"]["estimate"], log(4); atol=0.5) @@ -249,6 +263,87 @@ end end end + @testset "glmGamPoi matches a direct negative-binomial GLM with fixed dispersion" begin + if !HAVE_GLMGAMPOI + @test get(ENV, "CI", "false") != "true" # the pinned dependency is required 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 (MASS::negative.binomial " * + "with dispersion trend from glmGamPoi::glm_gp, log link) with a " * + "TSS size-factor offset; Wald test of the group coefficient; " * + "Benjamini-Hochberg adjustment over the fitted taxa; the trend " * + "is treated as fixed, without propagating its estimation uncertainty; " * + "this is not glmGamPoi's quasi-likelihood test" + byt = Dict(r["taxon"] => r for r in res["rows"]) + off = log.(tss_factors(effect, effect_samples)) + g = effect_groups + phi = rcopy(R""" + local({ + gp <- glmGamPoi::glm_gp(t($effect), + design = model.matrix(~ factor($g, levels = c("A", "B"))), + size_factors = exp($off), overdispersion_shrinkage = TRUE) + as.numeric(gp$overdispersion_shrinkage_list$dispersion_trend) + }) + """) + for (j, t) in enumerate(effect_taxa) + y = effect[:, j] + theta = 1 / phi[j] + direct = rcopy(R""" + local({ + fit <- glm($y ~ factor($g, levels = c("A", "B")) + offset($off), + family = MASS::negative.binomial(theta = $theta), + control = glm.control(maxit = 100)) + unname(summary(fit, dispersion = 1)$coefficients[2, ]) + }) + """) + r = byt[t] + @test r["status"] == "ok" + @test r["estimate"] ≈ direct[1] rtol=1e-6 + @test r["standard_error"] ≈ direct[2] rtol=1e-6 + @test r["statistic"] ≈ direct[3] rtol=1e-6 + @test r["pvalue"] ≈ direct[4] rtol=1e-6 + @test r["dispersion_theta"] ≈ theta rtol=1e-6 + end + end + end + + @testset "glmGamPoi invalid trends are isolated per taxon" begin + if !HAVE_GLMGAMPOI + @test_skip false + else + MetaManifold.RRuntime.with_r_lock() do + reval(raw""" + test_original_gp <- get("glm_gp", asNamespace("glmGamPoi")) + test_phi <- 0 + assignInNamespace("glm_gp", function(...) { + fit <- test_original_gp(...) + fit$overdispersion_shrinkage_list$dispersion_trend[1] <- test_phi + fit + }, ns="glmGamPoi") + """) + try + off = log.(tss_factors(effect, effect_samples)) + for (phi, status) in ((0.0, "boundary"), (-1.0, "failed"), (NaN, "failed")) + RCall.globalEnv[:test_phi] = phi + rows = fit_nb(effect, effect_groups, off; levels=("A", "B"), dispersion="glmGamPoi") + @test rows.status[1] == status + @test isnothing(rows.pvalue[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_phi) + """) + end + end + end + end + @testset "Degenerate and filtered taxa are reported, not dropped" begin if !HAVE_MASS @test_skip false @@ -353,6 +448,205 @@ 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")) + for dispersion in ("ml", "glmGamPoi", " ML ", " GLMGAMPOI ") + @test isnothing(SV._config_value_error("analysis.differential.dispersion", dispersion)) + end + for dispersion in ("shrinkage", "", nothing, 1, true, ["ml"]) + @test !isnothing(SV._config_value_error("analysis.differential.dispersion", dispersion)) + end + end + + @testset "Config rejects incompatible effective dispersion pairs" begin + SV = MetaManifold.Server + keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" + @test SV._differential_pair_error(Dict(keym=>"clr_welch"), Dict(keyd=>"glmGamPoi")) !== nothing + @test SV._differential_pair_error(Dict(keyd=>"glmGamPoi"), Dict(keym=>"clr_welch")) !== nothing + @test SV._differential_pair_error(Dict(keym=>(; value="clr_welch",source="study")), Dict(keyd=>" GLMGAMPOI ")) !== nothing + @test SV._differential_pair_error(Dict(keym=>"clr_welch"), Dict(keym=>"nb_glm",keyd=>"glmGamPoi")) === nothing + @test SV._differential_pair_error(Dict(keyd=>"glmGamPoi"), Dict(keym=>"clr_welch",keyd=>"ml")) === nothing + _da_fixture() do + root = SV.ServerState._root[] + path = joinpath(root,"data","studyD","pipeline.yml") + write(path,"analysis:\n differential:\n method: clr_welch\n") + before = read(path,String) + response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", "/api/v1/studies/studyD/config", ["Content-Type"=>"application/json"], JSON3.write(Dict(keyd=>"glmGamPoi")))) + @test response.status == 400 + @test read(path,String) == before + # A batch can change both options together, irrespective of key order. + batch = Dict(keyd=>"glmGamPoi", keym=>"nb_glm") + response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", "/api/v1/studies/studyD/config", ["Content-Type"=>"application/json"], JSON3.write(batch))) + @test response.status == 200 + before = read(path,String) + response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", "/api/v1/studies/studyD/config", ["Content-Type"=>"application/json"], JSON3.write(Dict(keym=>"clr_welch")))) + @test response.status == 400 + @test read(path,String) == before + end + end + + @testset "Config PATCH publishes canonical, complete documents" begin + SV = MetaManifold.Server + keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" + _da_fixture() do + root = SV.ServerState._root[] + group_dir = joinpath(root, "data", "studyD", "groupG") + run_dir = joinpath(group_dir, "runG") + mkpath(run_dir) + touch(joinpath(group_dir, ".group")) + touch(joinpath(run_dir, "sample_R1.fastq.gz")) + touch(joinpath(run_dir, "sample_R2.fastq.gz")) + scopes = [("/api/v1/config", joinpath(root, "config", "pipeline.yml")), + ("/api/v1/studies/studyD/config", joinpath(root, "data", "studyD", "pipeline.yml")), + ("/api/v1/studies/studyD/groups/groupG/config", joinpath(group_dir, "pipeline.yml")), + ("/api/v1/studies/studyD/runs/runG/config?group=groupG", joinpath(run_dir, "pipeline.yml"))] + for (url, path) in scopes + for (submitted, expected) in ((" ML ", "ml"), (" GLMGAMPOI ", "glmGamPoi")) + response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", url, + ["Content-Type"=>"application/json"], JSON3.write(Dict(keym=>"nb_glm", keyd=>submitted, "seed"=>456)))) + @test response.status == 200 + @test JSON3.read(String(response.body))[keyd].value == expected + @test YAML.load_file(path)["analysis"]["differential"]["dispersion"] == expected + response = SV.Oxygen.internalrequest(HTTP.Request("GET", url)) + @test response.status == 200 + @test JSON3.read(String(response.body))[keyd].value == expected + @test YAML.load_file(path)["seed"] == 456 + end + end + # Both transitions must preserve unrelated fields. Explicit ordering + # exercises the formerly invalid intermediate clr_welch/glmGamPoi pair. + path = scopes[2][2] + for changes in (OrderedDict(keyd=>"ml", keym=>"clr_welch"), + OrderedDict(keyd=>"glmGamPoi", keym=>"nb_glm"), + OrderedDict(keym=>"clr_welch", keyd=>"ml")) + @test SV._write_config_updates(path, changes) === nothing + saved = YAML.load_file(path) + @test saved["seed"] == 456 + @test saved["analysis"]["differential"]["method"] == changes[keym] + @test saved["analysis"]["differential"]["dispersion"] == changes[keyd] + end + before = read(path, String) + bad = SV._write_config_updates(path, Dict(keyd=>"glmGamPoi")) + @test bad.status == 400 + @test read(path, String) == before + + run_url, run_path = scopes[4] + before_run = read(run_path, String) + response = SV.Oxygen.internalrequest(HTTP.Request("PATCH", run_url, + ["Content-Type"=>"application/json"], JSON3.write(Dict(keym=>"clr_welch")))) + @test response.status == 400 + @test read(run_path, String) == before_run + + # The first valid key must stay unpublished when a later setter fails. + YAML.write_file(path, Dict("phylogeny"=>"not-a-section", "seed"=>456)) + before = read(path, String) + @test_throws ErrorException SV._write_config_updates(path, + OrderedDict("seed"=>789, "phylogeny.reference.trim.method"=>"none")) + @test read(path, String) == before + end + end + + @testset "Parent PATCH validates affected descendants before publication" begin + SV = MetaManifold.Server + keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" + _da_fixture() do + root = SV.ServerState._root[] + global_path = joinpath(root, "config", "pipeline.yml") + study_path = joinpath(root, "data", "studyD", "pipeline.yml") + group_dir = joinpath(dirname(study_path), "groupG") + run_dir = joinpath(group_dir, "runG") + mkpath(run_dir); touch(joinpath(group_dir, ".group")) + touch(joinpath(run_dir, "sample_R1.fastq.gz")) + group_path, run_path = joinpath(group_dir, "pipeline.yml"), joinpath(run_dir, "pipeline.yml") + save(path, values) = YAML.write_file(path, Dict("analysis"=>Dict("differential"=>values))) + patch(url, changes) = SV.Oxygen.internalrequest(HTTP.Request("PATCH", url, + ["Content-Type"=>"application/json"], JSON3.write(changes))) + save(global_path, Dict("method"=>"nb_glm", "dispersion"=>"ml")) + save(study_path, Dict("dispersion"=>"glmGamPoi")) + before = read(global_path, String) + response = patch("/api/v1/config", Dict(keym=>"clr_welch", "seed"=>789)) + @test response.status == 400 + @test occursin("studyD", JSON3.read(String(response.body)).message) + @test read(global_path, String) == before + # A deeper complete pair shields that subtree from a parent edit. + save(study_path, Dict("method"=>"nb_glm", "dispersion"=>"glmGamPoi")) + @test patch("/api/v1/config", Dict(keym=>"clr_welch", keyd=>"ml")).status == 200 + + save(global_path, Dict("method"=>"nb_glm", "dispersion"=>"ml")) + save(study_path, Dict("method"=>"nb_glm", "dispersion"=>"ml")) + save(group_path, Dict("dispersion"=>"glmGamPoi")) + before = read(study_path, String) + @test patch("/api/v1/studies/studyD/config", Dict(keym=>"clr_welch")).status == 400 + @test read(study_path, String) == before + save(group_path, Dict("method"=>"nb_glm", "dispersion"=>"ml")) + save(run_path, Dict("dispersion"=>"glmGamPoi")) + before = read(group_path, String) + group_url = "/api/v1/studies/studyD/groups/groupG/config" + response = patch(group_url, Dict(keym=>"clr_welch")) + @test response.status == 400 + @test occursin("runG", JSON3.read(String(response.body)).message) + @test read(group_path, String) == before + save(run_path, Dict("method"=>"nb_glm", "dispersion"=>"glmGamPoi")) + @test patch(group_url, Dict(keym=>"clr_welch")).status == 200 + + # A run without an override file still inherits a validated cascade. + run_a = joinpath(root, "data", "studyD", "runA", "pipeline.yml") + rm(run_a) + @test patch("/api/v1/studies/studyD/config", Dict(keyd=>"glmGamPoi")).status == 200 + response = SV.Oxygen.internalrequest(HTTP.Request("GET", "/api/v1/studies/studyD/runs/runA/config")) + @test response.status == 200 + @test JSON3.read(String(response.body))[keyd].value == "glmGamPoi" + + # A malformed sibling is outside a run edit's affected subtree. + write(joinpath(root,"data","studyD","runB","pipeline.yml"), "analysis: [broken") + @test SV._write_config_updates(run_a, Dict(keyd=>"ml")) === nothing + end + end + + @testset "DELETE validates inherited values at every config scope" begin + SV = MetaManifold.Server + keym, keyd = "analysis.differential.method", "analysis.differential.dispersion" + for scope in 1:4, deleted in (keym, keyd, "valid") + _da_fixture() do + root = SV.ServerState._root[] + group_dir = joinpath(root, "data", "studyD", "groupG") + run_dir = joinpath(group_dir, "runG") + mkpath(run_dir); touch(joinpath(group_dir, ".group")) + touch(joinpath(run_dir, "sample_R1.fastq.gz")) + factory = joinpath(root, "config", "defaults", "pipeline.yml") + mkpath(dirname(factory)) + paths = [joinpath(root,"config","pipeline.yml"), + joinpath(root,"data","studyD","pipeline.yml"), + joinpath(group_dir,"pipeline.yml"), joinpath(run_dir,"pipeline.yml")] + urls = ["/api/v1/config", "/api/v1/studies/studyD/config", + "/api/v1/studies/studyD/groups/groupG/config", + "/api/v1/studies/studyD/runs/runG/config"] + query = scope == 4 ? "?group=groupG" : "" + method, dispersion = deleted == keym ? ("clr_welch", "ml") : + deleted == keyd ? ("nb_glm", "glmGamPoi") : ("nb_glm", "ml") + doc(m, d, seed) = Dict("seed"=>seed, "analysis"=>Dict("differential"=>Dict("method"=>m, "dispersion"=>d))) + for path in [factory; paths] + YAML.write_file(path, doc(method, dispersion, 123)) + end + target_method, target_dispersion = deleted == keym ? ("nb_glm", "glmGamPoi") : + deleted == keyd ? ("clr_welch", "ml") : ("nb_glm", "ml") + YAML.write_file(paths[scope], doc(target_method, target_dispersion, 456)) + before = Dict(p=>read(p,String) for p in [factory; paths]) + if deleted != "valid" + response = SV.Oxygen.internalrequest(HTTP.Request("DELETE", "$(urls[scope])/$deleted$query")) + @test response.status == 400 + @test JSON3.read(String(response.body)).error == "invalid_config_value" + @test all(read(p,String) == text for (p,text) in before) + else + for (key, expected) in ((keym,"nb_glm"), (keyd,"ml"), ("seed",123)) + response = SV.Oxygen.internalrequest(HTTP.Request("DELETE", "$(urls[scope])/$key$query")) + @test response.status == 200 + @test JSON3.read(String(response.body))[key].value == expected + end + unchanged = read(paths[scope], String) + @test SV.Oxygen.internalrequest(HTTP.Request("DELETE", "$(urls[scope])/seed$query")).status == 200 + @test read(paths[scope], String) == unchanged + end + end + end end @testset "CLR transform" begin @@ -720,4 +1014,59 @@ 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 # no MASS means the ml 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)) + _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 must still fit normally with it in place. The + # injected-effect dataset, not the minimal matrix above: + # MEASURED with MASS::glm.nb, with 4 samples of small counts the + # alternating fit runs theta to its upper bound (status boundary) + # and one taxon's fit even produced NaNs (status failed), so that + # data never reaches status ok on any machine; the dataset + # "Negative-binomial fit recovers an injected effect" already + # proves fits. + 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