install.packages("ggplot2")
library("ggplot2")
install.packages("scales")
library(scales) 

#Working directory needs to be set to folder containing data file

data<- read.csv("Data_final.csv")
data$media_stoich<- as.factor(data$media_stoich)

#Subset values by Media C:P
attach(data) 
CP_100 <- data[ which(media_stoich=='100'),] 

CP_500 <- data[ which(media_stoich=='500'),]

CP_1000 <- data[ which(media_stoich=='1000'),]

detach(data)

#Combine CP_500 and CP_1000

CP_other <- rbind(CP_500, CP_1000)

# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

#Figure 1: Residual Glucose concentrations as a percetnage of the residual DOC pool
figure_1 <- ggplot(data, aes(x=media_stoich,y=glucose_percent,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  ylab("Residual Glucose (%)")+
  xlab("Media C:P")+ 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))

figure_1
ggsave("figure_1_updated.pdf", width = 8, height = 6, units = c("in"))

#Figure 2: Concentration of DOC per unit biomass (µM per mg biomass) produced by flexible and infelxible strains at each media level

figure_2 <- ggplot(data, aes(x=media_stoich,y=doc_produced_biomass,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  ylab("DOC produced (µM per mg biomass)")+
  xlab("Media C:P")+
  scale_y_log10(limits = c(100,10000))+
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))

figure_2
ggsave("figure_2_updated.pdf", width = 8, height = 6, units = c("in"))

#Figure 3: SRP removal efficiency by flexible and inflexible bacterial strains grown at three different media C:P ratios.
#Create Box Plot of SRP Removed @ C:P 100
figure_3a <- ggplot(CP_100, aes(x=media_stoich,y=percent_srp_removed,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  ylab("SRP Removed (%)")+
  xlab("Media C:P")+ 
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),legend.position="none",
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))

#Create Box Plot of C:P of residual media at 500 & 1000 and color by Flexibility status
figure_3b <- ggplot(CP_other, aes(x=media_stoich,y=percent_srp_removed,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  ylab("SRP Removed (%)")+
  xlab("Media C:P")+ 
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),legend.position=c(.75, .1),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))

multiplot(figure_3a,figure_3b,cols=2)


#Figure 4: Concentration of DOP per unit biomass (µM per mg biomass) produced by flexible and infelxible strains at each media level
figure_4 <- ggplot(data, aes(x=media_stoich,y=min_dop_produced_biomass,fill=flex_qual))+geom_boxplot()+
  scale_fill_manual(values=c("grey40", "grey80")) + ylab("DOP (µM per mg biomass)") +xlab("Media C:P")+
  scale_y_continuous(trans=log10_trans())+
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))
figure_4
ggsave("figure_4_updated.pdf", width = 8, height = 6, units = c("in"))

#Figure 5: Percentage DOP production (log scale) relative to initial media phosphate concentration.
figure_5 <- ggplot(data, aes(x=media_stoich,y=percent_dop_produced,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  scale_y_continuous(trans=log10_trans(), labels = scales::percent_format(accuracy = .1))+
  ylab("DOP Produced (%)")+
  xlab("Media C:P")+ 
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))
figure_5
ggsave("figure_5_updated.pdf", width = 8, height = 6, units = c("in"))

#Figure 6: Specific-UV-Absorbance at wavelength 254 (SUVA254, L * mg-1 * L-1) for flexible and inflexible strains grown at three different media C:P ratios
figure_6 <- ggplot(data, aes(x=media_stoich,y=suva_produced,fill=flex_qual))+geom_boxplot()+
  scale_fill_manual(values=c("grey40", "grey80")) + ylab(expression(SUVA[254])) +xlab("Media C:P")+
  ylim(0,3.2)+
  geom_hline(yintercept = 0.136, size=1.25)+
  annotate("text",x=2.25,y=6000 , label = "")+
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))
figure_6
ggsave("figure_6_updated.pdf", width = 8, height = 6, units = c("in"))

#Figure 7: C:P values in residual media for inflexible and flexible strains across all three media C:P levels. 
figure_7 <- ggplot(data, aes(x=media_stoich,y=cp_produced_biomass,fill=flex_qual))+geom_boxplot()+
  scale_fill_manual(values=c("grey40", "grey80")) + ylab("DOC:DOP Produced") +xlab("Media C:P")+
  scale_y_continuous(trans=log10_trans())+
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))
figure_7
ggsave("figure_7_updated.pdf", width = 8, height = 6, units = c("in"))


#Figure 8: Biomass production by flexible and inflexibly strains across all three media C:P levels. 
figure_8 <- ggplot(data, aes(x=media_stoich,y=mean_biomass,fill=flex_qual))+geom_boxplot()+ 
  scale_fill_manual(values=c("grey40", "grey80"))+ 
  ylab("Biomass (mg)")+
  xlab("Media C:P")+ 
  theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
        panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
        axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))
figure_8
ggsave("figure_8_updated.pdf", width = 8, height = 6, units = c("in"))


#Supplement Figure 1: Scatterplot of SUVA and Produced SUVA with a 1:1 reference line to examine how much the glucose correction changed the SUVA vlaues. Turns out not very much. 
Supplement_Figure_1 <- ggplot(data, aes(x=suva_produced,y=suva))+
geom_point(shape=3, size=4)+
geom_abline(slope= 1, intercept = 0, color="Black", linetype="dashed")+
ylab(expression(SUVA[254])) +
xlab(expression(Produced_SUVA[254]))+
theme(panel.grid.major = element_blank(), axis.line = element_line(colour = "black"),
panel.background = element_blank(),legend.title=element_blank(), legend.text=element_text(size=18), 
axis.title.x = element_text(size=20),axis.title.y = element_text(size=20), axis.text.y = element_text(size=14),axis.text.x = element_text(size=14))

Supplement_Figure_1
ggsave("sup_figure_1.pdf", width = 6, height = 4, units = c("in"))


#Figure 1 analysis Kruskal-Wallis Test for residual glucose 
kruskal.test(glucose_um ~  media_stoich, data = data)
kruskal.test(glucose_um ~  flex_qual, data = data)
wilcox.test(glucose_um ~  flex_qual, data = CP_100)
wilcox.test(glucose_um ~  flex_qual, data = CP_500)
wilcox.test(glucose_um ~  flex_qual, data = CP_1000)

#Figure 2 analysis Kruskal-Wallis Test for DOC produced for unit biomass
kruskal.test(doc_produced_biomass ~  media_stoich, data = data)
kruskal.test(doc_produced_biomass ~  flex_qual, data = data)
wilcox.test(doc_produced_biomass ~  flex_qual, data = CP_100)
wilcox.test(doc_produced_biomass ~  flex_qual, data = CP_500)
wilcox.test(doc_produced_biomass ~  flex_qual, data = CP_1000)

#Figure 3 analysis Kruskal-Wallis Test for %SRP Removed
kruskal.test(percent_srp_removed ~  media_stoich, data = data)
kruskal.test(percent_srp_removed ~  flex_qual, data = data)
wilcox.test(percent_srp_removed ~  flex_qual, data = CP_100)
wilcox.test(percent_srp_removed ~  flex_qual, data = CP_500)
wilcox.test(percent_srp_removed ~  flex_qual, data = CP_1000)

#Figure 4 analysis Kruskal-Wallis Test for minimum DOP production per unit biomass
kruskal.test(min_dop_produced_biomass ~  media_stoich, data = data)
kruskal.test(min_dop_produced_biomass ~  flex_qual, data = data)
wilcox.test(min_dop_produced_biomass ~  flex_qual, data = CP_100)
wilcox.test(min_dop_produced_biomass ~  flex_qual, data = CP_500)
wilcox.test(min_dop_produced_biomass ~  flex_qual, data = CP_1000)

#Figure 5 analysis Kruskal-Wallis Test for percent DOP production 
kruskal.test(percent_dop_produced ~  media_stoich, data = data)
kruskal.test(percent_dop_produced ~  flex_qual, data = data)
wilcox.test(percent_dop_produced ~  flex_qual, data = CP_100)
wilcox.test(percent_dop_produced ~  flex_qual, data = CP_500)
wilcox.test(percent_dop_produced ~  flex_qual, data = CP_1000)

#Figure 6 analysis Kruskal-Wallis Test for percent SUVA production 
kruskal.test(suva_produced ~  media_stoich, data = data)
kruskal.test(suva_produced ~  flex_qual, data = data)
wilcox.test(suva_produced ~  flex_qual, data = CP_100)
wilcox.test(suva_produced ~  flex_qual, data = CP_500)
wilcox.test(suva_produced ~  flex_qual, data = CP_1000)

#Figure 7 analysis Kruskal-Wallis Test for DOC:DOP Produced
kruskal.test(cp_produced_biomass ~  media_stoich, data = data)
kruskal.test(cp_produced_biomass ~  flex_qual, data = data)
wilcox.test(cp_produced_biomass ~  flex_qual, data = CP_100)
wilcox.test(cp_produced_biomass ~  flex_qual, data = CP_500)
wilcox.test(cp_produced_biomass ~  flex_qual, data = CP_1000)

#Figure 8 analysis Kruskal-Wallis Test for Biomass
kruskal.test(mean_biomass ~  media_stoich, data = data)
kruskal.test(mean_biomass ~  flex_qual, data = data)
wilcox.test(mean_biomass ~  flex_qual, data = CP_100)
wilcox.test(mean_biomass ~  flex_qual, data = CP_500)
wilcox.test(mean_biomass ~  flex_qual, data = CP_1000)