inst/doc/MicrobiotaProcess-basics.R

In MicrobiotaProcess: an R package for analysis, visualization and biomarker discovery of microbiome

## ---- 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_colour_manual(values=c("#2874C5", "#EABF00")) +
           scale_fill_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_colour_manual(values=c("#2874C5", "#EABF00")) +
            scale_fill_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()