CFUs Competition

# loading packages and functions
library(tidyverse)
library(kableExtra)
library(ggtext)
library(scales)
library(lme4)
library(lmerTest)
library(emmeans)
library(ggh4x)
library(broom)
library(ggpubr)

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 files
dataframeID <- "DatCol"
#Monocolonization CFU controls associated with competition experiments
CFUs_HNO_Controls <- readRDS(file.path(dataframes_folder, paste0(dataframeID, "_CFUs_HNO_Controls.rds")))

#Cocolonization CFU data
CFUs_HNO_Competition <- readRDS(file.path(dataframes_folder, paste0(dataframeID, "_CFUs_HNO_Competition.rds")))

# Filter the data for JE2 vs other parent strains
CFUs_HNO_JE2 <- CFUs_HNO_Competition %>%
  filter(str_detect(culture_pair, "_vs_JE2$"))

CFUs_HNO_Other <- CFUs_HNO_Competition %>%
  filter(!str_detect(culture_pair, "_vs_JE2$"))

CFUs

# Function to create a ggplot with jittered points, faceted by a variable, and customized aesthetics for fill and shape.
plot_comp_function <- function(data, ymin_exp, ymax_exp, title_x = " ", title_y = " ",
                          x_var, y_var, x_labels, fill_var, shape_var, facet_var, facet_labels,
                          fill_colors, shape_leves, facet_colors, fill_title, shape_title) {
  plot <- ggplot(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") +
    
    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) +
    scale_x_discrete(labels = x_labels) +
    
    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 = "black"),
      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)
}

Plot CFUs for JE2

plot_CFUs_Comp_JE2 <- plot_comp_function(
  data        = CFUs_HNO_JE2,
  x_var       = culture_pair, title_x = " ",
  y_var       = CFU_HNO_byStrain, title_y = "CFUs/HNO-ALI", ymin_exp = NA, ymax_exp = 9,
  fill_var    = strain, fill_colors = c("#5F05AF", "#2FD612", "#03C0FF", "#FC8EAC"), fill_title = "Strain",
  shape_var   = line, shape_leves = c(21, 22), shape_title = "HNO line",
  facet_var   = time, facet_colors = c("grey"), facet_labels = c("0" = "inoculum", "24" = "24 h on HNO-ALI"),
  x_labels = c(
    "dat_vs_JE2" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>",
    "alr1::Tn_vs_JE2" = "<b><span style='color:#03C0FF'><i>alr1</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>",
    "fmta::Tn_vs_JE2" = "<b><span style='color:#FC8EAC'><i>fmtA</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>"
  )
)

plot_CFUs_Comp_JE2

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

# Calculate mean CFUs for each strain and culture pair at 24h, then create a data frame with all pairwise comparisons of these means, including fold change and log2 fold change.
means_df_JE2 <- CFUs_HNO_JE2 %>%
  filter(time != 0) %>%
  group_by(strain, culture_pair) %>%
  summarise(mean_CFU = mean(CFU_HNO_byStrain, na.rm = TRUE), .groups = "drop")

pairs_df_JE2 <- crossing(
  means_df_JE2 %>% rename(strain1 = strain, culture_pair1 = culture_pair, mean1 = mean_CFU),
  means_df_JE2 %>% rename(strain2 = strain, culture_pair2 = culture_pair, mean2 = mean_CFU)
) 

pairs_df_JE2 <- pairs_df_JE2 %>%
  filter(!(strain1 == strain2 & culture_pair1 == culture_pair2)) %>%
  filter(
    # Keep all comparisons where JE2 is NOT involved
    !(strain1 == "JE2" | strain2 == "JE2") |
      
      # If JE2 IS involved, only keep comparisons within the same culture_pair
      (strain1 == "JE2" | strain2 == "JE2") & (culture_pair1 == culture_pair2)
  )

pairs_df_JE2 <- pairs_df_JE2 %>%
  mutate(
    FC = mean1 / mean2,
    log2FC = log2(FC)
  )

kable(
  pairs_df_JE2 %>%
    mutate(across(where(is.numeric), ~ format(.x, scientific = FALSE)))
)

strain1	culture_pair1	mean1	strain2	culture_pair2	mean2	FC	log2FC
JE2	dat_vs_JE2	64582000	dat	dat_vs_JE2	43000	1501.9069767442	10.552580
JE2	alr1::Tn_vs_JE2	52825000	alr1::Tn	alr1::Tn_vs_JE2	1550000	34.0806451613	5.090881
JE2	fmta::Tn_vs_JE2	51693750	fmta	fmta::Tn_vs_JE2	10181250	5.0773480663	2.344075
dat	dat_vs_JE2	43000	JE2	dat_vs_JE2	64582000	0.0006658202	-10.552580
dat	dat_vs_JE2	43000	alr1::Tn	alr1::Tn_vs_JE2	1550000	0.0277419355	-5.171788
dat	dat_vs_JE2	43000	fmta	fmta::Tn_vs_JE2	10181250	0.0042234500	-7.887362
alr1::Tn	alr1::Tn_vs_JE2	1550000	JE2	alr1::Tn_vs_JE2	52825000	0.0293421675	-5.090881
alr1::Tn	alr1::Tn_vs_JE2	1550000	dat	dat_vs_JE2	43000	36.0465116279	5.171788
alr1::Tn	alr1::Tn_vs_JE2	1550000	fmta	fmta::Tn_vs_JE2	10181250	0.1522406384	-2.715575
fmta	fmta::Tn_vs_JE2	10181250	JE2	fmta::Tn_vs_JE2	51693750	0.1969532100	-2.344075
fmta	fmta::Tn_vs_JE2	10181250	dat	dat_vs_JE2	43000	236.7732558140	7.887362
fmta	fmta::Tn_vs_JE2	10181250	alr1::Tn	alr1::Tn_vs_JE2	1550000	6.5685483871	2.715575

Plot CFUs for other parent strains

plot_CFUs_Comp_Other <- plot_comp_function(
  data        = CFUs_HNO_Other,
  x_var       = culture_pair, title_x = " ",
  y_var       = CFU_HNO_byStrain, title_y = "CFUs/HNO-ALI", ymin_exp = NA, ymax_exp = 9,
  fill_var    = strain, fill_colors = c("#5F05AF", "#2FD612", "#03C0FF","#5F05AF", "#2FD612", "#5F05AF", "#2FD612"), fill_title = "Strain",
  shape_var   = line, shape_leves = c(21, 22), shape_title = "HNO line",
  facet_var   = time, facet_colors = c("grey"), facet_labels = c("0" = "inoculum", "16" = "16 h on HNO-ALI", "24" = "24 h"),
  x_labels = c(
    "dat::Tn_vs_502a" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>502A</span> (CC8)</b>",
    "alr1::Tn_vs_502a" = "<b><span style='color:#03C0FF'><i>alr1</i>::Tn</span> vs. <span style='color:#5F05AF'>502A</span> (CC8)</b>",
    "dat::Tn_vs_KPL2715" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>KPL2715</span> (CC30)</b>",
    "dat::Tn_vs_KPL4559" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>AH2623</span> (CC5)</b>"
  )
)

plot_CFUs_Comp_Other

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

# Calculate mean CFUs for each strain and culture pair at 16/24h, then create a data frame with all pairwise comparisons of these means, including fold change and log2 fold change.
means_df_Other <- CFUs_HNO_Other %>%
  filter(time != 0) %>%
  group_by(strain, culture_pair) %>%
  summarise(mean_CFU = mean(CFU_HNO_byStrain, na.rm = TRUE), .groups = "drop")

pairs_df_Other <- crossing(
  means_df_Other %>% rename(strain1 = strain, culture_pair1 = culture_pair, mean1 = mean_CFU),
  means_df_Other %>% rename(strain2 = strain, culture_pair2 = culture_pair, mean2 = mean_CFU)
)

pairs_df_Other <- pairs_df_Other %>%
  mutate(
    fold_change = mean1 / mean2,
    log2FC = log2(fold_change)
  )

kable(
  pairs_df_Other %>%
    mutate(across(where(is.numeric), ~ format(.x, scientific = FALSE)))
)

strain1	culture_pair1	mean1	strain2	culture_pair2	mean2	fold_change	log2FC
502a	dat::Tn_vs_502a	73995750	502a	dat::Tn_vs_502a	73995750	1.00000000000	0.0000000
502a	dat::Tn_vs_502a	73995750	502a	alr1::Tn_vs_502a	65380000	1.13177959621	0.1785930
502a	dat::Tn_vs_502a	73995750	502a_dat	dat::Tn_vs_502a	4250	17410.76470588235	14.0876919
502a	dat::Tn_vs_502a	73995750	502a_alr1	alr1::Tn_vs_502a	370000	199.98851351351	7.6437733
502a	dat::Tn_vs_502a	73995750	KPL4559	dat::Tn_vs_KPL4559	23228000	3.18562726020	1.6715775
502a	dat::Tn_vs_502a	73995750	KPL4559_dat	dat::Tn_vs_KPL4559	22000	3363.44318181818	11.7157232
502a	dat::Tn_vs_502a	73995750	KPL2715	dat::Tn_vs_KPL2715	25734425	2.87536053360	1.5237429
502a	dat::Tn_vs_502a	73995750	KPL2715_dat	dat::Tn_vs_KPL2715	15575	4750.93097913323	12.2139945
502a	alr1::Tn_vs_502a	65380000	502a	dat::Tn_vs_502a	73995750	0.88356425876	-0.1785930
502a	alr1::Tn_vs_502a	65380000	502a	alr1::Tn_vs_502a	65380000	1.00000000000	0.0000000
502a	alr1::Tn_vs_502a	65380000	502a_dat	dat::Tn_vs_502a	4250	15383.52941176471	13.9090989
502a	alr1::Tn_vs_502a	65380000	502a_alr1	alr1::Tn_vs_502a	370000	176.70270270270	7.4651803
502a	alr1::Tn_vs_502a	65380000	KPL4559	dat::Tn_vs_KPL4559	23228000	2.81470638884	1.4929844
502a	alr1::Tn_vs_502a	65380000	KPL4559_dat	dat::Tn_vs_KPL4559	22000	2971.81818181818	11.5371301
502a	alr1::Tn_vs_502a	65380000	KPL2715	dat::Tn_vs_KPL2715	25734425	2.54056579854	1.3451498
502a	alr1::Tn_vs_502a	65380000	KPL2715_dat	dat::Tn_vs_KPL2715	15575	4197.75280898876	12.0354015
502a_dat	dat::Tn_vs_502a	4250	502a	dat::Tn_vs_502a	73995750	0.00005743573	-14.0876919
502a_dat	dat::Tn_vs_502a	4250	502a	alr1::Tn_vs_502a	65380000	0.00006500459	-13.9090989
502a_dat	dat::Tn_vs_502a	4250	502a_dat	dat::Tn_vs_502a	4250	1.00000000000	0.0000000
502a_dat	dat::Tn_vs_502a	4250	502a_alr1	alr1::Tn_vs_502a	370000	0.01148648649	-6.4439186
502a_dat	dat::Tn_vs_502a	4250	KPL4559	dat::Tn_vs_KPL4559	23228000	0.00018296883	-12.4161145
502a_dat	dat::Tn_vs_502a	4250	KPL4559_dat	dat::Tn_vs_KPL4559	22000	0.19318181818	-2.3719688
502a_dat	dat::Tn_vs_502a	4250	KPL2715	dat::Tn_vs_KPL2715	25734425	0.00016514843	-12.5639491
502a_dat	dat::Tn_vs_502a	4250	KPL2715_dat	dat::Tn_vs_KPL2715	15575	0.27287319422	-1.8736974
502a_alr1	alr1::Tn_vs_502a	370000	502a	dat::Tn_vs_502a	73995750	0.00500028718	-7.6437733
502a_alr1	alr1::Tn_vs_502a	370000	502a	alr1::Tn_vs_502a	65380000	0.00565922300	-7.4651803
502a_alr1	alr1::Tn_vs_502a	370000	502a_dat	dat::Tn_vs_502a	4250	87.05882352941	6.4439186
502a_alr1	alr1::Tn_vs_502a	370000	502a_alr1	alr1::Tn_vs_502a	370000	1.00000000000	0.0000000
502a_alr1	alr1::Tn_vs_502a	370000	KPL4559	dat::Tn_vs_KPL4559	23228000	0.01592905115	-5.9721959
502a_alr1	alr1::Tn_vs_502a	370000	KPL4559_dat	dat::Tn_vs_KPL4559	22000	16.81818181818	4.0719498
502a_alr1	alr1::Tn_vs_502a	370000	KPL2715	dat::Tn_vs_KPL2715	25734425	0.01437762841	-6.1200305
502a_alr1	alr1::Tn_vs_502a	370000	KPL2715_dat	dat::Tn_vs_KPL2715	15575	23.75601926164	4.5702212
KPL4559	dat::Tn_vs_KPL4559	23228000	502a	dat::Tn_vs_502a	73995750	0.31390992050	-1.6715775
KPL4559	dat::Tn_vs_KPL4559	23228000	502a	alr1::Tn_vs_502a	65380000	0.35527684307	-1.4929844
KPL4559	dat::Tn_vs_KPL4559	23228000	502a_dat	dat::Tn_vs_502a	4250	5465.41176470588	12.4161145
KPL4559	dat::Tn_vs_KPL4559	23228000	502a_alr1	alr1::Tn_vs_502a	370000	62.77837837838	5.9721959
KPL4559	dat::Tn_vs_KPL4559	23228000	KPL4559	dat::Tn_vs_KPL4559	23228000	1.00000000000	0.0000000
KPL4559	dat::Tn_vs_KPL4559	23228000	KPL4559_dat	dat::Tn_vs_KPL4559	22000	1055.81818181818	10.0441457
KPL4559	dat::Tn_vs_KPL4559	23228000	KPL2715	dat::Tn_vs_KPL2715	25734425	0.90260419652	-0.1478346
KPL4559	dat::Tn_vs_KPL4559	23228000	KPL2715_dat	dat::Tn_vs_KPL2715	15575	1491.36436597111	10.5424171
KPL4559_dat	dat::Tn_vs_KPL4559	22000	502a	dat::Tn_vs_502a	73995750	0.00029731437	-11.7157232
KPL4559_dat	dat::Tn_vs_KPL4559	22000	502a	alr1::Tn_vs_502a	65380000	0.00033649434	-11.5371301
KPL4559_dat	dat::Tn_vs_KPL4559	22000	502a_dat	dat::Tn_vs_502a	4250	5.17647058824	2.3719688
KPL4559_dat	dat::Tn_vs_KPL4559	22000	502a_alr1	alr1::Tn_vs_502a	370000	0.05945945946	-4.0719498
KPL4559_dat	dat::Tn_vs_KPL4559	22000	KPL4559	dat::Tn_vs_KPL4559	23228000	0.00094713277	-10.0441457
KPL4559_dat	dat::Tn_vs_KPL4559	22000	KPL4559_dat	dat::Tn_vs_KPL4559	22000	1.00000000000	0.0000000
KPL4559_dat	dat::Tn_vs_KPL4559	22000	KPL2715	dat::Tn_vs_KPL2715	25734425	0.00085488601	-10.1919803
KPL4559_dat	dat::Tn_vs_KPL4559	22000	KPL2715_dat	dat::Tn_vs_KPL2715	15575	1.41252006421	0.4982714
KPL2715	dat::Tn_vs_KPL2715	25734425	502a	dat::Tn_vs_502a	73995750	0.34778247399	-1.5237429
KPL2715	dat::Tn_vs_KPL2715	25734425	502a	alr1::Tn_vs_502a	65380000	0.39361310798	-1.3451498
KPL2715	dat::Tn_vs_KPL2715	25734425	502a_dat	dat::Tn_vs_502a	4250	6055.15882352941	12.5639491
KPL2715	dat::Tn_vs_KPL2715	25734425	502a_alr1	alr1::Tn_vs_502a	370000	69.55250000000	6.1200305
KPL2715	dat::Tn_vs_KPL2715	25734425	KPL4559	dat::Tn_vs_KPL4559	23228000	1.10790532977	0.1478346
KPL2715	dat::Tn_vs_KPL2715	25734425	KPL4559_dat	dat::Tn_vs_KPL4559	22000	1169.74659090909	10.1919803
KPL2715	dat::Tn_vs_KPL2715	25734425	KPL2715	dat::Tn_vs_KPL2715	25734425	1.00000000000	0.0000000
KPL2715	dat::Tn_vs_KPL2715	25734425	KPL2715_dat	dat::Tn_vs_KPL2715	15575	1652.29052969502	10.6902517
KPL2715_dat	dat::Tn_vs_KPL2715	15575	502a	dat::Tn_vs_502a	73995750	0.00021048506	-12.2139945
KPL2715_dat	dat::Tn_vs_KPL2715	15575	502a	alr1::Tn_vs_502a	65380000	0.00023822270	-12.0354015
KPL2715_dat	dat::Tn_vs_KPL2715	15575	502a_dat	dat::Tn_vs_502a	4250	3.66470588235	1.8736974
KPL2715_dat	dat::Tn_vs_KPL2715	15575	502a_alr1	alr1::Tn_vs_502a	370000	0.04209459459	-4.5702212
KPL2715_dat	dat::Tn_vs_KPL2715	15575	KPL4559	dat::Tn_vs_KPL4559	23228000	0.00067052695	-10.5424171
KPL2715_dat	dat::Tn_vs_KPL2715	15575	KPL4559_dat	dat::Tn_vs_KPL4559	22000	0.70795454545	-0.4982714
KPL2715_dat	dat::Tn_vs_KPL2715	15575	KPL2715	dat::Tn_vs_KPL2715	25734425	0.00060522044	-10.6902517
KPL2715_dat	dat::Tn_vs_KPL2715	15575	KPL2715_dat	dat::Tn_vs_KPL2715	15575	1.00000000000	0.0000000

Competition Index (CI)

Competition indices (CI = (mutant output / parent output) / (mutant input / parent input)) were calculated for each culture pair, where output is the CFU/HNO-ALI at the endpoint and input is the CFU/HNO-ALI at the inoculum.

# Function to create a ggplot with jittered points, emmeans predictions with error bars, and a dashed line at y=1.
plot_CI_function <- function(formated_stats_results, ymin_exp, ymax_exp, title_x = " ", title_y = " ",
                             x_var, y_var, x_labels, fill_var, shape_var, facet_var,
                             fill_colors, shape_leves, fill_title, shape_title) {
  
  annot <- formated_stats_results$main %>%
    filter(sign == "*") %>%
    mutate(sign = "")

  plot <- ggplot(data = formated_stats_results$input_data) +
    geom_hline(yintercept = 1, linetype = "dashed", color = "#5F05AF", linewidth = 0.8) +
    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) +
    
    scale_fill_manual(values = fill_colors, name = fill_title) +
    scale_shape_manual(values = shape_leves, name = shape_title) +
    scale_x_discrete(labels = x_labels) +
    
    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(
      breaks = 10^(ymin_exp:ymax_exp),
      labels = trans_format("log10", math_format(10^.x))) +
    
    labs(x = title_x, y = title_y) +
    theme_bw() +
    theme(
      panel.grid    = element_blank(),
      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")
    )
  
  if (nrow(annot) > 0) {
    plot <- plot +
      stat_pvalue_manual(annot, label = "sign", y.position = ymax_exp, 
                         step.increase = 0.04, tip.length = 0.01, bracket.shorten = 0.2, vjust = 0.6, bracket.size = 0.5) 
  }
  
  return(plot)
}

CI for JE2

Log-transformed CI were compared using a LMM to determine statistical significance with culture_pair as fixed effect and experimental date as a random effect (line was singular). All pairwise contrasts were performed and adjusted for multiple comparisons using the Holm method. A sign annotation column was created 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.

# Calculate Competition Index (CI) for JE2 
CFUs_HNO_JE2_CI <- CFUs_HNO_JE2 %>%
  mutate(type = if_else(strain == "JE2", "parent", "mutant")) %>%
  select(date, line, CC,culture_pair, type, time, CFU_HNO_byStrain) %>%
  group_by(date, culture_pair) %>%
  
  # pivot within each date × culture_pair group
  pivot_wider(
    names_from = c(type, time),
    values_from = CFU_HNO_byStrain
  ) %>%
  
  # compute CI 
  mutate(CI = (mutant_24 / parent_24) / (mutant_0 / parent_0)) %>%
  ungroup() 

# Run the mixed model with emmeans contrasts for the JE2 dataset.
stats_results_JE2 <- mixed_model_with_emmeans(
  input_data = CFUs_HNO_JE2_CI, response = "CI", log_response = TRUE, adjust = "holm",
  fixed_effects = c("culture_pair"), random_effects1 = "line", random_effects2 = "date"
)
stats_formated_JE2 <- format_stats_results(stats_results_JE2)

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

Model Formula: log(CI) ~ culture_pair + (1 | date)

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

Random Effects

grp	var1	var2	vcov	sdcor	proportion	singular
line:date	(Intercept)	NA	0.0000000	0.0000112	0.0	TRUE
line	(Intercept)	NA	0.2478644	0.4978598	20.9	TRUE
Residual	NA	NA	0.9380644	0.9685373	79.1	TRUE

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

ANOVA Table

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
culture_pair	63.83041	31.91521	2	8.826361	34.89034	6.44e-05

kable(stats_formated_JE2$main %>%
    mutate(
      p.value.anova = formatC(p.value.anova, format = "e", digits = 2),
      p.value = formatC(p.value, format = "e", digits = 2)
    ) %>%
    select(effect_name, p.value.anova,  condition1, condition2, predicted_mean.1, predicted_mean.2, FC, log2FC, p.value, 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	sign
culture_pair	6.44e-05	dat_vs_JE2	alr1::Tn_vs_JE2	0.0007727	0.0264973	0.0291629	-5.099724	1.99e-03	*
culture_pair	6.44e-05	dat_vs_JE2	fmta::Tn_vs_JE2	0.0007727	0.1200855	0.0064349	-7.279866	1.44e-05	*
culture_pair	6.44e-05	alr1::Tn_vs_JE2	fmta::Tn_vs_JE2	0.0264973	0.1200855	0.2206540	-2.180142	4.24e-02	*

plot_CI_JE2 <- plot_CI_function(
  formated_stats_results = stats_formated_JE2,
  x_var       = culture_pair, title_x = " ",
  y_var       = CI, title_y = "Competition Index (CI)", ymin_exp = -4, ymax_exp = 0.1,
  fill_var    = culture_pair, fill_colors = c("#2FD612", "#03C0FF", "#FC8EAC"), fill_title = "Strain",
  shape_var   = line, shape_leves = c(21, 22), shape_title = "HNO line",
  x_labels = c(
    "dat_vs_JE2" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>",
    "alr1::Tn_vs_JE2" = "<b><span style='color:#03C0FF'><i>alr1</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>",
    "fmta::Tn_vs_JE2" = "<b><span style='color:#FC8EAC'><i>fmtA</i>::Tn</span> vs. <span style='color:#5F05AF'>parent</span></b>"
  )
)

plot_CI_JE2

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

A one-sample t-test was used to determine whether the CI for each mutant was significantly lower than 1, where 1 indicates the same fitness as the parental strain (purple dashed line). Holm adjustment was used to adjust p-values for multiple tests.

# Run one sample T test comparing to 1 for each culture pair
CI.tTest_JE2 <- CFUs_HNO_JE2_CI %>%
  group_by(culture_pair) %>%
  summarise(
    test = list(tidy(t.test(log(CI), mu = 1)))
  ) %>%
  unnest(test) %>%
  mutate(p.adj = p.adjust(p.value, method = "holm"),
         sign = case_when(p.adj < 0.05 ~ "*", TRUE ~ ""))

kable(CI.tTest_JE2 %>%
        mutate(p.value = formatC(p.value, format = "e", digits = 2),
               p.adj = formatC(p.adj, format = "e", digits = 2))
)

culture_pair	estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative	p.adj	sign
dat_vs_JE2	-7.165570	-15.023814	2.37e-05	5	-8.562703	-5.7684373	One Sample t-test	two.sided	7.10e-05	*
alr1::Tn_vs_JE2	-3.630711	-13.983704	3.36e-05	5	-4.481960	-2.7794614	One Sample t-test	two.sided	7.10e-05	*
fmta::Tn_vs_JE2	-2.119552	-7.835838	4.33e-03	3	-3.386526	-0.8525773	One Sample t-test	two.sided	4.33e-03	*

# Save the stats results
stats_formated_JE2$tTest <- CI.tTest_JE2
saveRDS(stats_formated_JE2, file.path(outputs_folder, "stats_results_HNO_CI_JE2.rds"))

CI for other parent strains

Log-transformed CI were compared using a LMM to determine statistical significance with culture_pair as fixed effect and experimental date and line as random effects (line was singular). We analyzed pairwise contrasts comparing only between mutants from the same parent 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.

# Calculate Competition Index (CI) for other parent strains
CFUs_HNO_Other_CI <- CFUs_HNO_Other %>%
  mutate(type = if_else(!str_detect(strain, "_"), "parent", "mutant")) %>%
  select(date, line, culture_pair, type, CC, time, CFU_HNO_byStrain) %>%
  group_by(date, culture_pair) %>%
  
  # pivot within each date × culture_pair group
  pivot_wider(
    names_from = c(type, time),
    values_from = CFU_HNO_byStrain
  ) %>%
  # compute CI 
  mutate(
    CI = coalesce(
      (mutant_24 / parent_24) / (mutant_0 / parent_0),
      (mutant_16 / parent_16) / (mutant_0 / parent_0)
    )
  ) %>%
  ungroup()

# Run the mixed model with emmeans contrasts for the JE2 dataset.
stats_results_Other <- mixed_model_with_emmeans(
  input_data = CFUs_HNO_Other_CI, response = "CI", log_response = TRUE, adjust = "none",
  fixed_effects = c("culture_pair"), random_effects1 = "line", random_effects2 = "date"
)
stats_formated_Other <- format_stats_results(stats_results_Other)

# Filter the contrasts to only include those comparing between the same parent strain.
# Pvalues are adjusted using the Holm method only for the contrasts that are being kept and the significance annotation is added accordingly.
stats_formated_Other$main <- stats_formated_Other$main %>%
  filter(sub(".*_vs_", "", condition1) == sub(".*_vs_", "", condition2)) %>%
  mutate(
    p.value.adj = p.adjust(p.value, method = "holm"),
    sign = case_when(p.value.adj < 0.05 ~ "*", TRUE ~ "")
  )

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

Model Formula: log(CI) ~ culture_pair + (1 | date)

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

Random Effects

grp	var1	var2	vcov	sdcor	proportion	singular
line:date	(Intercept)	NA	0.0000000	0.0000000	0.00	TRUE
line	(Intercept)	NA	0.0009540	0.0308869	0.11	TRUE
Residual	NA	NA	0.8582918	0.9264404	99.89	TRUE

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

ANOVA Table

	Sum Sq	Mean Sq	NumDF	DenDF	F value	Pr(>F)
culture_pair	27.58891	9.196304	3	8	10.71466	0.0035578

kable(stats_formated_Other$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
culture_pair	3.56e-03	dat::Tn_vs_502a	alr1::Tn_vs_502a	7.14e-05	0.0042141	0.0169428	-5.883184	8.31e-03	8.31e-03	*

plot_CI_Other <- plot_CI_function(
  formated_stats_results = stats_formated_Other,
  x_var       = culture_pair, title_x = " ",
  y_var       = CI, title_y = "Competition Index (CI)", ymin_exp = -5, ymax_exp = 0.1,
  fill_var    = culture_pair, fill_colors = c("#2FD612", "#03C0FF", "#2FD612", "#2FD612"), fill_title = "Strain",
  shape_var   = line, shape_leves = c(21, 22), shape_title = "HNO line",
  x_labels = c(
    "dat::Tn_vs_502a" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>502A</span></b>",
    "alr1::Tn_vs_502a" = "<b><span style='color:#03C0FF'><i>alr1</i>::Tn</span> vs. <span style='color:#5F05AF'>502A</span></b>",
    "dat::Tn_vs_KPL2715" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>KPL2715</span></b>",
    "dat::Tn_vs_KPL4559" = "<b><span style='color:#2FD612'><i>dat</i>::Tn</span> vs. <span style='color:#5F05AF'>AH2623</span></b>"
  )
)

plot_CI_Other <- plot_CI_Other +
  geom_vline(xintercept = c(2.5, 3.5), color = "grey70", linewidth = 0.5) +
  annotate("text", x = 1.5, y = Inf, label = "CC8", vjust = -0.5, fontface = "bold") +
  annotate("text", x = 3.0, y = Inf, label = "CC5", vjust = -0.5, fontface = "bold") +
  annotate("text", x = 4, y = Inf, label = "CC30", vjust = -0.5, fontface = "bold") +
  coord_cartesian(clip = "off") +   
  theme(plot.margin = margin(t = 30))

plot_CI_Other

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

A one-sample t-test was used to determine whether the CI for each mutant was significantly lower than 1, where 1 indicates the same fitness as the parental strain (purple dashed line). Holm adjustment was used to adjust p-values for multiple tests.

# Run one sample T test comparing to 1 for each culture pair
CI.tTest_Other <- CFUs_HNO_Other_CI %>%
  group_by(culture_pair) %>%
  summarise(
    test = list(tidy(t.test(log(CI), mu = 1)))
  ) %>%
  unnest(test) %>%
  mutate(p.adj = p.adjust(p.value, method = "holm"),
         sign = case_when(p.adj < 0.05 ~ "*", TRUE ~ ""))

kable(CI.tTest_Other %>%
        mutate(p.value = formatC(p.value, format = "e", digits = 2),
               p.adj = formatC(p.adj, format = "e", digits = 2))
)

culture_pair	estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative	p.adj	sign
dat::Tn_vs_502a	-9.547233	-31.32453	1.02e-03	2	-10.995973	-8.098493	One Sample t-test	two.sided	4.07e-03	*
alr1::Tn_vs_502a	-5.469321	-21.68466	2.12e-03	2	-6.752959	-4.185683	One Sample t-test	two.sided	6.36e-03	*
dat::Tn_vs_KPL4559	-6.750516	-17.81165	3.14e-03	2	-8.622762	-4.878270	One Sample t-test	two.sided	6.36e-03	*
dat::Tn_vs_KPL2715	-8.071223	-10.45597	9.02e-03	2	-11.804049	-4.338396	One Sample t-test	two.sided	9.02e-03	*

# Save the stats results
stats_formated_Other$tTest <- CI.tTest_Other
saveRDS(stats_formated_Other, file.path(outputs_folder, "stats_results_HNO_CI_Other.rds"))