
## ----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, 25)


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


## ----run-deseq2----------------------------------------------------------
diagdds = phyloseq_to_deseq2(kostic, ~ DIAGNOSIS)
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())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
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(x=Genus, y=log2FoldChange, color=Phylum)) + geom_point(size=6) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))


