## ----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", empty.rm = TRUE)
altExp(tse,"species") <- agglomerateByRank(tse, rank = "Species", empty.rm = FALSE)
dim(x1)
dim(altExp(tse,"species"))

## -----------------------------------------------------------------------------
tse <- agglomerateByRanks(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 <- transformAssay(enterotype, method = "log10", pseudocount = 1)
assayNames(enterotype)

## -----------------------------------------------------------------------------
data(GlobalPatterns, package = "mia")

tse.subsampled <- rarefyAssay(GlobalPatterns, 
                                  sample = 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 <- addAlpha(tse, index = "shannon")
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, 
                              assay.type = "counts", 
                              keep_dist = TRUE)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(phyloseq)
data(esophagus, package = "phyloseq")

## -----------------------------------------------------------------------------
esophagus
esophagus <- convertFromPhyloseq(esophagus)
esophagus

## -----------------------------------------------------------------------------
# Specific taxa
assay(tse['522457',], "counts") |> head()
# Specific sample
assay(tse[,'CC1'], "counts") |> head()


## -----------------------------------------------------------------------------
data(esophagus, package = "mia")
top_taxa <- getTop(esophagus,
                    method = "mean",
                    top = 5,
                    assay.type = "counts")
top_taxa

## -----------------------------------------------------------------------------
molten_data <- meltAssay(tse,
                         assay.type = "counts",
                         add.row = TRUE,
                         add.col = TRUE
)
molten_data

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