Skip to content

Commit ec2b45d

Browse files
docs(analysis): clarify differential abundance dispersion and error contracts
1 parent 7c4e484 commit ec2b45d

1 file changed

Lines changed: 37 additions & 19 deletions

File tree

‎src/analysis/differential.jl‎

Lines changed: 37 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -61,13 +61,12 @@ have at least one read to be tested. `replacement_delta`, in (0, 1), is the
6161
above its taxon's smallest observed proportion is lowered to this fraction of it.
6262
`dispersion` names the per-taxon dispersion `nb_glm` uses (`"ml"` or
6363
`"glmGamPoi"`): `"ml"` takes the maximum-likelihood theta from `MASS::glm.nb`,
64-
`"glmGamPoi"` takes the shrunken dispersion from `glmGamPoi::glm_gp` -- fitted
65-
across all taxa in one call, which is where its shared information helps low-count
66-
taxa -- and fits each taxon through `MASS::negative.binomial(theta = ...)` so the
67-
Wald test and its checks are unchanged. It is `nb_glm` only.
68-
Method and offset names are stripped and lowercased. Throws `ArgumentError`
69-
for unsupported names or out-of-range fractions, including settings unused by
70-
the selected method.
64+
`"glmGamPoi"` uses the `overdispersions` returned by `glmGamPoi::glm_gp` across
65+
the tested taxa, with `theta = 1/dispersion` fixed in each taxon's model.
66+
Method and offset names are stripped and lowercased; dispersion names are
67+
stripped and matched without regard to case. Throws `ArgumentError` for
68+
unsupported names, `"glmGamPoi"` with `"clr_welch"`, or out-of-range fractions,
69+
including settings unused by the selected method.
7170
"""
7271
struct DifferentialConfig
7372
method :: String
@@ -147,6 +146,12 @@ Raised when the R package glmGamPoi, which provides the shrunken dispersion
147146
"""
148147
struct GlmGamPoiUnavailable <: Exception end
149148

149+
"""
150+
Base.showerror(io::IO, error::GlmGamPoiUnavailable)
151+
152+
Write the missing-package error to `io`, with installation instructions and the
153+
`ml` dispersion alternative.
154+
"""
150155
Base.showerror(io::IO, ::GlmGamPoiUnavailable) = print(io,
151156
"the R package glmGamPoi is not available, so the shrunken dispersion cannot be " *
152157
"fitted; install it from Bioconductor (Rscript -e \"BiocManager::install('glmGamPoi')\") " *
@@ -425,19 +430,27 @@ _num(x) = (ismissing(x) || isnothing(x) || !isfinite(x)) ? nothing : Float64(x)
425430
"""
426431
fit_nb(counts, groups, offset; levels, dispersion="ml") -> DataFrame
427432
428-
Fit `MASS::glm.nb(y ~ group + offset(offset))` to each column of the samples x
429-
taxa integer matrix `counts`. `levels` is `(reference, contrast)`, so the
433+
Fit a negative-binomial model to each column of the samples x taxa integer
434+
matrix `counts`. `offset` contains one natural-log size factor per sample;
435+
`groups` labels the same rows. `levels` is `(reference, contrast)`, so the
430436
estimate is the natural-log fold change of `contrast` over `reference`.
431437
`dispersion` names where each taxon's theta comes from: `"ml"` (the default)
432-
takes `glm.nb`'s own maximum-likelihood estimate, `"glmGamPoi"` takes the
433-
shrunken dispersion `glmGamPoi::glm_gp` fits across all taxa in one call and
434-
runs each per-taxon fit through `MASS::negative.binomial(theta = 1/phi)`, so
435-
the Wald test and its checks are unchanged.
438+
takes `MASS::glm.nb`'s maximum-likelihood estimate; `"glmGamPoi"` uses the
439+
`overdispersions` returned by `glmGamPoi::glm_gp` across all supplied taxa and
440+
fits each taxon with `glm` and `MASS::negative.binomial(theta = 1/phi)`.
436441
Returns one row per taxon with `status` (`ok`, `boundary` or `failed`), `note`,
437442
`estimate`, `se`, `statistic`, `pvalue` and `theta`; a failed taxon has
438-
`nothing` for every statistic. Throws `MASSUnavailable` when MASS cannot be
439-
loaded, `GlmGamPoiUnavailable` when `dispersion` is `"glmGamPoi"` and glmGamPoi
440-
cannot be loaded, and `RBusyError` when the R runtime stays busy.
443+
`nothing` for every statistic. Caught model-fitting errors become failed rows;
444+
a `glm_gp` fit error or an invalid dispersion vector fails every taxon without
445+
falling back to `ml`. Boundary rows retain their statistics when theta is at or
446+
outside [1e-8, 1e7], or non-finite (reported as `nothing`).
447+
448+
Throws `ArgumentError` for an unsupported dispersion name or mismatched sample,
449+
group and offset lengths. Dispersion names are stripped and matched without
450+
regard to case. Throws `MASSUnavailable` when MASS cannot be loaded,
451+
`GlmGamPoiUnavailable` when `dispersion` is `"glmGamPoi"` and glmGamPoi cannot
452+
be loaded, and `RBusyError` when the R runtime stays busy. R evaluation errors
453+
outside the caught model-fitting calls propagate to the caller.
441454
"""
442455
function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:AbstractString},
443456
offset::AbstractVector{<:Real}; levels::NTuple{2,String},
@@ -707,7 +720,10 @@ labelled `reference` in `groups`, by `config.method`.
707720
708721
`nb_glm`: size factors are computed on the full matrix by `config.offset` and
709722
enter each negative-binomial model as `log(factor)`; the estimate is a
710-
natural-log fold change.
723+
natural-log fold change. `config.dispersion` selects per-taxon maximum-likelihood
724+
theta (`"ml"`) or fixed theta from `glmGamPoi::glm_gp` (`"glmGamPoi"`), fitted
725+
jointly over the taxa passing `config.min_prevalence`. The selected dispersion
726+
is echoed in the result's `config`.
711727
712728
`clr_welch`, in this order: (1) `min_prevalence` chooses the taxa, and taxa with
713729
no reads in any sample are always left out; (2) the zeros of those taxa only
@@ -730,8 +746,10 @@ effect size and its label.
730746
731747
Throws `ArgumentError` for unusable input (including a zero replacement that
732748
cannot be made), `ScalingRefusal` when the offset cannot be computed,
733-
`MASSUnavailable` or `ZCompositionsUnavailable` when the R package a method
734-
needs is missing, and an `ErrorException` when no taxon could be fitted at all.
749+
`MASSUnavailable`, `GlmGamPoiUnavailable` or `ZCompositionsUnavailable` when the
750+
R package a method needs is missing, `RBusyError` when the R runtime stays busy,
751+
and an `ErrorException` when no taxon could be fitted at all. R evaluation
752+
errors not converted to failed rows by the fitting functions also propagate.
735753
"""
736754
function differential_abundance(counts::AbstractMatrix{<:Real},
737755
samples::AbstractVector{<:AbstractString},

0 commit comments

Comments
 (0)