# Set up global options for nice reports and keeping figures:
knitr::opts_chunk$set(fig.width=14, fig.height=8, fig.align="center",
                      warning=FALSE, message=FALSE)

library(phyloseq)
library(vegan)
library(ggplot2)
library(kableExtra)
library(tidyr)

nice_colors = c("#999999", "#E69F00", "#56B4E9","#e98756","#c08160","#5800e6", "#CDDC49", "#C475D3", 
                "#E94B30", "#233F57", "#FEE659", "#A1CFDD", "#F4755E", "#D6F6F7","#EB6D58", "#6898BF")

load(file=file.path("rdata_objects", "initial_rooted.Rdata"))

df = data.frame(ASVs=rowSums(otu_table(ps)>0), reads=sample_sums(ps), sample_data(ps))

ggplot(df, aes(x=reads)) + geom_histogram(bins=50, color='black', fill='grey') + 
  theme_bw() +  geom_vline(xintercept=10000, color= "red", linetype='dashed') +
  labs(title="Histogram: Reads per Sample") + xlab("Read Count") + ylab("Sample Count")

ggplot(df, aes(x = Group, y = reads, color = ASVs)) +
  geom_boxplot(color="black") + theme_bw() +
  geom_jitter(width=.2, height=0)  + 
  theme(axis.text.x = element_text(angle = 90)) +
  geom_hline(yintercept=10000, color= "purple", linetype='dashed') +
  geom_hline(yintercept=1000, color= "red", linetype='dashed') +
  ggtitle("Reads sequenced by SampleGroup")

out = rarecurve(otu_table(ps), col=nice_colors, step=100 , lwd=2, ylab="ASVs", label=F,
              main="Rarefaction Curve for all samples")

# We use the rarefaction curve data produce by vegan above
names(out) = rownames(otu_table(ps))

# Coerce data into "long" form.
protox <- mapply(FUN = function(x, y) {
  mydf <- as.data.frame(x)
  colnames(mydf) <- "value"
  mydf$SampleID <- y
  mydf$subsample <- attr(x, "Subsample")
  mydf
}, x = out, y = as.list(names(out)), SIMPLIFY = FALSE)

xy <- do.call(rbind, protox)
rownames(xy) <- NULL  # pretty
xy = data.frame(xy, 
  sample_data(ps)[match(xy$SampleID, rownames(sample_data(ps))), ])
# Plot Rarefaction curve
ggplot(xy, aes(x = subsample, y = value, color = SampleID)) +
  theme_bw() + 
  scale_color_discrete(guide = FALSE) +  # turn legend on or off
  geom_line() +
  geom_vline(xintercept=10000, color= "red", linetype='dashed') + 
  labs(title="Rarefaction curves") + xlab("Sequenced Reads") + ylab('ASVs Detected')

# Plot Rarefaction curve
ggplot(xy, aes(x = subsample, y = value, color = SampleID)) +
  theme_bw() + 
  scale_color_discrete(guide = FALSE) +  # turn legend on or off
  geom_line() +
  facet_wrap(~Group) + 
  geom_vline(xintercept=10000, color= "red", linetype='dashed') + 
  labs(title="Rarefaction curves by Group") + xlab("Sequenced Reads") + ylab('ASVs Detected')

# Plot Rarefaction curve, faceting by group
ggplot(xy, aes(x = subsample, y = value, color = SampleID)) +
  theme_bw() + 
  scale_color_discrete(guide = FALSE) +  # turn legend on or off
  geom_line() +
  facet_wrap(~Group) + 
  geom_vline(xintercept=10000, color= "red", linetype='dashed') + 
  labs(title="Rarefaction curves by Group") + xlab("Sequenced Reads") + ylab('ASVs Detected')

df = data.frame(names=names(refseq(ps)), ASV.Length = Biostrings::width(refseq(ps)))
ggplot(df, aes(x=ASV.Length)) + geom_histogram(binwidth=1, color='black', fill='grey') + 
  theme_bw() + labs(title="Histogram: Amplicon Lengths")

readsPerSample = rowSums(otu_table(ps))
fractionReadsAssigned = sapply(colnames(tax_table(ps)), function(x){
  rowSums(otu_table(ps)[, !is.na(tax_table(ps))[,x]]) / readsPerSample
})

fractionReadsAssigned = data.frame(SampleID = rownames(fractionReadsAssigned), fractionReadsAssigned)
fractionReadsAssigned.L = pivot_longer(fractionReadsAssigned, 
    cols=colnames(tax_table(ps)), names_to="taxlevel", values_to="fractionReadsAssigned")
fractionReadsAssigned.L = data.frame(fractionReadsAssigned.L, sample_data(ps)[fractionReadsAssigned.L$SampleID,])

fractionReadsAssigned.L$taxlevelf = factor(fractionReadsAssigned.L$taxlevel, levels=c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))

# Boxplot, fraction assigned by SampleType
ggplot(fractionReadsAssigned.L, aes(y = fractionReadsAssigned, x = taxlevelf, fill=Group, color=Replicate)) +
  theme_bw() +
  geom_boxplot(color='black', outlier.shape=NA) +  
  ggtitle("Fraction of reads identified by taxonomic level") +
  geom_point(position=position_jitterdodge(jitter.width=0, jitter.height=0), size=.5)

download.file("https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2021-May-Microbial-Community-Analysis/master/data_analysis/mca_part3.Rmd", "mca_part3.Rmd")

sessionInfo()