## ---- 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")
se <- GlobalPatterns
se

## -----------------------------------------------------------------------------
# print the available taxonomic ranks
colnames(rowData(se))
taxonomyRanks(se)
# subset to taxonomic data only
rowData(se)[,taxonomyRanks(se)]

## -----------------------------------------------------------------------------
table(taxonomyRankEmpty(se, "Species"))
head(getTaxonomyLabels(se))

## -----------------------------------------------------------------------------
# agglomerate at the Family taxonomic rank
x1 <- agglomerateByRank(se, rank = "Family")
## How many taxa before/after agglomeration?
nrow(se)
nrow(x1)

## -----------------------------------------------------------------------------
# with agglomeration of the tree
x2 <- agglomerateByRank(se, 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(se, "family") <- x2

## -----------------------------------------------------------------------------
x1 <- agglomerateByRank(se, rank = "Species", na.rm = TRUE)
altExp(se,"species") <- agglomerateByRank(se, rank = "Species", na.rm = FALSE)
dim(x1)
dim(altExp(se,"species"))

## -----------------------------------------------------------------------------
altExps(se) <- splitByRanks(se)
se
altExpNames(se)

## -----------------------------------------------------------------------------
taxa <- rowData(altExp(se,"Species"))[,taxonomyRanks(se)]
taxa_res <- resolveLoop(as.data.frame(taxa))
taxa_tree <- toTree(data = taxa_res)
taxa_tree$tip.label <- getTaxonomyLabels(altExp(se,"Species"))
rowNodeLab <- getTaxonomyLabels(altExp(se,"Species"), make_unique = FALSE)
altExp(se,"Species") <- changeTree(altExp(se,"Species"),
                                   rowTree = taxa_tree,
                                   rowNodeLab = rowNodeLab)

## -----------------------------------------------------------------------------
assayNames(enterotype)
anterotype <- transformCounts(enterotype, method = "log10", pseudocount = 1)
assayNames(enterotype)

## -----------------------------------------------------------------------------
se <- estimateDiversity(se)
colnames(colData(se))[8:ncol(colData(se))]

## -----------------------------------------------------------------------------
library(scater)
altExp(se,"Genus") <- runMDS(altExp(se,"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(se, "CC1", abund_values = "counts")
all(abund == assay(se, "counts")[,"CC1"])
abund <- getAbundanceFeature(se, "522457", abund_values = "counts")
all(abund == assay(se, "counts")["522457",])

## -----------------------------------------------------------------------------
top_taxa <- getTopTaxa(se,
                       method = "mean",
                       top = 5,
                       abund_values = "counts")
top_taxa

## -----------------------------------------------------------------------------
molten_data <- meltAssay(se,
                         abund_values = "counts",
                         add_row_data = TRUE,
                         add_col_data = TRUE
                         )
molten_data

## -----------------------------------------------------------------------------
sessionInfo()