manuscript.qmd

---
title: Orchestrating Microbiome Analysis with Bioconductor
format:
  nature-pdf:
    # journal.cite-style is included in the tex file but ignored by pandoc if
    # cite-method is not `natbib`.
    journal:
      cite-style: sn-nature
    # `citeproc` is the pandoc default. Set `cite-method: natbib` if required
    # to use the bst styles from the upstream template.
    cite-method: natbib
    keep-tex: true
    equal-margins: false
    # These settings below enable word calculation
    filters:
      - at: pre-quarto
        path: _extensions/andrewheiss/wordcount/citeproc.lua
      - at: pre-quarto
        path: _extensions/andrewheiss/wordcount/wordcount.lua
affiliations:
  - id: 1
    name: University of Turku
    department: Department of Computing
    address: Turku, Finland

  - id: 2
    name: University of Turku
    department: Department of Life Technologies
    address: Turku, Finland

  - id: 3
    name: University of Oxford
    department: Department of Biology
    address: Oxford, United Kingdom

  - id: 4
    name: University of Oxford
    department: Queen’s College
    address: Oxford, United Kingdom

  - id: 5
    name: University College Cork
    department: Department of Anatomy and Neuroscience
    address: Cork, Ireland

  - id: 6
    name: University College Cork
    department: APC Microbiome Ireland
    address: Cork, Ireland

  - id: 7
    name: LMU München
    department: Department of Statistics
    address: Munich, Germany

  - id: 8
    name: Cornell University
    department: Department of Population Health Sciences, Weill Cornell Medicine
    address: New York, USA

  - id: 9
    name: Cornell University
    department: Department of Biomedical Informatics & Data Science
    address: New York, USA

  - id: 10
    name: Cornell University
    department: Department of Biostatistics
    address: New York, USA

  - id: 11
    name: Yale University
    department: Department of Biomedical Informatics & Data Science
    address: New Haven, USA

  - id: 12
    name: Yale University
    department: Department of Biostatistics
    address: New Haven, USA

  - id: 13
    name: University of Helsinki
    department: Department of Microbiology
    address: Helsinki, Finland

  - id: 14
    name: Helmholtz Munich
    address: Munich, Germany

  - id: 15
    name: Flatiron Institute
    address: New York, USA

  - id: 16
    name: Institute for Molecular Medicine Finland, FIMM-HiLIFE
    address: Helsinki, Finland

  - id: 17
    name: Danone Research and Innovation
    address: Utrecht, Netherlands

  - id: 18
    name: CUNY Graduate School of Public Health and Health Policy
    department: Department of Epidemiology and Biostatistics
    address: New York, USA

  - id: 19
    name: No affiliation

  - id: 20
    name: Amsterdam University Medical Centres
    department: Department of Psychiatry
    address: Amsterdam, Netherlands

  - id: 21
    name: Amsterdam University Medical Centres
    department: Amsterdam Public Health
    address: Amsterdam, Netherlands

  - id: 22
    name: Cornell University
    department: Division of Gastroenterology and Hepatology, Weill Cornell Medicine
    address: New York, USA

  - id: 23
    name: Cornell University
    department: Department of Statistics and Data Science
    address: New York, USA

author:
  - name: Tuomas Borman
    orcid: 0000-0002-8563-8884
    email: tuomas.v.borman@utu.fi
    affiliations:
      - ref: 1
    attributes:
      corresponding: true

  - name: Giulio Benedetti
    orcid: 0000-0002-8732-7692
    affiliations:
      - ref: 1
    attributes:
      equal-contributor: true

  - name: Geraldson Muluh
    orcid: 0000-0002-9721-5192
    affiliations:
      - ref: 1
    attributes:
      equal-contributor: true

  - name: Aura Raulo
    orcid: 0000-0003-4860-7840
    affiliations:
      - ref: 3
      - ref: 4

  - name: Benjamin Valderrama
    orcid: 0000-0002-4740-8965
    affiliations:
      - ref: 5
      - ref: 6

  - name: Artur Sannikov
    orcid: 0000-0001-7765-123X
    affiliations:
      - ref: 2

  - name: Stefanie Peschel
    orcid: 0000-0002-7936-7627
    affiliations:
      - ref: 7

  - name: Yihan Liu
    orcid: 0009-0008-0462-7501
    affiliations:
      - ref: 8
      - ref: 9
      - ref: 10
      - ref: 11
      - ref: 12

  - name: Rasmus Hindström
    orcid: 0009-0004-5731-178X
    affiliations:
      - ref: 1

  - name: OMA consortium
    affiliations:
      - ref: 19

  - name: Katariina Pärnänen
    orcid: 0000-0003-2220-1738
    affiliations:
      - ref: 13

  - name: Christian L. Müller
    orcid: 0000-0002-3821-7083
    affiliations:
      - ref: 7
      - ref: 14
      - ref: 15

  - name: Aki S. Havulinna
    orcid: 0000-0002-4787-8959
    affiliations:
      - ref: 1
      - ref: 16

  - name: Sudarshan Shetty
    orcid: 0000-0001-7280-9915
    affiliations:
      - ref: 17

  - name: Marcel Ramos
    orcid: 0000-0002-3242-0582
    affiliations:
      - ref: 18

  - name: Domenick J. Braccia
    orcid: 0000-0002-4606-6254
    affiliations:
      - ref: 19

  - name: Héctor Corrada Bravo
    orcid: 0000-0002-1255-4444
    affiliations:
      - ref: 19

  - name: Felix M. Ernst
    orcid: 0000-0001-5064-0928
    affiliations:
      - ref: 19

  - name: Levi Waldron
    orcid: 0000-0003-2725-0694
    affiliations:
      - ref: 18

  - name: Thomaz F. S. Bastiaanssen
    orcid: 0000-0001-6891-734X
    affiliations:
      - ref: 20
      - ref: 21
    attributes:
      equal-contributor: true

  - name: Himel Mallick
    orcid: 0000-0003-4956-2429
    affiliations:
      - ref: 8
      - ref: 22
      - ref: 23
    attributes:
      equal-contributor: true

  - name: Leo Lahti
    orcid: 0000-0001-5537-637X
    email: leo.lahti@utu.fi
    affiliations:
      - ref: 1
    attributes:
      corresponding: true
abstract: |
  The expansion of microbiome research has led to the accumulation of
  interlinked datasets encompassing versatile taxonomic and functional
  assays. The analysis of increasingly large and heterogeneous multi-modal
  microbiome data would benefit from unified approaches supporting the design
  of modular data science workflows through interoperable methods. Building upon
  the optimized statistical programming framework for multi-assay data
  integration recently developed by the Bioconductor project,
  we introduce a community-developed open source framework
  for microbiome data science. In contrast to the previous alternatives,
  the methodology is specifically designed to support joint analysis
  of hierarchical, interlinked, and heterogeneous multi-table datasets that
  are increasingly common in modern microbiome research. This
  framework encompasses open data, methods, tutorials, and an active online
  community, thereby supporting standardized and reproducible data
  wrangling, joint analysis, and reporting. We have detailed the functionality
  and usage in the OMA book
  ([https://bioconductor.org/books/release/OMA](https://bioconductor.org/books/release/OMA)),
  which offers guidance for prospective users and contributors.

keywords: [microbiome, Bioconductor, data science, data integration]
bibliography: bibliography.bib
---

# Introduction {#sec-intro}

Modern data science applications critically rely on open-source methods created
by the research community [@barker2022; @pearson2025]. Open
software ecosystems have become central innovation hubs, driving
technological and cultural transformation [@hocquet2024]. The
Bioconductor project has become a significant distribution
channel for open data science methods in the life sciences.
Its vast developer community maintains over 2,300 quality-controlled R
packages for various fields of life science informatics [@gentleman2004].
Additionally, Bioconductor has evolved into a global community that supports
data science collaboration and training [@drnevich2025], playing an increasingly
important role in advancing computational skills that are essential for modern
microbiology education [@timmis2024].

Microbiome research focuses on how microbial communities vary, function and
interact with their environment and a potential host [@moreno-indias2021]. The
technologies used in this field generate vast amounts of high-dimensional and
hierarchically structured data, commonly in the form of microbial DNA sequences
detected in a set of samples. Following the rapid expansion of microbiome
research and technological advances, metagenomic sequencing is increasingly
complemented by other types of high-throughput measurements. This has
led to the accumulation of multi-modal datasets—collections
that combine different types of molecular data. These interlinked datasets
encompass diverse taxonomic and functional assays, ranging from microbial
abundances and phylogenetic relationships to functional profiles.

Microbiome data scientists
are thus increasingly facing the practical challenges of integrating
heterogeneous data across multiple data modes into reproducible workflows.
The unique statistical properties and complexity of microbiome data
necessitate specialized approaches in their analysis
[@moreno-indias2021; @jeganathan2021; @willis2019], yet the diversity of
proposed solutions and inconsistent outcomes
[@nearing2022; @pelto2025; @gamboa-tuz2025] can make it
difficult to identify the optimal methods and build interoperable workflows
[@wen2023; @shetty2019]. Accordingly, a need for standardization has been
highlighted [@moreno-indias2021].
Open microbiome data science _frameworks_ have emerged to provide sets of
interoperable methodologies in specific computational environments
[@eren2021; @schloss2009; @bolyen2019; @mcmurdie2012], and a broader ecosystem
of digital resources, engaged user communities and collaboration culture that
extend the capacities far beyond the underlying technical framework [@katz2019].
However, frameworks focused solely on microbiome data science fall short on
addressing the need to integrate data across multiple modalities, such as
metagenomics, metabolomics, and host genomics.
The Bioconductor project integrates the methods of microbiome data science
into a wider statistical programming framework, where
interoperable methods and
workflows for single-cell omics [@amezquita2020], spatial transcriptomics [@crowell2025],
metabolomics [@gatto2025], proteomics [@gatto2025; @anwar2025], and other
fields are concurrently developed.

Here, we introduce a stable and optimized data science ecosystem supporting the
integration and analysis of multi-modal datasets in microbiome research.
It provides harmonized data import and representation, and fast,
robust methods for transforming abundance tables into biological insights. These
methods have been extensively tested and distributed through R/Bioconductor
packages
([https://bioconductor.org/packages](https://bioconductor.org/packages)),
supporting the community-driven development of the framework and integration
within the broader Bioconductor ecosystem. To promote
evidence-based best practices, we introduce the online book _Orchestrating
Microbiome Analysis with Bioconductor_ (OMA), which represents a comprehensive
compilation of best practices and current trends in microbiome data science using the
R/Bioconductor ecosystem.

# Data analytical frameworks {#sec-data_frameworks}

Several data analytical frameworks are available for microbiome data science,
each distinguished by its underlying paradigms, user communities, and scope.
For example, tools such as QIIME2 [@bolyen2019], Anvi'o
[@eren2021], and Mothur [@schloss2009] emphasize standardized, modular
approaches tailored to microbiome analysis. These frameworks provide
validated and user-friendly solutions, enabling development of efficient
end-to-end workflows for microbiome research. However, they are typically
designed around specific data types, which limits flexibility
when integrating heterogeneous data modalities within a unified analytical
framework. In contrast, Bioconductor embraces a more exploratory, statistically
oriented paradigm spanning multiple areas of bioinformatics beyond
microbiome research.

This data analytical paradigm of Bioconductor is similar to scikit-bio's
[@aton2025]. While scikit-bio as a Python-based ecosystem benefits from rapidly
evolving machine learning and artificial intelligence tools, Bioconductor’s
strengths lie in its vast statistical ecosystem dedicated to omics research and
strong interoperability with other computational environments.
Another key distinction is community structure: scikit-bio is more
centralized, whereas Bioconductor packages are developed and maintained by an
open community, resulting in a broad toolkit and collaboration enabling easy
reuse of common solutions.
Although this approach can introduce some challenges, such as dependency management,
package compatibility, and differences in package stability, these are actively
addressed through several coordinated mechanisms. First, submitted packages undergo
peer review, ensuring relevance and adherence to standards such as documentation
and working examples. Secondly, Bioconductor follows synchronized, semiannual
release cycle, during which packages are frozen and mutually compatible. Thirdly,
all packages are continuously tested within Bioconductor's automated build system,
which ensures that breaking changes are rapidly detected before release. Each package has an active
maintainer who is notified of failures and resolves them in a timely manner.
Finally, community guidelines and continuous integration practices, together
with shared data structures such as (Tree)SummarizedExperiment, further improve
interoperability between packages.

# Data infrastructure {#sec-infra}

Data scientists routinely deal with interlinked assays, hierarchical data
(*e.g.*, ontologies, trees), relational data (*e.g.*, networks), and supporting
side information. This side information can include
sample attributes such as treatment group, collection site, or time point. The
ability to integrate heterogeneous data elements is essential in biological
research, where the individual components of a system can seldom be studied in
isolation. However, this growing complexity introduces additional overhead, as
researchers must track samples and features across multiple tables and manage an
increasing number of data elements and non-trivial links between them. Moreover,
input data formats vary across different methods, leading to technical
complications and time lost on data wrangling. Bioconductor is addressing these
challenges through standardized data structures that are supported by an
ecosystem of quality-controlled, interoperable methods. This allows the analyst
to invest more time in core analytical tasks.

Bioconductor data structures support statistical analysis of complex
data combinations [@huber2015] and are widely adopted by various
research communities. They have been implemented through specialized
object-oriented classes that integrate multiple data elements into
a single structure known as a *data container*
[@chambers2014]. The SummarizedExperiment (SE) class is
the primary data container underlying the Bioconductor ecosystem
[@huber2015; @morgan2025]. Its position among the top 1\% most-downloaded
packages in Bioconductor reflects its widespread recognition and
increasing adoption among developers (@fig-bioc_packages). SE provides a
standardized solution to store tabular data by linking numeric abundance
matrices with side information on features or rows  (*e.g.*, taxonomy),
and samples or columns (*e.g.*, origin, collection time, host health).

The interoperability of SE with other Bioconductor data structures
enables the transfer of methods across different fields of omics research,
thereby promoting the development of modular and scalable data science methods
and ultimately yielding improved overall quality control and user support [@drnevich2025].
In particular, the SE framework has been extended to
similar data integration tasks, notably in single-cell
and spatial omics through the SingleCellExperiment (SCE) [@amezquita2020] and
SpatialExperiment [@crowell2025] classes.

```{r}
#| label: bioc_packages
#| eval: false
#| echo: false
#| message: false

if (!require("BiocManager")) {
    install("BiocManager")
    library("BiocManager")
}
pkgs <- c("BiocPkgTools", "lubridate", "tidyverse")
temp <- sapply(pkgs, function(pkg) {
    if (!require(pkg, character.only = TRUE)) {
        install(pkg)
        library(pkg, character.only = TRUE)
    }
})

# Get packages and info based on their BiocViews
pkgs <- biocPkgList()
# Get how many times these packages are loaded
df <- biocDownloadStats()

# Take only packages that depend on SE
pkgs[["SE"]] <- sapply(pkgs[["Depends"]], function(x) any(grepl("summarizedexperiment", x, ignore.case = TRUE)))
pkgs <- pkgs[pkgs[["SE"]], ]

# Extract the package topic based on biocView
fields <- sapply(pkgs[["biocViews"]], function(x) {
    field <- "Other"
    if (any(grepl("genomic|genetic|gene|genom|DNA", x, ignore.case = TRUE))) field <- "Genomics"
    if (any(grepl("proteomics|protein", x, ignore.case = TRUE))) field <- "Proteomics"
    if (any(grepl("metabolomics|metabolome|metabolite|Lipidom|massspectro", x, ignore.case = TRUE))) field <- "Metabolomics"
    if (any(grepl("transcripto|RNA-seq|RNA", x, ignore.case = TRUE))) field <- "Transcriptomics"
    if (any(grepl("immuno", x, ignore.case = TRUE))) field <- "Immunology"
    if (any(grepl("cytom", x, ignore.case = TRUE))) field <- "Cytometrics"
    if (any(grepl("microarray|chip", x, ignore.case = TRUE))) field <- "Microarray"
    if (any(grepl("single-cell|singlecell", x, ignore.case = TRUE))) field <- "Single-cell"
    if (any(grepl("metagenom|microbiome|16S|microbiota|amplicon|shotgun|microb|metatranscript|metametabolo|metaproteo", x, ignore.case = TRUE))) field <- "Microbiome"
    return(field)
})
pkgs[["Field"]] <- fields

# Subset downloads table and add field to it
df <- df[df[["Package"]] %in% pkgs[["Package"]], ]
df[["Field"]] <- pkgs[match(df[["Package"]], pkgs[["Package"]]), "Field"][[1]]

# Get information on when the package became available in Bioconductor
dates <- df[["Date"]] |> unique()
dates <- dates[dates < floor_date(Sys.time(), "month")]
available <- lapply(dates, function(date) {
    temp <- df[df[["Date"]] == date, ]
    temp <- temp[match(unique(df[["Package"]]), temp[["Package"]]), ]
    temp <- temp[["Nb_of_distinct_IPs"]] > 0 | temp[["Nb_of_downloads"]] > 0
    temp[is.na(temp)] <- FALSE
    as.numeric(temp)
})
available <- do.call(cbind, available) |> as.data.frame()
colnames(available) <- dates
rownames(available) <- unique(df[["Package"]])

# Create a table showing the number of packages by field through time
ind <- pkgs[match(unique(df[["Package"]]), pkgs[["Package"]]), "Field"][[1]]
pkgs_date <- rowsum(available, group = ind)
pkgs_date <- pkgs_date |>
    rownames_to_column("Field") |>
    pivot_longer(
        cols = -Field,
        names_to = "Date",
        values_to = "N"
    )
pkgs_date[["Date"]] <- as.Date(pkgs_date[["Date"]])

# Save results
pkgs_date |>
    write.csv(file = file.path("data", "bioc_packages_stats.csv"), row.names = FALSE)
```

```{r}
#| label: fig-bioc_packages
#| eval: true
#| echo: false
#| message: false
#| warning: false
#| fig-height: 3
#| fig-width: 6
#| fig-cap: "**Adoption of the SummarizedExperiment (SE) data container among
#|   Bioconductor developers.** The growth in the number of Bioconductor
#|   packages supporting SE over time is shown, categorized by their respective
#|   fields. Packages are counted by the date they were added to Bioconductor.
#|   Some packages added the SE support after their original release."

if (!require("BiocManager")) {
    install("BiocManager")
    library("BiocManager")
}
pkgs <- c("tidyverse", "see")
temp <- sapply(pkgs, function(pkg) {
    if (!require(pkg, character.only = TRUE)) {
        install(pkg)
        library(pkg, character.only = TRUE)
    }
})

# Import data
pkgs_date <- read.csv(file.path("data", "bioc_packages_stats.csv"))
pkgs_date[["Date"]] <- pkgs_date[["Date"]] |> as.Date()

# Create the plot
p <- ggplot(pkgs_date, aes(x = Date, y = N, fill = Field)) +
    geom_area() +
    theme_classic() +
    scale_fill_okabeito() +
    labs(x = "Year", y = "Number of packages") +
    scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
    scale_x_date(expand = c(0, 0), limits = c(as.Date("2015-01-01"), max(pkgs_date[["Date"]])))

p
```

**Hierarchical data structures.** Microbiome data scientists frequently need to
integrate information on feature phylogenies and sample hierarchies
to accurately reflect fine-scale variability in microbiome data.
While maintaining interoperability with
the broader Bioconductor ecosystem, the TreeSummarizedExperiment
(TreeSE) class [@huang2021] extends the SE and SCE data structures to
hierarchical datasets by incorporating both feature and sample organizations
as tree structures (@fig-datacontainers a). Moreover, TreeSE supports the
storage of other common elements of microbiome data, for example, reference
sequences, and incorporates results from common analytical procedures such as
dimensionality reduction, statistical transformations, and agglomeration by
taxonomic, functional or other information. The interoperability with other
Bioconductor classes facilitates scalable analysis and the design of modular
yet flexible workflows, as TreeSE inherits full compatibility with methods
designed for SE and SCE. In particular, the direct interoperability with the
widely adopted SE data science ecosystem distinguishes our TreeSE
framework from phyloseq, which is another popular Bioconductor data container
developed initially with a focus on 16S amplicon data analysis [@mcmurdie2012].
TreeSE provides a more general framework than phyloseq, as the former not only
accommodates the typical components of a phyloseq object
(taxa abundances, sample metadata, phylogenetic tree, and reference sequences),
but also extends to multi-modal datasets encompassing
parallel taxonomic levels, functional predictions, and resistome, metabolome or
other profiles and alternative feature representations. It also offers
substantial gains in computing speed and memory efficiency over phyloseq
(@fig-benchmarking). While speedyseq is faster in certain operations, TreeSE
tends to scale more efficiently in terms of execution time (t) and allocated
memory (m) as sample size (n) increases. This benefits from the built-in support
for sparse and delayed matrices in SE-based data containers, which can improve
performance [@huber2015]. Despite these advantages, further optimizing time and
memory consumption represents a key area for further development in the (Tree)SE
framework. To facilitate migration from phyloseq, the OMA book provides
conversion functions between the two data container formats together with a
mapping table of commonly used functions. These resources enable users to
incrementally migrate existing phyloseq workflows to the (Tree)SE framework while
preserving familiar coding patterns and minimizing changes to established
analysis pipelines. Despite the advantages of the (Tree)SE framework,
further optimization of time and memory consumption remains an area for future
development.

```{r}
#| label: benchmark_analysis
#| eval: false
#| echo: false
#| message: false

# Adapted from https://github.com/microbiome/benchmarking/tree/main/article

# Import libraries
if (!require("BiocManager")) {
    install("BiocManager")
    library("BiocManager")
}
pkgs <- c(
    "bench", "DelayedArray", "mia", "microbiome", "microbiomeDataSets", "philr",
    "phyloseq", "picante", "tidyverse", "speedyseq"
)
temp <- sapply(pkgs, function(pkg) {
    if (!require(pkg, character.only = TRUE)) {
        install(pkg)
        library(pkg, character.only = TRUE)
    }
})

# Set seed for reproducibility
set.seed(123)

# Set benchmarking hyperparameters
n_iter <- 10
memory_threshold <- 10000
beta_threshold <- c(1000, 10000)
N <- c(10, 100, 1000, 10000)

# Define class names
classes <- c(tse = "TreeSE", pseq = "phyloseq", spseq = "speedyseq")

# Define method names
methods <- c(
    alpha = "Faith diversity",
    beta = "UniFrac dissimilarity",
    melt = "Melting",
    trans = "PhILR transformation",
    agg = "Family agglomeration"
)

# Import dataset
GTSD <- GrieneisenTSData()
# Agglomerate by Genus to reduce size
GTSD <- mia::agglomerateByRank(GTSD, rank = "Genus")

# Run benchmark for each sample size
benchmark_out <- bench::press(
    N = N,
    {
        # Select a random subset of samples
        tse <- GTSD[, sample(ncol(GTSD), N)]
        # Convert TreeSE to phyloseq
        pseq <- mia::convertToPhyloseq(tse)

        # List expressions to benchmark
        expressions <- list(
            # Estimate faith from phyloseq
            alpha_pseq = quote(picante::pd(samp = t(otu_table(pseq)), tree = phy_tree(pseq))[, 1]),
            # Estimate faith from TreeSE
            alpha_tse = quote(mia::getAlpha(tse, index = "faith_diversity")),
            # philr transform phyloseq
            trans_pseq = quote(philr::philr(t(otu_table(pseq)), tree = phy_tree(pseq), pseudocount = 1)),
            # philr transform TreeSE
            trans_tse = quote(mia::transformAssay(tse, method = "philr", MARGIN = 1L, pseudocount = 1)),
            # Melt phyloseq
            melt_pseq = quote(phyloseq::psmelt(pseq)),
            # Melt speedyseq
            melt_spseq = quote(speedyseq::psmelt(pseq)),
            # Melt TreeSE
            melt_tse = quote(mia::meltSE(tse, add.row = TRUE, add.col = TRUE)),
            # Agglomerate phyloseq
            agg_pseq = quote(phyloseq::tax_glom(pseq, taxrank = "Family")),
            # Agglomerate speedyseq
            agg_spseq = quote(speedyseq::tax_glom(pseq, taxrank = "Family")),
            # Agglomerate TreeSE
            agg_tse = quote(mia::agglomerateByRank(tse, rank = "Family"))
        )

        if (N <= beta_threshold[[1]]) {
            # Estimate unifrac from phyloseq
            expressions[["beta_pseq"]] <- quote(
                phyloseq::UniFrac(pseq)
            )
        } else {
            # Use dummy expression because all results must have the same length
            expressions[["beta_pseq"]] <- quote(1 + 1)
        }

        if (N <= beta_threshold[[2]]) {
            # Estimate unifrac from TreeSE
            expressions[["beta_tse"]] <- quote(
                mia::getDissimilarity(tse, method = "unifrac")
            )
        } else {
            # Use dummy expression because all results must have the same length
            expressions[["beta_tse"]] <- quote(1 + 1)
        }

        # Run benchmark
        bench::mark(
            iterations = n_iter,
            memory = if (N <= memory_threshold) TRUE else FALSE,
            check = FALSE,
            exprs = expressions
        )
    }
)

# Retrieve benchmarking results for each experiment and iteration
benchmark_df <- benchmark_out |>
    unnest(c(time, gc)) |>
    dplyr::select(expression, N, time, gc, mem_alloc) |>
    separate_wider_delim(
        cols = expression, delim = "_",
        names = c("method", "object")
    ) |>
    filter(method != "beta" |
        (object == "pseq" & N <= beta_threshold[[1]]) |
        (object == "tse" & N <= beta_threshold[[2]]))

# Summarise benchmarking results with mean time and standard deviation
benchmark_df <- benchmark_df |>
    group_by(method, object, N) |>
    summarise(
        Time = as.numeric(mean(time)), Memory = as.numeric(mean(mem_alloc)),
        TimeSD = sd(time), TimeSE = TimeSD / sqrt(n_iter),
        NoGC = sum(gc == "none"), .groups = "drop"
    ) |>
    mutate(
        method = factor(method, levels = names(methods)),
        object = factor(object, levels = names(classes))
    )

# Write to file
benchmark_df |>
    mutate(Time = as.numeric(Time), Memory = as.numeric(Memory)) |>
    write.csv(file = file.path("data", "benchmark_results.csv"), row.names = FALSE)
```

```{r}
#| label: fig-benchmarking
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "**Execution time and memory consumption for microbiome data
#|   containers.** The Bioconductor's TreeSE, phyloseq and speedyseq are
#|   commonly used data containers in microbiome analysis. We benchmarked them
#|   for five common operations applied to random subsets of samples from a
#|   large study on wild baboons [@grieneisen2021] (accessible through the
#|   microbiomeDataSets package in Bioconductor). Execution time (s) was
#|   measured as the mean value over 10 iterations, whereas allocated memory
#|   (MB) was estimated based on a single iteration. The benchmarking was
#|   conducted in R using the bench library with 8 CPU cores and 32 GB RAM.
#|   All operations, except for UniFrac dissimilarity with phyloseq, could
#|   also be performed on a regular laptop (*e.g.*, 4 CPU cores and 16 GB RAM).
#|   The source code for the benchmarking experiments is available through the
#|   OMA book."
#| fig-width: 12
#| fig-height: 6

# Adapted from https://github.com/microbiome/benchmarking/tree/main/article

pkgs <- c("patchwork", "tidyverse")
temp <- sapply(pkgs, function(pkg) {
    if (!require(pkg, character.only = TRUE)) {
        install(pkg)
        library(pkg, character.only = TRUE)
    }
})

classes <- c(tse = "TreeSE", pseq = "phyloseq", spseq = "speedyseq")
methods <- c(
    alpha = "Faith diversity",
    beta = "UniFrac dissimilarity",
    melt = "Melting",
    trans = "PhILR transformation",
    agg = "Family agglomeration"
)

# Import results
benchmark_df <- read.csv(file.path("data", "benchmark_results.csv")) %>%
    mutate(method = factor(method, levels = names(methods)),
           object = factor(object, levels = names(classes)))

# Select only points from 10, 100, 1000, 10000
benchmark_df <- benchmark_df[benchmark_df[["N"]] %in% c(10, 100, 1000, 10000), ]

# Specify plot layouts
scientific_10 <- function(y) {
    sapply(y, function(z) {
        if (is.na(z)) {
            return(NA)
        } else if (z == 1) {
            return("1")
        } else if (z == 10) {
            return("10")
        } else {
            return(parse(text = paste0("10^", log10(z))))
        }
    })
}
breaks <- benchmark_df[["N"]][log10(benchmark_df[["N"]]) %% 1 == 0] |> unique()

# Visualise benchmarking results: time
p1 <- ggplot(benchmark_df, aes(x = N, y = Time, colour = object)) +
    geom_errorbar(aes(ymin = Time - TimeSE, ymax = Time + TimeSE), width = 0) +
    geom_line() +
    geom_point() +
    scale_x_log10(breaks = breaks, limits = benchmark_df[["N"]] |> range(), labels = scientific_10) +
    scale_y_log10(labels = scientific_10) +
    scale_colour_manual(
        labels = classes,
        values = c("black", "darkgrey", "lightgrey")
    ) +
    facet_grid(. ~ method, labeller = labeller(method = methods)) +
    labs(x = "# Samples", y = "Execution time (s)", colour = "Object") +
    theme_bw() +
    theme(
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 15),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.ticks.x = element_blank(),
        strip.text = element_text(size = 12),
        strip.background = element_blank()
    )
# Visualise benchmarking results: memory
p2 <- ggplot(benchmark_df, aes(x = N, y = Memory / 1e6, colour = object)) +
    geom_line() +
    geom_point() +
    scale_x_log10(breaks = breaks, limits = benchmark_df[["N"]] |> range(), labels = scientific_10) +
    scale_y_log10(labels = scientific_10) +
    scale_colour_manual(
        labels = classes,
        values = c("black", "darkgrey", "lightgrey")
    ) +
    facet_grid(. ~ method,
        labeller = labeller(method = methods)
    ) +
    labs(x = "Sample size (n)", y = "Allocated memory (MB)", colour = "Object") +
    theme_bw() +
    theme(
        axis.title = element_text(size = 15),
        axis.text = element_text(size = 12),
        strip.text = element_blank()
    ) +
    guides(colour = "none")
# Combine results
p <- (p1 / p2) +
    plot_layout(guides = "collect") &
    theme(
        legend.position = "bottom",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 15),
        legend.key.size = unit(1.2, "cm")
    )

p
```

**Multi-assay data structures.** Contemporary microbiome research frequently
involves data across multiple measurement modalities for enhanced functional
and mechanistic insights. This brings the need to find an appropriate
representation for each data type and map the links between them. Whereas the
TreeSE data container is broadly applicable to different modalities such as
taxonomic, resistome, transcriptome, proteome, and metabolome profiling data,
dedicated data containers have been developed for certain modalities by the
Bioconductor community to address their specific characteristics; examples
include single cell sequencing [@amezquita2020] and host genomics
[@lawrence2013]. Linking data across modalities presents an additional technical
challenge, as multi-omics datasets are often sparse and only a subset of the
samples may be shared across multiple modalities [@sherwani2025]. In other
cases, a single sample in one modality may link to several samples in
another modality [@husso2023] (*e.g.*, due to technical or spatial replication).
These scenarios require the ability to integrate multiple
data modalities into a single object and create flexible sample mappings across them.
The MultiAssayExperiment (MAE) data container [@ramos2017] addresses these
challenges by introducing a formal mapping layer (sampleMap) that links
biological samples to their corresponding measurements in each modality
(@fig-datacontainers b). This design supports one-to-one, one-to-many, and
partially overlapping sample relationships across experiments. Importantly, this
mapping is also accessible for downstream operations. For instance when
subsetting a MAE object
(*e.g.*, selecting samples with both microbiome and metabolome data), the
framework automatically resolves the intersection of available samples and
ensures consistent alignment across all modalities. Such structured handling of
sample relationships is particularly important for integrative methods, as it
reduces the risk of sample mismatches, a common source of irreproducibility in
ad hoc multi-table analyses. For instance, diagonal integration methods
[@xu2022; @argelaguet2021computational]&mdash;where shared features across
data types may not be explicitly defined&mdash;can benefit from
such a systematic framework that enforces anchor alignment. This reflects lessons
learned from single-cell data integration, where shared anchors and hierarchical
mapping have proven essential for interpretability. To summarize, flexible mapping of samples
across several data objects is key for efficient data integration, method sharing,
and overall interoperability.

@tbl-containers summarizes the integrative data containers for microbiome
analysis, whereas @tbl-elements introduces the available data elements in each
container.

| Data container                                | Target use                                                  |
|-----------------------------------------------|-------------------------------------------------------------|
| SummarizedExperiment (SE)[@huber2015]         | generic container linking multiple data tables (assays)     |
| TreeSummarizedExperiment (TreeSE)[@huang2021] | extension incorporating hierarchical data structures        |
| MultiAssayExperiment (MAE)[@ramos2017]        | binding omic-specific containers in multi-omics experiments |

: Integrative data containers for microbiome analysis.
{#tbl-containers tbl-colwidths="[50,50]"}

| Slot          | Content                                      | SE | TreeSE | MAE$^a$ |
|---------------|----------------------------------------------|:--:|:------:|:-------:|
| **Tables**    |                                              |    |        |         |
| assays        | tables of features (rows) and samples (cols) | X  | X      | X       |
| altExps$^b$   | experiments with variable number of features |    | X      | X       |
| reducedDims   | results of dimensionality reduction          |    | X      | X       |
| **Rows**      |                                              |    |        |         |
| rowData       | side information on features                 | X  | X      | X       |
| rowPairs      | linkage between pairs of associated features |    | X      | X       |
| rowTrees      | hierarchical organization of features        |    | X      | X       |
| referenceSeq  | reference sequences of features              |    | X      | X       |
| **Columns**   |                                              |    |        |         |
| colData       | side information on samples                  | X  | X      | X       |
| colTrees      | hierarchical organization of samples         |    | X      | X       |
| sampleMap$^c$ | flexible sample mapping across experiments   |    |        | X       |
| **Other**     |                                              |    |        |         |
| metadata      | additional information on the experiment     | X  | X      | X       |

: Key data elements in the integrative data containers.
{#tbl-elements tbl-colwidths="[17,63,3,9,8]"}

$^a$ While MAE does not contain conventional slots (with the exception
of colData), its experiments can include any number of data containers that
provide those slots. $^b$ altExp samples must match assay samples in a 1-to-1
mapping. $^c$ The sampleMap supports diverse mappings between samples from
different experiments (1-to-1, 1-to-many, many-to-1 and many-to-many).

![**Bioconductor data containers for microbiome data.** **(a)
TreeSummarizedExperiment (TreeSE)** supports hierarchical, multi-table
microbiome analysis. A data object can accommodate: i) “assays” that describe
the abundances of features in each sample (*e.g.* taxonomic units, resistance
genes, or predicted functions), ii) “colData” with tabular metadata on each
sample, iii) “rowData” with tabular metadata on each feature (*e.g.* taxonomic
mappings), and iv) tree information describing feature hierarchies, such as
phylogenetic relations between the taxonomic units (rowTree) or host species
(colTree). TreeSE can thus link multiple numeric abundance tables and their
derived versions (*e.g.* statistical transformations, different taxonomic
levels, dimensionality reduction) with tabular and hierarchical side
information on the features (rows) and samples (columns). Additional slots for
reference sequences and experiment metadata are available. (Figure adapted from
[@huang2021]) **(b) MultiAssayExperiment (MAE)** facilitates the integration of
interlinked datasets that may represent different omics modalities
linked by potentially non-trivial sample mappings. (Figure adapted from
[@ramos2017])](figures/data_containers.png){#fig-datacontainers}

**Data resources.** Several open microbiome data resources are supported,
enabling direct data import into the aforementioned data containers for further
downstream analysis within the Bioconductor methods ecosystem
(@tbl-dataresources). Examples of the available data resources include
curatedMetagenomicData [@pasolli2017], HoloFood [@rogers2024],
microbiomeDataSets [@microbiomedatasets], and the EBI/MGnify database
[@gurbich2023]. In addition, Bioconductor packages often
include built-in demonstration datasets. Collectively, these resources
encompass taxonomic and functional profiles from thousands of
published microbiome studies and tens of thousands samples that can
be accessed with the available tools for research and educational purposes.

| Data resource                               | Package                                                                  | # studies or datasets | # samples |
|---------------------------------------------|--------------------------------------------------------------------------|-----------------------|-----------|
| Curated metagenomic data [@pasolli2017]     | curatedMetagenomicData [@pasolli2017]                                    | 93                    | 22,588    |
| HoloFood [@rogers2024]                      | HoloFoodR [@holofoodr]                                                   | -                     | 9,990     |
| MGnify [@gurbich2023]                       | MGnifyR [@mgnifyr]                                                       | 5,129                 | 629,821   |
| microbiomeDataSets [@microbiomedatasets]    | microbiomeDataSets [@microbiomedatasets]                                 | 6                     | 19,100    |
| Microbiome Benchmark Data [@gamboa-tuz2025] | MicrobiomeBenchmarkData [@gamboa-tuz2025]                                | 6                     | 1,125     |
| Human Microbiome Project 2 [@hmp2data]      | HMP2data [@hmp2data]                                                     | 3                     | 11,493    |
| A compilation of open microbiome datasets   | [https://github.com/microbiome/data](https://github.com/microbiome/data) | 8                     | 18,777    |

: Supporting open microbiome data resources (accessed 29/09/2025).
{#tbl-dataresources tbl-colwidths="[30,40,15,15]"}

# Data preparation {#sec-process}

While several tools for microbiome data analysis are available [@jeganathan2021; @willis2019;
@marcoszambrano2023], their implementations typically differ in terms of
expected input and output formats. Common data standards promote
modular data science workflows, speeding up methods development,
benchmarking, and replication. Building upon these principles, we describe a set
of interoperable software libraries that support microbiome data integration and
analysis based on the two integrative data containers, TreeSE and MAE
(Figure \ref{fig:workflow}).

This framework specifically supports the downstream analysis of
high-throughput sequencing data from microbiome experiments after summarizing
the original measurements into abundance
tables (*e.g.*, taxonomic or functional profiles). Taxonomic
abundance tables are one of the most commonly encountered data types.
Common measurement platforms include amplicon-based marker gene
studies (*e.g.*, 16S rRNA) and metagenomic sequencing, but alternative profiling
techniques are available including phylogenetic microarrays
[@rajilicstojanovic2009] and high-density optical mapping
[@abakumov2024]. The interoperability of the data containers
support modular workflow design for integrating different
omics data types into the broader Bioconductor ecosystem. This
mitigates reproducibility challenges posed by rapidly
evolving reference databases and profiling tools
(*e.g.*, MetaPhlAn2 vs. MetaPhlAn4) by enabling structured coexistence and
direct comparison within a unified, coherent framework
[@blanco2023; @truong2015].

**Data import.** While the TreeSE data container can be populated with microbiome
data from common tabular data types, certain non-standard file formats may
require additional data wrangling. To streamline data
import, we provide importers for a variety of standard microbiome file
formats including BIOM [@mcdonald2012], MetaPhlAn/HUMAnN (MetaPhlAn
v. 2&ndash;4, HUMAnN v. 3) [@mciver2017], Mothur [@schloss2009], QIIME2
[@bolyen2019] and the TAXonomic Profile Aggregation and
standardization format (taxpasta) [@beber2023]. These importers bridge the gap
between the raw sequence processing and downstream statistical analysis in R,
allowing researchers to focus on interpreting biological patterns rather than
wrangling file formats. Moreover, converters between alternative data
structures in R are available supporting seamless conversions between TreeSE
and BIOM, DADA2 [@callahan2016] and phyloseq (the mia::convert\* functions).
Similarly, microbiome data stored in TreeSE objects can also be
exported to external services
such as MicrobiomeDB [@oliveira2018] and other languages
(*e.g.*, Julia). These converters facilitate widespread access to
microbiome analysis methods. Imported datasets from different omics modalities
can then be interlinked within the MAE data container with its flexible sample
mappings as described above. This integrated data object can then be analyzed
with dedicated downstream methods.

**Quality control and exploration.** After importing the data into R,
the recommended first step involves basic data exploration and quality
control [@zuur2010]. In addition to methods for identifying
contaminant sequences and calculating summary statistics, the
framework includes a versatile set of custom visualization methods to
assess data variability, outliers, homogeneity of variance and other
aspects that may need to be considered in downstream data
analysis. The mia and miaViz packages provide
convenient access to many operations for the integrative
data containers regarding common tasks in microbiome
data exploration and visualization (@fig-downstream_analyses).

**Filtering and agglomeration.** Microbiome datasets can contain tens of
thousands of unique features depending on the sequencing resolution, used
reference databases, and the scope of analysis. Moreover, the integration of
multiple omics modalities adds another level of challenges. To summarize data
into biologically meaningful groups or to reduce multiple comparisons the data
is often subsetted, filtered or agglomerated to higher-level taxonomic or
functional categories or stratified by sample groups. The package ecosystem,
and the mia package in particular, feature various tools for such tasks.
For instance, the prevalence-based methods can be
helpful to remove rare features that are likely sequencing artifacts
or present otherwise limited value for the analysis. The mechanism of
alternative experiments (altExp) allows agglomeration
by taxonomic levels or other features such as
pathogenicity or functional potential (*e.g.*, butyrate production
[@kullberg2024], gut-brain modules [@vallescolomer2019], and other
host-associated signatures [@geistlinger2024]). The agglomerated data can be
incorporated as a new alternative experiment within the original data object.
The automated tracking of feature and sample links ensures that
changes are propagated across all relevant elements of the data object,
such as transformed abundance tables, phylogenetic trees, or interlinked
omics modalities. This helps to avoid redundant processing and enables
efficient management of interlinked multi-modal datasets.

**Transformation.** Statistical transformations commonly used in microbiome
research include relative abundance [@gloor2017],
tree-aware phylogenetic ILR (phILR) [@silverman2017], centered log-ratio (CLR),
and robust CLR (rCLR). The latter three aim to mitigate
compositionality bias [@martino2019]. A related concept that controls for
unequal sampling depth is data rarefaction. Rarefaction is
the process of repeatedly subsampling data to equal read depth and averaging
computed metrics [@schloss2024b]. Subsampling data to equal read depth
(*rarefying*) without iteration is not recommended
[@mcmurdie2014; @schloss2024; @willis2019b]. Our framework supports these and a
variety of other transformations commonly applied to microbiome data.

**Data enrichment.** Incorporating additional data, such as new versions of
feature abundance data resulting from data transformations, is a
frequent operation in applied data analysis, and it often takes place
iteratively. The integrative data containers
support gradual and organized accumulation of interlinked data elements
including transformations, aggregated data, dimension reduction, and derived
data on samples and features. We generally recommend
initializing the data objects at the beginning of the
workflow, where they can then be gradually enriched as the project advances.
Data can be enriched by importing new external information as well as by
deriving new information from the data itself into the same data container as
described in the next section. By systematically extending the upstream data
preparation step&mdash;instead of ad-hoc calculations scattered amidst the
workflow&mdash;one can achieve a well-organized and transparent
provenance of the integrated data object.

::: {.figure}
\begin{center}
\includegraphics[height=\textheight]{figures/workflow.png}
\end{center}
\captionof{figure}{\textbf{Multi-modal data integration workflow.} Following
data generation and preprocessing (top two panels), the integrative framework
covers five steps of microbiome data integration and analysis: 1) Data is
wrapped into a TreeSE data object (Figure~\ref{fig-datacontainers}),
containing information on the samples (tube icons) and their features
(\textit{e.g.} taxonomic, functional, or resistome profiles; colored circles).
TreeSE objects can accommodate multiple assay tables as well as additional
hierarchical structures. 2) The TreeSE data object can then be processed to get
cleaner data. This can include filtering out low-quality samples and features,
transforming assay data for example by normalizing it to enhance the
comparability of sample-specific abundances, or agglomerating the feature data
to specific taxonomic levels or functional categories. 3) The cleaned data can
then be analyzed using statistical models to uncover relations among samples
(\textit{e.g.}, with ordination methods), analyze drivers of variation in community diversity
and composition or test for differential abundance (\textit{e.g.} between
sampling time or location). 4) The framework supports the integration of
various types of omics modalities through MAE and 
omic-specific data containers provided by Bioconductor. Here,
for example, the taxonomic and metabolomic profiles of samples can be summarized
and compared with joint analysis among features. 5) Finally, Quarto notebooks
can aid reproducible analysis and reporting.}
\label{fig:workflow}
:::

# Downstream data analysis {#sec-analysis}

This framework is designed to improve the management of integrative analysis
workflows. It facilitates joint operations on interlinked data elements,
helping to shift the focus towards multi-modal analyses and
methods development while also supporting standard data analyses in
microbiome research.

Microbiome data analysis often relies on specialized statistical approaches,
since data sparsity, zero-inflation, compositionality, and heteroscedasticity
pose challenges for standard methods [@gloor2017; @jeganathan2021].
Furthermore, microbiome data is hierarchically
structured across taxonomic, ecological, temporal, and spatial levels
[@moreno-indias2021; @ruuskanen2021]. Alternative processing pipelines
introduce different biases and generate varying data formats; for example, the
same sample can yield differing taxonomic profiles depending on the
computational method used [@quince2017]. These and other unique characteristics
have wide-reaching implications for microbiome analysis, necessitating dedicated
analysis methods [@moreno-indias2021]. The shared data containers
support community-driven development and reproducible benchmarking of methods.
The following sections provide an overview of some of the currently
available tools, encompassing standard methods for taxonomic and other
abundance tables as well as integrative analyses between modalities.

**Community indicators.** Various general indicators are available
to quantify the microbial community state or other sample
properties. Community diversity represents a common measure in
microbial ecology [@bastiaanssen_bugs1_2023]
(@fig-downstream_analyses). Alpha diversity metrics generally quantify
the observed or estimated number of distinct taxa (_richness_) and
their distribution (_evenness_). Many diversity indices (*e.g.*,
Shannon index), combine these two aspects. Phylogenetic diversity
metrics additionally account the breadth of phylogenetic relatedness among
community members [@faith1992], but are more computationally demanding than
abundance-based metrics. However, a new implementation of the Stacked Faith's
Phylogenetic Diversity improves computational efficiency by 1&ndash;2
orders of magnitude [@armstrong2021]. The framework supports analyses
with rarefaction through the mia package [@schloss2024b]. Other
global indices of community state include measures of dominance,
rarity, skewness, kurtosis, modality, library size (number of reads),
or absolute abundance quantification. Sample information in the TreeSE object
(colData) can be enriched with these derived quantities. By default,
the mia::addAlpha method returns a set of complementary alpha
diversity measures [@cassol2025]. Alpha diversity indices have been frequently
applied to measure not only taxonomic diversity but also resistome
diversity [@salehi2025], and potentially other modalities.

**Community similarity** is commonly evaluated with beta diversity
measures, which quantify multivariate dissimilarity between feature communities
across all pairs of samples and can be used to
compare community composition quantitatively. The
framework supports common dissimilarity metrics in microbial ecology,
including Bray-Curtis, Jaccard, Aitchison
[@bastiaanssen_bugs1_2023; @martino2019],
and the phylogeny-aware UniFrac [@lozupone2005]. Community
dissimilarity is commonly visualized with unsupervised ordination
methods that project the data onto a lower-dimensional space
[@goodrich2014]. Among unsupervised techniques, those most frequently applied
include Principal Component Analysis (PCA), Principal Coordinates Analysis
(PCoA; also known as Multi-Dimensional Scaling, MDS), and Non-metric MDS
(@fig-downstream_analyses). Statistical tests
(*e.g.*, PERMANOVA) and supervised
ordination methods, such as
distance-based Redundancy Analysis (dbRDA), can complement the
unsupervised analyses, and help uncover and quantify covariation
between community composition and sample-specific covariates. This type of analysis can benefit
from compositionally robust measures (*e.g.*, robust
Aitchison distance [@martino2019]), or phylogeny-aware dissimilarity
measures specifically developed for microbiome research (*e.g.*, UniFrac
[@lozupone2005]). Sample dissimilarity provides the basis for many
specialized microbiome techniques, including microbial community typing
(clustering) with Dirichlet multinomial mixtures [@holmes2012]
and community state types [@digiulio2015], or factorization into
community signatures [@frioux2023]. These and other multivariate
methods are provided by mia, miaViz, and additional single-cell packages
(scater, scuttle [@amezquita2020]). Quantification of sample dissimilarities is
a common task, and many omics modalities can benefit from shared,
standardized methods provided by the integrative (Tree)SE framework.

**Differential abundance and mediation analysis.** Differential Abundance
Analysis (DAA)
and Differential Prevalence Analysis (DPA) methods are used to
identify indicator microbial features whose abundances vary across study groups
or continuous gradients, for example in terms of spatial or temporal distance between
samples (@fig-downstream_analyses). Most DAA methods used in microbiome research
support the (Tree)SE framework, such as
ANCOM-BC2 [@lin2024], LimROTS [@anwar2025], MaAsLin3
[@nickols2024], and radEmu [@clausen2024]. Leveraging the interoperability of
this framework, several of these methods have also found
applications with other omics modalities, including transcriptomics
[@love2014], metabolomics [@nickols2024], and proteomics [@anwar2025].
Importantly, recent studies have indicated that elementary methods, such as log-transformed
Total Sum Scaled (TSS) counts coupled with a linear model, or
presence/absence data with logistic regression, often yield the most consistent
results [@pelto2025; @gamboa-tuz2025]. Standard DAA typically focuses on
individual features, ignoring interactions or other covariation among feature
abundances. Thus, we support covariation analyses. While agglomerating collinear
features as a preprocessing step can reduce multiple testing, increase
statistical power, and improve interpretability, it is also prone to subjective
parameterization and limited generalizability. Feature can be grouped
based on clustering results from
abundance profiles or on external information drawn from phylogenetic
relations, taxonomic classifications or functional properties
(*e.g.*, aggregating based on metabolite production [@kullberg2024]). Beyond these,
early implementations of advanced methods, such as causal mediation
analysis tailored to microbiome multi-omics data, have been implemented for
SE [@jiang2025], highlighting its readiness to support the adoption of
constantly developing methodology.

**Function- and sequence-level analyses.** Taxonomic analysis is
frequently complemented by functional analyses, either based on
functional predictions derived from amplicon or metagenomic
sequences, or from parallel metabolomic, metatranscriptomic or other functional
assays. Our framework not only supports the incorporation of functional
abundance tables within the same data objects, but many statistical methods implemented in the ecosystem are generic and can therefore be applied across
different omics data types. Where appropriate, specialized methods are also
available; for instance, regarding metabolome
analyses, the notame package [@klavus2020] has recently added SE support,
and the R4MS project provides many interoperable tools [@gatto2025].
The framework has found recent applications in extended analyses of metagenomic
features, for instance, antimicrobial resistome composition [@salehi2025].
We have also converted a large-scale compilation of resistome profiles in 8,972
metagenome samples [@lee2023] into the TreeSE format to support
research and training activities in this area of microbiome analysis.

**Spatio-temporal, ecological, and epidemiological models.** Microbiome time
series [@berg2020] encompass a wide range of study designs, ranging
from short to long follow-up times, sparse to dense data, univariate
to multivariate and multi-omic monitoring [@velten2022], and real,
pseudo, or simulated time series. Moreover, an increasing number of
studies incorporate spatial dimension in microbiome analysis [@ruuskanen2021].
While spatio-temporal analyses require specialized methods, the miaTime R
package offers a set of data wrangling methods to support longitudinal,
multi-modal microbiome analysis. In addition, Bioconductor contains applicable
methods from other research domains, namely single-cell omics
[@amezquita2020] and spatial transcriptomics [@crowell2025], which use closely
related data containers. Spatial and temporal information can be
incorporated as a covariate (*e.g.*, location or time) or by estimating
temporally or spatially resolved metrics, such as divergence (dissimilarity
between consecutive time points), autocorrelation (similarity of closely
related samples), and feature modality (whether a feature exhibits multiple
abundance states). Another type of temporal analysis
includes the time-to-event, or survival models, which are encountered
regularly in epidemiological and population cohort studies; for microbiome
data, these analyses require specialized approaches, as proposed in previous
studies [@salosensaari2021; @calle2023]. We also provide methods to simulate
time series from standard models of microbial community dynamics including
Hubbell's neutral model, generalized Lotka-Volterra, and the
consumer-resource model through the miaSim package [@gao2023]. Our framework
provides a robust basis for developing extended methods for longitudinal
multi-omics, which is currently an active area of investigation [@sherwani2025].

**Multi-modal integration.** A key benefit of the ecosystem lies in its ability
to support the integration of microbial abundance data with other omics
modalities, such as resistomes [@salehi2025], metabolomes [@hintikka2021],
host genomes [@qin2022], transcriptomes, or proteomes
[@nyholm2020; @moreno-indias2021] (@fig-downstream_analyses). As we have
demonstrated, the integrative data containers support complex combinations of
omics datasets [@huber2015; @ramos2017; @amezquita2020]. By leveraging
standardized multi-assay data structures, users can apply a growing
number of integrative methods for multi-omics analysis. This eliminates the
need for manual data wrangling, reduces the risk of errors (*e.g.*, sample
mismatches), and enables the creation of modular, efficient, and reproducible
workflows. Such integrative approaches in microbiome
research can be broadly categorized into three groups: Machine Learning (ML)
prediction, association-based approaches, and
latent variable modeling [@jeganathan2021; @sherwani2025;
@moreno-indias2021; @mallick2024; @bastiaanssen2023;
@argelaguet2020; @singh2019; @mangnier2025].
ML approaches aim to integrate multiple omics modalities into a
single predictive model. However, naive concatenation of feature tables can
obscure modality-specific signals and reduce predictive performance. Methods
such as IntegratedLearner [@mallick2024] address this limitation by modeling
each modality separately while leveraging their complementary predictive
information within a unified framework.
Association-based methods aim to quantify relationships between datasets.
Whereas the cross-correlation analysis can be used to capture direct
associations of samples or features between two datasets, anansi can supervise
such analyses  by incorporating known functional relations. While these
association-based approaches are often easier to interpret, they do not
explicitly account for the fact that microbes operate within complex
communities. In contrast, latent variable methods such
as multi-omics factor analysis [@argelaguet2020],
joint robust PCA [@martino2025] and multiblock partial least
squares discriminant analysis [@singh2019] model shared latent
structure across modalities, enabling the discovery and visualization of
coordinated variation at the community level.
The available methodologies are rapidly increasing, and we aim to
accelerate the development, benchmarking, and standardization of
multi-modal extensions for microbiome analysis through a systematic integrative
framework proposed here.

**Other methods and programming environments.** Several other methods have
been employed in microbiome data science. Examples include network analysis
[@peschel2021], ML-based prediction and classification [@topcuoglu2021], and
microbe-set enrichment
analysis [@kou2020; @pham2025]. The recent adoption of interoperable
data containers and shared methodological approaches in related research fields
can significantly expand the transfer of existing toolkits to adjacent or emerging
disciplines. Importantly, this framework enables integration not only within the Bioconductor ecosystem
but also with external pipelines, including bioBakery [@navarro2024; @mciver2017], and foreign platforms through language-specific interfaces,
namely `Rcpp` for C+, `reticulate` for Python, and `JuliaCall` for Julia.
These connections allow users to combine the
most suitable tools across multiple frameworks, while ensuring
cross-language interoperability of structured data through standard containers such as TreeSE and MAE.

```{r}
#| label: downstream_analyses
#| eval: false
#| echo: false
#| message: false

# Load packages
pkgs <- c("maaslin3", "mia", "miaTime", "microbiomeDataSets")
temp <- sapply(pkgs, function(x) {
    suppressMessages(require(x, character.only = TRUE))
})

# Load data
tse <- SprockettTHData()
# Subset by taking only fecal samples from infants
tse <- tse[, tse[["Sample_Type"]] %in% c("Feces")]
tse <- tse[, tse[["Age_Class"]] %in% c("Infant")]

# Agglomerate and apply transformations
tse <- agglomerateByRanks(tse)
tse <- transformAssay(tse, method = "relabundance", altexp = altExpNames(tse))
tse <- transformAssay(
    tse,
    method = "clr", pseudocount = TRUE, altexp = altExpNames(tse)
)
altExp(tse, "prev_genus") <- subsetByPrevalent(
    altExp(tse, "Genus"),
    prevalence = 0.5, detection = 0.001
)
rownames(altExp(tse, "prev_genus")) <- gsub(
    "/| ", "_", rownames(altExp(tse, "prev_genus"))
)

# Calculate alpha diversity
tse <- addAlpha(tse, index = "faith")

# Perform PCoA
tse <- addMDS(tse, assay.type = "counts", method = "unifrac", niter = 10L)

# Perform dbRDA
tse <- addRDA(
    tse,
    assay.type = "clr", method = "euclidean",
    formula = x ~ Age_Months + Country, na.action = na.exclude
)
# Remove the RDA object as it takes a lot of space
attr(reducedDim(tse, "RDA"), "obj") <- NULL

# Perform DAA
res_full <- maaslin3(
    input_data = altExp(tse, "prev_genus"),
    output = "DAA_example",
    formula = "~ Age_Months + Sex + Country",
    normalization = "TSS",
    transform = "LOG",
    verbosity = "ERROR"
)
res <- res_full[["fit_data_abundance"]][["results"]]
res <- res[res[["value"]] == "Age_Months", ]
signif_feat <- res[
    !is.na(res[["qval_individual"]]) & res[["qval_individual"]] < 0.05,
    "feature"
]
rowData(altExp(tse, "prev_genus"))[["signif"]] <- rownames(
    altExp(tse, "prev_genus")
) %in% signif_feat

# Calculate time divergence
tse <- addStepwiseDivergence(
    tse,
    time.col = "Age_Months",
    group = "Subject_ID", method = "bray",
    tree = rowTree(tse)
)

# Save object so that it can be used in visualization step without running the
# analyses
saveRDS(tse, file.path("data", "tse.rds"))
```

```{r}
#| label: fig-downstream_analyses
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Common visualization methods applied to time series data on infants
#|   from three different countries (Finland, Bangladesh, and Bolivia) based on
#|   public datasets - SprockettTHData()
#|   [@sprockett2020; @subramanian2014; @vatanen2016] available via
#|   microbiomeDataSets [@microbiomedatasets]. **a**, Phylogenetic
#|   relationships. **b**, Phyla abundances of the samples ordered by age,
#|   stratified by country and gender. Actinobacteria is more abundant in
#|   Bolivia and Bangladesh compared to Finland.  **c&ndash;d**,
#|   Multidimensional Scaling (MDS) using UniFrac distances and distance-based
#|   Redundancy Analysis (dbRDA) based on Aitchison distance, respectively. Both
#|   ordination methods shows clear pattern separating microbial profiles in
#|   different countries. dbRDA shows an additional association between
#|   microbial profile and age. **e&ndash;f**, Faith's phylogenetic diversity by
#|   country and age respectively. **g**, Differential Abundance Analysis (DAA)
#|   results. Two genera exhibit an association with age. **h**, UniFrac
#|   divergence between consecutive time points. Results indicate slow
#|   stabilization of microbial communities."
#| fig-width: 19
#| fig-height: 10

# Load packages
pkgs <- c("miaViz", "patchwork", "scater")
temp <- sapply(pkgs, function(x) {
    suppressMessages(require(x, character.only = TRUE))
})

# Load data
file_path <- file.path("data", "tse.rds")
if (!file.exists(file_path)) {
    stop("Run manually the chunk above. '", file_path, "' cannot found.")
}
tse <- readRDS(file_path)

# Get phylum names
phylum_names <- sort(unique(rowData(altExp(tse, "Phylum"))[, "Phylum"]))
# Named vector of phylum colors
phylum_colors <- setNames(
    scater:::.get_palette("tableau20")[seq_along(phylum_names)],
    phylum_names
)

# Define a font size
title_size <- 16
text_size <- 12

# Visualize phylogeny
p_tree <- plotRowTree(
    altExp(tse, "Family"),
    edge.colour.by = "Phylum", layout = "ape"
) +
    scale_color_manual(values = phylum_colors) + # Manually set colors
    guides(color = "none")

# Visualize abundances with barplot
p_bar <- plotAbundance(
    altExp(tse, "Phylum"),
    assay.type = "relabundance",
    col.var = c("Country", "Sex"), facet.cols = TRUE, scales = "free",
    order.col.by = "Age_Months"
) +
    guides(fill = guide_legend(title = "Phylum")) +
    scale_color_manual(values = phylum_colors) +
    theme(
        axis.title.y = element_text(size = title_size),
        axis.text.y = element_text(size = text_size)
    )

# Visualize alpha diversity
p_alpha <- plotColData(
    tse,
    y = "faith", x = "Country", colour_by = "Country",
    show_boxplot = TRUE, show_violin = FALSE, point_alpha = 0.2
) +
    labs(y = "Faith", x = "Country") +
    theme(
        legend.position = "none",
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size, angle = 20, hjust = 0.6, vjust = 0.8),
        axis.text.y = element_text(size = text_size)
    )

p_alpha2 <- plotColData(
    tse,
    y = "faith", x = "Age_Months"
) +
    geom_smooth(method = "lm") +
    labs(y = "Faith", x = "Age in months") +
    theme(
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size),
        axis.text.y = element_text(size = text_size)
    )

# Visualize PCoA
p_mds <- plotMDS(tse, colour_by = "Country") +
    theme(
        legend.position = "inside",
        legend.justification = c("right", "top"),
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size),
        axis.text.y = element_text(size = text_size),
        legend.title = element_text(size = title_size),
        legend.text = element_text(size = text_size)
      )

# Visualize dbRDA
p_rda <- plotRDA(tse, "RDA", colour.by = "Country", add.significance = FALSE) +
    theme(
        legend.position = "none",
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size),
        axis.text.y = element_text(size = text_size)
    )

# Visualize abundances vs age
p_daa <- plotExpression(
    altExp(tse, "prev_genus"),
    features = rownames(altExp(tse, "prev_genus"))[rowData(altExp(tse, "prev_genus"))[["signif"]]],
    assay.type = "relabundance", x = "Age_Months", ncol = 3L,
) +
    scale_y_log10() +
    geom_smooth() +
    labs(x = "Age in months", y = "Relative abundance") +
    theme(
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size),
        axis.text.y = element_text(size = text_size)
    )

# Visualize divergence
tse_sub <- tse[, !is.na(tse[["time_diff"]]) & tse[["time_diff"]] >= 0.8 & tse[["time_diff"]] <= 1.2]
p_divergence <- plotColData(tse_sub, y = "divergence", x = "Age_Months") +
    geom_smooth(method = "lm") +
    labs(x = "Age in months", y = expression(Delta ~ "UniFrac" ~ "(N, N-1)")) +
    theme(
        legend.position = "none",
        axis.title.x = element_text(size = title_size),
        axis.title.y = element_text(size = title_size),
        axis.text.x = element_text(size = text_size),
        axis.text.y = element_text(size = text_size)
    )

# Combine plots
p1 <- (p_tree / p_bar) +
    plot_layout(guides = "collect")
p2 <- (p_mds / p_rda) +
    plot_layout(heights = c(1, 1))
p3 <- ((p_alpha | p_alpha2) / p_daa / p_divergence) +
    plot_layout(heights = c(1, 1, 1))
p <- (p1 | p2 | p3) +
    plot_layout(widths = c(0.33, 0.4, 0.32)) +
    plot_annotation(tag_levels = "a") &
    theme(plot.tag = element_text(size = 22, face = "bold"))
p
```

# Accessible analysis and community {#sec-community}

Modern science is embodied by a critical community capable of
assessing discoveries and replicating results [@wootton2016].
Open source developer communities present a timely manifestation of this
principle, with technical, legal, sociological and epistemic consequences
[@hocquet2024]. FAIR principles for software support inclusive
access to well-documented and user-friendly methods [@barker2022],
and are facilitated by the extensive quality control and
documentation standards of Bioconductor [@huber2015].
In addition to these technical
demands, the development and quality control of research software
significantly relies on a strong community engagement, which
stimulates continuous software improvement and promotes the emergence and
dissemination of evidence-based best practices [@berg2020;
@barker2022; @hocquet2024]. We recommend starting with simple
methods [@pelto2025], calculating derived quantities such
as diversity, prevalence, and other such measures early in
the workflow, and complementing analyses with visual exploration. These
practices can improve the transparency and interpretability of the
analysis and its results. Expressive code and meaningful comments support
reproducibility and verification.

Harmonized data structures and access to validated methods support standardized
and transparent data wrangling across studies while facilitating method sharing,
benchmarking, and reproducible comparisons. These capabilities help establish
community best practices and build consensus while allowing users to flexibly
select and evaluate methods appropriate for their specific data and research
questions. The Bioconductor framework naturally aligns with these goals through
its open ecosystem of data resources, interoperable software, and educational
materials. Community building is actively supported through online training
resources and communication channels, making it easier for both developers and
users to contribute to the continued development and standardization of the
ecosystem. Together with interactive applications, these resources offer
accessible entry points for developers and users, regardless of their
programming experience.

**The OMA book.** _Orchestrating Microbiome Analysis with Bioconductor_ (OMA)
is an online book published at [https://bioconductor.org/books/release/OMA](https://bioconductor.org/books/release/OMA)
that provides a comprehensive guide to the (Tree)SE framework for microbiome
data science in Bioconductor, including practical guidelines,
reproducible workflows, and training material for integrative microbiome
analyses. The book is versioned according to Bioconductor's biannual software
release cycles. The work incorporates rich feedback from the user community and
microbiome training courses. The OMA book aligns to a broader Bioconductor convention, where
prominent data structures [@huber2015; @ramos2017; @huang2021]
designed for multi-table data integration have been described across several
domains, such as single-cell omics
(OSCA book [@amezquita2020]) and spatial transciptomics (OSTA book [@crowell2025]).
These resources
demonstrate how standard data containers can support statistical
programming on interlinked multi-table datasets.

**Interactive applications.** Graphical user interfaces provide
accessible entry points to the framework. The iSEEtree and miaDash
packages provide interactive analysis and visualization tools for
microbiome data [@benedetti2025]. The former extends iSEE to provide support for
TreeSE objects, while the latter builds on iSEEtree to enable analyses without
requiring any programming. Their interface supports automatic generation of
source code and replicable workflows as well as analysis templates for users
without a strong programming expertise. Moreover, the tidy R paradigm, available
through tidyomics, simplifies the syntax and streamlines analytical
workflows [@hutchison2024]. These tools allow users to analyze and
visualize complex data combinations with ease.

**Community.** Its widespread adoption has made R one of the most cited software
resources of our time [@pearson2025] and the introduction of the
shared data structures facilitates the broad interoperability of the
methods ecosystem. Bioconductor has a strong tradition in fostering the
R user community through shared data structures and training activities
[@drnevich2025]. An active community of developers and users thrives online,
supporting the software sustainability and promoting good practices in
scientific computing [@wilson2017]. The community interacts via several online
communication channels, including Bioconductor Zulip, GitHub platform,
and an email list. The development of the framework is influenced by
feedback from the user community through these channels, regular
meetings, and training workshops.

# Outlook {#sec-discussion}

```{r}
#| label: bioc_views
#| eval: false
#| echo: false
#| message: false

if (!require("BiocPkgTools")) {
    BiocManager::install("BiocPkgTools")
    library("BiocPkgTools")
}
# Get packages and info based on their BiocViews
pkgs <- biocPkgList()
#
fields <- sapply(pkgs[["biocViews"]], function(x) {
    field <- "Other"
    if (any(grepl("genomic|genetic|gene|genom|DNA", x, ignore.case = TRUE))) {
        field <- "Genomics"
    }
    if (any(grepl("proteomics|protein", x, ignore.case = TRUE))) {
        field <- "Proteomics"
    }
    if (any(grepl("metabolomics|metabolome|metabolite|Lipidom|massspectro", x, ignore.case = TRUE))) {
        field <- "Metabolomics"
    }
    if (any(grepl("transcripto|RNA-seq|RNA", x, ignore.case = TRUE))) {
        field <- "Transcriptomics"
    }
    if (any(grepl("immuno", x, ignore.case = TRUE))) {
        field <- "Immunology"
    }
    if (any(grepl("cytom", x, ignore.case = TRUE))) {
        field <- "Cytometrics"
    }
    if (any(grepl("microarray|chip", x, ignore.case = TRUE))) {
        field <- "Microarray"
    }
    if (any(grepl("single-cell|singlecell", x, ignore.case = TRUE))) {
        field <- "Single-cell"
    }
    if (any(grepl("metagenom|microbiome|16S|microbiota|amplicon|shotgun|microb|metatranscript|metametabolo|metaproteo", x, ignore.case = TRUE))) {
        field <- "Microbiome"
    }
    return(field)
})
pkgs[["Field"]] <- fields

pkgs[, c("Package", "Field")] |>
    write.csv(file = file.path("data", "packages.csv"), row.names = FALSE)
```

```{r}
#| label: bioc_views2
#| eval: true
#| echo: false
pkgs <- read.csv(file.path("data", "packages.csv"))
```

To enjoy the benefits of automation while also
maintaining flexibility, data scientists often recombine open-source
software components into custom workflows. However, the lack of commonly adopted
data representation standards can hinder the development, sharing,
and application of modular, interoperable statistical workflows. Here, we present
an R/Bioconductor ecosystem for multi-modal microbiome data science. Building upon
previously established data structures [@huang2021; @ramos2017],
we have introduced a comprehensive, interoperable ecosystem
developed as an extensive community effort. We also provide exhaustive guidance in
the OMA book, thereby supporting independent learning and reproducible analysis in
microbiome research and beyond.

**The Bioconductor community** is widely recognized for
developing open research software for life sciences. Its
strengths include minimal documentation standards,
portability across operating platforms, quality
control (*e.g.* peer-review, unit testing, and continuous integration),
usability (APIs), robustness (dependency testing and deprecation mechanisms),
and a biannual release cycle with semantic versioning. Over the past two decades
Bioconductor has grown into a repository of more than 2,300 actively maintained
software packages, of which `r sum(pkgs$Field == "Microbiome")` are now
classified under the _Microbiome_ Task View. At the core of the Bioconductor's
data science strategy has been the idea of shared data representations
[@gentleman2004; @lawrence2013]. This has enhanced interoperability, mitigated
redundant efforts and facilitated the translation of methods between research
domains [@delia2023; @barker2022]. In addition to these technical advances, the
shared standards have played an important role in the formation of broad
developer and user communities [@drnevich2025].

**Microbiome research** has shifted towards integrative analyses
linking various omics towards enhanced mechanistic and causal insights
[@joos2024; @bastiaanssen_bugs1_2023]. Technological advances in long-read
sequencing, functional assays, spatial profiling, and single-cell analysis
have created an increasing demand for integrative analysis
techniques that can link microbiome data with other
omics along varying spatial and temporal dimensions. In concert with parallel
developments in other omics, our framework specifically addresses the need to
integrate and jointly analyze interlinked abundance tables and supporting
hierarchical and tabular side information across different omics. Whereas our
framework supports the basic data wrangling and microbiome data science
operations, it also provides the structured input necessary for probabilistic
inference and scalable computation. This is essential as the recent
developments in microbiome and multi-omics
workflows increasingly rely on advanced statistical and ML
techniques such as Bayesian modeling and Gaussian processes to accommodate
uncertainty and temporal or spatial structures [@cheng2019]. Building upon these
statistical foundations, the rapid emergence of foundation models and deep
learning approaches in biomedicine calls for robust, scalable data
representations [@awais2025] such as those we have adopted.
Applications include transfer learning across studies,
interpretable embeddings of microbiome multi-omic profiles,
computational prediction of microbial function and metabolites
[@mallick2019; @weiss2016], and integration with AI agents for hypothesis
generation and automated analysis, among others [@zou2025; @gao2024].

**Open data science infrastructures** have become critical components
of reproducible research and a vast body of applications critically
rely on open-source methods and workflows [@berg2020; @barker2022;
@hocquet2024]. Standardized data structures can enhance the transparency,
reliability, and explainability of analysis workflows. We have
documented the open data science ecosystem for multi-modal microbiome research
built on widely adopted Bioconductor data containers, software libraries,
data resources, tutorials, and a growing community of developers.
Further extensions could include methods for temporal, spatial, and multi-omic
analyses, improved support for community metadata standards and FAIR data
practices
(*e.g.*, Minimum Information about any (x) Sequence, MIxS [@field2008]),
along with an end-to-end case study demonstrating how an example
analysis can be conducted in practice from raw data to results.
Finally, supporting interoperability between R and other
programming languages could facilitate the application of other ML
frameworks. As the field evolves, the developer communities will continue
to integrate new tools and approaches to improve integrative analyses of
microbiome data.

# References
<!-- the pattern below controls the placement of the references -->
:::{#refs}
:::

## Consortia {#sec-consortium}

OMA consortium: Bastiaanssen, Thomaz; Beber, Moritz E.; Bektanov, Aituar;
Benedetti, Giulio; Borman, Tuomas; Braccia, Domenick J.; Calgaro, Matteo;
Callan, Danielle; Chaudhary, Jiya; Corrada Bravo, Héctor; Courbayre Dussau,
Basil; Das, Sneha; Eckermann, Henrik; Erawijantari, Pande; Ernst, Felix G.M.;
Gaudron-Parry, Elliot; Gopalakrishnan, Shyam; Hakula, Maija;
Hanhineva, Kati; Havulinna, Aki; Hindström, Rasmus; Jeba, Akewak;
Lahti, Leo; Limborg, Morten; Lindgren, Himmi; Liu, Yihan; Mallick, Himel;
Mongad, Dattatray; Müller, Christian L.; Muluh, Geraldson; Nynäs, Signe;
Pasanen, Jesse; Paulo, Joao; Pelto, Juho; Peschel, Stefanie; Potbhare, Renuka;
Pralas, Théotime; Ramos, Marcel; Raulo, Aura; Sannikov, Artur; Seraidarian, Ely;
Serizay, Jacques; Shetty, Sudarshan; Shigdel, Rajesh; Suksi, Vilhelm;
Tammi, Eetu; Valderamma, Benjamin; Waldron, Levi; Yang, Lu.

## Contributions {#sec-contributions}

LL and TB initiated and coordinated the work. All authors contributed
material, participated in manuscript preparation, and approved the final
version.

## Acknowledgments {#sec-acknowledgments}

This project received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 952914. TB
was additionally supported by the University of Turku Graduate School. We would
like to thank the following individuals for their support for the project:
Richa Ashma; Nitin Bayal; Chouaib Benchraka; Yang Cao; Axel Dagnaud;
Teo Dallier; Noah de Gunst; Karoline Faust; Robert D. Finn; Samuel Hillman; Ruizhu
Huang; Hervé Pagés; Lorna Richardson; Daniel Rios Garza; Mark Robinson; Matti
Ruuskanen; Mahkameh Salehi. We would like to thank Christopher Quince and
Kihyun Lee for providing the resistome profiling data,
Justin Silverman for philr support, and Ben Allen for MGnifyR support.
We also extend our thanks to all contributors of the open-source resources used
in this work and acknowledge CSC – IT Center for Science, Finland,
for computational resources.