ASLeafOverview.rmd

---
title: "ASLeaf example"
author: "JReceveur"
output: 
  html_document:
    code_folding: hide
    toc: true
    toc_depth: 3
    toc_float: true
editor_options: 
  chunk_output_type: console
---
  
  ```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width=14, fig.height=10)
knitr::opts_chunk$set(fig.align="center")
```




# Import{.tabset}

## Sample Overview

```{r,warning=F,include=F}
#From https://gist.github.com/grabear/018e86413b19b62a6bb8e72a9adba349
#Parse Silva function

parse_taxonomy_silva_138 <- function(char.vec){
  # Use default to assign names to elements in case problem with greengenes prefix
  char.vec = parse_taxonomy_default(char.vec)
  # Check for unassigned taxa
  if (char.vec["Rank1"] == "Unassigned") {
    char.vec <- c(Rank1="d__Unassigned", Rank2="d__Unassigned", Rank3="d__Unassigned", Rank4="d__Unassigned",
                  Rank5="d__Unassigned", Rank6="d__Unassigned", Rank7="d__Unassigned")
  }
  # Define the meaning of each prefix according to SILVA taxonomy
  Tranks = c(Rank1="Kingdom", Rank2="Phylum", Rank3="Class", Rank4="Order", Rank5="Family", Rank6="Genus", Rank7="Species")
  # Check for prefix using regexp, warn if there were none. trim indices, ti
  ti = grep("[[:alpha:]]\\_\\_", char.vec)
  if( length(ti) == 0L ){
    warning(
      "No silva prefixes were found. \n",
      "Consider using parse_taxonomy_delfault() instead if true for all OTUs. \n",
      "Dummy ranks may be included among taxonomic ranks now."
    )
    # Will want to return without further modifying char.vec
    taxvec = char.vec
    # Replace names of taxvec according to prefix, if any present...
  } else {
    # Format character vectors for Ambiguous taxa
    if( length(ti) < 7 ){
      for (key in names(char.vec)) {
        if ( char.vec[key] == "Ambiguous_taxa" ) {
          tax_no <- (as.numeric(substr(key, 5, 5)) - 1)
          char.vec[key] = sprintf("d__Ambiguous_taxa", tax_no)
        }
      }
      # Reset the trimmed indicies if Ambiguous taxa
      ti = grep("[[:alpha:]]\\_\\_", char.vec)
    }
    # Remove prefix using sub-"" regexp, call result taxvec
    taxvec = gsub("[[:alpha:]]\\_\\_", "", char.vec)
    # Define the ranks that will be replaced
    repranks = Tranks[substr(char.vec[ti], 1, 3)]
    # Replace, being sure to avoid prefixes notK present in Tranks
    names(taxvec)[ti[!is.na(repranks)]] = repranks[!is.na(repranks)]
  }
  return(taxvec)
}




library(vegan)
library(ggplot2)
library(plyr)
library(dplyr)
library(phyloseq)
library(randomForest)
library(knitr)




set.seed(10)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#000000","#CC79A7") #Build a user defined colorblind palette



theme_set(theme_bw(base_size = 18)+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())) #Set font size and change default plot theme

biom=import_biom("ASLeafDec2020_7k.biom",parseFunction=parse_taxonomy_silva_138) #Import biom file which contains taxa names and abundances 
colnames(tax_table(biom))<-c(Rank1="Kingdom", Rank2="Phylum", Rank3="Class", Rank4="Order", Rank5="Family", Rank6="Genus", Rank7="Species")




metadata=read.table("LauraMetadata9Dec2020.txt",header = TRUE)

metadata$Date = factor(metadata$Date, levels = c("3_1","24_1","14_2","7_3","29_3","19_4")) #fixes x-axis labels
metadata$Paired<-paste0(metadata$Site,metadata$Date)

tree=read_tree("ASLeafTreeDec2020.nwk")

sampdat=sample_data(metadata)
sample_names(sampdat)=metadata$ID # If using fungal data change ID to FungalID
physeq=merge_phyloseq(biom,sampdat,tree)
physeq<-filter_taxa(physeq, function (x) {sum(x > 0) > 1}, prune=TRUE)#Remove singletons and zeros
physeq<-subset_samples(physeq, Site!="<NA>") # Remove Control
physeq
#physeq 



GPr  = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
PhylumAll=tax_glom(GPr, "Phylum")

PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 3e-2, TRUE) #filter out any taxa lower than 3%
FamilyAll=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(FamilyAll, function(x) mean(x) > 3e-2, TRUE) #filter out any taxa lower than 3%
GenusAll=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GenusAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower than 2%

```

```{r}
physeq
```

## Metadata
```{r}
metadata
```



```{r}
#Converting and summarizing c + Q mass for metadata
# Chestnut<-read.csv2("Chestnut_Decomposition.csv",header=T)
# Oak<-read.csv2("Oak_Decomposition.csv")
# head(Chestnut)
# head(Oak)
# #write.csv(Chestnut, "ChestnutDecompConverted.csv")
# #write.csv(Oak,"OakDecompConverted.csv")
# Oak<-read.csv("OakDecompConverted.csv",header=T)
# 
# OakSummary <- ddply(Oak, c("Site","Date"), summarise,
#                  N    = length(MassLoss.),
#                  meanMassRemaining = mean(MassLoss.),
#                  sd   = sd(MassLoss.),
#                  se   = sd / sqrt(N)
# )
# OakSummary
# 
# Chestnut<-read.csv("ChestnutDecompConverted.csv",header=T)
# 
# ChestnutSummary <- ddply(Chestnut, c("Site","Date"), summarise,
#                  N    = length(MassLoss.),
#                  meanMassRemaining = mean(MassLoss.),
#                  sd   = sd(MassLoss.),
#                  se   = sd / sqrt(N)
# )
# ChestnutSummary
# 
# 
# ggplot(Oak,aes(x=Date,y=MassLoss.,color=Date))+geom_point()+facet_wrap(~Site)
# 

```


# Alpha Diversity

## Shannon Diversity

```{r,warning=F}
Shannon<-plot_richness(physeq, x="Date",color="RiverType1",shape="Treatment", measures=c("Shannon"))+ylab("Shannon Diversity")+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)+facet_wrap(~Site)
ShannonData<-Shannon$data
shapiro.test(ShannonData$value)


ShannonAOV<-aov(value~RiverType1*Treatment*Date,data=ShannonData)
summary(ShannonAOV)
# plot(ShannonAOV) #plots ok at best, no major bias by fitted values
hist(resid(ShannonAOV))
shapiro.test(resid(ShannonAOV))
library(ggpubr)
ggpaired(ShannonData,x="Treatment",y="value",id="Paired",color="Treatment",line.color = "gray", line.size = 0.4,)+ylab("Bacterial Shannon Diversity")+xlab("Leaf Type")+stat_compare_means(paired=T)


CData<-subset(ShannonData,Treatment=="C")
QData<-subset(ShannonData,Treatment=="Q")


CShannonAOV<-aov(value~RiverType1*Date,data=CData)
summary(CShannonAOV)


QShannonAOV<-aov(value~RiverType1*Date,data=QData)
summary(QShannonAOV)



kruskal.test(CData$value, g=CData$RiverType1)
kruskal.test(QData$value, g=QData$RiverType1)



compare_means(value ~ RiverType1, data = CData, group.by = "Date", p.adjust.method = "fdr",method="kruskal.test") # ns
compare_means(value ~ RiverType1, data = QData, group.by = "Date", p.adjust.method = "fdr",method="kruskal.test") # ns



kruskal.test(CData$value, g= CData$Date)
kruskal.test(QData$value, g= QData$Date)


CDataPer<-subset(CData, RiverType1=="Perennial")
kruskal.test(CDataPer$value, g=CDataPer$Date)
QDataPer<-subset(QData, RiverType1=="Perennial")
kruskal.test(QDataPer$value, g=QDataPer$Date)



# compare_means(
#   value~Treatment,
#   ShannonData,
#   method = "wilcox.test",
#   paired = T,
#   id = "Paired",
#   ref.group = "C",
#   p.adjust.method = "fdr"
# )
# 

CombinedDifference<-data.frame()
for (i in unique(ShannonData$Paired)){
  Subset<-subset(ShannonData, Paired==i)
  LeafQ<-subset(Subset,Treatment=="Q")
  LeafC<-subset(Subset,Treatment=="C")
  Difference<-LeafQ$value-LeafC$value
  CombinedDifference<-rbind(CombinedDifference,Difference)
}
colnames(CombinedDifference)<-"Difference"
hist(CombinedDifference$Difference)
mean(CombinedDifference$Difference)
sd(CombinedDifference$Difference)/sqrt(length(CombinedDifference$Difference))



ShannonSummary <- ddply(ShannonData, c("Treatment","Days","RiverType1"), summarise,
                 N    = length(value),
                 mean = mean(value),
                 sd   = sd(value),
                 se   = sd / sqrt(N)
)

ggplot(ShannonSummary,aes(x=Days,y=mean,group=RiverType1,color=RiverType1,shape=RiverType1))+geom_line(size=2)+facet_wrap(~Treatment)+ylab("Bacterial Shannon Diversity (+/- SEM)")+
  geom_point(data=ShannonData,aes(x=Days,y=value))+ geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_jitter()
```


```{r}

ggplot(ShannonData,aes(x=value,y=LeafMassRemaining,color=Treatment))+geom_point()+facet_wrap(~RiverType1)

ggplot(ShannonData,aes(x=value,y=LeafMassRemaining,color=Treatment))+geom_point()+facet_wrap(~RiverType1)


summary(aov(LeafMassRemaining~value+RiverType1,data=ShannonData))
library(nlme)

Diversity = lme(value~Days+Treatment+RiverType1,random=~1|Site,data=ShannonData)


summary(Diversity)
hist(resid(Diversity))
```





Shannon Summary (Crissolo vs Revello)

```{r}
CrisVRevello<-subset(ShannonData,Site=="CRI"|Site=="RE")
CrisVRevello


ggplot(CrisVRevello,aes(x=LeafMassRemaining,y=value,color=Treatment))+geom_point()+facet_wrap(~Site)


ggplot(CrisVRevello,aes(x=Site,y=value))+geom_boxplot()+facet_wrap(~Treatment)


summary(aov(value~Site,data=CrisVRevello))



ggplot(ShannonData,aes(x=RiverType1,y=value,color=RiverType1))+geom_boxplot()+facet_wrap(~Treatment)


ShannonOak<-subset(ShannonData,Treatment=="Q")
ShannonC<-subset(ShannonData,Treatment=="C")
summary(aov(value~RiverType1+Date+Site,data=ShannonOak))
summary(aov(value~RiverType1+Date+Site,data=ShannonC))


ggplot(ShannonC,aes(x=Date,y=value,color=RiverType1))+geom_boxplot()


ggplot(ShannonOak,aes(x=Date,y=value,color=RiverType1))+geom_boxplot()


```


## Observed richness

```{r}
plot_richness(physeq, x="Date",color="RiverType2",shape="Treatment", measures=c("Observed"))+ylab("Observed Species")+
  facet_wrap(~SampleLoc)+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)+facet_wrap(~Site)


Observed<-plot_richness(physeq, x="Date",color="RiverType1",shape="Treatment", measures=c("Observed"))+ylab("Observed Species")+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)+facet_wrap(~Site)
ObservedData<-Observed$data

ObservedAOV<-aov(value~RiverType1*Treatment*Date,data=ObservedData)
summary(ObservedAOV)


ggpaired(ObservedData,x="Treatment",y="value",id="Paired",color="Treatment",line.color = "gray", line.size = 0.4,)+ylab("Bacterial (ASV) Richness")+xlab("Leaf Type")+stat_compare_means(paired=T)


compare_means(
  value~Treatment,
  ObservedData,
  method = "wilcox.test",
  paired = T,
  id = "Paired",
  ref.group = "C",
  p.adjust.method = "fdr"
)


CombinedDifference<-data.frame()
for (i in unique(ObservedData$Paired)){
  Subset<-subset(ObservedData, Paired==i)
  LeafQ<-subset(Subset,Treatment=="Q")
  LeafC<-subset(Subset,Treatment=="C")
  Difference<-LeafQ$value-LeafC$value
  CombinedDifference<-rbind(CombinedDifference,Difference)
}
hist(CombinedDifference$X104)
mean(CombinedDifference$X104)
sd(CombinedDifference$X104)/sqrt(length(CombinedDifference$X104))






ObservedSummary <- ddply(ObservedData, c("Treatment","Date","RiverType1"), summarise,
                 N    = length(value),
                 mean = mean(value),
                 sd   = sd(value),
                 se   = sd / sqrt(N)
)

ggplot(ObservedSummary,aes(x=Date,y=mean,group=RiverType1,color=RiverType1,shape=RiverType1))+geom_line(size=2)+facet_wrap(~Treatment)+ylab(" Observed Bacterial Richness (+/- SEM)")+theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  geom_point(data=ObservedData,aes(x=Date,y=value))+ geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_jitter()


CData<-subset(ObservedData,Treatment=="C")
QData<-subset(ObservedData,Treatment=="Q")


CRichnessAOV<-aov(value~RiverType1*Date,data=CData)
summary(CRichnessAOV)


QRichnessAOV<-aov(value~RiverType1*Date,data=QData)
summary(QRichnessAOV)



CDataPer<-subset(CData, RiverType1=="Perennial")
kruskal.test(CDataPer$value, g=CDataPer$Date)
QDataPer<-subset(QData, RiverType1=="Perennial")
kruskal.test(QDataPer$value, g=QDataPer$Date)







```

# Taxa plots

```{r,warning=F}
###########
#Phylum level plot
###########
df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum","Treatment","Date","RiverType1"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
PhylumPlotTreatment=ggplot(Trtdata, aes(x=Date,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+
  facet_grid(RiverType1~Treatment)+xlab("Date")+ylab("Relative Abundance (> 3%, SEM)")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
PhylumPlotTreatment

dfCP<-subset(df,Treatment=="C"&RiverType1=="Perennial")
compare_means(Abundance ~ Date, data = dfCP, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # ns

dfCI<-subset(df,Treatment=="C"&RiverType1=="Intermittent")
compare_means(Abundance ~ Date, data = dfCI, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # ns



dfQP<-subset(df,Treatment=="Q"&RiverType1=="Perennial")
compare_means(Abundance ~ Date, data = dfQP, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # Actinobacteriota

dfQPActinobacteria<-subset(dfQP,Phylum=="Actinobacteriota")
compare_means(Abundance ~ Date, data = dfQPActinobacteria, group.by = "Phylum", p.adjust.method = "fdr",method="wilcox.test") # Actinobacteriota

#Actinobacteriota plot
TrtdataActino <- ddply(dfQPActinobacteria, c("Phylum","Date"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
ggplot(TrtdataActino, aes(x=Date,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("Date")+ylab("Relative Abundance Actinobacteriota (%, SEM)")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+theme(legend.position="none")












dfQI<-subset(df,Treatment=="Q"&RiverType1=="Intermittent")
compare_means(Abundance ~ Date, data = dfQI, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # ns


dfC<-subset(df,Treatment=="C")
compare_means(Abundance ~ Date, data = dfC, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # ns # Actinobacteria close


dfQ<-subset(df,Treatment=="Q")
compare_means(Abundance ~ Date, data = dfQ, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # Actinobacteriota



dfC<-subset(df,Treatment=="C")
compare_means(Abundance ~ RiverType1, data = dfC, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")


dfQ<-subset(df,Treatment=="Q")
compare_means(Abundance ~ RiverType1, data = dfQ, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")


dfPer<-subset(df,RiverType1=="Perennial")
compare_means(Abundance ~ Treatment, data = dfPer, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # Actinobacteriota significant, Bacteroidota nearly sig

ddply(dfPer, c("Phylum","Treatment","RiverType1"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N))



dfInt<-subset(df,RiverType1=="Intermittent")
compare_means(Abundance ~ Treatment, data = dfInt, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test") # Actinobacteriota significant, Bacteroidota sig

ddply(dfInt, c("Phylum","Treatment","RiverType1"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N))



```


```{r}
# Note: Family level code not updated recently!!

df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family","Treatment","Date","Site"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
FamilyPlotTreatment=ggplot(Trtdata, aes(x=Date,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+
  facet_grid(Site~Treatment)+xlab("Date")+ylab("Relative Abundance (> 3%, SEM)")
FamilyPlotTreatment

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "Treatment","Date"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
#Trtdata
cdataplot=ggplot(Trtdata, aes(x=Date,y=mean))+geom_bar(aes(fill = Treatment),colour="black", stat="identity")+ facet_grid(Phylum~Treatment)+xlab("Date")+ylab("Relative Abundance (%, SEM)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
```


# Beta diversity

```{r,warning=F}
#Multivariate dispersions
set.seed(12345)
GPdist=phyloseq::distance(physeq, "bray")
beta=betadisper(GPdist, sample_data(physeq)$Site)
permutest(beta)
boxplot(beta)

#Multivariate dispersions
set.seed(12345)

GPdist=phyloseq::distance(physeq, "bray")
beta=betadisper(GPdist, sample_data(physeq)$Treatment)
permutest(beta)
boxplot(beta)

set.seed(12345)

GPdist=phyloseq::distance(physeq, "bray")
beta=betadisper(GPdist, sample_data(physeq)$Date)
permutest(beta)
boxplot(beta)

set.seed(12345)
GPdist=phyloseq::distance(physeq, "bray")
beta=betadisper(GPdist, sample_data(physeq)$RiverType1)
permutest(beta)
boxplot(beta)
DistanceFromCenter<-data.frame("Distance"=beta$distances,"Group"=beta$group)
ddply(DistanceFromCenter, c("Group"), summarise,
                 N    = length(Distance),
                 mean = mean(Distance),
                 sd   = sd(Distance),
                 se   = sd / sqrt(N)
)

adonis(GPdist ~ RiverType1*Treatment*Site+Date, as(sample_data(physeq), "data.frame"))

set.seed(12345)
ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="Site",shape="RiverType1")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+facet_wrap(~Treatment)+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))


ordplot2=plot_ordination(physeq, ord,"samples", color="Treatment",shape="Treatment")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot2<-ordplot2+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Treatment))+facet_wrap(~Date)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ordplot2


ordplot3=plot_ordination(physeq, ord,"samples", color="Treatment",shape="RiverType1")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot3<-ordplot3+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Treatment))+facet_wrap(~Date)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
ordplot3

set.seed(12345)
ord=ordinate(physeq,"PCoA", "bray")
ordplot4=plot_ordination(physeq, ord,"samples", color="RiverType1",shape="RiverType1")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot4+facet_wrap(~Treatment)+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))


```



# Random Forest{.tabset}

```{r,warning=F}
print("Treatment")
set.seed(155258)
GenusLevelRF = filter_taxa(GenusAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower than 1%

ForestData=GenusLevelRF#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Treatment)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000,importance=T)
print(MozzieForest)#returns overall Random Forest results

imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseAccuracy))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:10, ]#Eventually you'll probably need to do this step manually 

ggplot(imp.20, aes(x = predictors, y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  samples\n by Treatment")#\n in a string tells it to start a new line



destroyX = function(es) {
  f = es
  for (row in c(1:nrow(f))){ #for each column in dataframe
    if (startsWith(row.names(f)[row], "X") == TRUE)  { #if starts with 'X' ..
      row.names(f)[row] <- substr(row.names(f)[row], 2, 100) #get rid of it
    }
    assign(deparse(substitute(es)), f, inherits = TRUE)
  }
}
destroyX(imp.20)


otunames <- row.names(imp.20)
r <- rownames(tax_table(ForestData)) %in% otunames
otunames
PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification

GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)

df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100

dfSorted<- df[with(df, order(Genus)), ]
dfSorted$FamilyGenus<-paste0(dfSorted$Family,": \n",dfSorted$Genus)
Trtdata <- ddply(dfSorted, c("FamilyGenus", "Treatment"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
TrtdataSorted<- Trtdata[with(Trtdata, order(FamilyGenus)), ]
TrtdataSorted
cdataplot=ggplot(Trtdata, aes(x=Treatment,y=mean))+geom_bar(aes(fill = Treatment),colour="black", stat="identity")+xlab("Leaf Type ")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~FamilyGenus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

Comparisons<-compare_means(Abundance ~ Treatment, data = dfSorted, group.by = "FamilyGenus", p.adjust.method = "fdr",method="kruskal.test")
Comparisons
#write.csv(Comparisons, "GenusLevelRFIndicatorsLeafType.csv")



```

Indicators for River type within chestnut leaves

```{r}
print("River")
set.seed(155258)
GenusLevelRF = filter_taxa(GenusAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower than 1%
GenusLevelRFChestnut<-subset_samples(GenusLevelRF,Treatment=="C")
ForestData=GenusLevelRFChestnut#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$RiverType1)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000,importance=T)
print(MozzieForest)#returns overall Random Forest results

imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseAccuracy))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:10, ]#Eventually you'll probably need to do this step manually 

ggplot(imp.sort, aes(x = predictors, y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  chestnut samples\n by River")#\n in a string tells it to start a new line



destroyX = function(es) {
  f = es
  for (row in c(1:nrow(f))){ #for each column in dataframe
    if (startsWith(row.names(f)[row], "X") == TRUE)  { #if starts with 'X' ..
      row.names(f)[row] <- substr(row.names(f)[row], 2, 100) #get rid of it
    }
    assign(deparse(substitute(es)), f, inherits = TRUE)
  }
}
destroyX(imp.20)


otunames <- row.names(imp.20)
r <- rownames(tax_table(ForestData)) %in% otunames


GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)


df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100

dfSorted<- df[with(df, order(Genus)), ]

dfSorted$FamilyGenus<-paste0(dfSorted$Family,": \n",dfSorted$Genus)
ggplot(dfSorted, aes(x=RiverType1,y=Abundance,color=RiverType1))+geom_boxplot()+geom_point()+facet_wrap(~FamilyGenus,scales = "free_y")


Trtdata <- ddply(dfSorted, c("FamilyGenus", "RiverType1"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)

TrtdataSorted<- Trtdata[with(Trtdata, order(FamilyGenus)), ]
TrtdataSorted
cdataplot=ggplot(Trtdata, aes(x=RiverType1,y=mean))+geom_bar(aes(fill = RiverType1),colour="black", stat="identity")+xlab("Leaf Type ")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~FamilyGenus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

Comparisons<-compare_means(Abundance ~ Treatment, data = dfSorted, group.by = "FamilyGenus", p.adjust.method = "fdr",method="kruskal.test")
Comparisons
#write.csv(Comparisons, "GenusLevelRFIndicatorsRiverChestnut.csv")

```

Indicators for River type within Oak leaves

```{r}
print("River")
set.seed(155258)
GenusLevelRF = filter_taxa(GenusAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower than 1%
GenusLevelRFOak<-subset_samples(GenusLevelRF,Treatment=="Q")
GenusLevelRFOak
ForestData=GenusLevelRFOak#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$RiverType1)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000,importance=T)
print(MozzieForest)#returns overall Random Forest results

imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseAccuracy))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:10, ]#Eventually you'll probably need to do this step manually 

ggplot(imp.sort, aes(x = predictors, y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  oak samples\n by River")#\n in a string tells it to start a new line



destroyX = function(es) {
  f = es
  for (row in c(1:nrow(f))){ #for each column in dataframe
    if (startsWith(row.names(f)[row], "X") == TRUE)  { #if starts with 'X' ..
      row.names(f)[row] <- substr(row.names(f)[row], 2, 100) #get rid of it
    }
    assign(deparse(substitute(es)), f, inherits = TRUE)
  }
}
destroyX(imp.20)


otunames <- row.names(imp.20)
r <- rownames(tax_table(ForestData)) %in% otunames


GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)


df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100

dfSorted<- df[with(df, order(Genus)), ]

dfSorted$FamilyGenus<-paste0(dfSorted$Family,": \n",dfSorted$Genus)
Trtdata <- ddply(dfSorted, c("FamilyGenus", "RiverType1"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)

TrtdataSorted<- Trtdata[with(Trtdata, order(FamilyGenus)), ]
TrtdataSorted
cdataplot=ggplot(Trtdata, aes(x=RiverType1,y=mean))+geom_bar(aes(fill = RiverType1),colour="black", stat="identity")+xlab("River Type")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~FamilyGenus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot
ggplot(dfSorted, aes(x=RiverType1,y=Abundance,color=RiverType1))+geom_boxplot()+geom_point()+facet_wrap(~FamilyGenus,scales = "free_y")


Comparisons<-compare_means(Abundance ~ Treatment, data = dfSorted, group.by = "FamilyGenus", p.adjust.method = "fdr",method="kruskal.test")
Comparisons
#write.csv(Comparisons, "GenusLevelRFIndicatorsRiverOak.csv")

```

Comparing OST to other  sites

```{r}
OSTComparison=merge_phyloseq(biom,sampdat,tree)
OSTComparison<-subset_samples(OSTComparison, RiverType1=="Perennial"|RiverType1=="Intermittent")

ShannonOST<-plot_richness(OSTComparison, x="Date",color="Treatment",shape="Treatment", measures=c("Shannon"))+ylab("Shannon Diversity")+geom_point()+scale_fill_manual(cbPalette)+geom_point(size=5)+facet_wrap(~Site)
ShannonOST



set.seed(12345)
ord=ordinate(OSTComparison,"PCoA", "bray")
ordplot=plot_ordination(OSTComparison, ord,"samples", color="Site",shape="Treatment")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+facet_wrap(RiverType1~Treatment)+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

set.seed(12345)
ord=ordinate(OSTComparison,"PCoA", "bray")
ordplot=plot_ordination(OSTComparison, ord,"samples", color="Date",shape="Treatment")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+facet_wrap(RiverType1~Treatment)+
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))




GPr  = transform_sample_counts(OSTComparison, function(x) x / sum(x) ) #transform samples based on relative abundance
PhylumAll=tax_glom(GPr, "Phylum")

PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 3e-2, TRUE) #filter out any taxa lower than 3%

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum","Treatment","Site","RiverType1","Date"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
TrtdataSubset<-subset(Trtdata, Treatment=="C")
PhylumPlotTreatment=ggplot(TrtdataSubset, aes(x=Site,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+
  facet_grid(RiverType1~Date,scale="free_x")+xlab("Date (Only C shown)")+ylab("Relative Abundance (> 3%, SEM)")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
PhylumPlotTreatment



```



Package info

```{r}
sessionInfo()

```