inst/doc/phyloseq-mixture-models.R

In phyloseq: Handling and analysis of high-throughput microbiome census data

## ----load-phyloseq, message=FALSE, warning=FALSE------------------------------
library("phyloseq"); packageVersion("phyloseq")

## ----filepath-----------------------------------------------------------------
filepath = system.file("extdata", "study_1457_split_library_seqs_and_mapping.zip", package="phyloseq")
kostic = microbio_me_qiime(filepath)

## ----example-path-local, eval=FALSE-------------------------------------------
#  filepath = "~/Downloads/study_1457_split_library_seqs_and_mapping.zip"
#  kostic = microbio_me_qiime(filepath)

## ----example-path-remote, eval=FALSE------------------------------------------
#  kostic = microbio_me_qiime(1457)

## ----show-variables-----------------------------------------------------------
kostic
head(sample_data(kostic)$DIAGNOSIS, 10)

## ----deseq2, message=FALSE, warning=FALSE-------------------------------------
library("DESeq2"); packageVersion("DESeq2")

## ----rm-bad-samples-----------------------------------------------------------
kostic <- subset_samples(kostic, DIAGNOSIS != "None")
kostic <- prune_samples(sample_sums(kostic) > 500, kostic)
kostic

## ----run-deseq2---------------------------------------------------------------
diagdds = phyloseq_to_deseq2(kostic, ~ DIAGNOSIS)
# calculate geometric means prior to estimate size factors
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(diagdds), 1, gm_mean)
diagdds = estimateSizeFactors(diagdds, geoMeans = geoMeans)
diagdds = DESeq(diagdds, fitType="local")

## ----grab-results-process-table-----------------------------------------------
res = results(diagdds)
res = res[order(res$padj, na.last=NA), ]
alpha = 0.01
sigtab = res[(res$padj < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(kostic)[rownames(sigtab), ], "matrix"))
head(sigtab)

## ----table-prelim-------------------------------------------------------------
posigtab = sigtab[sigtab[, "log2FoldChange"] > 0, ]
posigtab = posigtab[, c("baseMean", "log2FoldChange", "lfcSE", "padj", "Phylum", "Class", "Family", "Genus")]

## ----make-markdown-table, echo=FALSE, results='asis'--------------------------
# Make a markdown table
posigtab = data.frame(OTU=rownames(posigtab), posigtab)
cat(paste(colnames(posigtab), collapse=" | "), fill=TRUE)
cat(paste(rep("---", times=ncol(posigtab)), collapse=" | "), fill=TRUE)
dummy = apply(posigtab, 1, function(x){
  cat(paste(x, collapse=" | "), fill=TRUE)
})

## ----bar-plot-----------------------------------------------------------------
library("ggplot2")
theme_set(theme_bw())
sigtabgen = subset(sigtab, !is.na(Genus))
# Phylum order
x = tapply(sigtabgen$log2FoldChange, sigtabgen$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Phylum = factor(as.character(sigtabgen$Phylum), levels=names(x))
# Genus order
x = tapply(sigtabgen$log2FoldChange, sigtabgen$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Genus = factor(as.character(sigtabgen$Genus), levels=names(x))
ggplot(sigtabgen, aes(y=Genus, x=log2FoldChange, color=Phylum)) + 
  geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
  geom_point(size=6) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))