## ----setup, echo=FALSE--------------------------------------------------------
knitr::opts_chunk$set(fig.width=7, fig.height=5, message=FALSE, warning=FALSE)

# mia and TreeSummarizedExperiment are in Suggests, so show the code but
# do not evaluate it when they are missing.
tse.available <- requireNamespace("mia", quietly=TRUE) &&
    requireNamespace("TreeSummarizedExperiment", quietly=TRUE) &&
    requireNamespace("SummarizedExperiment", quietly=TRUE)
knitr::opts_chunk$set(eval=tse.available)

## ----loading------------------------------------------------------------------
library(microbiome)

## ----masking, eval=FALSE------------------------------------------------------
# microbiome::transform(tse, "compositional")
# microbiome::coverage(tse)

## ----convert------------------------------------------------------------------
library(mia)

data(dietswap)
tse <- convertFromPhyloseq(dietswap)
tse

## ----accessors----------------------------------------------------------------
a <- abundances(tse)
dim(a)
a[seq_len(4), seq_len(3)]

head(meta(tse), 3)

head(taxa(tse), 5)

## ----dims---------------------------------------------------------------------
c(features=nrow(tse), samples=ncol(tse))

## ----assaynames---------------------------------------------------------------
assayNames(tse)

## ----assaytype----------------------------------------------------------------
identical(abundances(tse), abundances(tse, assay.type="counts"))

## ----assaytype-error, error=TRUE----------------------------------------------
try({
abundances(tse, assay.type="nonexistent")
})

## ----transform----------------------------------------------------------------
tse <- microbiome::transform(tse, "compositional")
assayNames(tse)

## ----transform-read-----------------------------------------------------------
abundances(tse, assay.type="compositional")[seq_len(4), seq_len(3)]

## ----transform-counts---------------------------------------------------------
abundances(tse)[seq_len(4), seq_len(3)]

## ----transform-name-----------------------------------------------------------
tse <- microbiome::transform(tse, "compositional", name="relabundance")
assayNames(tse)

## ----transform-chain----------------------------------------------------------
tse <- microbiome::transform(tse, "clr")
tse <- microbiome::transform(tse, "Z")
assayNames(tse)

## ----transform-abundances-----------------------------------------------------
clr <- abundances(tse, transform="clr")
clr[seq_len(4), seq_len(3)]

## ----alpha--------------------------------------------------------------------
head(alpha(tse, index=c("shannon", "gini_simpson")), 3)

## ----indices------------------------------------------------------------------
head(richness(tse), 3)
head(microbiome::diversity(tse, index="shannon"), 3)
head(evenness(tse, index="pielou"), 3)
head(dominance(tse, index="absolute"), 3)
head(rarity(tse, index="log_modulo_skewness"), 3)

## ----dominant-----------------------------------------------------------------
head(dominant(tse), 3)
top_taxa(tse, n=5)

## ----divergence---------------------------------------------------------------
ref <- apply(abundances(tse), 1, median)
head(divergence(tse, ref, method="bray"), 3)

## ----summaries----------------------------------------------------------------
head(readcount(tse), 3)
head(microbiome::coverage(tse), 3)
head(low_abundance(tse), 3)

## ----prevalence---------------------------------------------------------------
head(prevalence(tse, detection=1, sort=TRUE), 5)

## ----core-members-------------------------------------------------------------
core_members(tse, detection=1, prevalence=50/100)

## ----core---------------------------------------------------------------------
tse.core <- core(tse, detection=1, prevalence=50/100)
c(features=nrow(tse.core), samples=ncol(tse.core))

## ----core-class---------------------------------------------------------------
class(tse.core)

## ----rare---------------------------------------------------------------------
head(core_abundance(tse, detection=1, prevalence=50/100), 3)
length(rare_members(tse, detection=1, prevalence=50/100))

## ----remove-------------------------------------------------------------------
drop <- taxa(tse)[seq_len(5)]
tse.sub <- remove_taxa(drop, tse)
c(before=nrow(tse), after=nrow(tse.sub))

## ----remove-samples-----------------------------------------------------------
drop.s <- colnames(tse)[seq_len(5)]
tse.sub <- remove_samples(drop.s, tse)
c(before=ncol(tse), after=ncol(tse.sub))

## ----plot-core, fig.cap="Core microbiota heatmap."----------------------------
plot_core(tse, prevalences=seq(.1, 1, .2), detections=3)

## ----plot-boxplot, fig.cap="Abundance of a single taxon across groups."-------
boxplot_abundance(tse, x="nationality", y=taxa(tse)[[1]])

## ----plot-alpha, fig.cap="Alpha diversity across groups."---------------------
boxplot_alpha(tse, x_var="nationality", index="shannon")

## ----plot-density, fig.cap="Abundance density for a single taxon."------------
plot_density(tse, variable=taxa(tse)[[1]])

## ----time---------------------------------------------------------------------
tse.t <- tse
tse.t$time <- tse.t$timepoint

length(timesplit(tse.t))
head(time_sort(tse.t)$time, 3)
class(time_normalize(tse.t))

## ----sessioninfo, eval=TRUE---------------------------------------------------
sessionInfo()