mia 1.4.0
mia implements tools for microbiome analysis based on the
SummarizedExperiment (Morgan et al. 2020), SingleCellExperiment (Amezquita et al. 2020) and
TreeSummarizedExperiment (Huang 2021) infrastructure. Data wrangling and analysis
are the main scope of this package.
To install mia, install BiocManager first, if it is not installed.
Afterwards use the install function from BiocManager.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("mia")library("mia")TreeSummarizedExperiment objectA few example datasets are available via mia. For this vignette the
GlobalPatterns dataset is loaded first.
data(GlobalPatterns, package = "mia")
tse <- GlobalPatterns
tse## class: TreeSummarizedExperiment
## dim: 19216 26
## metadata(0):
## assays(1): counts
## rownames(19216): 549322 522457 ... 200359 271582
## rowData names(7): Kingdom Phylum ... Genus Species
## colnames(26): CL3 CC1 ... Even2 Even3
## colData names(7): X.SampleID Primer ... SampleType Description
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (19216 rows)
## rowTree: 1 phylo tree(s) (19216 leaves)
## colLinks: NULL
## colTree: NULLOne of the main topics for analysing microbiome data is the application of taxonomic data to describe features measured. The interest lies in the connection between individual bacterial species and their relation to each other.
mia does not rely on a specific object type to hold taxonomic data, but
uses specific columns in the rowData of a TreeSummarizedExperiment object.
taxonomyRanks can be used to construct a character vector of available
taxonomic levels. This can be used, for example, for subsetting.
# print the available taxonomic ranks
colnames(rowData(tse))## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"taxonomyRanks(tse)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"# subset to taxonomic data only
rowData(tse)[,taxonomyRanks(tse)]## DataFrame with 19216 rows and 7 columns
## Kingdom Phylum Class Order Family
## <character> <character> <character> <character> <character>
## 549322 Archaea Crenarchaeota Thermoprotei NA NA
## 522457 Archaea Crenarchaeota Thermoprotei NA NA
## 951 Archaea Crenarchaeota Thermoprotei Sulfolobales Sulfolobaceae
## 244423 Archaea Crenarchaeota Sd-NA NA NA
## 586076 Archaea Crenarchaeota Sd-NA NA NA
## ... ... ... ... ... ...
## 278222 Bacteria SR1 NA NA NA
## 463590 Bacteria SR1 NA NA NA
## 535321 Bacteria SR1 NA NA NA
## 200359 Bacteria SR1 NA NA NA
## 271582 Bacteria SR1 NA NA NA
## Genus Species
## <character> <character>
## 549322 NA NA
## 522457 NA NA
## 951 Sulfolobus Sulfolobusacidocalda..
## 244423 NA NA
## 586076 NA NA
## ... ... ...
## 278222 NA NA
## 463590 NA NA
## 535321 NA NA
## 200359 NA NA
## 271582 NA NAThe columns are recognized case insensitive. Additional functions are available to check for validity of taxonomic information or generate labels based on the taxonomic information.
table(taxonomyRankEmpty(tse, "Species"))##
## FALSE TRUE
## 1413 17803head(getTaxonomyLabels(tse))## [1] "Class:Thermoprotei" "Class:Thermoprotei_1"
## [3] "Species:Sulfolobusacidocaldarius" "Class:Sd-NA"
## [5] "Class:Sd-NA_1" "Class:Sd-NA_2"For more details see the man page ?taxonomyRanks.
Agglomeration of data based on these taxonomic descriptors can be performed
using functions implemented in mia. In addition to the aggValue functions
provide by TreeSummarizedExperiment agglomerateByRank is available.
agglomerateByRank does not require tree data to be present.
agglomerateByRank constructs a factor to guide merging from the available
taxonomic information. For more information on merging have a look at the man
page via ?mergeRows.
# agglomerate at the Family taxonomic rank
x1 <- agglomerateByRank(tse, rank = "Family")
## How many taxa before/after agglomeration?
nrow(tse)## [1] 19216nrow(x1)## [1] 603Tree data can also be shrunk alongside agglomeration, but this is turned of by default.
# with agglomeration of the tree
x2 <- agglomerateByRank(tse, rank = "Family",
agglomerateTree = TRUE)
nrow(x2) # same number of rows, but## [1] 603rowTree(x1) # ... different##
## Phylogenetic tree with 19216 tips and 19215 internal nodes.
##
## Tip labels:
## 549322, 522457, 951, 244423, 586076, 246140, ...
## Node labels:
## , 0.858.4, 1.000.154, 0.764.3, 0.995.2, 1.000.2, ...
##
## Rooted; includes branch lengths.rowTree(x2) # ... tree##
## Phylogenetic tree with 496 tips and 260 internal nodes.
##
## Tip labels:
## Family:Cenarchaeaceae, Family:SAGMA-X, Family:Nitrososphaeraceae, Family:Sulfolobaceae, Family:Halobacteriaceae, Family:MSBL1, ...
## Node labels:
## root:ALL, Kingdom:Archaea, Phylum:Crenarchaeota, Class:C2, Class:Sd-NA, Class:Thaumarchaeota, ...
##
## Rooted; includes branch lengths.For agglomerateByRank to work, taxonomic data must be present. Even though
only one rank is available for the enterotype dataset, agglomeration can be
performed effectively de-duplicating entries for the genus level.
data(enterotype, package = "mia")
taxonomyRanks(enterotype)## [1] "Genus"agglomerateByRank(enterotype)## class: TreeSummarizedExperiment
## dim: 552 280
## metadata(1): agglomerated_by_rank
## assays(1): counts
## rownames(552): NA Prosthecochloris ... Syntrophococcus Mogibacterium
## rowData names(1): Genus
## colnames(280): AM.AD.1 AM.AD.2 ... TS98_V2 TS99.2_V2
## colData names(9): Enterotype Sample_ID ... Age ClinicalStatus
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULLTo keep data tidy, the agglomerated data can be stored as an alternative experiment in the object of origin. With this synchronized sample subsetting becomes very easy.
altExp(tse, "family") <- x2Keep in mind, that if you set na.rm = TRUE, rows with NA or similar value
(defined via the empty.fields argument) will be removed. Depending on these
settings different number of rows will be returned.
x1 <- agglomerateByRank(tse, rank = "Species", na.rm = TRUE)
altExp(tse,"species") <- agglomerateByRank(tse, rank = "Species", na.rm = FALSE)
dim(x1)## [1] 944 26dim(altExp(tse,"species"))## [1] 2307 26For convenience the function splitByRanks is available, which agglomerates
data on all ranks selected. By default all available ranks will be used.
The output is compatible to be stored as alternative experiments.
altExps(tse) <- splitByRanks(tse)
tse## class: TreeSummarizedExperiment
## dim: 19216 26
## metadata(0):
## assays(1): counts
## rownames(19216): 549322 522457 ... 200359 271582
## rowData names(7): Kingdom Phylum ... Genus Species
## colnames(26): CL3 CC1 ... Even2 Even3
## colData names(7): X.SampleID Primer ... SampleType Description
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(7): Kingdom Phylum ... Genus Species
## rowLinks: a LinkDataFrame (19216 rows)
## rowTree: 1 phylo tree(s) (19216 leaves)
## colLinks: NULL
## colTree: NULLaltExpNames(tse)## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species"Constructing a taxonomic tree from taxonomic data stored in rowData is quite
straightforward and uses mostly functions implemented in
TreeSummarizedExperiment.
taxa <- rowData(altExp(tse,"Species"))[,taxonomyRanks(tse)]
taxa_res <- resolveLoop(as.data.frame(taxa))
taxa_tree <- toTree(data = taxa_res)
taxa_tree$tip.label <- getTaxonomyLabels(altExp(tse,"Species"))
rowNodeLab <- getTaxonomyLabels(altExp(tse,"Species"), make_unique = FALSE)
altExp(tse,"Species") <- changeTree(altExp(tse,"Species"),
rowTree = taxa_tree,
rowNodeLab = rowNodeLab)Transformation of count data stored in assays is also a main task when work
with microbiome data. transformCounts can be used for this and offers a few
choices of available transformations. A modified object is returned and the
transformed counts are stored in a new assay.
assayNames(enterotype)## [1] "counts"anterotype <- transformCounts(enterotype, method = "log10", pseudocount = 1)
assayNames(enterotype)## [1] "counts"For more details have a look at the man page ?transformCounts.
Sub-sampling to equal number of counts per sample. Also known as rarefying.
data(GlobalPatterns, package = "mia")
tse.subsampled <- subsampleCounts(GlobalPatterns,
min_size = 60000,
name = "subsampled",
replace = TRUE,
seed = 1938)
tse.subsampled## class: TreeSummarizedExperiment
## dim: 12403 25
## metadata(1): subsampleCounts_min_size
## assays(2): counts subsampled
## rownames(12403): 549322 522457 ... 535321 200359
## rowData names(7): Kingdom Phylum ... Genus Species
## colnames(25): CL3 CC1 ... Even2 Even3
## colData names(7): X.SampleID Primer ... SampleType Description
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (12403 rows)
## rowTree: 1 phylo tree(s) (19216 leaves)
## colLinks: NULL
## colTree: NULLAlternatively, one can save both original TreeSE and subsampled TreeSE within a MultiAssayExperiment object.
library(MultiAssayExperiment)
mae <- MultiAssayExperiment(c("originalTreeSE" = GlobalPatterns,
"subsampledTreeSE" = tse.subsampled))
mae## A MultiAssayExperiment object of 2 listed
## experiments with user-defined names and respective classes.
## Containing an ExperimentList class object of length 2:
## [1] originalTreeSE: TreeSummarizedExperiment with 19216 rows and 26 columns
## [2] subsampledTreeSE: TreeSummarizedExperiment with 12403 rows and 25 columns
## Functionality:
## experiments() - obtain the ExperimentList instance
## colData() - the primary/phenotype DataFrame
## sampleMap() - the sample coordination DataFrame
## `$`, `[`, `[[` - extract colData columns, subset, or experiment
## *Format() - convert into a long or wide DataFrame
## assays() - convert ExperimentList to a SimpleList of matrices
## exportClass() - save data to flat files# To extract specifically the subsampled TreeSE
experiments(mae)$subsampledTreeSE## class: TreeSummarizedExperiment
## dim: 12403 25
## metadata(1): subsampleCounts_min_size
## assays(2): counts subsampled
## rownames(12403): 549322 522457 ... 535321 200359
## rowData names(7): Kingdom Phylum ... Genus Species
## colnames(25): CL3 CC1 ... Even2 Even3
## colData names(7): X.SampleID Primer ... SampleType Description
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (12403 rows)
## rowTree: 1 phylo tree(s) (19216 leaves)
## colLinks: NULL
## colTree: NULLIn the field of microbiome ecology several indices to describe samples and community of samples are available. In this vignette we just want to give a very brief introduction.
Functions for calculating alpha and beta diversity indices are available.
Using estimateDiversity multiple diversity indices are calculated by default
and results are stored automatically in colData. Selected indices can be
calculated individually by setting index = "shannon" for example.
tse <- estimateDiversity(tse)
colnames(colData(tse))[8:ncol(colData(tse))]## [1] "coverage" "fisher" "gini_simpson"
## [4] "inverse_simpson" "log_modulo_skewness" "shannon"
## [7] "faith"Beta diversity indices are used to describe inter-sample connections.
Technically they are calculated as dist object and reduced dimensions can
be extracted using cmdscale. This is wrapped up in the runMDS function
of the scater package, but can be easily used to calculated beta diversity
indices using the established functions from the vegan package or any other
package using comparable inputs.
library(scater)
altExp(tse,"Genus") <- runMDS(altExp(tse,"Genus"),
FUN = vegan::vegdist,
method = "bray",
name = "BrayCurtis",
ncomponents = 5,
exprs_values = "counts",
keep_dist = TRUE)JSD and UniFrac are implemented in mia as well. calculateJSD can be used
as a drop-in replacement in the example above (omit the method argument as
well) to calculate the JSD. For calculating the UniFrac distance via
calculateUniFrac either a TreeSummarizedExperiment must be used or a tree
supplied via the tree argument. For more details see ?calculateJSD,
?calculateUniFrac or ?calculateDPCoA.
runMDS performs the decomposition. Alternatively runNMDS can also be used.
estimateDominance and estimateEvenness implement other sample-wise indices.
The function behave equivalently to estimateDiversity. For more information
see the corresponding man pages.
To make migration and adoption as easy as possible several utility functions are available.
Functions to load data from biom files, QIIME2 output, DADA2 objects
(Callahan et al. 2016) or phyloseq objects are available.
data("esophagus", package = "phyloseq")esophagus## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 58 taxa and 3 samples ]
## phy_tree() Phylogenetic Tree: [ 58 tips and 57 internal nodes ]esophagus <- makeTreeSummarizedExperimentFromPhyloseq(esophagus)
esophagus## class: TreeSummarizedExperiment
## dim: 58 3
## metadata(0):
## assays(1): counts
## rownames(58): 59_8_22 59_5_13 ... 65_9_9 59_2_6
## rowData names(0):
## colnames(3): B C D
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (58 rows)
## rowTree: 1 phylo tree(s) (58 leaves)
## colLinks: NULL
## colTree: NULLFor more details have a look at the man page, for examples
?makeTreeSummarizedExperimentFromPhyloseq.
getAbundanceFeature and getAbundanceSample are wrappers on row-wise or
column-wise assay data subsetting.
abund <- getAbundanceSample(tse, "CC1", abund_values = "counts")
all(abund == assay(tse, "counts")[,"CC1"])## [1] TRUEabund <- getAbundanceFeature(tse, "522457", abund_values = "counts")
all(abund == assay(tse, "counts")["522457",])## [1] TRUEgetTopTaxa returns a vector of the most top abundant feature IDs.
data("esophagus", package = "mia")
top_taxa <- getTopTaxa(esophagus,
method = "mean",
top = 5,
abund_values = "counts")
top_taxa## [1] "59_2_6" "59_7_6" "59_8_22" "59_5_19" "65_6_2"To generate tidy data as used and required in most of the tidyverse,
meltAssay can be used. A data.frame in the long format will be returned.
molten_data <- meltAssay(tse,
abund_values = "counts",
add_row_data = TRUE,
add_col_data = TRUE
)
molten_data## # A tibble: 499,616 × 24
## FeatureID SampleID counts Kingdom Phylum Class Order Family Genus Species
## <fct> <fct> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 549322 CL3 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 2 549322 CC1 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 3 549322 SV1 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 4 549322 M31Fcsw 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 5 549322 M11Fcsw 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 6 549322 M31Plmr 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 7 549322 M11Plmr 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 8 549322 F21Plmr 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 9 549322 M31Tong 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## 10 549322 M11Tong 0 Archaea Crenarcha… Ther… <NA> <NA> <NA> <NA>
## # … with 499,606 more rows, and 14 more variables: X.SampleID <fct>,
## # Primer <fct>, Final_Barcode <fct>, Barcode_truncated_plus_T <fct>,
## # Barcode_full_length <fct>, SampleType <fct>, Description <fct>,
## # coverage <int>, fisher <dbl>, gini_simpson <dbl>, inverse_simpson <dbl>,
## # log_modulo_skewness <dbl>, shannon <dbl>, faith <dbl>sessionInfo()## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
##
## 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
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] phyloseq_1.40.0 scater_1.24.0
## [3] ggplot2_3.3.5 scuttle_1.6.0
## [5] mia_1.4.0 MultiAssayExperiment_1.22.0
## [7] TreeSummarizedExperiment_2.4.0 Biostrings_2.64.0
## [9] XVector_0.36.0 SingleCellExperiment_1.18.0
## [11] SummarizedExperiment_1.26.0 Biobase_2.56.0
## [13] GenomicRanges_1.48.0 GenomeInfoDb_1.32.0
## [15] IRanges_2.30.0 S4Vectors_0.34.0
## [17] BiocGenerics_0.42.0 MatrixGenerics_1.8.0
## [19] matrixStats_0.62.0 BiocStyle_2.24.0
##
## loaded via a namespace (and not attached):
## [1] ggbeeswarm_0.6.0 colorspace_2.0-3
## [3] ellipsis_0.3.2 BiocNeighbors_1.14.0
## [5] ggrepel_0.9.1 bit64_4.0.5
## [7] fansi_1.0.3 decontam_1.16.0
## [9] codetools_0.2-18 splines_4.2.0
## [11] sparseMatrixStats_1.8.0 cachem_1.0.6
## [13] knitr_1.38 ade4_1.7-19
## [15] jsonlite_1.8.0 cluster_2.1.3
## [17] BiocManager_1.30.17 compiler_4.2.0
## [19] assertthat_0.2.1 Matrix_1.4-1
## [21] fastmap_1.1.0 lazyeval_0.2.2
## [23] cli_3.3.0 BiocSingular_1.12.0
## [25] htmltools_0.5.2 tools_4.2.0
## [27] igraph_1.3.1 rsvd_1.0.5
## [29] gtable_0.3.0 glue_1.6.2
## [31] GenomeInfoDbData_1.2.8 reshape2_1.4.4
## [33] dplyr_1.0.8 Rcpp_1.0.8.3
## [35] jquerylib_0.1.4 rhdf5filters_1.8.0
## [37] vctrs_0.4.1 multtest_2.52.0
## [39] ape_5.6-2 nlme_3.1-157
## [41] DECIPHER_2.24.0 iterators_1.0.14
## [43] DelayedMatrixStats_1.18.0 xfun_0.30
## [45] stringr_1.4.0 beachmat_2.12.0
## [47] lifecycle_1.0.1 irlba_2.3.5
## [49] zlibbioc_1.42.0 MASS_7.3-57
## [51] scales_1.2.0 biomformat_1.24.0
## [53] parallel_4.2.0 rhdf5_2.40.0
## [55] yaml_2.3.5 memoise_2.0.1
## [57] gridExtra_2.3 yulab.utils_0.0.4
## [59] sass_0.4.1 stringi_1.7.6
## [61] RSQLite_2.2.12 foreach_1.5.2
## [63] ScaledMatrix_1.4.0 tidytree_0.3.9
## [65] permute_0.9-7 BiocParallel_1.30.0
## [67] rlang_1.0.2 pkgconfig_2.0.3
## [69] bitops_1.0-7 evaluate_0.15
## [71] lattice_0.20-45 Rhdf5lib_1.18.0
## [73] purrr_0.3.4 treeio_1.20.0
## [75] bit_4.0.4 tidyselect_1.1.2
## [77] plyr_1.8.7 magrittr_2.0.3
## [79] bookdown_0.26 R6_2.5.1
## [81] generics_0.1.2 DelayedArray_0.22.0
## [83] DBI_1.1.2 pillar_1.7.0
## [85] withr_2.5.0 mgcv_1.8-40
## [87] survival_3.3-1 RCurl_1.98-1.6
## [89] tibble_3.1.6 crayon_1.5.1
## [91] utf8_1.2.2 rmarkdown_2.14
## [93] viridis_0.6.2 grid_4.2.0
## [95] data.table_1.14.2 blob_1.2.3
## [97] vegan_2.6-2 digest_0.6.29
## [99] tidyr_1.2.0 munsell_0.5.0
## [101] DirichletMultinomial_1.38.0 beeswarm_0.4.0
## [103] viridisLite_0.4.0 vipor_0.4.5
## [105] bslib_0.3.1Amezquita, Robert, Aaron Lun, Etienne Becht, Vince Carey, Lindsay Carpp, Ludwig Geistlinger, Federico Marini, et al. 2020. “Orchestrating Single-Cell Analysis with Bioconductor.” Nature Methods 17: 137–45. https://www.nature.com/articles/s41592-019-0654-x.
Callahan, Benjamin J, Paul J McMurdie, Michael J Rosen, Andrew W Han, Amy Jo A Johnson, and Susan P Holmes. 2016. “DADA2: High-Resolution Sample Inference from Illumina Amplicon Data.” Nature Methods 13: 581–83. https://doi.org/10.1038/nmeth.3869.
Huang, Ruizhu. 2021. TreeSummarizedExperiment: TreeSummarizedExperiment: A S4 Class for Data with Tree Structures.
Morgan, Martin, Valerie Obenchain, Jim Hester, and Hervé Pagès. 2020. SummarizedExperiment: SummarizedExperiment Container. https://bioconductor.org/packages/SummarizedExperiment.