Manuscript_DatCol_HNO_Complementation.qmd

---
format:
  html:
    suppress-bibliography: true
---

# CFUs Complementation {.unnumbered}

```{r}
# loading packages and functions
library(tidyverse)
library(kableExtra)
library(ggtext)
library(scales)
library(lme4)
library(lmerTest)
library(emmeans)
library(ggh4x)

source("R/LMM_functions.R")
```

## Data Input

All data inputs and outputs in this repo are relative to **FOLDER_PATH**. Please, save the provided data frames in RDS format in the `dataframes` folder.

```{r}
# Load the BOX_PATH environment variable
FOLDER_PATH <- Sys.getenv("BOX_PATH") 

# Define subdirectories within the Box folder
dataframes_folder <- file.path(FOLDER_PATH, "dataframes")

outputs_folder <- file.path(FOLDER_PATH, "outputs/CFUs")
if (!dir.exists(outputs_folder)) {
  dir.create(outputs_folder, recursive = TRUE)
}
```

```{r}
# Read the provided RDS file
dataframeID <- "DatCol"
CFUs_HNO_Genetic <- readRDS(file.path(dataframes_folder, paste0(dataframeID, "_CFUs_HNO_Genetic.rds")))
CFUs_HNO_Chemical <- readRDS(file.path(dataframes_folder, paste0(dataframeID, "_CFUs_HNO_Chemical.rds")))

# Filter the data for 24h time point and drop unused factor levels
CFUs_24h_Gen <- CFUs_HNO_Genetic %>%
  filter(time == 24) %>%
  droplevels()

CFUs_24h_Chem <- CFUs_HNO_Chemical %>%
  filter(time == 24) %>%
  droplevels() 
```

```{r}
# Function to create a ggplot with jittered points, emmeans predictions with error bars, and faceting. 
plot_CFUs_function <- function(formated_stats_results, ymin_exp, ymax_exp, title_x = " ", title_y = " ",
                          x_var, y_var, fill_var, shape_var, facet_var,
                          fill_colors, shape_leves, facet_colors, facet_labels, fill_title, shape_title) {
  plot <- ggplot(data = formated_stats_results$input_data) +
    geom_jitter(
      aes(x = {{ x_var }}, y = {{ y_var }}, fill = {{ fill_var }}, shape = {{ shape_var }}),
      size  = 4, stroke = 0.01, width = 0.2, color = "white") +
    
    geom_point(data = formated_stats_results$predictions, aes(x = {{ x_var }}, y = predicted_mean), shape = 3, size = 1) +
    geom_errorbar(data = formated_stats_results$predictions, aes(x = {{ x_var }},
                                         y = predicted_mean,
                                         ymin = predicted_lower,
                                         ymax = predicted_upper),
                  width = .15) +
    
    facet_grid2(
      rows = NULL, cols = vars({{ facet_var }}),
      scales = "free_x", space = "free_x",
      strip = strip_themed(background_x = elem_list_rect(fill = facet_colors)),
      labeller = labeller(.cols = facet_labels)) +
    
    scale_fill_manual(values = fill_colors, name = fill_title) +
    scale_shape_manual(values = shape_leves, name = shape_title) +
    
    guides(
      fill  = guide_legend(override.aes = list(shape = 21, color = "white", size = 3, stroke = 2)),
      shape = guide_legend(override.aes = list(fill = "black", size = 3, stroke = 2))) +
    
    scale_y_log10(
      limits = c(10^ymin_exp, 10^ymax_exp),
      breaks = 10^(1:ymax_exp),
      labels = trans_format("log10", math_format(10^.x))) +
    labs(x = title_x, y = title_y) +
    theme_bw() +
    theme(
      panel.grid    = element_blank(),
      panel.spacing = unit(0, "lines"),
      strip.text.x = element_markdown(size = 14, face = "bold", color = "white"),
      text          = element_text(size = 16, color = "black"),
      axis.text.x   = element_markdown(angle = 45, hjust = 1, color = "black"),
      axis.text.y   = element_text(color = "black")
    )
  
  return(plot)
}
```

## Genetic Complementation

Log-transformed CFU/HNO-ALI between groups were compared using a LMM to determine statistical significance with `parent_strain` and `vector` as fixed effects and experimental date as a random effect (line was singular). 
We analyzed pairwise contrasts for the `parent_strain` and `vector` effects separately, comparing each condition only to the reference for that group (JE2 for parent_strain and pLL39 for vector).
We then adjusted the p-values for multiple comparisons using the Holm method and created a `sign` annotation column to indicate which comparisons are statistically significant with p<0.05.
Vertical brackets represent model-predicted mean values and confidence intervals (± 2 × SEM), back-transformed to the original response scale. 

```{r}
# Run the mixed model with emmeans contrasts
stats_results_Gen <- mixed_model_with_emmeans(
  input_data = CFUs_24h_Gen, response = "CFU_HNO", log_response = TRUE, adjust = "none",
  fixed_effects = c("vector", "parent_strain"), random_effects1 = "line", random_effects2 = "date"
)
stats_formated_Gen <- format_stats_results(stats_results_Gen, grouping_vars = c("condition", "label", "plot_label"))

# Filter the contrasts to only compare all conditions to the parent strain JE2 carrying the same vector and to compare each vector-bearing mutant to its isogenic empty-vector (pLL39) mutant strain. 
# Pvalues are adjusted using the Holm method only for the contrasts that are being kept and the significance annotation is added accordingly with "a" for the parent strain comparisons and "b" for the vector comparisons.
stats_formated_Gen$main <- stats_formated_Gen$main %>%
  filter(
    (effect_name == "parent_strain" & str_detect(condition1, "JE2")) |
    (effect_name == "vector" & str_detect(condition1, "pLL39"))
  ) %>% 
  filter(is.na(predicted_mean.2) == FALSE) %>% # Filter out comparisons that do not exist in the data
  mutate(
    p.value.adj = p.adjust(p.value, method = "holm"),
    sign = case_when(
      p.value.adj < 0.05 & effect_name == "parent_strain" ~ "a",
      p.value.adj < 0.05 & effect_name == "vector" ~ "b",
      TRUE ~ ""
    )
  )

# Save the stats results
saveRDS(stats_formated_Gen, file.path(outputs_folder, "stats_results_HNO_CFUs_Genetic.rds"))
```

```{r}
# Render the results
cat("Model Formula:", format(stats_results_Gen$model_formula))
kable(stats_formated_Gen$random_effects, caption = "Random Effects")
kable(stats_formated_Gen$anova, caption = "ANOVA Table")
kable(stats_formated_Gen$main %>%
    mutate(
      p.value.anova = formatC(p.value.anova, format = "e", digits = 2),
      p.value = formatC(p.value, format = "e", digits = 2),
            p.value.adj = formatC(p.value.adj, format = "e", digits = 2)
    ) %>%
    select(effect_name, p.value.anova,  condition1, condition2, predicted_mean.1, predicted_mean.2, FC, log2FC, p.value, p.value.adj, sign),
  caption = "Main Contrasts with Predicted Log2FC"
)
```

```{r}
# Plot
plot_HNO_CFUs_Gen <- plot_CFUs_function(
  formated_stats_results = stats_formated_Gen,
  x_var       = plot_label, title_x = " ",
  y_var       = CFU_HNO, title_y = "CFUs/HNO-ALI", ymin_exp = 1, ymax_exp = 10,
  fill_var    = vector, fill_colors = c("gray50", "#2FD612", "#8DFD65", "#03C0FF"), fill_title = "Vector",
  shape_var   = line, shape_leves = c(21, 22), shape_title = "HNO line",
  facet_var   = parent_strain, facet_colors = c("#5F05AF", "#2FD612", "#03C0FF"),
  facet_labels = c("JE2" = "parent", "dat" = "<i>dat</i>::Tn", "alr1" = "Δ<i> alr1</i>")
)
```

```{r, fig.width=10, fig.height=6, fig.align="center"}
# Add significance annotations
annots <- stats_formated_Gen$main %>%
  select(condition2, sign) %>%
  filter(sign != "") %>%
  left_join(
    stats_formated_Gen$predictions %>%
      select(condition, predicted_upper, plot_label, parent_strain), 
    by = c("condition2" = "condition")
  ) %>%
  mutate(y_pos = predicted_upper * 5) %>%
  arrange(sign) %>%
  group_by(condition2, plot_label, parent_strain, y_pos) %>%
  summarise(sign = paste(sign, collapse = " "), .groups = "drop")

plot_HNO_CFUs_Gen_annot <- plot_HNO_CFUs_Gen +
  geom_text(data = annots, aes(x = plot_label, y = y_pos, label = sign, group = NULL),color = "black", size = 4) 

plot_HNO_CFUs_Gen_annot

# Save the plot
ggsave(filename = file.path(outputs_folder, "HNO_CFUs_Genetic.png"), plot = plot_HNO_CFUs_Gen_annot, width = 10, height = 6, dpi = 300)
saveRDS(plot_HNO_CFUs_Gen_annot, file.path(outputs_folder, "HNO_CFUs_Genetic.rds"))
```

## Chemical Complementation

Log-transformed CFU/HNO-ALI between groups were compared using a LMM to determine statistical significance with `label` (strain/media condition combined varaible) as a fixed effect and experimental date as a random effect (line was singular).
We analyzed pairwise contrasts comparing all conditions to the unsupplemented parent strain (JE2) and all mutant (dat) conditions to the unsupplemented mutant strain.
We then adjusted the p-values for multiple comparisons using the Holm method and created a `sign` annotation column to indicate which comparisons are statistically significant with p<0.05.
Vertical brackets represent model-predicted mean values and confidence intervals (± 2 × SEM), back-transformed to the original response scale. 

```{r}
# Run the mixed model with emmeans contrasts
stats_results_Chem <- mixed_model_with_emmeans(
  input_data = CFUs_24h_Chem, response = "CFU_HNO", log_response = TRUE, adjust = "none",
  fixed_effects = c("label"), random_effects1 = "line", random_effects2 = "date"
)
stats_formated_Chem <- format_stats_results(stats_results_Chem, grouping_vars = c("condition", "plot_label", "parent_strain", "aa"))
```

```{r}
# Filter the contrasts to only include those comparing to the unsupplemented JE2, dat and alr1 strains.
# Pvalues are adjusted using the Holm method only for the contrasts that are being kept and the significance annotation is added accordingly.
stats_formated_Chem$main = stats_formated_Chem$main %>%
  filter(condition1 == "JE2" | condition1 == "dat" | condition1 == "alr1") %>%
  filter(!(grepl("dat", condition1) & grepl("JE2|alr1", condition2))) %>%
  filter(!(grepl("alr1", condition1) & grepl("JE2|dat", condition2))) %>%
  mutate(
    p.value.adj = p.adjust(p.value, method = "holm"),
    sign = case_when(
      p.value.adj < 0.05 & condition1 == "JE2" ~ "a",
      p.value.adj < 0.05 & condition1 == "dat" ~ "b",
      p.value.adj < 0.05 & condition1 == "alr1" ~ "c",
      TRUE ~ ""
    )
  )

# Save the stats results
saveRDS(stats_formated_Chem, file.path(outputs_folder, "stats_results_HNO_CFUs_Chemical.rds"))
```

```{r}
# Render the results
cat("Model Formula:", format(stats_results_Chem$model_formula))
kable(stats_formated_Chem$random_effects, caption = "Random Effects")
kable(stats_formated_Chem$anova, caption = "ANOVA Table")
kable(stats_formated_Chem$main %>%
    mutate(
      p.value.anova = formatC(p.value.anova, format = "e", digits = 2),
      p.value = formatC(p.value, format = "e", digits = 2),
      p.value.adj = formatC(p.value.adj, format = "e", digits = 2)
    ) %>%
    select(effect_name, p.value.anova, contrast, condition1, condition2, FC, log2FC, p.value, p.value.adj, sign),
  caption = "Main Contrasts with Predicted Log2FC"
)
```

```{r}
plot_HNO_CFUs_Chem <- plot_CFUs_function(
  formated_stats_results  = stats_formated_Chem,
  x_var       = plot_label, title_x = " ",
  y_var       = CFU_HNO, title_y = "CFUs/HNO-ALI", ymin_exp = 1, ymax_exp = 10,
  fill_var    = aa, fill_colors = c("gray80", "#09582A", "#A5DEA7", "#F50BEF", "#FFB0FD"), fill_title = "Vector",
  shape_var   = line, shape_leves = c(21, 22, 23), shape_title = "HNO line",
  facet_var   = parent_strain, facet_colors = c("#5F05AF", "#2FD612", "#03C0FF"),
  facet_labels = c("JE2" = "parent", "dat" = "<i>dat</i>::Tn", "alr1" = "Δ<i> alr1</i>")
)
```

```{r, fig.width=10, fig.height=6, fig.align="center"}
# Add significance annotations
annots <- stats_formated_Chem$main %>%
  select(condition2, sign) %>%
  filter(sign != "") %>%  
  left_join(
    stats_formated_Chem$predictions %>%
      select(condition, predicted_upper, plot_label, parent_strain), 
    by = c("condition2" = "condition")
  ) %>%
  mutate(y_pos = predicted_upper * 5) %>%
  arrange(sign) %>%
  group_by(condition2, plot_label, parent_strain, y_pos) %>%
  summarise(sign = paste(sign, collapse = " "), .groups = "drop")

plot_HNO_CFUs_Chem_annot <- plot_HNO_CFUs_Chem +
  geom_text(
    data = annots, 
    aes(x = plot_label, y = y_pos, label = sign, group = NULL),
    color = "black",
    size = 4,
    inherit.aes = FALSE
  ) 

plot_HNO_CFUs_Chem_annot

# Save the plot
ggsave(filename = file.path(outputs_folder, "HNO_CFUs_Chemical.png"), plot = plot_HNO_CFUs_Chem_annot, width = 10, height = 6, dpi = 300)
saveRDS(plot_HNO_CFUs_Chem_annot, file.path(outputs_folder, "HNO_CFUs_Chemical.rds"))
```