## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) knitr::opts_chunk$set(dev = "png", dev.args = list(type = "cairo-png")) ## ----------------------------------------------------------------------------- if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("miaViz") ## ----setup, message=FALSE----------------------------------------------------- library(miaViz) data(GlobalPatterns, package = "mia") ## ----------------------------------------------------------------------------- plotAbundance(GlobalPatterns, rank = NULL, features = "549322", assay.type = "counts") ## ----------------------------------------------------------------------------- GlobalPatterns <- transformAssay(GlobalPatterns, method = "relabundance") ## ----------------------------------------------------------------------------- plotAbundance(GlobalPatterns, rank = "Kingdom", assay.type = "relabundance") ## ----------------------------------------------------------------------------- prev_phylum <- getPrevalentTaxa(GlobalPatterns, rank = "Phylum", detection = 0.01, onRankOnly = TRUE) ## ----------------------------------------------------------------------------- plotAbundance(GlobalPatterns[rowData(GlobalPatterns)$Phylum %in% prev_phylum], rank = "Phylum", assay.type = "relabundance") ## ----------------------------------------------------------------------------- library(patchwork) plots <- plotAbundance(GlobalPatterns[rowData(GlobalPatterns)$Phylum %in% prev_phylum], features = "SampleType", rank = "Phylum", assay.type = "relabundance") plots$abundance / plots$SampleType + plot_layout(heights = c(9, 1)) ## ----------------------------------------------------------------------------- plotFeaturePrevalence(GlobalPatterns, rank = "Phylum", detections = c(0, 0.001, 0.01, 0.1, 0.2)) ## ----------------------------------------------------------------------------- plotPrevalentAbundance(GlobalPatterns, rank = "Family", colour_by = "Phylum") + scale_x_log10() ## ----------------------------------------------------------------------------- plotPrevalence(GlobalPatterns, rank = "Phylum", detections = c(0.01, 0.1, 1, 2, 5, 10, 20)/100, prevalences = seq(0.1, 1, 0.1)) ## ----message=FALSE------------------------------------------------------------ library(scater) library(mia) ## ----------------------------------------------------------------------------- altExp(GlobalPatterns,"Genus") <- agglomerateByRank(GlobalPatterns,"Genus") altExp(GlobalPatterns,"Genus") <- addPerFeatureQC(altExp(GlobalPatterns,"Genus")) rowData(altExp(GlobalPatterns,"Genus"))$log_mean <- log(rowData(altExp(GlobalPatterns,"Genus"))$mean) rowData(altExp(GlobalPatterns,"Genus"))$detected <- rowData(altExp(GlobalPatterns,"Genus"))$detected / 100 top_taxa <- getTopFeatures(altExp(GlobalPatterns,"Genus"), method="mean", top=100L, assay.type="counts") ## ----plot1, fig.cap="Tree plot using ggtree with tip labels decorated by mean abundance (colour) and prevalence (size)"---- plotRowTree(altExp(GlobalPatterns,"Genus")[top_taxa,], tip_colour_by = "log_mean", tip_size_by = "detected") ## ----plot2, fig.cap="Tree plot using ggtree with tip labels decorated by mean abundance (colour) and prevalence (size). Tip labels of the tree are shown as well."---- plotRowTree(altExp(GlobalPatterns,"Genus")[top_taxa,], tip_colour_by = "log_mean", tip_size_by = "detected", show_label = TRUE) ## ----plot3, fig.cap="Tree plot using ggtree with tip labels decorated by mean abundance (colour) and prevalence (size). Selected node and tip labels are shown."---- labels <- c("Genus:Providencia", "Genus:Morganella", "0.961.60") plotRowTree(altExp(GlobalPatterns,"Genus")[top_taxa,], tip_colour_by = "log_mean", tip_size_by = "detected", show_label = labels, layout="rectangular") ## ----plot4, fig.cap="Tree plot using ggtree with tip labels decorated by mean abundance (colour) and edges labeled Kingdom (colour) and prevalence (size)"---- plotRowTree(altExp(GlobalPatterns,"Genus")[top_taxa,], edge_colour_by = "Phylum", tip_colour_by = "log_mean") ## ----------------------------------------------------------------------------- data(col_graph) ## ----------------------------------------------------------------------------- plotColGraph(col_graph, altExp(GlobalPatterns,"Genus"), colour_by = "SampleType", edge_colour_by = "weight", edge_width_by = "weight", show_label = TRUE) ## ----------------------------------------------------------------------------- metadata(altExp(GlobalPatterns,"Genus"))$graph <- col_graph ## ----include=FALSE------------------------------------------------------------ plotColGraph(altExp(GlobalPatterns,"Genus"), name = "graph", colour_by = "SampleType", edge_colour_by = "weight", edge_width_by = "weight", show_label = TRUE) ## ----eval=FALSE--------------------------------------------------------------- # if(!requireNamespace("devtools", quietly = TRUE)){ # BiocManager::install("devtools") # } # if(!requireNamespace("miaTime", quietly = TRUE)){ # devtools::install_github("microbiome/miaTime", upgrade = "never") # } ## ----eval=FALSE--------------------------------------------------------------- # # Load data from miaTime package # library("miaTime") # data(SilvermanAGutData, package="miaTime") # tse <- SilvermanAGutData # tse <- transformAssay(tse, method = "relabundance") # taxa <- getTopFeatures(tse, 2) ## ----eval=FALSE--------------------------------------------------------------- # plotSeries(tse, # x = "DAY_ORDER", # y = taxa, # colour_by = "Family") ## ----eval=FALSE--------------------------------------------------------------- # plotSeries(tse[taxa,], # x = "DAY_ORDER", # colour_by = "Family", # linetype_by = "Phylum", # assay.type = "relabundance") ## ----eval=FALSE--------------------------------------------------------------- # plotSeries(tse, # x = "DAY_ORDER", # y = getTopFeatures(tse, 5), # colour_by = "Family", # linetype_by = "Phylum", # assay.type = "counts") ## ----------------------------------------------------------------------------- data(GlobalPatterns, package="mia") se <- GlobalPatterns plotColTile(se,"SampleType","Primer") + theme(axis.text.x.top = element_text(angle = 45, hjust = 0)) ## ----------------------------------------------------------------------------- data(dmn_se, package = "mia") names(metadata(dmn_se)) # plot the fit plotDMNFit(dmn_se, type = "laplace") ## ----MDS, eval=FALSE---------------------------------------------------------- # library(miaTime) # data(hitchip1006, package = "miaTime") # tse <- hitchip1006 # tse <- transformAssay(tse, method = "relabundance") # ## Ordination with PCoA with Bray-Curtis dissimilarity # tse <- runMDS(tse, FUN = vegan::vegdist, method = "bray", name = "PCoA_BC", # assay.type = "relabundance", na.rm = TRUE) # # plot # p <- plotReducedDim(tse, dimred = "PCoA_BC") # p ## ----eval=FALSE--------------------------------------------------------------- # library(dplyr) # # # List subjects with two time points # selected.subjects <- names(which(table(tse$subject)==2)) # # # Subjects counts per number of time points available in the data # table(table(tse$subject)) %>% as.data.frame() %>% # rename(Timepoints=Var1, Subjects=Freq) ## ----eval=FALSE--------------------------------------------------------------- # # plot # p + geom_path(aes(x=X1, y=X2, group=subject), # arrow=arrow(length = unit(0.1, "inches")), # # combining ordination data and metadata then selecting the subjects # # Note, scuttle::makePerCellDF could also be used for the purpose. # data = subset(data.frame(reducedDim(tse), colData(tse)), # subject %in% selected.subjects) %>% arrange(time))+ # labs(title = "All trajectories with two time points")+ # theme(plot.title = element_text(hjust = 0.5)) ## ----eval=FALSE--------------------------------------------------------------- # library(miaTime) # # calculating step wise divergence based on the microbial profiles # tse <- getStepwiseDivergence(tse, group = "subject", time_field = "time") # # retrieving the top 10% divergent subjects having two time points # top.selected.subjects <- subset(data.frame(reducedDim(tse), colData(tse)), # subject %in% selected.subjects) %>% # top_frac(0.1, time_divergence) %>% select(subject) %>% .[[1]] # # plot # p + geom_path(aes(x=X1, y=X2, # color=time_divergence, group=subject), # # the data is sorted in descending order in terms of time # # since geom_path will use the first occurring observation # # to color the corresponding segment. Without the sorting # # geom_path will pick up NA values (corresponding to initial time # # points); breaking the example. # data = subset(data.frame(reducedDim(tse), colData(tse)), # subject %in% top.selected.subjects) %>% # arrange(desc(time)), # # arrow end is reversed, due to the earlier sorting. # arrow=arrow(length = unit(0.1, "inches"), ends = "first"))+ # labs(title = "Top 10% divergent trajectories from time point one to two")+ # scale_color_gradient2(low="white", high="red")+ # theme(plot.title = element_text(hjust = 0.5)) ## ----eval=FALSE--------------------------------------------------------------- # # Get subject with the maximum total divergence # selected.subject <- data.frame(reducedDim(tse), colData(tse)) %>% # group_by(subject) %>% # summarise(total_divergence = sum(time_divergence, na.rm = TRUE)) %>% # filter(total_divergence==max(total_divergence)) %>% select(subject) %>% .[[1]] # # plot # p + geom_path(aes(x=X1, y=X2, group=subject), # data = subset(data.frame(reducedDim(tse), colData(tse)), # subject %in% selected.subject) %>% arrange(time), # arrow=arrow(length = unit(0.1, "inches")))+ # labs(title = "Longest trajectory by divergence")+ # theme(plot.title = element_text(hjust = 0.5)) ## ----------------------------------------------------------------------------- sessionInfo()