DAssemble

DAssemble is an R package for ensemble differential-abundance / differential-expression (DA/DE) analysis across multiple omics domains. It wraps a collection of existing DA/DE methods as core methods, combines them via Cauchy Combination Tests (CCT), and augments them with lightweight enhancers.

---

Installation

1. Install the DAssemble package

You can install DAssemble directly from GitHub:

# install.packages("remotes")  # if not already installed
remotes::install_github("himelmallick/DAssemble")

After installation, load the package with:

library(DAssemble)

---

2. Install and load all method dependencies

DAssemble depends on packages from CRAN and Bioconductor for its core DA methods and enhancers. The script below will automatically install any missing packages and then load them.

## ===========================================
## Install + load all required DAssemble deps
## ===========================================

# All required packages
req_pkgs <- c(
  "MAST", "DESeq2", "edgeR", "limma", "metagenomeSeq",
  "dearseq", "SummarizedExperiment", "ALDEx2", "MicrobiomeStat",
  "LOCOM2", "Maaslin2", "maaslin3", "MultiAssayExperiment",
  "cplm", "dfadjust", "glmmTMB",
  "logging", "MASS", "pbapply", "preprocessCore",
  "ANCOMBC", "TreeSummarizedExperiment", "S4Vectors"
)

# Install required managers
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

# Install + load loop
for (pkg in req_pkgs) {

  # Install if missing
  if (!requireNamespace(pkg, quietly = TRUE)) {
    message("Installing missing package: ", pkg)

    tryCatch({
      install.packages(pkg)
    }, error = function(e) {
      message("  -> CRAN failed, installing via Bioconductor: ", pkg)
      BiocManager::install(pkg, ask = FALSE)
    })
  }

  # Load package
  message("Loading package: ", pkg)
  ok <- suppressPackageStartupMessages(
    require(pkg, character.only = TRUE)
  )
  if (!ok)
    stop("Failed to load package: ", pkg, call. = FALSE)
}

cat("All required packages installed and loaded successfully.\n")

---

Core Methods Supported

Method	Package(s)
DESeq2	DESeq2
edgeR	edgeR
limma-voom	limma + edgeR
metagenomeSeq	metagenomeSeq
MAST	MAST
dearseq	dearseq + SummarizedExperiment
ALDEx2	ALDEx2
LinDA	MicrobiomeStat
LOCOM	LOCOM2
Maaslin2	Maaslin2
Maaslin3	maaslin3
Tweedieverse	Internal DAssemble implementation using cplm/glmmTMB
Robseq	Internal DAssemble implementation using MASS/dfadjust
ANCOM-BC2	ANCOMBC + TreeSummarizedExperiment + S4Vectors

Note on singleton core methods

Some methods already combine two models internally (e.g., hurdle or two‑part models) and/or directly output a combined p‑value. These should normally be used alone in an ensemble, i.e. as the only core method. Combining a hurdle model with other core methods can double‑count the same underlying signal and lead to inflated type I error. The following methods fall into this category:

• Maaslin3 – fits both a prevalence and abundance model under the hood;
• MAST – implements a hurdle model for single‑cell expression data;
• metagenomeSeq – fits a zero‑inflated log‑normal model.

You may still add enhancers such as WLX or KS to these models, but LR is not recommended here because it can reintroduce evidence. Combining them with additional core methods is also not recommended.

---

Compatible combinations

DAssemble can be run in a number of configurations depending on whether you include a core method, enhancers, or both. The current framework allows you to combine one core method with up to three enhancers, or to run enhancer‑only analyses when core_method = NULL. The table below lists the valid combinations of enhancers by the number of enhancers used. The currently available enhancers are WLX (Wilcoxon rank‑sum), LR logistic regression) and KS (Kolmogorov–Smirnov).

# enhancers	With core	Example combinations
0	Yes	core_method = "DESeq2", enhancers = NULL
1	Yes / No	c("WLX"), c("LR"), c("KS")
2	Yes / No	c("WLX","LR"), c("WLX","KS"), c("LR","KS")
3	Yes / No	c("WLX","LR","KS")

For example, you might call:

# core + two enhancers
DAssemble(
  X, metadata,
  core_method = "DESeq2",
  enhancers = c("WLX", "LR"),
  expVar = "group",
  coVars = c("batch", "age")
)

# enhancer‑only, all three
DAssemble(
  X, metadata,
  core_method = NULL,
  enhancers = c("WLX", "LR", "KS"),
  expVar = "group"
)

When return_subensembles = TRUE the returned $ensembles element will include every sub‑combination of the requested methods. For instance, with three enhancers and no core, the sub‑ensembles will report CCT results for each single test ("WLX", "LR", "KS"), each pair ("WLX+LR", etc.) and the full trio ("WLX+LR+KS").

---

Enhancers

DAssemble currently implements three enhancer methods:

• DA_fit_enhancer_WLX() – Wilcoxon rank-sum test (WLX)
• DA_fit_enhancer_LR() – presence–absence logistic regression (LR)
• DA_fit_enhancer_KS() – Kolmogorov–Smirnov test (KS)

All enhancers return:

feature     # character
pval_<TAG>  # numeric, where <TAG> is WLX / LR / KS

Within DAssemble(), you usually specify enhancers using their short tags:

enhancers = c("WLX", "LR", "KS")

DAssemble then routes to the corresponding DA_fit_enhancer_*() functions.

---

DAssemble main function

The main entry point is:

DAssemble(
  features,
  metadata      = NULL,
  core_method   = NULL,
  enhancers     = NULL,
  expVar        = "group",
  coVars        = NULL,
  random_effects = NULL,
  method_args   = NULL,
  assay_name    = NULL,
  p_adj         = "BY",
  enhancer_norm = "TSS",
  return_components   = TRUE,
  return_subensembles = FALSE
)

Key arguments:

• features – a MultiAssayExperiment, or a data frame with samples in rows and features in columns
• metadata – sample metadata data frame; leave as NULL when features is a MultiAssayExperiment
• expVar – name (or column) of the primary exposure variable (binary, 2 levels)
• coVars – optional character vector of adjustment covariates for compatible methods
• random_effects – optional character vector of grouping variables for longitudinal-compatible methods
• method_args – optional named list of method-specific control arguments passed through to the selected core method and/or enhancers
• assay_name – experiment name to use when features is a MultiAssayExperiment; required when multiple experiments are present
• core_method – one of the supported core method names, or NULL / "none" to run an enhancer-only analysis
• enhancers – NULL or a subset of c("WLX", "LR", "KS")
• p_adj – multiple testing correction method (passed to p.adjust)
• enhancer_norm – normalization used by enhancers, one of "TSS", "CLR", "TMM", or "SCRAN"
• return_components – if TRUE, return per-method results in $components
• return_subensembles – if TRUE, compute CCT sub-ensembles and return them in $ensembles

The main output includes:

• $res – data frame of features ranked by joint CCT p-value and adjusted q-value
• $components – named list of individual core/enhancer result tables (returned by default)
• $ensembles – (optional) table of sub-ensemble CCT results

For example, with core_method = "DESeq2" and enhancers = c("WLX", "LR"), the individual results are available as result$components$DESeq2, result$components$WLX, and result$components$LR. These tables do not alter the joint ranking in result$res.

Passing through method-specific options

DAssemble() supports an open-ended method_args control list so users can pass additional arguments to underlying methods without the package needing to enumerate every possible option in the top-level API. When an argument overlaps with a DAssemble wrapper default, the value in method_args takes precedence.

The expected pattern is:

• method_args$core for arguments shared by the selected core method
• method_args$enhancer for arguments shared by the requested enhancers
• method_args$<MethodName> for method-specific overrides such as method_args$Maaslin3 or method_args$LR

For the LR enhancer specifically, the wrapper also recognizes method_args$LR$separation_method. The default is "augment", which follows the MaAsLin3-style augmented logistic fit. For cross-sectional presence/absence logistic models, users can instead request "firth" to use bias-reduced logistic regression via brglm2.

For the Maaslin2 core, the wrapper also recognizes method_args$Maaslin2$median_comparison = TRUE. When requested, DAssemble uses raw fit$results, applies the internal median_comparison_tweedie() adjustment there, and only then reduces the output back to the standardized exposure-specific result table.

For wrappers with multiple internal stages, method_args$<MethodName> can also target subcalls. For example, the DESeq2 wrapper can receive entries such as DESeq, results, or estimateSizeFactors, while the edgeR wrapper can receive entries such as glmFit or glmLRT.

Example:

res <- DAssemble(
  features = features,
  metadata = metadata,
  enhancers = "LR",
  expVar = "group",
  method_args = list(
    LR = list(separation_method = "firth")
  )
)

For example, MaAsLin3's prevalence median comparison can be enabled while retaining the other DAssemble defaults:

res <- DAssemble(
  features = features,
  metadata = metadata,
  core_method = "Maaslin3",
  expVar = "group",
  method_args = list(
    Maaslin3 = list(median_comparison_prevalence = TRUE)
  )
)

---

Sub-ensembles and core_method = NULL

If return_subensembles = TRUE and the helper DA_build_subensembles() is available, DAssemble computes sub-ensembles in addition to the main joint CCT p-values.

• When core_method is not NULL, sub-ensembles typically correspond to:

  • the core alone,
  • each core + single-enhancer combination,
  • and (optionally) the full core + all-enhancers ensemble.

• When core_method = NULL (enhancer-only mode), the sub-ensembles are defined purely in terms of the enhancers.
  For example, with enhancers = c("WLX", "LR", "KS"), the sub-ensembles can include:

  • WLX alone, LR alone, KS alone
  • WLX + LR, WLX + KS, LR + KS
  • WLX + LR + KS

In the resulting $ensembles object, the sub_ensemble (or similarly named) column labels each combination (e.g., "WLX", "WLX+LR", "LR+KS", "WLX+LR+KS").

This is useful for sensitivity analysis: you can see how conclusions change as you move from individual tests to different CCT combinations, including the enhancer-only setting when no core method is used.

---

Real data examples

To illustrate how to use DAssemble on real datasets, this section provides two short examples: one for a bulk RNA‑Seq experiment and one for a microbiome study. These examples rely on publicly available data from Bioconductor packages; citations are provided for a brief description of each dataset.

Bulk RNA‑Seq: Airway dexamethasone experiment

The airway dataset from the Bioconductor package airway is a small RNA‑Seq experiment in which four airway smooth muscle cell lines are treated with the asthma medication dexamethasone. It is often used as an introductory example for differential expression analysis; the dataset comprises counts for ~63 k genes measured in eight samples (4 × 2 design). A Bioconductor vignette notes that the airway dataset provides a typical small‑scale RNA‑Seq experiment where four ASM cell lines are treated with dexamethasone.

To run DAssemble on this dataset using the DESeq2 core and two enhancers (Wilcoxon and logistic regression):

library(DESeq2)
library(airway)
library(MultiAssayExperiment)
library(DAssemble)

# load counts and metadata
data("airway")
counts   <- assay(airway, "counts")
metadata <- as.data.frame(colData(airway))
metadata <- metadata[colnames(counts), , drop = FALSE]

mae <- MultiAssayExperiment(
  experiments = list(rnaseq = SummarizedExperiment::SummarizedExperiment(
    assays = list(counts = counts)
  )),
  colData = S4Vectors::DataFrame(metadata)
)

# the exposure variable is dex (treated vs untreated)
res <- DAssemble(
  features    = mae,
  assay_name  = "rnaseq",
  core_method = "DESeq2",
  enhancers   = c("WLX", "LR"),
  expVar      = "dex",
  coVars      = c("cell", "albut"),
  p_adj       = "BH",
  enhancer_norm = "tmm",
  method_args = list(
    DESeq2 = list(
      DESeq = list(fitType = "local")
    ),
    LR = list(
      control = glm.control(maxit = 100)
    )
  ),
  return_components = TRUE,
  return_subensembles = TRUE
)

# inspect the top differential genes
head(res$res)
# inspect per‑method p‑values
names(res$components)

This call normalizes the counts using the TMM scheme (enhancer_norm = "tmm"), adjusts for cell and albut, combines DESeq2 with the Wilcoxon and logistic regression enhancers, and returns the full set of sub‑ensembles. You can access individual method outputs through the $components element of the result. The method_args list illustrates how to forward extra options to the underlying DESeq2 and logistic-regression wrappers.

Microbiome: Global Patterns study

The Global Patterns dataset, available via the phyloseq package, contains 16S rRNA profiles for 25 environmental samples and three synthetic “mock communities,” representing nine sample types in total at an average sequencing depth of 3.1 million reads per sample【156744987103578†L128-L133】. It has been used to explore diversity patterns across a variety of ecosystems.

To run an enhancer‑only DAssemble analysis on a microbiome count table extracted from GlobalPatterns:

library(phyloseq)
library(MultiAssayExperiment)
library(DAssemble)

# load Global Patterns as a phyloseq object
data("GlobalPatterns")
gp <- GlobalPatterns

# extract count matrix and sample metadata
otu  <- as(otu_table(gp), "matrix")
meta <- as.data.frame(sample_data(gp))

# here we compare human fecal samples against soil samples as an example
keep <- meta$SampleType %in% c("Feces", "Soil")
X    <- t(otu[, keep])
meta <- meta[keep, , drop = FALSE]
meta$group <- droplevels(factor(meta$SampleType))
stopifnot(nlevels(meta$group) == 2L)

mae <- MultiAssayExperiment(
  experiments = list(microbiome = SummarizedExperiment::SummarizedExperiment(
    assays = list(counts = t(as.matrix(X)))
  )),
  colData = S4Vectors::DataFrame(meta)
)

res <- DAssemble(
  features    = mae,
  assay_name  = "microbiome",
  core_method = NULL,
  enhancers   = c("WLX", "LR", "KS"),
  expVar      = "group",
  enhancer_norm = "clr",
  return_components = TRUE,
  return_subensembles = TRUE
)

head(res$res)

In this example the core is set to NULL, so the analysis combines three nonparametric enhancers only. Because microbiome data are compositional, we use centered log‑ratio (CLR) normalization (enhancer_norm = "clr"). The $components element contains the individual Wilcoxon, logistic regression and KS results, while $ensembles summarises every sub‑combination of the three enhancers.

---

Covariate and longitudinal support

The main DAssemble() entry point now supports both multiple covariates (coVars) and repeated-measures / longitudinal designs (random_effects) for the subset of methods whose wrappers can express those models directly.

Multiple covariates

The following methods support coVars in the current implementation:

Method	Multiple covariates
DESeq2	Yes
edgeR	Yes
limmaVOOM	Yes
metagenomeSeq	Yes
MAST	Yes
dearseq	Yes
ALDEx2	Yes
LinDA	Yes
Maaslin2	Yes
Maaslin3	Yes
Tweedieverse	Yes
Robseq	Yes
ANCOMBC2	Yes
LR	Yes
WLX	No
KS	No
LOCOM	No

Longitudinal / repeated-measures

When random_effects is supplied, DAssemble() currently supports the following longitudinal-compatible methods:

Method	Longitudinal support
Maaslin2	Yes
Maaslin3	Yes
Tweedieverse	Yes
LR	Yes

All other current methods are treated as non-longitudinal in DAssemble(). If random_effects is provided with an unsupported core or enhancer, DAssemble will stop with a clear error instead of fitting a mis-specified model.

For the LR enhancer, the default fitting path uses augmentation in both cross-sectional and longitudinal settings. The alternative method_args$LR$separation_method = "firth" is available only for the non-longitudinal model.

Longitudinal example

Maaslin3 is treated through its standard package pathway in DAssemble().

res <- DAssemble(
  features = features,
  metadata = metadata,
  core_method = "Maaslin3",
  enhancers = "LR",
  expVar = "group",
  coVars = c("age", "sex", "batch"),
  random_effects = "subject",
  method_args = list(
    Maaslin3 = list(max_significance = 0.25),
    LR = list(control = glm.control(maxit = 100))
  ),
  p_adj = "BH"
)