CFUs Complementation

# 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.

# 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)
}

# 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() 

# 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.

# 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"))

# Render the results
cat("Model Formula:", format(stats_results_Gen$model_formula))

Model Formula: log(CFU_HNO) ~ vector * parent_strain + (1 | date)

kable(stats_formated_Gen$random_effects, caption = "Random Effects")

Random Effects

grp	var1	var2	vcov	sdcor	proportion	singular
line:date	(Intercept)	NA	0.1647199	0.4058570	16.11	TRUE
line	(Intercept)	NA	0.0000000	0.0000000	0.00	TRUE
Residual	NA	NA	0.8576640	0.9261015	83.89	TRUE

kable(stats_formated_Gen$anova, caption = "ANOVA Table")

ANOVA Table

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
vector	177.5813	59.19378	3	45	69.01745	0
parent_strain	281.3469	140.67343	2	45	164.01928	0
vector:parent_strain	209.4649	52.36622	4	45	61.05680	0

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"
)

Main Contrasts with Predicted Log2FC

effect_name	p.value.anova	condition1	condition2	predicted_mean.1	predicted_mean.2	FC	log2FC	p.value	p.value.adj	sign
vector	7.22e-17	pLL39.JE2	pAB004.JE2	92586244.85	52809854.40	1.7532002	0.8099907	2.99e-01	1.00e+00
vector	7.22e-17	pLL39.JE2	pAB009.JE2	92586244.85	61468267.35	1.5062446	0.5909561	4.48e-01	1.00e+00
vector	7.22e-17	pLL39.JE2	pAB001.JE2	92586244.85	54553082.86	1.6971771	0.7631371	3.28e-01	1.00e+00
vector	7.22e-17	pLL39.dat	pAB004.dat	13216.98	9269384.11	0.0014259	-9.4539371	6.14e-16	6.14e-15	b
vector	7.22e-17	pLL39.dat	pAB009.dat	13216.98	91898573.97	0.0001438	-12.7634342	9.41e-21	1.18e-19	b
vector	7.22e-17	pLL39.dat	pAB001.dat	13216.98	24455.36	0.5404534	-0.8877579	2.56e-01	1.00e+00
vector	7.22e-17	pLL39.alr1	pAB001.alr1	1709665.80	69078978.55	0.0247494	-5.3364605	1.35e-08	1.08e-07	b
parent_strain	2.15e-21	pLL39.JE2	pLL39.dat	92586244.85	13216.98	7005.0983266	12.7741896	9.10e-21	1.18e-19	a
parent_strain	2.15e-21	pLL39.JE2	pLL39.alr1	92586244.85	1709665.80	54.1545867	5.7590116	2.09e-09	1.88e-08	a
parent_strain	2.15e-21	pAB004.JE2	pAB004.dat	52809854.40	9269384.11	5.6972344	2.5102618	2.16e-03	1.51e-02	a
parent_strain	2.15e-21	pAB009.JE2	pAB009.dat	61468267.35	91898573.97	0.6688707	-0.5802007	4.56e-01	1.00e+00
parent_strain	2.15e-21	pAB001.JE2	pAB001.dat	54553082.86	24455.36	2230.7213208	11.1232946	1.75e-18	1.93e-17	a
parent_strain	2.15e-21	pAB001.JE2	pAB001.alr1	54553082.86	69078978.55	0.7897205	-0.3405860	6.61e-01	1.00e+00

# 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>")
)

# 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.

# 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"))

# 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"))

# Render the results
cat("Model Formula:", format(stats_results_Chem$model_formula))

Model Formula: log(CFU_HNO) ~ label + (1 | date)

kable(stats_formated_Chem$random_effects, caption = "Random Effects")

Random Effects

grp	var1	var2	vcov	sdcor	proportion	singular
line:date	(Intercept)	NA	0.1304513	0.3611804	10.25	TRUE
line	(Intercept)	NA	0.0000000	0.0000000	0.00	TRUE
Residual	NA	NA	1.1426579	1.0689518	89.75	TRUE

kable(stats_formated_Chem$anova, caption = "ANOVA Table")

ANOVA Table

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
label	1802.136	106.008	17	71.26213	92.7724	0

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"
)

Main Contrasts with Predicted Log2FC

effect_name	p.value.anova	contrast	condition1	condition2	FC	log2FC	p.value	p.value.adj	sign
label	1.28e-41	JE2-dat	JE2	dat	7.208314e+03	12.8154461	1.30e-26	3.38e-25	a
label	1.28e-41	JE2-alr1	JE2	alr1	1.940502e+01	4.2783582	2.22e-04	3.33e-03	a
label	1.28e-41	JE2-JE2_37.5mM_Dala	JE2	JE2_37.5mM_Dala	6.767993e-01	-0.5632001	5.04e-01	1.00e+00
label	1.28e-41	JE2-JE2_37.5mM_Lala	JE2	JE2_37.5mM_Lala	5.880842e-01	-0.7659053	3.64e-01	1.00e+00
label	1.28e-41	JE2-JE2_37.5mM_Dglu	JE2	JE2_37.5mM_Dglu	1.724813e+00	0.7864398	3.52e-01	1.00e+00
label	1.28e-41	JE2-JE2_37.5mM_Lglu	JE2	JE2_37.5mM_Lglu	6.492691e-01	-0.6231115	4.59e-01	1.00e+00
label	1.28e-41	JE2-dat_0.375mM_Dala	JE2	dat_0.375mM_Dala	2.305015e+03	11.1705604	1.80e-21	3.78e-20	a
label	1.28e-41	JE2-dat_3.75mM_Dala	JE2	dat_3.75mM_Dala	6.116778e+00	2.6127718	2.61e-03	3.65e-02	a
label	1.28e-41	JE2-dat_37.5mM_Dala	JE2	dat_37.5mM_Dala	5.161223e-01	-0.9542151	2.59e-01	1.00e+00
label	1.28e-41	JE2-dat_37.5mM_Lala	JE2	dat_37.5mM_Lala	4.558183e-01	-1.1334693	1.80e-01	1.00e+00
label	1.28e-41	JE2-dat_0.375mM_Dglu	JE2	dat_0.375mM_Dglu	1.267696e+04	13.6299214	2.59e-26	6.22e-25	a
label	1.28e-41	JE2-dat_3.75mM_Dglu	JE2	dat_3.75mM_Dglu	4.832199e+03	12.2384641	1.24e-23	2.86e-22	a
label	1.28e-41	JE2-dat_37.5mM_Dglu	JE2	dat_37.5mM_Dglu	3.794741e+04	15.2117137	3.52e-29	9.86e-28	a
label	1.28e-41	JE2-dat_37.5mM_Lglu	JE2	dat_37.5mM_Lglu	2.425389e+00	1.2782163	1.31e-01	1.00e+00
label	1.28e-41	JE2-alr1_3.75mM_Dala	JE2	alr1_3.75mM_Dala	4.187689e+01	5.3880825	5.67e-06	9.64e-05	a
label	1.28e-41	JE2-alr1_37.5mM_Dala	JE2	alr1_37.5mM_Dala	8.278199e+03	13.0151013	9.72e-19	1.85e-17	a
label	1.28e-41	JE2-alr1_37.5mM_Lala	JE2	alr1_37.5mM_Lala	6.162314e-01	-0.6984558	5.28e-01	1.00e+00
label	1.28e-41	dat-dat_0.375mM_Dala	dat	dat_0.375mM_Dala	3.197717e-01	-1.6448857	5.36e-02	6.44e-01
label	1.28e-41	dat-dat_3.75mM_Dala	dat	dat_3.75mM_Dala	8.486000e-04	-10.2026743	1.94e-19	3.89e-18	b
label	1.28e-41	dat-dat_37.5mM_Dala	dat	dat_37.5mM_Dala	7.160000e-05	-13.7696612	1.42e-26	3.55e-25	b
label	1.28e-41	dat-dat_37.5mM_Lala	dat	dat_37.5mM_Lala	6.320000e-05	-13.9489154	7.50e-27	2.03e-25	b
label	1.28e-41	dat-dat_0.375mM_Dglu	dat	dat_0.375mM_Dglu	1.758658e+00	0.8144753	3.35e-01	1.00e+00
label	1.28e-41	dat-dat_3.75mM_Dglu	dat	dat_3.75mM_Dglu	6.703647e-01	-0.5769820	4.94e-01	1.00e+00
label	1.28e-41	dat-dat_37.5mM_Dglu	dat	dat_37.5mM_Dglu	5.264395e+00	2.3962676	5.54e-03	7.20e-02
label	1.28e-41	dat-dat_37.5mM_Lglu	dat	dat_37.5mM_Lglu	3.365000e-04	-11.5372298	3.44e-22	7.57e-21	b
label	1.28e-41	alr1-alr1_3.75mM_Dala	alr1	alr1_3.75mM_Dala	2.158044e+00	1.1097243	3.81e-01	1.00e+00
label	1.28e-41	alr1-alr1_37.5mM_Dala	alr1	alr1_37.5mM_Dala	4.266009e+02	8.7367431	1.31e-09	2.36e-08	c
label	1.28e-41	alr1-alr1_37.5mM_Lala	alr1	alr1_37.5mM_Lala	3.175630e-02	-4.9768140	1.76e-04	2.82e-03	c

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>")
)

# 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"))