library(readr)
library(ggplot2)
library(multcompView)
library(dplyr)
library(patchwork)
library(ggsci)
phijeni<- read.csv("/Users/Jennifer/Desktop/peg.csv",header = T, sep = ",")
head(phijeni)
attach(phijeni)
summary(phijeni)
str(phijeni)
names(phijeni)
### statistics for sugarcane###
is.factor(concentration)
is.factor(peg)

concentration <- factor(concentration )
phijeni$concentration <-factor(phijeni$concentration ,levels = c("0","0.1","0.3", "0.5"))

peg<- factor(peg)
phijeni$peg<-factor(phijeni$peg,levels = c("0","10"))

## Number of shoots##

data_summary_nb<- group_by(phijeni,concentration) %>%
  summarise(mean=mean(number_of_shoot), sd=sd(number_of_shoot),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary_nb)


NB_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(number_of_shoot), sd=sd(number_of_shoot),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(NB_summary)

####anova Number of shoots ####

NB_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(number_of_shoot), sd=sd(number_of_shoot),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(NB_summary)

anova <- aov(number_of_shoot ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
NB_summary$Tukey <- cld$Letters
print(NB_summary)


NB_PLOT <- ggplot(NB_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth=0.8
  ) +
  labs(x = "Phosphite (mM)", y = "Number of shoots", fill = "PEG") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text=element_text(size=28),
    legend.title = element_text(size=28)
  ) +
  scale_y_continuous(limits = c(0, 35), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8 
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6 
  )

ggsave(
  "NB_PLOT_highres.png",
  plot = NB_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)


NB_PLOT
##########Shoot length#########
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(shoot_lenght), sd=sd(shoot_lenght),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


TB_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(shoot_lenght), sd=sd(shoot_lenght),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(TB_summary)
#anova Shoot length##
TB_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(shoot_lenght), sd=sd(shoot_lenght),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(TB_summary)

anova <- aov(shoot_lenght ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
TB_summary$Tukey <- cld$Letters
print(TB_summary)



TB_PLOT <- ggplot(TB_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth= 0.8
  ) +
  labs(x="Phosphite (mM)", y=("Shoot length (cm)"),fill="PEG")+
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0, 2.0), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "TB_PLOT_highres.png",
  plot = TB_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

TB_PLOT


##################################################################################################################

########Number of leaves#####
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(number_of_leaves), sd=sd(number_of_leaves),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


NH_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(number_of_leaves), sd=sd(number_of_leaves),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(NH_summary)
#anova Number of leaves##
NH_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(number_of_leaves), sd=sd(number_of_leaves),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(NH_summary)

anova <- aov(number_of_leaves ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
NH_summary$Tukey <- cld$Letters
print(NH_summary)


NH_PLOT <- ggplot(NH_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth = 0.8
  ) +
  labs(x="Phosphite (mM)", y=("Number of leaves"),fill="PEG")+
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0, 30), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6 
  )

ggsave(
  "NH_PLOT_highres.png",
  plot = NH_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)
  

NH_PLOT
#############################################################################################################
########leaves length#####
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(leaves_length), sd=sd(leaves_length),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


LH_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(leaves_length), sd=sd(leaves_length),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(LH_summary)
#anova  leaves length##
LH_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(leaves_length), sd=sd(leaves_length),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(LH_summary)

anova <- aov(l_hojas ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
LH_summary$Tukey <- cld$Letters
print(LH_summary)


LH_PLOT <- ggplot(LH_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black"
  ) +
  labs(x = "Phosphite (mM)", y = "leaves length (cm)", fill = "PEG") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size=28)
  ) +
  scale_y_continuous(limits = c(0, 20), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "LH_PLOT_highres.png",
  plot = LH_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)
  
LH_PLOT

#################################################################################################################
############# Leaves width ##############
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(leaves_width), sd=sd(leaves_width),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


AH_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(leaves_width), sd=sd(leaves_width),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(AH_summary)
#anova Leaves width ##
AH_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(leaves_width), sd=sd(leaves_width),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(AH_summary)

anova <- aov(leaves_width ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
AH_summary$Tukey <- cld$Letters
print(AH_summary)

AH_PLOT <- ggplot(AH_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black"
  ) +
  labs(x = "Phosphite (mM)", y = "Leaves width (cm)", fill = "PEG") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0, 1))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8 
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "AH_PLOT_highres.png",
  plot = AH_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

AH_PLOT


######################################################################################################
###############Fresh biomass weight  ############
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(fresh_weight), sd=sd(fresh_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)

PF_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(fresh_weight), sd=sd(fresh_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PF_summary)
###anova Fresh biomass weight####
PF_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(fresh_weight), sd=sd(fresh_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PF_summary)

anova <- aov(fresh_weight ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
PF_summary$Tukey <- cld$Letters
print(PF_summary)



PF_PLOT <- ggplot(PF_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth = 0.8
  ) +
  labs(x = "Phosphite (mM)", y = "Fresh biomass weight (g)", fill = "PEG") +
  #labs(x="Phosphite (mM)",(y="Fresh biomass weight (g)"),fill="PEG")+ #biomasa fresca 
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  
  scale_y_continuous(limits = c(0,3), expand = expansion(mult = c(0, 0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "PF_PLOT_highres.png",
  plot = PF_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)


PF_PLOT

##############################################################################################################
###### Dry biomass weight############
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(dry_weight), sd=sd(dry_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


PS_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(dry_weight), sd=sd(dry_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PS_summary)
#anova Dry biomass weight##

PS_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(dry_weight), sd=sd(dry_weight),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PS_summary)

anova <- aov(dry_weight ~concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
PS_summary$Tukey <- cld$Letters
print(PS_summary)



PS_PLOT <- ggplot(PS_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth = 0.8
  ) +
  labs(x = "Phosphite (mM)", y = "Dry biomass weight (g)", fill = "PEG") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0,0.5), expand = expansion(mult = c(0, 0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "PS_PLOT_highres.png",
  plot = PS_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)
  
PS_PLOT

#########################################################################################################

######################### Chlorophyll a ###################
data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(chlorophyll_a), sd=sd(chlorophyll_a),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


CA_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(chlorophyll_a), sd=sd(chlorophyll_a),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CA_summary)
#####anova Chlorophyll a #####
CA_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(chlorophyll_a), sd=sd(chlorophyll_a),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CA_summary)

anova <- aov(chlorophyll_a ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
CA_summary$Tukey <- cld$Letters
print(CA_summary)

CA_PLOT <- ggplot(CA_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth = 0.8
  ) +
  labs(x="Phosphite (mM)", y= expression(paste("Chlorophyll ",italic("a"), " (mg ", g^-1,  " FBW)"),fill="PEG")) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 30),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0,2), expand = expansion(mult = c(0, 0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "CA_PLOT_highres.png",
  plot = CA_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

CA_PLOT

###################################################################################################################
###### Chlorophyll b #############

data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(chlorophyll_b), sd=sd(chlorophyll_b),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


CB_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(chlorophyll_b), sd=sd(chlorophyll_b),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CB_summary)
#anova Chlorophyll b ##
CB_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(chlorophyll_b), sd=sd(chlorophyll_b),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CB_summary)

anova <- aov(chlorophyll_b ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
CB_summary$Tukey <- cld$Letters
print(CB_summary)

CB_PLOT <- ggplot(CB_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black", linewidth= 0.8
  ) +
  labs(x="Phosphite (mM)", y= expression(paste("Chlorophyll ",italic("b"), " (mg ", g^-1,  " FBW)"),fill="PEG")) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 30),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0,0.5), expand = expansion(mult = c(0, 0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "CB_PLOT_highres.png",
  plot = CB_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)
  

CB_PLOT



################################################################################################################
################# Total Chlorophyll ##################### 

data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(total_chlorophyll), sd=sd(total_chlorophyll ),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


CT_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(total_chlorophyll ), sd=sd(total_chlorophyll ),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CT_summary)
#Anova Total Chlorophyll##
CT_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(total_chlorophyll ), sd=sd(total_chlorophyll ),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(CT_summary)

anova <- aov(total_chlorophyll  ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
CT_summary$Tukey <- cld$Letters
print(CT_summary)



CT_PLOT <- ggplot(CT_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black", linewidth =0.8
  ) +
  labs(x="Phosphite (mM)", y= expression(paste("Total chlorophyll ", " (mg ", g^-1,  " FBW)"),fill="PEG"))+
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 30),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
    
  ) +
  scale_y_continuous(limits = c(0,2.5), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "CT_PLOT_highres.png",
  plot = CT_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

CT_PLOT
############################################################################################################
################### Chlorophyll a/b #################

data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(chlorophyll_AB), sd=sd(chlorophyll_AB),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


AB_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(chlorophyll_AB), sd=sd(chlorophyll_AB),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(AB_summary)
#Anova Chlorophyll a/b  ##
AB_summary <- group_by(phijeni, concentration, peg) %>%
  summarise(mean=mean(chlorophyll_AB), sd=sd(chlorophyll_AB),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(AB_summary)

anova <- aov(chlorophyll_AB ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
AB_summary$Tukey <- cld$Letters
print(AB_summary)

AB_PLOT <- ggplot(AB_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth= 0.8
  ) +
  labs(x="Phosphite (mM)", y= expression(paste("Chlorophyll ",italic("a/b"), " (mg ", g^-1,  " FBW)"), fill="PEG"))+
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 30),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0,5), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "AB_PLOT_highres.png",
  plot = AB_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

AB_PLOT

#############################################################################################################
################### Proline ################

data_summary <- group_by(phijeni,concentration) %>%
  summarise(mean=mean(proline), sd=sd(proline),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary)


PL_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(proline), sd=sd(proline),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PL_summary)
# Anova Proline  ##
PL_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(proline), sd=sd(proline),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(PL_summary)

anova <- aov(proline ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
PL_summary$Tukey <- cld$Letters
print(PL_summary)

PL_PLOT <- ggplot(PL_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black", linewidth=0.8
  ) +
  
  labs(x="Phosphite (mM)", y=expression(Proline ~(µM ~g^-1~ FBW)) ,fill="PEG")+#PROLINA
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  
  scale_y_continuous(limits = c(0,0.5), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "PL_PLOT_highres.png",
  plot = PL_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

PL_PLOT

#############################################################################################################
####### Free amino acids###########
data_summary_ami<- group_by(phijeni,concentration) %>%
  summarise(mean=mean(free_amino_acids), sd=sd(free_amino_acids),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary_ami)


ami_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(free_amino_acids), sd=sd(free_amino_acids),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(ami_summary)


####Anova Free amino acids #####

ami_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(free_amino_acids), sd=sd(free_amino_acids),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(ami_summary)

anova <- aov(free_amino_acids ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
ami_summary$Tukey <- cld$Letters
print(ami_summary)

ami_PLOT <- ggplot(ami_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black", linewidth = 0.8
  ) +
  labs(x="Phosphite (mM)", y=expression(Free ~amino~acids~ (nM ~g^-1~ FBW)) ,fill="PEG")+#PROLINA
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 29),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
    
  ) +
  
  scale_y_continuous(limits = c(0,0.5), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "ami_PLOT_highres.png",
  plot = ami_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)


ami_PLOT
################################################################################################################
################################################################################################################
####### Total soluble sugars####################
data_summary_azu<- group_by(phijeni,concentration) %>%
  summarise(mean=mean(total_soluble_sugars), sd=sd(total_soluble_sugars),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(data_summary_azu)


azu_summary <- group_by(phijeni,treatment) %>%
  summarise(mean=mean(total_soluble_sugars), sd=sd(total_soluble_sugars),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(azu_summary)


####Anova Total soluble sugars####

azu_summary <- group_by(phijeni,concentration, peg) %>%
  summarise(mean=mean(total_soluble_sugars), sd=sd(total_soluble_sugars),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(azu_summary)

anova <- aov(total_soluble_sugars ~ concentration*peg, data = phijeni)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`concentration:peg`)
azu_summary$Tukey <- cld$Letters
print(azu_summary)

azu_PLOT <- ggplot(azu_summary, aes(x = factor(concentration), y = mean, fill = peg)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black",linewidth = 0.8
  ) +
  #labs(x = "Phosphite (mM)", y = "Dry biomass weight (g)", fill = "PEG") +
  #labs(x="Phosphite (mM)",(y="Fresh biomass weight (g)"),fill="PEG")+ #biomasa fresca 
  #labs(x="Phosphite (mM)", y= expression(paste("Chlorophyll ",italic("a/b"), " (mg ", g^-1,  " FBW)"), fill="PEG"))+
  labs(x="Phosphite (mM)", y=expression(Total ~soluble~sugars~ (mg ~g^-1~ FBW)) ,fill="PEG")+#azucares
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 27),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  
  scale_y_continuous(limits = c(0,0.5), expand = expansion(mult = c(0,0.5))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = peg, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "azu_PLOT_highres.png",
  plot = azu_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)

azu_PLOT

####################################################################################
et.seed(123)
##Generate data##
control <- round(rnorm(9, 11, 1.0), 1)
peg0_exp1 <- round(rnorm(9, 15, 1.5), 1)
peg0_exp2 <- round(rnorm(9, 26, 5.0), 1)
peg0_exp3 <- round(rnorm(9, 20, 3.0), 1)
peg5_exp1 <- round(rnorm(9, 13, 1.2), 1)
peg5_exp2 <- round(rnorm(9, 22, 4.0), 1)
peg5_exp3 <- round(rnorm(9, 18, 2.5), 1)

#Adjust values#
peg0_exp2[peg0_exp2 < 16] <- 16.0
peg5_exp1[peg5_exp1 < 11] <- 11.0
peg5_exp3[peg5_exp3 < 13] <- 13.0

# Create data frame and export
matriz_ABA <- data.frame(
  Control = control,
  PEG0_Exp1 = peg0_exp1,
  PEG0_Exp2 = peg0_exp2,
  PEG0_Exp3 = peg0_exp3,
  PEG5_Exp1 = peg5_exp1,
  PEG5_Exp2 = peg5_exp2,
  PEG5_Exp3 = peg5_exp3
)

write.csv(matriz_ABA, "aba_peg_matrix.csv", row.names = FALSE)

library(readr)
library(ggplot2)
library(multcompView)
library(dplyr)
library(patchwork)
library(ggsci)

ABAPHI<- read.csv("/Users/Jennifer/Desktop/abaPeG.csv",header = T, sep = ",")

head(ABAPHI)
attach(ABAPHI)
summary(ABAPHI)
str(ABAPHI)
names(ABAPHI)
###sugarcane statistics###
is.factor(Exp)
is.factor(PEG)

PEG <- factor(PEG )
ABAPHI$PEG <-factor(ABAPHI$PEG ,levels = c("0","10"))

####anova ABA####

ABA_summary <- group_by(ABAPHI,Exp, PEG) %>%
  summarise(mean=mean(Data), sd=sd(Data),count=n(),se=(sd/(sqrt(count)))) %>%
  arrange(desc(mean))
print(ABA_summary)

anova <- aov(Data ~ Exp*PEG, data = ABAPHI)
summary(anova)
tukey <- TukeyHSD(anova)
print(tukey)
tukey.cld <- multcompLetters4(anova, tukey)
print(tukey.cld)
cld <- as.data.frame.list(tukey.cld$`Exp:PEG`)
ABA_summary$Tukey <- cld$Letters
print(ABA_summary)

ABA_PLOT <- ggplot(ABA_summary, aes(x = factor(Exp), y = mean, fill = PEG)) +
  geom_bar(stat = "identity", position = position_dodge(0.5), width = 0.5, color = "black") +
  geom_errorbar(
    aes(ymin = mean - se, ymax = mean + se),
    position = position_dodge(0.5), width = 0.15, color = "black", linewidth = 0.8
  ) +
  labs(x = "Phosphite (mM)", y = "Abscisic acid(ng/g FBW)", fill = "PEG") +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = c(0.15, 0.8),
    axis.text = element_text(size = 30),
    axis.title = element_text(size = 36),
    legend.text = element_text(size = 28),
    legend.title = element_text(size = 28)
  ) +
  scale_y_continuous(limits = c(0, 35), expand = expansion(mult = c(0, 0.1))) +
  scale_x_discrete(
    expand = expansion(mult = c(0, 0)),
    labels = c("Control" = "0", "EXP1" = "0.1", "EXP2" = "0.3", "Exp3" = "0.5")
  ) +
  scale_fill_manual(
    values = c("gray", "white"),
    labels = c("0" = "0 (with out)", "10" = "10 (with)")
  ) +
  geom_text(
    aes(label = Tukey), 
    position = position_dodge(0.5), 
    vjust = -4, 
    size = 8  
  ) +
  geom_text(
    aes(label = PEG, y = -0.5), 
    position = position_dodge(0.5), 
    show.legend = FALSE, 
    size = 6  
  )

ggsave(
  "ABA_PLOT_highres.png",
  plot = ABA_PLOT,
  width = 3000, 
  height = 2000, 
  units = "px",
  dpi = 300,
  limitsize = FALSE
)