ALDEx3

Overview

ALDEx3 is the successor to ALDEx2, providing improved computational efficiency and support for both fixed- and mixed-effects models. It enables statistically principled inference on both relative and absolute abundances from sequencing count data (e.g., RNA-seq or 16S rRNA gene sequencing).

Like ALDEx2, ALDEx3 models sequencing counts using a Dirichlet–multinomial framework, which accounts for sampling variability and naturally accommodates zero counts through the multinomial likelihood (e.g., the model does not require treating zeros as true structural absences).

High-throughput sequencing data contain only relative information: the overall biological scale is unobserved (unidentifiable). Here, scale refers to quantities such as total microbial load, total cellular RNA content, or total biomass. Standard normalization-based workflows implicitly fix this unmeasured scale using strong and often untestable assumptions. When these assumptions are violated, normalization can induce substantial bias, inflate false discovery rates, or mask true biological effects.

ALDEx3 addresses this limitation by explicitly modeling scale uncertainty. Rather than fixing scale through normalization, ALDEx3 uses scale models—probability distributions that represent uncertainty in the unobserved scale and propagate that uncertainty into downstream inference. This approach generalizes normalization while making its assumptions explicit and testable.

Scale models can be specified in several ways, including: 1. As uncertainty induced by normalization procedures (e.g., total-sum scaling or log-ratio–based approaches), 2. Using external measurements of scale (e.g., flow cytometry, qPCR, or spike-ins), or 3. Through weak biological prior information (e.g., assuming an antibiotic reduces total microbial load by 20–80%).

By treating scale as a latent random variable rather than a fixed normalization constant, ALDEx3 enables coherent inference for both relative and absolute effects and supports complex experimental designs via mixed-effects modeling.

Before using ALDEx3, users are strongly encouraged to read:

Nixon, Gloor, and Silverman (2025).
Incorporating scale uncertainty in microbiome and gene expression analysis as an extension of normalization.
Genome Biology.
https://link.springer.com/article/10.1186/s13059-025-03609-3

This article provides the theoretical motivation for scale models, explains why normalization implicitly fixes unmeasured scale, and demonstrates how explicitly modeling scale uncertainty improves robustness and reproducibility in differential abundance and expression analyses.

Installation

if (!require("devtools", quietly = TRUE)) {
  install.packages("devtools")
}
devtools::install_github("jsilve24/ALDEx3")

Quickstart (Inference with Relative Abundances)

Note: the following examples use the CLR scale model often with zero scale uncertainty gamma=0. This is generally a bad idea but is shown here just for backwards comparability with the defaults in ALDEx2. In general we strongly recommend carefully thinking about the choice of scale model (see https://link.springer.com/article/10.1186/s13059-025-03609-3 for a background on the subject). In the worst case scenario where the user is unwilling to carefully design a scale model, we recommend choosing scale=tss.sm and gamma=0.5as a default.

set.seed(42)
# Load Crohn's Disease vs. Control Data
data(gut_crohns_data)

# NOTE: We recommend a ~75% sparsity filter of counts!
# This is a pragmatic default for speed and stability in very sparse tables;
# adjust based on your depth and scientific question. Always keep an aggregated
# ‘other’ category so each sample’s total count is preserved under the count
# model.
Y <- gut_crohns_data$counts
keep_names <- row.names(Y[((rowSums(Y==0))/ncol(Y))<=0.75,])
other <- colSums(Y[((rowSums(Y==0))/ncol(Y))>0.75,])
Y <- Y[keep_names,]
Y <- rbind(Y, other)

# Metadata Samples x Metadata Values
X <- gut_crohns_data$metadata
X$Health.status <- factor(X$Health.status, levels=c("Control", "CD"))

# Relative Abundances CD vs. Control
aldex.gut.raw <- aldex(Y,
                       ~Health.status,
                       X,
                       nsample=2000,
                       scale=clr.sm,   # CLR assumption
                       gamma=0)        # Gamma=0 no scale uncertainty
aldex.gut.summary <- summary(aldex.gut.raw)
head(aldex.gut.summary)
#>         parameter          entity   estimate std.error   p.val.adj
#> 1 Health.statusCD    Agathobacter  3.4427036  1.500150 0.146352623
#> 2 Health.statusCD     Akkermansia -0.9509523  1.119066 0.797674468
#> 3 Health.statusCD       Alistipes -0.8547089  1.188096 0.868896117
#> 4 Health.statusCD Anaerobutyricum  2.0096506  1.209952 0.411322714
#> 5 Health.statusCD    Anaerostipes -1.0223508  1.010964 0.729975305
#> 6 Health.statusCD     Bacteroides  5.3739786  1.362744 0.003758113

Mixed Effects Models

ALDEx3 implements the published SR-MEM method (Scale Reliant Mixed Effects Models) through the exact lme4 and nlme engines and the approximate blmm engine. The quickstart no longer walks through those models inline. See the dedicated mixed-effects vignette vignettes/ALDEx3-mixed-effects.Rmd for the model setup, the blmm mathematical formulation, and a practical comparison between exact lme4 and approximate blmm results on the oral_mouthwash_data dataset.

Inference with Absolute Abundances: Defining Scale Models

The Quickstart considered the relative abundances of taxa in the guts of patients with Crohn’s Disease (CD) vs. Control. Here absolute abundances are modeled using three different approaches to scale models.

Weak Prior Biological Knowledge

Consider prior studies suggest individuals with CD have on average 2 log2-fold lower total gut microbial load. If θ^tot_CD is the average log2-fold change in gut microbial load in CD vs. Control, we assume θ^tot_CD=-2. However, to account for potential uncertainty in this assumption, we define a scale model as a Normal distribution with standard deviation of \(0.5\): θ^tot_CD ~ N (-2, 0.5²). In ALDEx3 this is done with the coef.sm scale model and the parameters c.mu and c.cor:

set.seed(42)
# Load Data
data(gut_crohns_data)

# NOTE: We recommend a ~75% sparsity filter of counts!
# Genera w/ sparsity > 75% amalgamated into 'other' row
# Keep other counts for multinomial-Dirichlet accuracy
Y <- gut_crohns_data$counts
keep_names <- row.names(Y[((rowSums(Y==0))/ncol(Y))<=0.75,])
other <- colSums(Y[((rowSums(Y==0))/ncol(Y))>0.75,])
Y <- Y[keep_names,]
Y <- rbind(Y, other)

# Metadata Samples x Metadata Values
X <- gut_crohns_data$metadata
X$Health.status <- factor(X$Health.status,
                           levels=c("Control", "CD"))
# Fit Before vs. After teeth brushing
# relative abundance changes
aldex.gut.abs.raw <- aldex(Y,
                       ~Health.status,
                       X,
                       nsample=2000,
                       scale=coefficient.sm,   # coef scale model
                       c.mu=c(0, -2),   # c(intercept, CD)
                       c.cor=rbind(c(0, 0), c(0, 0.5**2)))

aldex.gut.abs.summary <- summary(aldex.gut.abs.raw)
head(aldex.gut.abs.summary)
#>         parameter          entity   estimate std.error    p.val.adj
#> 1 Health.statusCD    Agathobacter -2.3196415  1.585377 2.110150e-01
#> 2 Health.statusCD     Akkermansia -6.7132973  1.137395 5.479350e-06
#> 3 Health.statusCD       Alistipes -6.6170540  1.243634 1.377402e-05
#> 4 Health.statusCD Anaerobutyricum -3.7526945  1.311741 1.652680e-02
#> 5 Health.statusCD    Anaerostipes -6.7846958  1.098533 1.643941e-06
#> 6 Health.statusCD     Bacteroides -0.3883664  1.412088 7.682775e-01

External Scale Measurements (Ideal)

In the Crohn’s dataset we have flow cytometry measurements of fecal microbial load for each sample. If w^tot_n and q_n are the scale and flow cytometry measurement for sample n, respectively, we define the scale model as w^tot_n ~ N (log2 q_n, 0.5²). In ALDEx3, here we use the sample.sm scale model which takes s.mu and s.var as arguments:

set.seed(42)
# Load Data
data(gut_crohns_data)

# NOTE: We recommend a ~75% sparsity filter of counts!
# Genera w/ sparsity > 75% amalgamated into 'other' row
# Keep other counts for multinomial-Dirichlet accuracy
Y <- gut_crohns_data$counts
keep_names <- row.names(Y[((rowSums(Y==0))/ncol(Y))<=0.75,])
other <- colSums(Y[((rowSums(Y==0))/ncol(Y))>0.75,])
Y <- Y[keep_names,]
Y <- rbind(Y, other)

# Metadata Samples x Metadata Values
X <- gut_crohns_data$metadata
X$Health.status <- factor(X$Health.status,
                           levels=c("Control", "CD"))
# Log2 Scale Measurements
s.mu <- log2(X$Average.cell.count..per.gram.of.frozen.feces.)
# Account For Potential Technical Error in Measurements
s.var <- rep(0.25, length(s.mu))
# Estimate Absolute Abundance Changes
aldex.gut.abs.raw <- aldex(Y,
                       ~Health.status,
                       X,
                       nsample=2000,
                       scale=sample.sm,         # coef scale model
                       s.mu=s.mu,
                       s.var=s.var)
aldex.gut.abs.summary <- summary(aldex.gut.abs.raw)
head(aldex.gut.abs.summary)
#>         parameter          entity   estimate std.error    p.val.adj
#> 1 Health.statusCD    Agathobacter -1.8036338  1.704191 3.585030e-01
#> 2 Health.statusCD     Akkermansia -6.1972896  1.171886 1.341654e-05
#> 3 Health.statusCD       Alistipes -6.1010462  1.306783 4.760054e-05
#> 4 Health.statusCD Anaerobutyricum -3.2366868  1.403536 4.306112e-02
#> 5 Health.statusCD    Anaerostipes -6.2686881  1.171873 8.778170e-06
#> 6 Health.statusCD     Bacteroides  0.1276413  1.499862 9.689928e-01

Normalization Uncertainty

The Quickstart used the CLR to infer changes in relative abundances. While not the most powerful approach, we can define a scale model to account for uncertainty in how closely CLR transformed abundances match absolute abundances. Nixon et al. show for this exact dataset that this approach can still substantially reduce false positives.

set.seed(42)
# Load Data
data(gut_crohns_data)

# NOTE: We recommend a ~75% sparsity filter of counts!
# Genera w/ sparsity > 75% amalgamated into 'other' row
# Keep other counts for multinomial-Dirichlet accuracy
Y <- gut_crohns_data$counts
keep_names <- row.names(Y[((rowSums(Y==0))/ncol(Y))<=0.75,])
other <- colSums(Y[((rowSums(Y==0))/ncol(Y))>0.75,])
Y <- Y[keep_names,]
Y <- rbind(Y, other)

# Metadata Samples x Metadata Values
X <- gut_crohns_data$metadata
X$Health.status <- factor(X$Health.status,
                           levels=c("Control", "CD"))
# Fit Before vs. After teeth brushing; absolute abundances
aldex.gut.raw <- aldex(Y,
                       ~Health.status,
                       X,
                       nsample=2000,
                       scale=clr.sm,  # CLR assumption
                       gamma=1)       # Gamma=1: variance in CLR assumption 
aldex.gut.summary <- summary(aldex.gut.raw)
head(aldex.gut.summary)
#>         parameter          entity   estimate std.error  p.val.adj
#> 1 Health.statusCD    Agathobacter  3.4571869  1.500150 0.24152676
#> 2 Health.statusCD     Akkermansia -0.9364690  1.119066 0.68923913
#> 3 Health.statusCD       Alistipes -0.8402256  1.188096 0.72889977
#> 4 Health.statusCD Anaerobutyricum  2.0241339  1.209952 0.45935038
#> 5 Health.statusCD    Anaerostipes -1.0078675  1.010964 0.63593939
#> 6 Health.statusCD     Bacteroides  5.3884619  1.362744 0.02308018

Model Outputs

The output of summary.aldex is typically most useful. It looks like:

        parameter          entity   estimate std.error    p.val.adj
1 Health.statusCD    Agathobacter -1.7951101  1.702948 3.665291e-01
2 Health.statusCD     Akkermansia -6.1827824  1.168809 9.867942e-06
3 Health.statusCD       Alistipes -6.1123999  1.307723 4.871389e-05
4 Health.statusCD Anaerobutyricum -3.2153272  1.397338 4.276630e-02
5 Health.statusCD    Anaerostipes -6.2694178  1.177097 9.477381e-06
6 Health.statusCD     Bacteroides  0.1164949  1.503918 9.687594e-01

The output is a data.frame which includes: - parameter: The name of the regression parameter for the result, typical of the summary output of lm, lmer, etc. - entity: The taxon or gene name for that result (corresponding to row.names of input counts Y) - estimate: The average regression coefficient across Monte Carlo draws - std.err: The average regression standard error across Monte Carlo draws - p.val.adj: The BH adjusted p-values combined from Monte Carlo draws using special procedure

Plotting pairwise analyses

ALDEx3 has simple plotting functions available for pairwise plotting of outputs. Plots include: a volcano plot that plots the model estimate against the -log10 adjusted p-value; an effect plot that uses aldex.effect() to plot the pooled within-group standard deviation against the model estimate; a Bland-Altman, or MA plot that plots reconstructed log abundance (logComp + logScale) against the model estimate; and a waterfall plot that displays the largest significant positive and negative estimates. The contrast name must identify one binary model coefficient. The effect plot and MA plot require the Monte Carlo logComp and logScale arrays, so they are unavailable when aldex() streams large jobs and omits those arrays.

The aldex.effect() helper reports ALDEx2-inspired effect diagnostics for ALDEx3 model contrasts. It returns the mean contrast estimate, mean pooled within-group standard deviation, mean Cohen’s d, and overlap. The overlap statistic is a directional uncertainty diagnostic based on the signs of the Monte Carlo Cohen’s d draws: values near 0 indicate consistent effect direction across Monte Carlo instances, while values near 0.5 indicate uncertainty about direction. It is not a p-value or a literal density-overlap area.

It is possible to change the following aspects of the plots:

plot type: plot: default is ’volcano
minimum fold-change: min.diff, default is 0.5 for the MA plot, this is used as a cutoff by the waterfall plot
Cohen’s d: cohen, default is 0.5 for effect plot
significant feature color: sig.col, default is rgb(1,0,0,0.5)
minimum FDR: threshold, default is 0.05
waterfall number of features: water.show, default is 5 per group
waterfall colors: water.col, default is c(“red”, “blue”)
waterfall names: water.names, logical default is TRUE


par(mfrow=c(2,2))
aldex.plot(aldex.gut.raw, contrast="Health.status", plot='vol', main="Volcano plot")
aldex.plot(aldex.gut.raw, contrast="Health.status", plot='eff', main="Effect plot")
aldex.plot(aldex.gut.raw, contrast="Health.status", plot='MA', main="MA plot")
aldex.plot(aldex.gut.raw, contrast="Health.status", plot="water", 
    water.col=c("orange", "purple"), main="Waterfall plot")

Volcano, effect, Bland-Altman (MA) and waterfall plots of ALDEx3 output for the first example dataset. In the volcano plot, the horizontal line marks an adjusted p-value of 0.05. In the effect plot, the grey lines mark the requested Cohen’s d threshold. In the Bland-Altman or MA plot, the grey lines mark the requested minimum model-estimate difference. Significant features are drawn with the requested significant-feature color. The waterfall plot sorts displayed features by adjusted p-value and estimate.

ALDEx3-Quickstart