## ---- echo=FALSE-------------------------------------------------------------- knitr::opts_chunk$set(cache = FALSE, fig.width = 9, message = FALSE, warning = FALSE) ## ---- eval=FALSE-------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("mia") ## ----load-packages, message=FALSE, warning=FALSE------------------------------ library("mia") ## ----------------------------------------------------------------------------- data(GlobalPatterns, package = "mia") tse <- GlobalPatterns tse ## ----------------------------------------------------------------------------- # print the available taxonomic ranks colnames(rowData(tse)) taxonomyRanks(tse) # subset to taxonomic data only rowData(tse)[,taxonomyRanks(tse)] ## ----------------------------------------------------------------------------- table(taxonomyRankEmpty(tse, "Species")) head(getTaxonomyLabels(tse)) ## ----------------------------------------------------------------------------- # agglomerate at the Family taxonomic rank x1 <- agglomerateByRank(tse, rank = "Family") ## How many taxa before/after agglomeration? nrow(tse) nrow(x1) ## ----------------------------------------------------------------------------- # with agglomeration of the tree x2 <- agglomerateByRank(tse, rank = "Family", agglomerateTree = TRUE) nrow(x2) # same number of rows, but rowTree(x1) # ... different rowTree(x2) # ... tree ## ----------------------------------------------------------------------------- data(enterotype, package = "mia") taxonomyRanks(enterotype) agglomerateByRank(enterotype) ## ----------------------------------------------------------------------------- altExp(tse, "family") <- x2 ## ----------------------------------------------------------------------------- x1 <- agglomerateByRank(tse, rank = "Species", na.rm = TRUE) altExp(tse,"species") <- agglomerateByRank(tse, rank = "Species", na.rm = FALSE) dim(x1) dim(altExp(tse,"species")) ## ----------------------------------------------------------------------------- altExps(tse) <- splitByRanks(tse) tse altExpNames(tse) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- assayNames(enterotype) anterotype <- transformCounts(enterotype, method = "log10", pseudocount = 1) assayNames(enterotype) ## ----------------------------------------------------------------------------- data(GlobalPatterns, package = "mia") tse.subsampled <- subsampleCounts(GlobalPatterns, min_size = 60000, name = "subsampled", replace = TRUE, seed = 1938) tse.subsampled ## ----------------------------------------------------------------------------- library(MultiAssayExperiment) mae <- MultiAssayExperiment(c("originalTreeSE" = GlobalPatterns, "subsampledTreeSE" = tse.subsampled)) mae ## ----------------------------------------------------------------------------- # To extract specifically the subsampled TreeSE experiments(mae)$subsampledTreeSE ## ----------------------------------------------------------------------------- tse <- estimateDiversity(tse) colnames(colData(tse))[8:ncol(colData(tse))] ## ----------------------------------------------------------------------------- library(scater) altExp(tse,"Genus") <- runMDS(altExp(tse,"Genus"), FUN = vegan::vegdist, method = "bray", name = "BrayCurtis", ncomponents = 5, exprs_values = "counts", keep_dist = TRUE) ## ---- message=FALSE, warning=FALSE-------------------------------------------- data("esophagus", package = "phyloseq") ## ----------------------------------------------------------------------------- esophagus esophagus <- makeTreeSummarizedExperimentFromPhyloseq(esophagus) esophagus ## ----------------------------------------------------------------------------- abund <- getAbundanceSample(tse, "CC1", abund_values = "counts") all(abund == assay(tse, "counts")[,"CC1"]) abund <- getAbundanceFeature(tse, "522457", abund_values = "counts") all(abund == assay(tse, "counts")["522457",]) ## ----------------------------------------------------------------------------- data("esophagus", package = "mia") top_taxa <- getTopTaxa(esophagus, method = "mean", top = 5, abund_values = "counts") top_taxa ## ----------------------------------------------------------------------------- molten_data <- meltAssay(tse, abund_values = "counts", add_row_data = TRUE, add_col_data = TRUE ) molten_data ## ----------------------------------------------------------------------------- sessionInfo()