CFUs Epithelial Lines

library(tidyverse)
library(ggpubr)
library(ggtext)
library(scales)
library(lme4)
library(afex)
library(emmeans)
library(gridExtra)

Data Input and Selection

File Paths

# Folder paths
input_path <- "data/input_data/CFUs/"
metadata_path <- "data/metadata/CFUs"

# Create subfolders for output files
dataframes_folder <- "data/dataframes"
if (!file.exists("data/dataframes")) {
  dir.create("data/dataframes", recursive = TRUE)
}

# Load data and metadata
input_data <- read_csv(file.path(input_path, "HNObac_6h.csv")) %>% mutate_if(is.character, factor)
input_data$Time <- as.factor(input_data$Time)
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_CFUs_6h.csv"))

Data clean-up

# Setting zero values to the limit of detection
CFU_data <- input_data %>%
  mutate(
    LOD = CFUs == 0,
    NewCFU = ifelse(LOD, 3.75, CFUs)
  )

# Factor Ordering and Styling
CFU_data <- CFU_data %>%
  filter(Time != 24) 
CFU_data <- merge(CFU_data, Bac_order, by = "bacteria") 
CFU_data$bacteria_label <- factor(CFU_data$bacteria_label, levels = Bac_order$bacteria_label)
CFU_data$Line <- fct_recode(CFU_data$Line, "HNO918" = "A", "HNO204" = "B", "HNO919" = "C") 
CFU_data$Model <- fct_recode(CFU_data$Model, "HNO" = "HNO", "Calu-3" = "Calu-3", "RPMI" = "RPMI 2650") %>%
  fct_relevel("HNO", "C", "R")

# Calculate log2 fold change between final time point vs. time 0
CFU_FC <- CFU_data %>%
  group_by(Date, Model, Line, bacteria, bacteria_label) %>%
  reframe(FC = log2(NewCFU[Time == 6]/NewCFU[Time == 0])) %>%
  mutate(Species.Model = interaction(bacteria, Model, lex.order = T, sep = ".")) %>%
  mutate(Species.Model_label = interaction(bacteria_label, Model, lex.order = T, sep = "<br>"))

Saving files

# Save data frame as CSV files in the output folder
write_csv(CFU_data, file.path(dataframes_folder, "CFU_values_6h.csv"))
write_csv(CFU_FC, file.path(dataframes_folder, "CFU_FC_6h.csv"))

# Save data frame as R objects in the output folder
saveRDS(CFU_data, file.path(dataframes_folder, "CFU_values_6h.rds"))
saveRDS(CFU_FC, file.path(dataframes_folder, "CFU_FC_6h.rds"))

# Cleaning-up all objects from the environment
rm(list = ls())

# Use this to read the final objects
CFU_data <- readRDS("data/dataframes/CFU_values_6h.rds")
CFU_FC <- readRDS("data/dataframes/CFU_FC_6h.rds")

Stats and Plots

File Paths

# Folder paths
dataframes_path <- "data/dataframes"
metadata_path <- "data/metadata/CFUs"

# Create subfolders for output files
figures_folder <- "data/outputs/CFUs/figures"
if (!file.exists("data/outputs/CFUs/figures")) {
  dir.create("data/outputs/CFUs/figures", recursive = TRUE)
}
stats_folder <- "data/outputs/CFUs/stats"
if (!file.exists("data/outputs/CFUs/stats")) {
  dir.create("data/outputs/CFUs/stats", recursive = TRUE)
}

# Load data and metadata
CFU_FC <- readRDS("data/dataframes/CFU_FC_6h.rds")

# Function to compare Epithelial Lines
analysis_celllines_function <- function(data, cutoff_pvalue, cutoff_FC) {
  
  # Subset the data if needed
  data_subset <- data
  
  # Mixed-effects model with random effects
  model <- lmer(FC ~ bacteria * Model
                + (1|Line) + (1|Line:Date), 
                data = data_subset)
  #Anova
  anova <- anova(model)
  anova_df <- as.data.frame(anova) %>%
    mutate(sign = case_when(
      `Pr(>F)` < cutoff_pvalue ~ "*",
      TRUE ~ "")) %>%
    mutate_if(is.numeric, ~ format(., digits = 2, scientific = TRUE))
  
  # Calculate Individual contrasts
  emmeans_model <- emmeans(model, ~ bacteria * Model)
  emmeans_Model <- pairs(emmeans_model, simple = "Model", adjust = "none")   
  
  # Extract random effects and convert to dataframe (if not singular)
  random_effects_df <- as.data.frame(VarCorr(model)) %>%
    mutate(proportion = round(100 * (vcov / sum(vcov)), 2)) 
  
  # Adds predictions based on fixed effects, averaged over random effects. It gives a population estimate
  data_subset <- cbind(data_subset, predval = predict(model,re.form = NA, se.fit = TRUE))
  data_summary_df <- data_subset %>%
    group_by(Species.Model_label, Species.Model, Model, bacteria) %>%
    summarize(mean.real = mean(FC),
              mean.predval = mean(predval.fit), 
              mean.predval.se = mean(predval.se.fit)) %>%
    mutate(max = mean.predval + 2*mean.predval.se,
           min = mean.predval - 2*mean.predval.se)
  
  # Convert contrasts to dataframe and adjust pvalues.
  contrasts_df <- as.data.frame(summary(emmeans_Model)) %>%
    mutate(p.adj.holm = p.adjust(p.value, method = "holm")) %>%
    mutate(sign = case_when(
      p.adj.holm < 0.05 ~ "*",
      TRUE ~ ""))
  
  # Edits to the contrast dataframe to include pvalue brackets in plot
  contrasts_df <- contrasts_df %>%
    separate(contrast, into = c("Model1", "Model2"), sep = " - ") %>%
    mutate(Model1 = sub("Model", "", Model1),
           Model2 = sub("Model", "", Model2)) %>%
    mutate(condition1 = paste(bacteria, Model1, sep = "."),
           condition2 = paste(bacteria, Model2, sep = "."))
  
  # Calculate fold-change values for each contrast
  contrasts_df <- contrasts_df %>%
    ungroup() %>%
    left_join(select(data_summary_df, Species.Model, mean.predval), by = join_by(condition1 == Species.Model)) %>%
    left_join(select(data_summary_df, Species.Model, mean.predval), by = join_by(condition2 == Species.Model), suffix = c(".1", ".2")) %>%
    mutate(FC = mean.predval.1 / mean.predval.2,
           FC = if_else(FC < 1, -1 / FC, FC),
           highlighted = case_when(
             FC <= -cutoff_FC ~ "+",
             FC >= cutoff_FC ~ "-",
             TRUE ~ "")) 
  
  # Plot with predicted means and standard errors of the estimates
  plot_2 <- ggplot() +
    geom_point(data = data_subset, 
               aes(x = Species.Model_label, y = FC, fill = bacteria_label, color = bacteria_label, shape = Line), 
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               size = 3, alpha = 0.75, stroke = 0.75) +
    
    scale_fill_manual(values = c("#800080","#1E90FF","#927ed1")) +
    scale_color_manual(values = c("#800080","#1E90FF","#927ed1")) +
    scale_shape_manual(values = c(21, 22, 24, 23, 25)) +
    
    geom_point(data = data_summary_df, aes(x = Species.Model_label, 
                                           y = mean.predval), shape = 3, size = 3) +
    geom_errorbar(data = data_summary_df, aes(x = Species.Model_label,
                                              y = mean.predval,
                                              ymax = max,
                                              ymin = min),
                  width = 0.4) +

    scale_y_continuous(expand = c(0.1,0)) +
    
    labs(x = "",
         y = "log2 (6h CFUs/inoculum CFUs)",
         fill = "Bacteria", color = "Bacteria", shape = "HNO/Cell Line") +
    
    theme_bw() +
    theme(panel.grid = element_blank(), 
          legend.text = element_markdown(),
          text = element_text(size = 20), 
          axis.text.x = element_markdown(), 
          axis.text.y = element_text(color = "black"))
  
  # Save files
  ggsave(plot_2, filename = paste0(figures_folder, "/plotCFU_", "HNOvsCells", ".png"), width = 7, height = 6)
  saveRDS(plot_2, file.path(figures_folder, paste0("plotCFU_", "HNOvsCells", ".rds")))
  write_csv(anova_df, file.path(stats_folder, paste0("anova_", "HNOvsCells", ".csv")))
  write_csv(random_effects_df, file.path(stats_folder, paste0("stats_random_effects", "HNOvsCells", ".csv")))
  write_csv(contrasts_df, file.path(stats_folder, paste0("stats_contrasts", "HNOvsCells", ".csv")))
  write_csv(data_summary_df, file.path(stats_folder, paste0("stats_summary", "HNOvsCells", ".csv")))
  
  return(list(
    anova = anova_df,
    random_effects = random_effects_df,
    contrasts = contrasts_df,
    data_summary = data_summary_df,
    plot_2 = plot_2
  ))
}

analysis_celllines_function(CFU_FC, cutoff_pvalue = 0.05, cutoff_FC = 1)

$anova
                Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F) sign
bacteria       5.1e+01 2.6e+01     2 1.8e+01 1.6e+01 1.1e-04    *
Model          2.1e+01 1.1e+01     2 9.0e+00 6.5e+00 1.8e-02    *
bacteria:Model 3.5e+01 8.9e+00     4 1.8e+01 5.5e+00 4.6e-03    *

$random_effects
        grp        var1 var2         vcov        sdcor proportion
1 Line:Date (Intercept) <NA> 8.239120e-01 9.076960e-01       33.7
2      Line (Intercept) <NA> 1.958177e-15 4.425129e-08        0.0
3  Residual        <NA> <NA> 1.620894e+00 1.273143e+00       66.3

$contrasts
    Model1   Model2 bacteria   estimate       SE       df    t.ratio    p.value
1      HNO (Calu-3)      Dpi -1.7689085 1.187922 1.224298 -1.4890777 0.34289499
2      HNO     RPMI      Dpi -0.6034053 1.270784 1.745515 -0.4748290 0.68759918
3 (Calu-3)     RPMI      Dpi  1.1655033 1.194210 1.059678  0.9759620 0.50031665
4      HNO (Calu-3)      Sau -1.0677213 1.187922 1.224298 -0.8988142 0.51055683
5      HNO     RPMI      Sau -1.6321132 1.270784 1.745515 -1.2843352 0.34326190
6 (Calu-3)     RPMI      Sau -0.5643918 1.194210 1.059678 -0.4726070 0.71536152
7      HNO (Calu-3)      Spn -3.4261805 1.187922 1.224298 -2.8841791 0.17417683
8      HNO     RPMI      Spn -6.2490460 1.270784 1.745515 -4.9174714 0.05077755
9 (Calu-3)     RPMI      Spn -2.8228655 1.194210 1.059678 -2.3637939 0.24338444
  p.adj.holm sign   condition1   condition2
1  1.0000000           Dpi.HNO Dpi.(Calu-3)
2  1.0000000           Dpi.HNO     Dpi.RPMI
3  1.0000000      Dpi.(Calu-3)     Dpi.RPMI
4  1.0000000           Sau.HNO Sau.(Calu-3)
5  1.0000000           Sau.HNO     Sau.RPMI
6  1.0000000      Sau.(Calu-3)     Sau.RPMI
7  1.0000000           Spn.HNO Spn.(Calu-3)
8  0.4569979           Spn.HNO     Spn.RPMI
9  1.0000000      Spn.(Calu-3)     Spn.RPMI
                                         Species.Model_label.1 Model.1
1 <b><i><span style='color:#1E90FF;'>Dpi</span></i></b><br>HNO     HNO
2 <b><i><span style='color:#1E90FF;'>Dpi</span></i></b><br>HNO     HNO
3                                                         <NA>    <NA>
4 <b><i><span style='color:#800080;'>Sau</span></i></b><br>HNO     HNO
5 <b><i><span style='color:#800080;'>Sau</span></i></b><br>HNO     HNO
6                                                         <NA>    <NA>
7 <b><i><span style='color:#927ED1;'>Spn</span></i></b><br>HNO     HNO
8 <b><i><span style='color:#927ED1;'>Spn</span></i></b><br>HNO     HNO
9                                                         <NA>    <NA>
  mean.predval.1                                         Species.Model_label.2
1     -1.6156653                                                          <NA>
2     -1.6156653 <b><i><span style='color:#1E90FF;'>Dpi</span></i></b><br>RPMI
3             NA <b><i><span style='color:#1E90FF;'>Dpi</span></i></b><br>RPMI
4      0.8918372                                                          <NA>
5      0.8918372 <b><i><span style='color:#800080;'>Sau</span></i></b><br>RPMI
6             NA <b><i><span style='color:#800080;'>Sau</span></i></b><br>RPMI
7     -3.9852415                                                          <NA>
8     -3.9852415 <b><i><span style='color:#927ED1;'>Spn</span></i></b><br>RPMI
9             NA <b><i><span style='color:#927ED1;'>Spn</span></i></b><br>RPMI
  Model.2 mean.predval.2        FC highlighted
1    <NA>             NA        NA            
2    RPMI      -1.012260  1.596097           -
3    RPMI      -1.012260        NA            
4    <NA>             NA        NA            
5    RPMI       2.523950 -2.830057           +
6    RPMI       2.523950        NA            
7    <NA>             NA        NA            
8    RPMI       2.263804  0.568047            
9    RPMI       2.263804        NA            

$data_summary
# A tibble: 9 × 9
# Groups:   Species.Model_label, Species.Model, Model [9]
  Species.Model_label        Species.Model Model bacteria mean.real mean.predval
  <fct>                      <fct>         <fct> <fct>        <dbl>        <dbl>
1 <b><i><span style='color:… Sau.HNO       HNO   Sau          0.892        0.892
2 <b><i><span style='color:… Sau.Calu-3    Calu… Sau          1.96         1.96 
3 <b><i><span style='color:… Sau.RPMI      RPMI  Sau          2.52         2.52 
4 <b><i><span style='color:… Dpi.HNO       HNO   Dpi         -1.62        -1.62 
5 <b><i><span style='color:… Dpi.Calu-3    Calu… Dpi          0.153        0.153
6 <b><i><span style='color:… Dpi.RPMI      RPMI  Dpi         -1.01        -1.01 
7 <b><i><span style='color:… Spn.HNO       HNO   Spn         -3.99        -3.99 
8 <b><i><span style='color:… Spn.Calu-3    Calu… Spn         -0.559       -0.559
9 <b><i><span style='color:… Spn.RPMI      RPMI  Spn          2.26         2.26 
# ℹ 3 more variables: mean.predval.se <dbl>, max <dbl>, min <dbl>

$plot_2