#Include libraries
library(ggplot2)
library(vegan)
library(dplyr)
library(phyloseq)
library(randomForest)
library(rfUtilities)
library(rfPermute)
library(knitr)
library("data.table"); packageVersion("data.table")

#Change plots theme
theme_set(theme_bw())

#Set Working Directory
setwd("C:/path/RcodeBurns")

#Create a Phyloseq object for the BIOM table
physeq <- import_biom('miseq.biom', 'miseq.tree.phy')

colnames(tax_table(physeq)) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
sample_variables(physeq)

#The next section applies if you have a separate metadata file for variables

#Import metadata dataset into R
metadata.springfall2018 <- read.delim("C:/path/RcodeBurns/springfall2018-metadata1.txt", header=FALSE, colClasses = "character", row.names=1, stringsAsFactor = TRUE, comment.char="#")
colnames(metadata.springfall2018) <- c("PCRProgram","Plot", "Treatment", "Tree", "SamplingTime", "Burn", "Random", "TreeNum", "SampNum", "PCRNum", "Sample", "Treat", "Rep", "MSB") #Name the variables in the metadata file
View(metadata.springfall2018)

sample_data(metadata.springfall2018) #creates a phyloseq object for the data
sam = sample_data(metadata.springfall2018) #prepares the data for merge

sample_names(physeq)

tree = phy_tree(physeq)
tax  = tax_table(physeq)
otu  = otu_table(physeq)
otutax = phyloseq(otu, tax)
otutax

myData = merge_phyloseq(otutax, sam, tree) #merge
print(myData)

sample_variables(myData) #check variables

#Remove mock community, controls, and other unwanted samples
myData2 = subset_samples(myData, sample_names(myData) != "036370")
myData2 = subset_samples(myData2, sample_names(myData2) != "036371")
myData2 = subset_samples(myData2, sample_names(myData2) != "036372")
myData2 = subset_samples(myData2, sample_names(myData2) != "036373")
myData2 = subset_samples(myData2, sample_names(myData2) != "036369")
myData2 = subset_samples(myData2, sample_names(myData2) != "036368")
myData2 = subset_samples(myData2, sample_names(myData2) != "036283")
myData2 = subset_samples(myData2, sample_names(myData2) != "039657")
myData2 = subset_samples(myData2, sample_names(myData2) != "039658-Fall2018-MockCommunity-Plot-MC-Tree-DropPCR")
myData2 = subset_samples(myData2, sample_names(myData2) != "041307")
myData2 = subset_samples(myData2, sample_names(myData2) != "041308")
myData2 = subset_samples(myData2, sample_names(myData2) != "041394")
myData2 = subset_samples(myData2, sample_names(myData2) != "036214") #Pine
myData2 = subset_samples(myData2, sample_names(myData2) != "036299") #Pine
print(myData2)

# prune taxa which were only present in removed samples
myData3 <- prune_taxa(taxa_sums(myData2) > 0, myData2)
print(myData3)

#Remove replicate samples
myData4 = subset_samples(myData3, Rep=="0")
myData4 <- prune_taxa(taxa_sums(myData4) > 0, myData4)  
print(myData4)

#Merge Samples that were replicated with different PCR programs
myData5 = merge_samples(myData4, "Sample")
print(myData5)

OTU = t( otu_table(myData5) ) #Transposes the OTU table so taxa_are_rows = T
tree = phy_tree(myData5)
tax  = tax_table(myData5)
otutax = phyloseq(OTU, tax)
sam = sample_data(myData5)
myData6 = merge_phyloseq(otutax, sam, tree)
print(myData6)

#View a historgram of the Sequencing Depth
sdt = data.table(as(sample_data(myData6), "data.frame"),
                 TotalReads = sample_sums(myData6), keep.rownames = TRUE)
setnames(sdt, "rn", "SampleID")
pSeqDepth = ggplot(sdt, aes(TotalReads)) + geom_histogram() + ggtitle("Sequencing Depth")+ scale_x_continuous(name="Number of Reads", limits=c(0, 120000))
pSeqDepth

#Prune samples with low number of reads
myData7 = prune_samples(names(which(sample_sums(myData6) >= 24000)), myData6)
print(myData7)

#Subset samples by the various treatments
burnConShelClear = subset_samples(myData7, Treat != "4")
burnConShelClear <- prune_taxa(taxa_sums(burnConShelClear) > 0, burnConShelClear) # prune taxa which were only present in removed samples 
burnControlShelt = subset_samples(burnConShelClear, Treat != "2")
burnControlShelt <- prune_taxa(taxa_sums(burnControlShelt) > 0, burnControlShelt)
burnControl = subset_samples(burnControlShelt, Treat != "5")
burnControl <- prune_taxa(taxa_sums(burnControl) > 0, burnControl)
Control = subset_samples(burnControl, Treat != "1")
print(burnConShelClear)
print(burnControlShelt)
print(burnControl)

#Revert burnControl back to character vector variables
sample_data(burnControl)$Burn <- as.character(sample_data(burnControl)$Burn)
sample_data(burnControl)$Plot <- as.character(sample_data(burnControl)$Plot)
sample_data(burnControl)$Random <- as.character(sample_data(burnControl)$Random)
sample_data(burnControl)$TreeNum <- as.character(sample_data(burnControl)$TreeNum)
sample_data(burnControl)$SampNum <- as.character(sample_data(burnControl)$SampNum)
sample_data(burnControl)$Treat <- as.character(sample_data(burnControl)$Treat)
sample_data(burnControl)$Rep <- as.character(sample_data(burnControl)$Rep)

#Subset random samples and tree sampling regimes
BC.norand = subset_samples(burnControl, Random == "0")
print(BC.norand)
BC.norand <- prune_taxa(taxa_sums(BC.norand) > 0, BC.norand)
print(BC.norand)

BC.rand = subset_samples(burnControl, Random == "1")
print(BC.rand)
BC.rand <- prune_taxa(taxa_sums(BC.rand) > 0, BC.rand)
print(BC.rand)

#Subset samples by burn and control
BC.noburn = subset_samples(burnControl, Burn == "0")
print(BC.noburn)
BC.noburn <- prune_taxa(taxa_sums(BC.noburn) > 0, BC.noburn)
print(BC.noburn)

BC.burn = subset_samples(burnControl, Burn == "1")
print(BC.burn)
BC.burn <- prune_taxa(taxa_sums(BC.burn) > 0, BC.burn)
print(BC.burn)

#Subset burns and control by time of year
BC.spring = subset_samples(burnControl, SampNum=="1")
BC.spring <- prune_taxa(taxa_sums(BC.spring) > 0, BC.spring)
BC.spring

BC.summer = subset_samples(burnControl, SampNum=="2")
BC.summer <- prune_taxa(taxa_sums(BC.summer) > 0, BC.summer)
BC.summer

BC.fall = subset_samples(burnControl, SampNum=="3")
BC.fall <- prune_taxa(taxa_sums(BC.fall) > 0, BC.fall)
BC.fall

#Subset burns by time of year
B.spring = subset_samples(BC.burn, SampNum =="1")
B.spring <- prune_taxa(taxa_sums(B.spring) > 0, B.spring)
B.spring

B.summer = subset_samples(BC.burn, SampNum =="2")
B.summer <- prune_taxa(taxa_sums(B.summer) > 0, B.summer)
B.summer

B.fall = subset_samples(BC.burn, SampNum =="3")
B.fall <- prune_taxa(taxa_sums(B.fall) > 0, B.fall)
B.fall

NB.spring = subset_samples(BC.noburn, SampNum =="1")
NB.spring <- prune_taxa(taxa_sums(NB.spring) > 0, NB.spring)
NB.spring

NB.summer = subset_samples(BC.noburn, SampNum =="2")
NB.summer <- prune_taxa(taxa_sums(NB.summer) > 0, NB.summer)
NB.summer

NB.fall = subset_samples(BC.noburn, SampNum =="3")
NB.fall <- prune_taxa(taxa_sums(NB.fall) > 0, NB.fall)
NB.fall


#Subset random plots vs. tree plots for Burn-Control Spring & Fall
BC.fall.random = subset_samples(BC.fall, Random == "1")
BC.fall.trees = subset_samples(BC.fall, Random == "0")
print(BC.fall.random)
print(BC.fall.trees)
BC.fall.random <- prune_taxa(taxa_sums(BC.fall.random) > 0, BC.fall.random)
print(BC.fall.random)
BC.fall.trees <- prune_taxa(taxa_sums(BC.fall.trees) > 0, BC.fall.trees)
print(BC.fall.trees)

BC.summer.random = subset_samples(BC.summer, Random == "1")
BC.summer.trees = subset_samples(BC.summer, Random == "0")
print(BC.summer.random)
print(BC.summer.trees)
BC.summer.random <- prune_taxa(taxa_sums(BC.summer.random) > 0, BC.summer.random)
print(BC.summer.random)
BC.summer.trees <- prune_taxa(taxa_sums(BC.summer.trees) > 0, BC.summer.trees)
print(BC.summer.trees)

# Subset burns by random and tree samples
B.fall.random = subset_samples(B.fall, Random == "1")
B.fall.trees = subset_samples(B.fall, Random == "0")
print(B.fall.random)
print(B.fall.trees)
B.fall.random <- prune_taxa(taxa_sums(B.fall.random) > 0, B.fall.random)
print(B.fall.random)
B.fall.trees <- prune_taxa(taxa_sums(B.fall.trees) > 0, B.fall.trees)
print(B.fall.trees)

B.summer.random = subset_samples(B.summer, Random == "1")
B.summer.trees = subset_samples(B.summer, Random == "0")
print(B.summer.random)
print(B.summer.trees)
B.summer.random <- prune_taxa(taxa_sums(B.summer.random) > 0, B.summer.random)
print(B.summer.random)
B.summer.trees <- prune_taxa(taxa_sums(B.summer.trees) > 0, B.summer.trees)
print(B.summer.trees)


### Historgram of the Sequencing Depth ###
sdt = data.table(as(sample_data(burnControl), "data.frame"),
                 TotalReads = sample_sums(burnControl), keep.rownames = TRUE)
setnames(sdt, "rn", "SampleID")
pSeqDepth = ggplot(sdt, aes(TotalReads)) + geom_histogram() + ggtitle("Sequencing Depth") + scale_x_continuous(name="Number of Reads", limits=c(24000, 120000)) + scale_y_continuous(name="Number of Samples", limits=c(0, 6))
pSeqDepth

burn.labels <- c("Control", "Burn")
names(burn.labels) <- c("3", "1")

season.labels <- c("Spring", "Summer", "Fall")
names(season.labels) <- c("1", "2", "3")

pSeqDepth + facet_wrap(~Treat, labeller = labeller(Treat = burn.labels))
pSeqDepth + facet_wrap(~SampNum, labeller = labeller(SampNum = season.labels))

data(sdt)
mean(sdt[["TotalReads"]])
median(sdt[["TotalReads"]])

### Rarefaction Curve for burns and controls ###
rarecurve(t(otu_table(burnControl)), step=50, cex=0.5, label=FALSE)
rarefyBC = rarefy_even_depth(burnControl, rngseed = 711)
rarecurve(t(otu_table(rarefyBC)), step=50, cex=0.5, label=FALSE, xlab = "Number of Reads", ylab = "Number of OTUs")

rarefyBC.burn = rarefy_even_depth(BC.burn, rngseed = 711)
print(rarefyBC)
rarefyBC.noburn = rarefy_even_depth(BC.noburn, rngseed = 711)
rarefyBC.nbspr = rarefy_even_depth(NB.spring, rngseed = 711)
rarefyBC.nbsum = rarefy_even_depth(NB.summer, rngseed = 711)
rarefyBC.nbfall = rarefy_even_depth(NB.fall, rngseed = 711)
rarefyBC.bspr = rarefy_even_depth(B.spring, rngseed = 711)
rarefyBC.bsum = rarefy_even_depth(B.summer, rngseed = 711)
rarefyBC.bfall = rarefy_even_depth(B.fall, rngseed = 711)

rarefyBC.rand = rarefy_even_depth(BC.rand, rngseed = 711)
print(rarefyBC.rand)
rarefyBC.norand = rarefy_even_depth(BC.norand, rngseed = 711)
print(rarefyBC.norand)
rarefyBC.sf = rarefy_even_depth(BC.springfall, rngseed = 711)
rarefyBC.f = rarefy_even_depth(BC.fall, rngseed = 711)
rarefyBC.spr = rarefy_even_depth(BC.spring, rngseed = 711)
rarefyBC.sum = rarefy_even_depth(BC.summer, rngseed = 711)
rarefyBC.f.r = rarefy_even_depth(BC.fall.random, rngseed = 711)
rarefyBC.f.t = rarefy_even_depth(BC.fall.trees, rngseed = 711)
rarefyBC.sum.r = rarefy_even_depth(BC.summer.random, rngseed = 711)
rarefyBC.sum.t = rarefy_even_depth(BC.summer.trees, rngseed = 711)

### Species accumulation curve for burns and controls ###

#Convert the phyloseq object to vegan-compatible 
# http://joey711.github.io/phyloseq-demo/phyloseq-demo.html

veganotu <- function(physeq) {
  require("vegan")
  OTU <- otu_table(physeq)
  if (taxa_are_rows(OTU)) {
    OTU <- t(OTU)
  }
  return(as(OTU, "matrix"))
}

#Plot a species accumulation curve for the total sampling
vegan_BC = veganotu(rarefyBC)
print(vegan_BC)

#Create a species accumulation curve - Total Burn & Control
sp1 <- specaccum(vegan_BC, "random")
sp1
summary(sp1)
plot(sp1, ci.type="poly", col="blue", lwd=2, ci.lty=0, ci.col="lightblue", ylab = "Number of OTUs", xlab = "Number of Samples") 



##### Alpha Diversity Measures #####

#Alpha diversity for all Burn & Control Samples
alpha_chart = c("Observed", "Chao1", "InvSimpson")
estimate_richness(rarefyBC, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.burn, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.noburn, split = FALSE, measures = alpha_chart)

print(rarefyBC.f.t)

#Alpha diveristy T tests used in paper

erich <- estimate_richness(rarefyBC, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC)$Random)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.sum, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.sum)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.f, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.f)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.rand, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.rand)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.norand, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.norand)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.sum.r, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.sum.r)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.f.r, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.f.r)$Burn)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.sf, measures = c("Observed", "Chao1", "InvSimpson", ))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.sf)$PCRProgram)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.spr, measures = c("Observed", "Chao1", "InvSimpson", ))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.spr)$PCRProgram)[c("estimate","p.value","statistic","conf.int")])))
ttest

erich <- estimate_richness(rarefyBC.f, measures = c("Observed", "Chao1", "InvSimpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.f)$PCRProgram)[c("estimate","p.value","statistic","conf.int")])))
ttest

#Alpha diveristy plots used in paper

alpha_graph = c("Observed", "InvSimpson")
(p <- plot_richness(rarefyBC, "Burn", measures=alpha_graph))
#(p <- plot_richness(rarefyBC, "SamplingTime", "Random", measures=alpha_graph)) #Parse by a second variable
p + geom_boxplot(data=p$data, aes(x=Burn, y=value, color=NULL, group=Burn), alpha=0.1) +
  labs(x="Burn Treatment", y = "Number of OTUs")

alpha_graph = c("Observed")
(p <- plot_richness(rarefyBC.spr, "Burn", measures=alpha_graph))
p + geom_boxplot(data=p$data, aes(x=Burn, y=value, color=NULL, group=Burn), alpha=0.1)+
  labs(x="Burn Treatment", y = "Number of OTUs")

alpha_graph = c("Observed")
(p <- plot_richness(rarefyBC.sum, "Burn", measures=alpha_graph))
p + geom_boxplot(data=p$data, aes(x=Burn, y=value, color=NULL, group=Burn), alpha=0.1)+
  labs(x="Burn Treatment", y = "Number of OTUs")

alpha_graph = c("Observed")
(p <- plot_richness(rarefyBC.f, "Burn", measures=alpha_graph))
p + geom_boxplot(data=p$data, aes(x=Burn, y=value, color=NULL, group=Burn), alpha=0.1)+
  labs(x="Burn Treatment", y = "Number of OTUs")

#Alpha diversity for tree vs. random samples

alpha_chart = c("Observed", "Chao1", "Simpson")
estimate_richness(rarefyBC, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.rand, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.norand, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.spr, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.sum.r, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.sum.t, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.f.r, split = FALSE, measures = alpha_chart)
estimate_richness(rarefyBC.f.t, split = FALSE, measures = alpha_chart)

erich <- estimate_richness(rarefyBC, measures = c("Observed", "Chao1", "Simpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC)$Rand)[c("estimate","p.value","statistic","conf.int")])))
ttest
erich <- estimate_richness(rarefyBC.sum, measures = c("Observed", "Chao1", "Simpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.sum)$Rand)[c("estimate","p.value","statistic","conf.int")])))
ttest
erich <- estimate_richness(rarefyBC.f, measures = c("Observed", "Chao1", "Simpson"))
ttest <- t(sapply(erich, function(x) unlist(t.test(x~sample_data(rarefyBC.f)$Rand)[c("estimate","p.value","statistic","conf.int")])))
ttest


##### ADONIS ###########

# Function to run adonis test on a phyloseq object and a variable from metadata
# Make sure OTU data is standardized/normalized before 
phyloseq_to_adonis <- function(physeq, distmat = NULL, dist = "bray", formula) {
  
  if(!is.null(distmat)) {
    phydist <- distmat
  } else {
    phydist <- phyloseq::distance(physeq, dist)
  }
  
  metadata <- as(sample_data(physeq), "data.frame")
  
  # Adonis test
  f <- reformulate(formula, response = "phydist")
  adonis.test <- adonis(f, data = metadata)
  print(adonis.test)
  
  # Run homogeneity of dispersion test if there is only 1 variable
  if (grepl("\\+", formula)) {
    l <- list(
      dist = phydist, 
      formula = f, 
      adonis = adonis.test
    )
  } else {
    group <- metadata[ ,formula]
    beta <- betadisper(phydist, group)
    disper.test = permutest(beta)
    print(disper.test)
    
    l <- list(
      dist = phydist, 
      formula = f, 
      adonis = adonis.test, 
      disper = disper.test
    )
  }
  return (l)
}



### Beta Diversity Measures ###

bray.nmds1 <- ordinate(physeq = rarefyBC, method = "NMDS", distance = "bray")

plot_ordination(rarefyBC, 
                bray.nmds1, 
                color="Burn")+
  stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
  geom_point(size=3)+
  labs(main="stress .10")+
  ggtitle("Total Fungal Community - Burn vs. Control")+
  theme(panel.grid=element_blank())

adonis.site <- phyloseq_to_adonis(
  physeq = rarefyBC, 
  dist = "bray", 
  formula = "Burn")

bray.nmds2 <- ordinate(physeq = rarefyBC.noburn, method = "NMDS", distance = "bray")

plot_ordination(rarefyBC.noburn, 
                bray.nmds2, 
                color="Plot")+
  stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
  geom_point(size=3)+
  labs(main="stress .10")+
  ggtitle("Control Areas")+
  theme(panel.grid=element_blank())

adonis.site <- phyloseq_to_adonis(
  physeq = rarefyBC.noburn, 
  dist = "bray", 
  formula = "Plot")

bray.nmds3 <- ordinate(physeq = rarefyBC, method = "NMDS", distance = "bray")

plot_ordination(rarefyBC, 
                bray.nmds3, 
                color="Burn",
                shape="Plot")+
  stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
  geom_point(size=3)+
  labs(main="stress .10")+
  ggtitle("Burn vs. Control")+
  theme(panel.grid=element_blank())

adonis.site <- phyloseq_to_adonis(
  physeq = rarefyBC, 
  dist = "bray", 
  formula = "Burn")



##### RANDOM FOREST MODELS #####

library(metagMisc)
burnControlPA = phyloseq_standardize_otu_abundance(burnControl, method = "pa")
print(burnControlPA)

# Set prunescale = Prune taxa with a relative abundance below:
# 0.0001 Reduces taxa for the final forest by about 90%. .000001 to around 1/3. 
prunescale = 0.00001
minlib = 24000

# Prune out rare OTUs by mean relative abundance set by prunescale
tax.mean1 <- taxa_sums(rarefyBC)/nsamples(rarefyBC)
treat.prune1 <- prune_taxa(tax.mean1 > prunescale*minlib, rarefyBC)
treat.prune1

treat.prune1PA = phyloseq_standardize_otu_abundance(treat.prune1, method = "pa") #Convert to a presence/absense table

tax.mean2 <- taxa_sums(rarefyBC.norand)/nsamples(rarefyBC.norand)
treat.prune2 <- prune_taxa(tax.mean2 > prunescale*minlib, rarefyBC.norand)
treat.prune2

treat.prune2PA = phyloseq_standardize_otu_abundance(treat.prune2, method = "pa") #Convert to a presence/absence table

tax.mean3 <- taxa_sums(rarefyBC.rand)/nsamples(rarefyBC.rand)
treat.prune3 <- prune_taxa(tax.mean3 > prunescale*minlib, rarefyBC.rand)
treat.prune3

treat.prune3PA = phyloseq_standardize_otu_abundance(treat.prune3, method = "pa")

# Make a dataframe of training data with OTUs as column and samples as rows
predictors1 <- t(otu_table(treat.prune1PA))
dim(predictors1)

predictors2 <- t(otu_table(treat.prune2PA))
dim(predictors2)

predictors3 <- t(otu_table(treat.prune3PA))
dim(predictors3)


# Make one column for our outcome/response variable 
response1 <- as.factor(sample_data(treat.prune1PA)$Burn) #Classification
print(response1)
response2 <- as.factor(sample_data(treat.prune2PA)$Burn) #Classification
print(response2)
response3 <- as.factor(sample_data(treat.prune3PA)$Burn) #Classification
print(response3)

response7 <- as.numeric(sample_data(treat.prune1PA)$Burn) #Regression
print(response7)
response8 <- as.numeric(sample_data(treat.prune2PA)$Burn) #Regression
print(response8)
response9 <- as.numeric(sample_data(treat.prune3PA)$Burn) #Regression
print(response9)



# Combine them into 1 data frame
rf.data1 <- data.frame(response1, predictors1)
rf.data2 <- data.frame(response2, predictors2)
rf.data3 <- data.frame(response3, predictors3)

rf.data7 <- data.frame(response7, predictors1)
rf.data8 <- data.frame(response8, predictors2)
rf.data9 <- data.frame(response9, predictors3)


#Run the model
set.seed(2)
burn1.classify <- randomForest(response1~., data = rf.data1, proximity = TRUE, importance = TRUE, ntree = 5001)
print(burn1.classify) #Burn vs. Control - All Samples
set.seed(2)
burn2.classify <- randomForest(response2~., data = rf.data2, proximity = TRUE, importance = TRUE, ntree = 5001)
print(burn2.classify) #Burn vs. Control - Tree Samples
set.seed(2)
burn3.classify <- randomForest(response3~., data = rf.data3, proximity = TRUE, importance = TRUE, ntree = 5001)
print(burn3.classify) #Burn vs. Control - Random Samples
vs. Control - Presence Absence - Tree Samples


#Random Forest Regression
set.seed(2)
burn7.regress <- randomForest(response7~., data = rf.data7, proximity = TRUE, ntree = 5001)
print(burn7.regress) #Burns - All Samples
set.seed(2)
burn8.regress <- randomForest(response8~., data = rf.data8, proximity = TRUE, ntree = 5001)
print(burn8.regress) #Burns - Tree Samples
set.seed(2)
burn9.regress <- randomForest(response9~., data = rf.data9, proximity = TRUE, ntree = 5001)
print(burn9.regress) #Burns - Random Samples


#Variables that are stored in the outputs
names(burn1.classify)
print(burn1.classify$proximity)

# Make a data frame with predictor names and their importance
imp1 <- importance(burn1.classify)
imp1 <- data.frame(predictors1 = rownames(imp1), imp1)


#RF Plots

proximityPlot(burn5.classify, dim.x = 1, dim.y = 2, legend.loc = c("right"), 
              point.size = 2, circle.size = 0,
              circle.border = 1, hull.alpha = 0.3, plot = TRUE)

confusionMatrix(burn1.classify, conf.level = 0.95, threshold = 0.8)
confusionMatrix(burn2.classify, conf.level = 0.95, threshold = 0.8)
confusionMatrix(burn3.classify, conf.level = 0.95, threshold = 0.8)

plotInbag(burn1.classify, sampsize = NULL, bins = 20, plot = TRUE)


### DESeq2 - Differential Abundance  ###

#BiocManager::install("DESeq2")
library("DESeq2")
sample_data(rarefyBC)$Burn <- as.factor(sample_data(rarefyBC)$Burn) #Converts variable to a factor for analysis

#add 1 to each sample count
otu_table(rarefyBC) <- otu_table(rarefyBC) + 1

BC = phyloseq_to_deseq2(rarefyBC, ~ Burn)
BC = DESeq(BC)

alpha = 0.01
res = results(BC, contrast=c("Burn", "0", "1"), alpha=alpha)
res = res[order(res$padj, na.last=NA), ]
res_sig = res[(res$padj < alpha), ]
res_sig

res_sig = cbind(as(res_sig, "data.frame"), as(tax_table(rarefyBC)[rownames(res_sig), ], "matrix"))

res_sig$predictors1<-rownames(res_sig)

allDESeq <- merge(res_sig, imp1, all.y = predictors1)

#Drop records that do not have an output from DESeq2
library(tidyverse)
allDESeq <- allDESeq %>% drop_na("log2FoldChange")

#Order by importance
imp2 <- arrange(allDESeq, desc(MeanDecreaseAccuracy))
imp2$predictors1 <- factor(imp2$predictors1, levels = imp2$predictors1)
# Select the top 20 predictors
diff <- imp2[1:20, ]

#Export DESeq-RF Indicator Species Analysis as Excel Document
write.csv(imp2, file = "indicator-species.csv")


#library("tidyverse")
diff$Genus <- str_remove(diff$Genus, "g__")
#diff$Family <- str_remove(diff$Family, "f__")

#Cleanup diff dataframe for charts - Updated genus for several OTUs manually using BLAST
levels(diff$Genus) <- c(levels(diff$Genus), "Neurospora")
diff["1", "Genus"] <- "Sordariaceae" #OTU902
diff["3", "Genus"] <- "Fungi"      #OTU1194
diff["4", "Genus"] <- "Ascomycota" #OTU822
diff["5", "Genus"] <- "Umbelopsidomycetes"   #OTU1131
diff["6", "Genus"] <- "Ascomycota"     #OTU579
diff["10", "Genus"] <- "Pezizomycotina"   #OTU18718
diff["11", "Genus"] <- "Basidiomycota"   #OTU1900
diff["13", "Genus"] <- "Cladophialophora"   #OTU18813
diff["14", "Genus"] <- "Cucurbitariaceae"   #OTU1350
diff["15", "Genus"] <- "Mortierella"   #OTU1092
diff["16", "Genus"] <- "Helotiales"   #OTU2453
diff["19", "Genus"] <- "Glomeraceae"   #OTU2361

#Variable for Positive vs. Negative
diff[["sign"]] = ifelse(diff[["log2FoldChange"]] >= 0, "positive", "negative")

diff<-diff[order(diff$MeanDecreaseAccuracy),]

#Plot 1 - Classification Accuracy

plot1 <- ggplot(diff, aes(x = MeanDecreaseAccuracy, y = reorder(predictors1, MeanDecreaseAccuracy))) +
  geom_point(size=3, width = 0.2) +
  scale_y_discrete(label=function(y){return(paste(diff$Genus, diff$predictors1))})+
  xlab("OTU Importance")

plot1.1 <- plot1 + theme(axis.title.x = element_text(color="black", size=8),
                         axis.title.y = element_blank())

#Plot 2 - Differential Abundance
plot2 <- ggplot(diff, aes(x=log2FoldChange, y=reorder(predictors1, MeanDecreaseAccuracy), color = sign)) +
  geom_point(size=3, width = 0.2) +
  xlab("  Burned                           Unburned")+
  geom_vline(xintercept=0, linetype="dotted")

plot2.1 <- plot2 + theme(
                         axis.title.x = element_text(color="black", size=8),
                         axis.title.y = element_blank(),
                         legend.position="none",
                         axis.text.y = element_blank())

library(cowplot)
plot_grid(plot1.1, plot2.1)



##Boxplot of Alpha diversity for months since burn

devtools::source_gist("8d0ca4206a66be7ff6d76fc4ab8e66c6")

BC.burn.t = subset_samples(BC.burn, Random == "0")
print(BC.burn.t)
rarefyBC.burn.t = rarefy_even_depth(BC.burn.t)
print(rarefyBC.burn.t)

alpha_graph = c("Observed", "InvSimpson")

#install.packages("ggpmisc")
library(ggpmisc)

formula <- y ~ x

p <- plot_richness(rarefyBC.burn.t, "MSB", measures = alpha_graph)
p <- p + geom_boxplot(data=p$data, 
                      aes(x = MSB, y=value, color=NULL, group=MSB), 
                      alpha=0.1)+
  labs(x="Months Since Prescribed Burn", y="Alpha Diversity Measure") + 
  geom_smooth(method = "lm", se=TRUE, color="red", aes(group=1))+
  stat_poly_eq(aes(label = paste(..rr.label..)),
               label.x = 0.8, 
               label.y = 0.78, 
               formula = formula, 
               parse = TRUE, size = 3) + 
  stat_fit_glance(method = 'lm', method.args = list(formula = formula),
                  label.x = 0.82,
                  label.y = 0.82,
                  aes(label = paste("P = ", 
                                    signif(..p.value.., digits = 2), sep = "")),
                  size = 3)
p


### Heatmap & Alpha Diversity Graph

#Heatmap showing phyla for months since burn
library(ampvis2)
reorder <- c("1", "0")
sample_data(treat.prune1PA)$Burn <- factor(sample_data(rarefyBC)$Burn, levels = reorder)
sample_data(treat.prune1PA)$MSB <- as.character(sample_data(rarefyBC)$MSB)

ampBC <- phyloseq_to_ampvis2(treat.prune1PA)
amp_heatmap(ampBC, group_by = "MSB", facet_by = "Burn", measure = "mean", normalise = TRUE, tax_aggregate = "Phylum", tax_show = 5)


#Boxplots showing phyla for months since burn
alpha_graph2 = c("Observed")

ampBC2 <- phyloseq_to_ampvis2(treat.prune1PA)
ampBC2 <- amp_subset_samples(ampBC2, Burn %in% c("1"))

p <- amp_boxplot(ampBC2,
            group_by = "MSB",
            tax_show = 5,
            tax_aggregate = "Phylum",
            plot_flip = TRUE
      )


p$data <- unique(p$data)
p$data$Abundance <- as.numeric(p$data$Abundance)
p$data$Group <- as.numeric(as.character(p$data$Group))

p$data$Display <- factor(p$data$Display,levels=c("Ascomycota","Basidiomycota","Mucoromycota","Mortierellomycota","Rozellomycota"))

p1 <- ggplot(data=p$data, aes(Group, Abundance)) + 
      geom_point()
p1 <- p1 + geom_boxplot(data=p$data, 
                      aes(x = Group, y=Abundance, color=NULL, group=Group), 
                      alpha=0.1)+
                      facet_grid(cols = vars(p$data$Display))+
           labs(x="Months Since Prescribed Burn", y="Abundance %") +
           ylim(0,70) + 
           geom_smooth(method = "lm", se=TRUE, color="red", aes(group=Display))+
           stat_poly_eq(aes(label = paste(..rr.label..)),
                        label.x = 0.05, 
                        label.y = 0.92, 
                        formula = y ~ x, 
                        parse = TRUE, size = 3)+
           stat_fit_glance(method = 'lm', 
                           method.args = list(formula = y ~ x),
                           aes(label = paste("P = ", 
                               signif(..p.value.., digits = 2),
                               sep = ""), npcx=0.4, npcy=1),
                               size = 3, parse = FALSE)

p1




#amp_venn(ampBC, group_by = "Random", cut_f = 1, cut_a = .001) #Not used in this paper