We recommend switching from phyloseq to the
TreeSummarizedExperiment based methods described in the
Orchestrating Microbiome Analysis
(OMA) online book, which
is where method development now takes place. The
mia package provides the
analysis methods for that container.
To ease the transition, most functions in this package accept a
TreeSummarizedExperiment (or any other SummarizedExperiment-derived
object) wherever they accept a phyloseq object. This vignette walks
through the standard operations on a TreeSummarizedExperiment. The
companion vignette, vignette("vignette", package="microbiome"), gives
the general introduction to the package.
Almost everything is unchanged: the same function called on a
TreeSummarizedExperiment returns the same value it would return for
the equivalent phyloseq object. There are two things to know before
starting, both covered in their own sections below:
transform() stores its result as a new named assay rather than
overwriting the counts, so the result is read back with
abundances(x, assay.type=).library(microbiome)
Attaching mia also attaches its dependencies, two of which export
functions that share a name with a microbiome function. Since they
are attached later, they take precedence:
IRanges::transform() masks microbiome::transform()Biostrings::coverage() masks microbiome::coverage()Neither raises an error when called on a microbiome object, so the
symptom is a silently wrong result rather than a failure. Call these
two with their package prefix whenever mia is attached, as this
vignette does throughout:
microbiome::transform(tse, "compositional")
microbiome::coverage(tse)
The rest of the package is unaffected.
mia converts an existing phyloseq object. Here we use the
dietswap data set (O’Keefe et al. 2015) that ships with this package.
library(mia)
data(dietswap)
tse <- convertFromPhyloseq(dietswap)
tse
## class: TreeSummarizedExperiment
## dim: 130 222
## metadata(0):
## assays(1): counts
## rownames(130): Actinomycetaceae Aerococcus ... Xanthomonadaceae
## Yersinia et rel.
## rowData names(3): Phylum Family Genus
## colnames(222): Sample-1 Sample-2 ... Sample-221 Sample-222
## colData names(8): subject sex ... timepoint.within.group bmi_group
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
mia also reads the common file formats directly, with
importBiom(), importMothur(), importQIIME2() and others, so a
conversion step is not needed for new analyses. See the OMA
book for those.
abundances(), meta() and taxa() are the three accessors the rest
of the package is built on. They work the same way on both containers,
returning the abundance matrix, the sample metadata and the feature
names.
a <- abundances(tse)
dim(a)
## [1] 130 222
a[seq_len(4), seq_len(3)]
## Sample-1 Sample-2 Sample-3
## Actinomycetaceae 0 1 0
## Aerococcus 0 0 0
## Aeromonas 0 0 0
## Akkermansia 18 97 67
head(meta(tse), 3)
## subject sex nationality group sample timepoint
## Sample-1 byn male AAM DI Sample-1 4
## Sample-2 nms male AFR HE Sample-2 2
## Sample-3 olt male AFR HE Sample-3 2
## timepoint.within.group bmi_group
## Sample-1 1 obese
## Sample-2 1 lean
## Sample-3 1 overweight
head(taxa(tse), 5)
## [1] "Actinomycetaceae" "Aerococcus"
## [3] "Aeromonas" "Akkermansia"
## [5] "Alcaligenes faecalis et rel."
For a TreeSummarizedExperiment these read the assay, colData and
rownames respectively. The number of features and samples comes from
the object itself.
c(features=nrow(tse), samples=ncol(tse))
## features samples
## 130 222
A SummarizedExperiment can hold several assays at once. abundances()
takes an assay.type argument to pick between them. It defaults to the
counts assay when one is present, and otherwise to the first assay.
assayNames(tse)
## [1] "counts"
identical(abundances(tse), abundances(tse, assay.type="counts"))
## [1] TRUE
Asking for an assay that does not exist is an error rather than a silent fallback:
abundances(tse, assay.type="nonexistent")
## Error in `.se_assay()`:
## ! assay.type 'nonexistent' not found. Available assays: counts
This is the one place where the result differs from the phyloseq
version, and it is deliberate. For a phyloseq object, transform()
overwrites the abundance table. For a SummarizedExperiment,
transforming is additive: the result is stored as a new assay named
after the transformation, and the original counts are left untouched.
This follows the mia::transformAssay() convention.
tse <- microbiome::transform(tse, "compositional")
assayNames(tse)
## [1] "counts" "compositional"
The transformed values are read back by naming the assay:
abundances(tse, assay.type="compositional")[seq_len(4), seq_len(3)]
## Sample-1 Sample-2 Sample-3
## Actinomycetaceae 0.000000000 4.742483e-05 0.000000000
## Aerococcus 0.000000000 0.000000e+00 0.000000000
## Aeromonas 0.000000000 0.000000e+00 0.000000000
## Akkermansia 0.002127911 4.600209e-03 0.002319704
The counts are still there, unchanged:
abundances(tse)[seq_len(4), seq_len(3)]
## Sample-1 Sample-2 Sample-3
## Actinomycetaceae 0 1 0
## Aerococcus 0 0 0
## Aeromonas 0 0 0
## Akkermansia 18 97 67
Use the name argument to control the assay name, for instance to
follow the mia naming convention:
tse <- microbiome::transform(tse, "compositional", name="relabundance")
assayNames(tse)
## [1] "counts" "compositional" "relabundance"
Transformations chain, since each one adds an assay:
tse <- microbiome::transform(tse, "clr")
tse <- microbiome::transform(tse, "Z")
assayNames(tse)
## [1] "counts" "compositional" "relabundance" "clr"
## [5] "Z"
When you only need the transformed matrix and not a modified object,
pass transform to abundances() directly. This never round-trips
through the object and works identically for both containers:
clr <- abundances(tse, transform="clr")
clr[seq_len(4), seq_len(3)]
## Sample-1 Sample-2 Sample-3
## Actinomycetaceae -2.6385281 -2.016809 -3.114983
## Aerococcus -2.6385281 -3.115421 -3.114983
## Aeromonas -2.6385281 -3.115421 -3.114983
## Akkermansia 0.9723898 2.157578 1.790291
head(prevalence(tse, detection=1, sort=TRUE), 5)
## Vibrio Uncultured Mollicutes
## 1 1
## Uncultured Clostridiales II Uncultured Clostridiales I
## 1 1
## Tannerella et rel.
## 1
core_members() lists the taxa passing a detection and prevalence
threshold, and core() returns the object filtered down to them.
core_members(tse, detection=1, prevalence=50/100)
## [1] "Akkermansia"
## [2] "Alcaligenes faecalis et rel."
## [3] "Allistipes et rel."
## [4] "Anaerostipes caccae et rel."
## [5] "Anaerotruncus colihominis et rel."
## [6] "Anaerovorax odorimutans et rel."
## [7] "Bacteroides fragilis et rel."
## [8] "Bacteroides intestinalis et rel."
## [9] "Bacteroides ovatus et rel."
## [10] "Bacteroides plebeius et rel."
## [11] "Bacteroides splachnicus et rel."
## [12] "Bacteroides stercoris et rel."
## [13] "Bacteroides uniformis et rel."
## [14] "Bacteroides vulgatus et rel."
## [15] "Bifidobacterium"
## [16] "Bryantella formatexigens et rel."
## [17] "Bulleidia moorei et rel."
## [18] "Butyrivibrio crossotus et rel."
## [19] "Campylobacter"
## [20] "Clostridium (sensu stricto)"
## [21] "Clostridium cellulosi et rel."
## [22] "Clostridium colinum et rel."
## [23] "Clostridium difficile et rel."
## [24] "Clostridium leptum et rel."
## [25] "Clostridium nexile et rel."
## [26] "Clostridium orbiscindens et rel."
## [27] "Clostridium ramosum et rel."
## [28] "Clostridium sphenoides et rel."
## [29] "Clostridium stercorarium et rel."
## [30] "Clostridium symbiosum et rel."
## [31] "Collinsella"
## [32] "Coprobacillus catenaformis et rel."
## [33] "Coprococcus eutactus et rel."
## [34] "Desulfovibrio et rel."
## [35] "Dialister"
## [36] "Dorea formicigenerans et rel."
## [37] "Eggerthella lenta et rel."
## [38] "Enterobacter aerogenes et rel."
## [39] "Enterococcus"
## [40] "Escherichia coli et rel."
## [41] "Eubacterium biforme et rel."
## [42] "Eubacterium cylindroides et rel."
## [43] "Eubacterium hallii et rel."
## [44] "Eubacterium rectale et rel."
## [45] "Eubacterium siraeum et rel."
## [46] "Eubacterium ventriosum et rel."
## [47] "Faecalibacterium prausnitzii et rel."
## [48] "Fusobacteria"
## [49] "Helicobacter"
## [50] "Klebisiella pneumoniae et rel."
## [51] "Lachnobacillus bovis et rel."
## [52] "Lachnospira pectinoschiza et rel."
## [53] "Lactobacillus gasseri et rel."
## [54] "Lactobacillus plantarum et rel."
## [55] "Megasphaera elsdenii et rel."
## [56] "Mitsuokella multiacida et rel."
## [57] "Oceanospirillum"
## [58] "Oscillospira guillermondii et rel."
## [59] "Outgrouping clostridium cluster XIVa"
## [60] "Oxalobacter formigenes et rel."
## [61] "Papillibacter cinnamivorans et rel."
## [62] "Parabacteroides distasonis et rel."
## [63] "Peptococcus niger et rel."
## [64] "Peptostreptococcus micros et rel."
## [65] "Phascolarctobacterium faecium et rel."
## [66] "Prevotella melaninogenica et rel."
## [67] "Prevotella oralis et rel."
## [68] "Prevotella ruminicola et rel."
## [69] "Prevotella tannerae et rel."
## [70] "Proteus et rel."
## [71] "Roseburia intestinalis et rel."
## [72] "Ruminococcus bromii et rel."
## [73] "Ruminococcus callidus et rel."
## [74] "Ruminococcus gnavus et rel."
## [75] "Ruminococcus lactaris et rel."
## [76] "Ruminococcus obeum et rel."
## [77] "Sporobacter termitidis et rel."
## [78] "Streptococcus bovis et rel."
## [79] "Streptococcus intermedius et rel."
## [80] "Streptococcus mitis et rel."
## [81] "Subdoligranulum variable at rel."
## [82] "Sutterella wadsworthia et rel."
## [83] "Tannerella et rel."
## [84] "Uncultured Clostridiales I"
## [85] "Uncultured Clostridiales II"
## [86] "Uncultured Mollicutes"
## [87] "Veillonella"
## [88] "Vibrio"
## [89] "Yersinia et rel."
tse.core <- core(tse, detection=1, prevalence=50/100)
c(features=nrow(tse.core), samples=ncol(tse.core))
## features samples
## 89 222
Subsetting preserves the class, so the result is a
TreeSummarizedExperiment that can be fed back into any of these
functions:
class(tse.core)
## [1] "TreeSummarizedExperiment"
## attr(,"package")
## [1] "TreeSummarizedExperiment"
rare() and rare_members() are the complements, and
core_abundance() gives the fraction of the community the core
accounts for:
head(core_abundance(tse, detection=1, prevalence=50/100), 3)
## [1] NA
length(rare_members(tse, detection=1, prevalence=50/100))
## [1] 41
remove_taxa() and remove_samples() drop features and samples by
name, again returning the same class they were given.
drop <- taxa(tse)[seq_len(5)]
tse.sub <- remove_taxa(drop, tse)
c(before=nrow(tse), after=nrow(tse.sub))
## before after
## 130 125
drop.s <- colnames(tse)[seq_len(5)]
tse.sub <- remove_samples(drop.s, tse)
c(before=ncol(tse), after=ncol(tse.sub))
## before after
## 222 217
The plotting functions that go through the accessors take a
TreeSummarizedExperiment and return a ggplot object as usual.
plot_core(tse, prevalences=seq(.1, 1, .2), detections=3)
Figure 1: Core microbiota heatmap
boxplot_abundance(tse, x="nationality", y=taxa(tse)[[1]])
Figure 2: Abundance of a single taxon across groups
boxplot_alpha(tse, x_var="nationality", index="shannon")
Figure 3: Alpha diversity across groups
plot_density(tse, variable=taxa(tse)[[1]])
Figure 4: Abundance density for a single taxon
plot_atlas() and plot_tipping() accept a
TreeSummarizedExperiment in the same way.
Given a time field in the sample metadata, the time series helpers
work as they do for phyloseq. time_normalize() writes back into the
object and so returns a TreeSummarizedExperiment.
tse.t <- tse
tse.t$time <- tse.t$timepoint
length(timesplit(tse.t))
## [1] 4
head(time_sort(tse.t)$time, 3)
## NULL
class(time_normalize(tse.t))
## [1] "TreeSummarizedExperiment"
## attr(,"package")
## [1] "TreeSummarizedExperiment"
baseline() and collapse_replicates() likewise accept and return the
same container.
The functions that build or manipulate the taxonomy table are
phyloseq only, and error on a TreeSummarizedExperiment:
aggregate_taxa(), aggregate_rare()plot_composition()map_levels()psmelt2(), otu_tibble(), tax_tibble(), sample_tibble()add_besthit(), add_refseq()read_*() familymia covers this ground for TreeSummarizedExperiment objects, and
does it better: use mia::agglomerateByRank() in place of
aggregate_taxa(), mia::meltSE() in place of psmelt2(), and the
miaViz package for
composition plots. The OMA
book documents these.
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mia_1.21.7 TreeSummarizedExperiment_2.21.0
## [3] Biostrings_2.81.9 XVector_0.53.0
## [5] SingleCellExperiment_1.35.2 MultiAssayExperiment_1.39.1
## [7] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [9] GenomicRanges_1.65.4 Seqinfo_1.3.2
## [11] IRanges_2.47.5 S4Vectors_0.51.9
## [13] BiocGenerics_0.59.12 generics_0.1.4
## [15] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [17] microbiome_1.35.1 ggplot2_4.0.3
## [19] phyloseq_1.57.0 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 gridExtra_2.3.1
## [3] permute_0.9-10 rlang_1.3.0
## [5] magrittr_2.0.5 ade4_1.7-24
## [7] scater_1.41.2 otel_0.2.0
## [9] compiler_4.6.1 mgcv_1.9-4
## [11] DelayedMatrixStats_1.35.0 vctrs_0.7.3
## [13] reshape2_1.4.5 stringr_1.6.0
## [15] pkgconfig_2.0.3 crayon_1.5.3
## [17] fastmap_1.2.0 magick_2.9.1
## [19] labeling_0.4.3 scuttle_1.23.2
## [21] rmarkdown_2.32 ggbeeswarm_0.7.3
## [23] DirichletMultinomial_1.55.0 tinytex_0.60
## [25] purrr_1.2.2 xfun_0.60
## [27] bluster_1.23.1 cachem_1.1.0
## [29] beachmat_2.29.2 jsonlite_2.0.0
## [31] biomformat_1.41.0 DelayedArray_0.39.6
## [33] BiocParallel_1.47.0 irlba_2.3.7
## [35] parallel_4.6.1 cluster_2.1.8.3
## [37] R6_2.6.1 bslib_0.12.0
## [39] stringi_1.8.9 RColorBrewer_1.1-3
## [41] jquerylib_0.1.4 iterators_1.0.14
## [43] Rcpp_1.1.2 bookdown_0.48
## [45] knitr_1.52 DECIPHER_3.9.4
## [47] Matrix_1.7-6 splines_4.6.1
## [49] igraph_2.3.3 tidyselect_1.2.1
## [51] dichromat_2.0-1 abind_1.4-8
## [53] yaml_2.3.12 viridis_0.6.5
## [55] vegan_2.7-6 codetools_0.2-20
## [57] lattice_0.23-1 tibble_3.3.1
## [59] plyr_1.8.9 withr_3.0.3
## [61] treeio_1.37.0 S7_0.2.2
## [63] Rtsne_0.17 evaluate_1.0.5
## [65] survival_3.8-12 pillar_1.11.1
## [67] BiocManager_1.30.27 foreach_1.5.2
## [69] sparseMatrixStats_1.25.0 scales_1.4.0
## [71] tidytree_0.4.8 glue_1.8.1
## [73] lazyeval_0.2.3 tools_4.6.1
## [75] data.table_1.18.6.1 BiocNeighbors_2.7.3
## [77] ScaledMatrix_1.21.0 fs_2.1.0
## [79] grid_4.6.1 tidyr_1.3.2
## [81] ape_5.8-1 nlme_3.1-171
## [83] beeswarm_0.4.0 BiocSingular_1.29.1
## [85] vipor_0.4.7 cli_3.6.6
## [87] rsvd_1.0.5 rappdirs_0.3.4
## [89] S4Arrays_1.13.0 viridisLite_0.4.3
## [91] dplyr_1.2.1 gtable_0.3.6
## [93] yulab.utils_0.2.5 sass_0.4.10
## [95] digest_0.6.39 SparseArray_1.13.2
## [97] ggrepel_0.9.8 decontam_1.33.0
## [99] farver_2.1.2 multtest_2.69.0
## [101] htmltools_0.5.9 lifecycle_1.0.5
## [103] MASS_7.3-66