## ---- echo=FALSE, message=FALSE-------------------------------------------- library(knitr) require(knitcitations) # cleanbib() # options("citation_format" = "pandoc") bib <- read.bibtex("bibliography.bib") # Can also write math expressions, e.g. # $Y = X\beta + \epsilon$, footnotes^[A footnote here.] # > "He who gives up [code] safety for [code] speed deserves neither." # ([via](https://twitter.com/hadleywickham/status/504368538874703872)) knitr::opts_chunk$set(fig.width=10, fig.height=10, message=FALSE, warning=FALSE) ## ----install, eval=FALSE--------------------------------------------------- # library(BiocInstaller) # source("http://www.bioconductor.org/biocLite.R") # useDevel() # biocLite("microbiome") ## ----loading--------------------------------------------------------------- library(microbiome) ## ----atlasdata------------------------------------------------------------- # Data from # http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html data(atlas1006) atlas1006 ## ----dietswapdata---------------------------------------------------------- data(dietswap) # Data from http://dx.doi.org/10.1038/ncomms7342 dietswap ## ----peerjdata------------------------------------------------------------- data(peerj32) # Data from https://peerj.com/articles/32/ ## ----compositional--------------------------------------------------------- # dietswap is a phyloseq object; see above dietswap.compositional <- transform(dietswap, "compositional") ## ----global, results="asis"------------------------------------------------ g <- global(atlas1006, index = "gini") ## ----reg, fig.width=8, fig.height=4, dev="CairoPNG"------------------------ # Estimate Shannon diversity and add it to the phyloseq object sample_data(atlas1006)$diversity <- global(atlas1006, index = "shannon")[,1] # Compare age and microbiome diversity plot_regression(diversity ~ age, meta(atlas1006)) ## ----core-prevalence2------------------------------------------------------ p <- prevalence(dietswap, detection = 0, sort = TRUE) ## ----core-data------------------------------------------------------------- # Taxa with over 50% prevance at .2% relative abundance dietswap.core <- core(dietswap.compositional, detection = .2/100, prevalence = 50/100) dietswap.core ## ----corevisu, fig.width=20, fig.heigth=10, out.width = '350px'------------ library(ggplot2, quiet = TRUE) p <- plot_core(transform(dietswap.core, "compositional"), plot.type = "heatmap", colours = gray(seq(0,1,length=5)), prevalences = seq(.05, 1, .05), detections = 10^seq(log10(1e-3), log10(.2), length = 10), horizontal = TRUE) + xlab("Detection Threshold (Relative Abundance (%))") print(p) ## ----compheat, fig.height=8, fig.width=12, out.width="300px"--------------- tmp <- plot_composition(dietswap.core, plot.type = "heatmap", transform = "Z", mar = c(6, 13, 1, 1), sample.sort = "nationality") ## ----compbar, fig.height=6, fig.width=12----------------------------------- plot_composition(transform(dietswap.core, "compositional"), plot.type = "barplot", sample.sort = "neatmap") ## ----landscape4, fig.width=8, fig.height=6, out.width="400px"-------------- # Project the samples with the given method and dissimilarity measure. # Ordinate the data; note that some ordinations are sensitive to random seed # "quiet" is used to suppress intermediate outputs set.seed(423542) # Get ordination x <- dietswap.core quiet(x.ord <- ordinate(x, method = "NMDS", distance = "bray")) # Pick the projected data (first two columns + metadata) quiet(proj <- phyloseq::plot_ordination(x, x.ord, justDF=TRUE)) # Rename the projection axes names(proj)[1:2] <- paste("Comp", 1:2, sep=".") p <- plot_landscape(proj[, 1:2], col = proj$nationality, legend = TRUE) print(p) ## ----heatmap-crosscorrelate4, fig.keep='none'------------------------------ # Define data sets to cross-correlate data(peerj32) x <- log10(peerj32$microbes) # OTU Log10 (44 samples x 130 genera) y <- as.matrix(peerj32$lipids) # Lipids (44 samples x 389 lipids) # Cross correlate data sets and return a table correlation.table <- associate(x, y, method = "spearman", mode = "table", p.adj.threshold = 0.05, n.signif = 1) kable(head(correlation.table)) ## ----hm3, fig.width=12, fig.height=10, out.width="400px"------------------- heat(correlation.table, "X1", "X2", fill = "Correlation", star = "p.adj", p.adj.threshold = 0.05) ## ----varpl2, fig.show='hold', out.width="250px"---------------------------- # Use relative abundances and plot subject variation # during the follow-up period pseq <- transform(atlas1006, "compositional") plot_tipping(pseq, "Dialister", tipping.point = 0.004) # Show the bimodal population distribution hotplot(pseq, "Dialister", tipping.point = 0.004) ## ---- echo=FALSE,results="asis"-------------------------------------------- bibliography()