library(ggplot2)
library(plyr)
library(dplyr)
library(phyloseq)

#Load the previously generated phyloseq object using `readRDS`
physeq = readRDS("otu_physeq.rds")
physeq

# For future referencing, create a dataframe of just the sample data from your phyloseq object
physeq.m = sample_data(physeq)
physeq.m

# We can use the `sample_sums()` function of phyloseq to add up total read counts for each sample.
sample_sum_df <- data.frame(sum = sample_sums(physeq))

# Histogram of sample read counts
ggplot(sample_sum_df, aes(x = sum)) + 
  geom_histogram(color = "black", fill = "indianred", binwidth = 2500) +
  ggtitle("Distribution of sample sequencing depth") + 
  xlab("Read counts") +
  theme(axis.title.y = element_blank())

# Originally, the "Age variable was read as numerical, while I want it to be categorical. I use the `as.factor()` function to do this."
sample_data(physeq)$Age <- as.factor(sample_data(physeq)$Age)

#Phyloseq contains the `plot_richness()` function to display multiple alpha diversity measures at once. You can modify the plots based on sample metadata as well.
plot_richness(physeq, x="Site", measures=c("Observed", "Shannon", "Simpson", "Chao1", "InvSimpson"), color="Age") +   theme_bw()

# First we will create a dataframe containing all diversity measures from above using the `estimate_richness()` function.
adiv <- estimate_richness(physeq, measures = c("Observed", "Shannon", "Simpson", "Chao1", "InvSimpson"))

# Then we can run an analysis of variance test using the `aov()` function.
summary(aov(adiv$Observed ~ physeq.m$Site))

# First, we normalize the sequence counts by converting from raw abundance to relative abundance. This removes any bias due to total sequence counts per sample.
pn = transform_sample_counts(physeq, function(x) 100 * x/sum(x))

# Next, we use the `distance()` function from phyloseq to generate a distance matrix from our phyloseq object. You can select from multiple methods; we will use Bray-Curtis
iDist <- distance(pn, method = "bray")

# Using the distance matrix created above, we now make an NMDS ordination using the `ordinate()` function.
pn.nmds = ordinate(pn, 
                method = "NMDS", 
                distance = iDist)

# Finally, we create an plot of the previous ordination using the `plot_ordination()` function. The `justDF` option is set to `TRUE`, which indicates that we only want the dataframe created with plot ordination to be returned, not a ggplot object.
plot.pn.nmds = plot_ordination(pn, pn.nmds, justDF = TRUE)

# Next we can display the figure using the `print(ggplot())` command, with any ggplot aesthetic arguments you want. 
print(ggplot(plot.pn.nmds, aes(x = NMDS1, y = NMDS2, color = Site, shape = Age)) + 
        geom_point( size = 4, alpha = 0.75))

library(vegan)
p.df = as(sample_data(pn), "data.frame")
p.d = distance(pn, method = "bray")
p.adonis = adonis(p.d ~ Site + Age + Site*Age, p.df)
p.adonis

# This normalization is redundant from before, but usefull if you want to just copy this chunk of code
pn = transform_sample_counts(physeq, function(x) 100 * x/sum(x))

# Here we use the `tax_glom()` function from phyloseq to collapse our OTU table at the Phylum level. This takes all OTUs under a given phylum and combines them, adding up the total sequence counts. 
phylum = tax_glom(pn, taxrank="Phylum")

# Next, we want to look at just the average community by Site, not each individual sample. We will use the `merge_samples()` function to do this. You can use any categorical metadata variable for ths merge. 
pm  = merge_samples(phylum, "Site")

# Phyloseq does some odd things with sample names when you merge, so this next line of code makes sure that everything will be labelled correctly for our figures.
sample_data(pm)$DeRep <- levels(sample_data(pm)$Sample)

# Again, we make a sample_data dataframe for future reference
pm.m = pm %>% sample_data

# Since each 'sample' in this new object is an agglomerate sample, we need to re-normalize our sequence counts for each OTU. 
pm = transform_sample_counts(pm, function(x) 100 * x/sum(x))

# Finally, we can use the `plot_bar()` function to display our barchart.
plot_bar(pm, "Sample", fill = "Phylum")

# Here, we will take a look at the top 5 most abundant phyla using the a faceted plot.

# First, create a list of the top 5 phyla
trimlist = names(sort(taxa_sums(pm), TRUE)[1:5])

# Then we will use the `prune_taxa()` function to select only those top 5 phyla
pm.5 = prune_taxa(trimlist, pm)
    
# Finally, we can plot the abundance of those phyla using the `facet_grid` argument for `plot_bar()`
plot_bar(pm.5, "Sample", fill = "Phylum", facet_grid = ~Phylum)

# Create a new `phyloseq` object with only the BH and MP sites.
mb = subset_samples(physeq, Site == "MP" | Site == "BH") %>%
                      filter_taxa(function(x) sum(x) > 0, TRUE)

# Next, we need to check the order of our Site levels. DESeq2 takes the first level as the 'Control' and the second level as the 'Treatment', and this is needed for downstream interpretation of results.
head(sample_data(mb)$Site)

# Now we can prepare the DESeq2 object and run our analysis.
library(DESeq2)

# First, we will convert the sequence count data from our phyloseq object into the proper format for DESeq2
mb.dds <- phyloseq_to_deseq2(mb, ~ Site)

# Next we need to estimate the size factors for our sequences
gm_mean = function(x, na.rm=TRUE){
    exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(mb.dds), 1, gm_mean)
mb.dds = estimateSizeFactors(mb.dds, geoMeans = geoMeans)

# Finally, we can run the core DESeq2 algorithm
mb.dds = DESeq(mb.dds, fitType="local")

# Create a new dataframe from the results of our DESeq2 run
mb.res = results(mb.dds)

# Reorder the sequences by their adjusted p-values
mb.res = mb.res[order(mb.res$padj, na.last=NA), ]

# Set your alpha for testing significance, and filter out non-significant results
alpha = 0.01
mb.sigtab = mb.res[(mb.res$padj < alpha), ]

# Add taxonomy information to each sequence
mb.sigtab = cbind(as(mb.sigtab, "data.frame"), as(tax_table(mb)[rownames(mb.sigtab), ], "matrix"))
head(mb.sigtab)

# Set ggplot2 options
theme_set(theme_bw(base_size = 20))
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
# Rearrange the order of our OTUs by Phylum.
x = tapply(mb.sigtab$log2FoldChange, mb.sigtab$Phylum, function(x) max(x))
x = sort(x, TRUE)
mb.sigtab$Phylum = factor(as.character(mb.sigtab$Phylum), levels=names(x))

# Create and display the plot
p = ggplot(mb.sigtab, aes(x=Phylum, y=log2FoldChange, color=Phylum)) + geom_jitter(size=4, width=0.25) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))
p