---

title: “Microbiome Differential Abundance Testing Tutorial” Dr. upasna srivastava output: html_notebook

#for this study we can use The 16S sequence data have been submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi) under accession number SRP000383.

#Import data with phyloseq, convert to DESeq2
library("phyloseq"); 
packageVersion("phyloseq")
#Defined file path, and import the published OTU count data into R.
filepath = system.file("extdata", "study_1457_split_library_seqs_and_mapping.zip", package="phyloseq")
colon = microbio_me_qiime(filepath)
#in case if you qiime id then you could directly download the data from the microbio.me/qiime server using this code
colon = microbio_me_qiime(1457)
colon
library("DESeq2"); 
packageVersion("DESeq2")
colon <- subset_samples(colon, DIAGNOSIS != "None")
colon <- prune_samples(sample_sums(colon) > 500, colon)
colon
diagdds = phyloseq_to_deseq2(colon, ~ 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")
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(colon)[rownames(sigtab), ], "matrix"))
head(sigtab)
#OTUs that were significantly enriched
posigtab = sigtab[sigtab[, "log2FoldChange"] > 0, ]
posigtab = posigtab[, c("baseMean", "log2FoldChange", "lfcSE", "padj", "Phylum", "Class", "Family", "Genus")]
## ggplot2 commands to show bar plot showing the log2-fold-change, showing Genus and Phylum. 
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))
  dev.off()