# Pathogen: Pe = P. entomophila; Mb = M. brunneum
#HL UC = hemolymph unchallenged
#HL SI = hemolymph 24h after sterile injury (TrisNacl injection)
#HL PP = hemolymph 24h after 'priming' with heat-killed pathogen matching the one being tested


library("data.table")
library("tidyverse")
library("ggplot2")
library("outliers")
library("lme4")
library("multcomp")
library("car")

data_HL_antimicrobial <- fread("Fig5AB_data.csv")

#PE PLOT

data_HL_antimicrobial$Treatment <- factor(data_HL_antimicrobial$Treatment, levels = c("TrisNacl", "HLUC", "HLSI", "HLPP"))


data_HL_antimicrobial %>%
  filter(Pathogen %in% "Pe") %>%
  ggplot() +
  aes(x = Treatment, y = Normalized_colony_count) +
  geom_boxplot(shape = "circle", fill = c("#c0c0c0", "#6080ab","#51A066"), outlier.shape = NA) +
  theme_minimal() +
  theme_bw() +
  scale_y_continuous(limits = function(x) c(0, 2.5)) +
  labs(y = bquote(italic("P. entomophila") ~ "viability"), x = "") +
  scale_x_discrete(labels = c("HLUC" = "Unchallenged", "HLSI" = "Sterile injury", "HLPP" = "Pe injury")) +  # Replace x-axis labels
  theme(
    text = element_text(size = 12),  # Set the overall text size to 12
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, color = "black"),  # Increase x-axis text size
    axis.text.y = element_text(size = 12, color = "black"),  # Increase y-axis text size
    axis.title = element_text(size = 12),  # Increase axis title size
    strip.text = element_text(size = 14)) +  # Increase facet label size
  geom_point(position = position_jitter(width = 0.1), size = 0.9, shape = 16) +
  geom_text(x = 1, y = 2.5, label = "a", size = 4) +
  geom_text(x = 2, y = 2.5, label = "a", size = 4) +
  geom_text(x = 3, y = 2.5, label = "b", size = 4) -> plotPe
print(plotPe)

#MB plot

data_HL_antimicrobial %>%
  filter(Pathogen %in% "Mb") %>%
  ggplot() +
  aes(x = Treatment, y = Normalized_colony_count) +
  geom_boxplot(shape = "circle", fill = c("#c0c0c0", "#6080ab","#9c8762"), outlier.shape = NA) +
  theme_minimal() +
  theme_bw() +
  scale_y_continuous(limits = function(x) c(0, 2.5)) +
  labs(y = bquote(italic("M. brunneum") ~ "viability"), x = "") +
  scale_x_discrete(labels = c("HLUC" = "Unchallenged", "HLSI" = "Sterile injury", "HLPP" = "Mb injury")) +  # Replace x-axis labels
  theme(
    text = element_text(size = 12),  # Set the overall text size to 12
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, color = "black"),  # Increase x-axis text size
    axis.text.y = element_text(size = 12, color = "black"),  # Increase y-axis text size
    axis.title = element_text(size = 12),  # Increase axis title size
    strip.text = element_text(size = 14)) +  # Increase facet label size
  geom_point(position = position_jitter(width = 0.1), size = 0.9, shape = 16) +
  geom_text(x = 1, y = 2.5, label = "a,b", size = 4) +
  geom_text(x = 2, y = 2.5, label = "a", size = 4) +
  geom_text(x = 3, y = 2.5, label = "b", size = 4) -> plotMb
print(plotMb)


#Model PE

data_HL_antimicrobial %>%
  filter(Pathogen %in% "Pe") -> dataPe

modelPe <- lmer((sqrt(Normalized_colony_count)) ~ Treatment + (1|Replicate), family="gaussian", data=dataPe)

residuals_modelPe <- residuals(modelPe)
shapiro.test(residuals_modelPe)

summary(modelPe)

Anova(modelPe, type="II") 

post_hocPe <- summary(glht(modelPe, linfct=mcp(Treatment="Tukey")), test=adjusted("BH"))
print(cld(post_hocPe))

#Model MB

data_HL_antimicrobial %>%
  filter(Pathogen %in% "Mb") -> dataMb

modelMb <- lmer(((Normalized_colony_count)) ~ Treatment + (1|Replicate), data=dataMb)

residuals_modelMb <- residuals(modelMb)
hist(residuals_modelMb)
shapiro.test(residuals_modelMb)

summary(modelMb)

Anova(modelMb, type="II")