## ---- echo=FALSE, results="asis", message=FALSE, KnitrSetUp------------------- knitr::opts_chunk$set(tidy=FALSE,warning=FALSE,message=FALSE) Biocpkg <- function (pkg){ sprintf("[%s](http://bioconductor.org/packages/%s)", pkg, pkg) } CRANpkg <- function(pkg){ cran <- "https://CRAN.R-project.org/package" fmt <- "[%s](%s=%s)" sprintf(fmt, pkg, cran, pkg) } ## ---- echo=FALSE, results="hide", message=FALSE, Loadpackages----------------- library(ggplot2) library(DT) library(tidyverse) library(phyloseq) library(ggtree) library(treeio) library(tidytree) library(MicrobiotaProcess) ## ---- error=FALSE, importFunction--------------------------------------------- # import data from dada2 pipeline. seqtabfile <- system.file("extdata", "seqtab.nochim.rds", package="MicrobiotaProcess") seqtab <- readRDS(seqtabfile) taxafile <- system.file("extdata", "taxa_tab.rds",package="MicrobiotaProcess") seqtab <- readRDS(seqtabfile) taxa <- readRDS(taxafile) # the seqtab and taxa are output of dada2 sampleda <- system.file("extdata", "mouse.time.dada2.txt", package="MicrobiotaProcess") ps_dada2 <- import_dada2(seqtab=seqtab, taxatab=taxa, sampleda=sampleda) ps_dada2 # import data from qiime2 pipeline otuqzafile <- system.file("extdata", "table.qza", package="MicrobiotaProcess") taxaqzafile <- system.file("extdata", "taxa.qza", package="MicrobiotaProcess") mapfile <- system.file("extdata", "metadata_qza.txt", package="MicrobiotaProcess") ps_qiime2 <- import_qiime2(otuqza=otuqzafile, taxaqza=taxaqzafile, mapfilename=mapfile) ps_qiime2 ## ---- error=FALSE, fig.align="center", fig.height=3.2, fig.width=7.5, rarefaction---- # for reproducibly random number set.seed(1024) p_rare <- ggrarecurve(obj=ps_dada2, indexNames=c("Observe","Chao1","ACE"), chunks=300) + theme(legend.spacing.y=unit(0.02,"cm"), legend.text=element_text(size=6)) p_rare ## ---- error=FALSE, fig.width=7.5, fig.height=3.2, fig.align="center", alphaindex---- alphaobj <- get_alphaindex(ps_dada2) head(as.data.frame(alphaobj)) p_alpha <- ggbox(alphaobj, geom="violin", factorNames="time") + scale_fill_manual(values=c("#2874C5", "#EABF00"))+ theme(strip.background = element_rect(colour=NA, fill="grey")) p_alpha ## ---- error=FALSE, fig.align="center", fig.height=4.5, fig.width=7, taxabar---- # relative abundance otubar <- ggbartax(obj=ps_dada2) + xlab(NULL) + ylab("relative abundance (%)") otubar ## ---- error=FALSE, fig.align="center", fig.height=4.5, fig.width=7, phylumAbundance---- phytax <- get_taxadf(obj=ps_dada2, taxlevel=2) phybar <- ggbartax(obj=phytax) + xlab(NULL) + ylab("relative abundance (%)") phybar ## ---- error=FALSE, fig.align="center", fig.height=4.5, fig.width=7, classAbundance---- phybar2 <- ggbartax(obj=phytax, facetNames="time", count=TRUE) + xlab(NULL) + ylab("abundance") phybar2 classtax <- get_taxadf(obj=ps_dada2, taxlevel=3) classbar <- ggbartax(obj=classtax, facetNames="time", count=FALSE) + xlab(NULL) + ylab("relative abundance (%)") classbar ## ---- error=FALSE, fig.align="center", fig.height=4.5, fig.width=6, ordanalysis---- pcares <- get_pca(obj=ps_dada2, method="hellinger") # Visulizing the result pcaplot <- ggordpoint(obj=pcares, biplot=TRUE, speciesannot=TRUE, factorNames=c("time"), ellipse=TRUE) + scale_color_manual(values=c("#2874C5", "#EABF00")) pcaplot pcoares <- get_pcoa(obj=ps_dada2, distmethod="euclidean", method="hellinger") # Visualizing the result pcoaplot <- ggordpoint(obj=pcoares, biplot=TRUE, speciesannot=TRUE, factorNames=c("time"), ellipse=TRUE) + scale_color_manual(values=c("#2874C5", "#EABF00")) pcoaplot ## ---- fig.align="center", fig.height=5, fig.width=6, error=FALSE, hclustAnalysis---- hcsample <- get_clust(obj=ps_dada2, distmethod="euclidean", method="hellinger", hclustmethod="average") # rectangular layout clustplot1 <- ggclust(obj=hcsample, layout = "rectangular", pointsize=1, fontsize=0, factorNames=c("time")) + scale_color_manual(values=c("#2874C5", "#EABF00")) + theme_tree2(legend.position="right", plot.title = element_text(face="bold", lineheight=25,hjust=0.5)) clustplot1 # circular layout clustplot2 <- ggclust(obj=hcsample, layout = "circular", pointsize=1, fontsize=2, factorNames=c("time")) + scale_color_manual(values=c("#2874C5", "#EABF00")) + theme(legend.position="right") clustplot2 ## ---- echo=FALSE-------------------------------------------------------------- sessionInfo()