Prepare isolated R glmGamPoi dispersion option for upstream

.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
           '

R/_renv_dependencies.R

@@ -3,6 +3,7 @@
 if (FALSE) {
     library(dada2)
     library(vegan)
+    library(glmGamPoi)
     library(MASS)
     library(zCompositions)
 }

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

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
   }

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'],
+  })
+})

frontend/src/components/PipelineStages.tsx

@@ -152,6 +152,7 @@ export const CONFIG_DESCRIPTIONS: Record<string, string> = {
   '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<string, ConfigType> = {
   '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' },

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)],