LDH Epithelial Lines

library(tidyverse)
library(readxl)
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/LDH/"
metadata_path <- "data/metadata/LDH"

# 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_excel(file.path(input_path, "HNObac_6h_LDH.xlsx")) %>% mutate_if(is.character, factor)
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_LDH.csv"))

Data clean-up

# Filter only samples matched to CFU wells
LDH_data <- input_data %>% 
  filter(CFUmatched == "T")

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

# Average technical replicates for each individual experiment
LDH_avg <- LDH_data %>%
  group_by(Date, Model, Line, Species, bacteria_label, Temp, Well_Endpoint, Collection_Time) %>% 
  summarise(avg_Value = mean(Value)) %>%
  # Add labels to distinguish final collection times from time 0
  mutate(Collection_Label = if_else(Collection_Time == 6 | Collection_Time == 6, "final", "initial"))

# Calculate fold change between final time point vs. time 0
LDH_FC <- LDH_avg %>%
  group_by(Date, Model, Line, Species, bacteria_label, Temp, Well_Endpoint) %>%
  reframe(FC = avg_Value[Collection_Label == "final"]/avg_Value[Collection_Label == "initial"]) %>%
  mutate(Species.Model = interaction(Species, Model, lex.order = T, sep = ".")) %>%
  mutate(Species.Model_label = interaction(bacteria_label, Model, lex.order = T, sep = "<br>"))

Saving files

# Save data frames as CSV files in the output folder
write_csv(LDH_data, file.path(dataframes_folder, "LDH_values_6h.csv"))
write_csv(LDH_FC, file.path(dataframes_folder, "LDH_FC_6h.csv"))

# Save data frames as R objects in the output folder
saveRDS(LDH_data, file.path(dataframes_folder, "LDH_values_6h.rds"))
saveRDS(LDH_FC, file.path(dataframes_folder, "LDH_FC_6h.rds"))

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

# Use this to read the final objects
LDH_data <- readRDS("data/dataframes/LDH_values_6h.rds")
LDH_FC <- readRDS("data/dataframes/LDH_FC_6h.rds")

Stats and Plots

File Paths

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

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

# Load data and metadata
LDH_FC <- readRDS("data/dataframes/LDH_FC_6h.rds")

Analysis Epithelial Lines

# 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 ~ Species * 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, ~ Species * 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_label) %>%
    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(Species, Model1, sep = "."),
           condition2 = paste(Species, 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("#5b5b5b","#800080","#1E90FF","#927ed1")) +
    scale_color_manual(values = c("#5b5b5b","#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 = "Fold change in LDH from -1 hour",
         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, "/plotLDH_", "HNOvsCells", ".png"), width = 7, height = 6)
  saveRDS(plot_2, file.path(figures_folder, paste0("plotLDH_", "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(LDH_FC, cutoff_pvalue = 0.05, cutoff_FC = 1)

$anova
               Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F) sign
Species       9.2e-02 3.1e-02     3 2.7e+01 5.2e+00 5.7e-03    *
Model         4.5e-02 2.3e-02     2 9.0e+00 3.8e+00 6.2e-02     
Species:Model 7.5e-02 1.2e-02     6 2.7e+01 2.1e+00 8.4e-02     

$random_effects
        grp        var1 var2         vcov        sdcor proportion
1 Line:Date (Intercept) <NA> 8.247084e-03 9.081346e-02       58.4
2      Line (Intercept) <NA> 1.411378e-11 3.756831e-06        0.0
3  Residual        <NA> <NA> 5.874784e-03 7.664714e-02       41.6

$contrasts
     Model1   Model2     Species    estimate         SE        df    t.ratio
1       HNO (Calu-3)         Dpi -0.25630257 0.09257284 0.8138099 -2.7686584
2       HNO     RPMI         Dpi -0.23022073 0.09872463 1.1654776 -2.3319482
3  (Calu-3)     RPMI         Dpi  0.02608184 0.09076209 0.7033925  0.2873650
4       HNO (Calu-3)         Sau -0.18639879 0.09257284 0.8138099 -2.0135365
5       HNO     RPMI         Sau -0.14811415 0.09872463 1.1654776 -1.5002755
6  (Calu-3)     RPMI         Sau  0.03828464 0.09076209 0.7033925  0.4218131
7       HNO (Calu-3)         Spn -0.11118587 0.09257284 0.8138099 -1.2010636
8       HNO     RPMI         Spn -0.22458378 0.09872463 1.1654776 -2.2748505
9  (Calu-3)     RPMI         Spn -0.11339791 0.09076209 0.7033925 -1.2493974
10      HNO (Calu-3) Uncolonized -0.07782702 0.09257284 0.8138099 -0.8407111
11      HNO     RPMI Uncolonized -0.04814995 0.09872463 1.1654776 -0.4877197
12 (Calu-3)     RPMI Uncolonized  0.02967707 0.09076209 0.7033925  0.3269765
     p.value p.adj.holm sign           condition1           condition2
1  0.2643732          1                   Dpi.HNO         Dpi.(Calu-3)
2  0.2285556          1                   Dpi.HNO             Dpi.RPMI
3  0.8357759          1              Dpi.(Calu-3)             Dpi.RPMI
4  0.3351479          1                   Sau.HNO         Sau.(Calu-3)
5  0.3484847          1                   Sau.HNO             Sau.RPMI
6  0.7668712          1              Sau.(Calu-3)             Sau.RPMI
7  0.4746167          1                   Spn.HNO         Spn.(Calu-3)
8  0.2343990          1                   Spn.HNO             Spn.RPMI
9  0.4882296          1              Spn.(Calu-3)             Spn.RPMI
10 0.5793265          1           Uncolonized.HNO Uncolonized.(Calu-3)
11 0.7015759          1           Uncolonized.HNO     Uncolonized.RPMI
12 0.8147842          1      Uncolonized.(Calu-3)     Uncolonized.RPMI
                                          Species.Model_label.1 Model.1
1  <i><b><span style='color:#1E90FF;'>Dpi</span></i></b><br>HNO     HNO
2  <i><b><span style='color:#1E90FF;'>Dpi</span></i></b><br>HNO     HNO
3                                                          <NA>    <NA>
4  <i><b><span style='color:#800080;'>Sau</span></i></b><br>HNO     HNO
5  <i><b><span style='color:#800080;'>Sau</span></i></b><br>HNO     HNO
6                                                          <NA>    <NA>
7  <i><b><span style='color:#927ED1;'>Spn</span></i></b><br>HNO     HNO
8  <i><b><span style='color:#927ED1;'>Spn</span></i></b><br>HNO     HNO
9                                                          <NA>    <NA>
10      <b><span style='color:#5b5b5b;'>Uncol</span></b><br>HNO     HNO
11      <b><span style='color:#5b5b5b;'>Uncol</span></b><br>HNO     HNO
12                                                         <NA>    <NA>
   mean.predval.1                                         Species.Model_label.2
1       0.8357285                                                          <NA>
2       0.8357285 <i><b><span style='color:#1E90FF;'>Dpi</span></i></b><br>RPMI
3              NA <i><b><span style='color:#1E90FF;'>Dpi</span></i></b><br>RPMI
4       0.8120503                                                          <NA>
5       0.8120503 <i><b><span style='color:#800080;'>Sau</span></i></b><br>RPMI
6              NA <i><b><span style='color:#800080;'>Sau</span></i></b><br>RPMI
7       0.8467198                                                          <NA>
8       0.8467198 <i><b><span style='color:#927ED1;'>Spn</span></i></b><br>RPMI
9              NA <i><b><span style='color:#927ED1;'>Spn</span></i></b><br>RPMI
10      0.8344994                                                          <NA>
11      0.8344994      <b><span style='color:#5b5b5b;'>Uncol</span></b><br>RPMI
12             NA      <b><span style='color:#5b5b5b;'>Uncol</span></b><br>RPMI
   Model.2 mean.predval.2        FC highlighted
1     <NA>             NA        NA            
2     RPMI      1.0659492 -1.275473           +
3     RPMI      1.0659492        NA            
4     <NA>             NA        NA            
5     RPMI      0.9601645 -1.182395           +
6     RPMI      0.9601645        NA            
7     <NA>             NA        NA            
8     RPMI      1.0713036 -1.265240           +
9     RPMI      1.0713036        NA            
10    <NA>             NA        NA            
11    RPMI      0.8826493 -1.057699           +
12    RPMI      0.8826493        NA            

$data_summary
# A tibble: 12 × 9
# Groups:   Species.Model_label, Species.Model, Model [12]
   Species.Model_label Species.Model Model bacteria_label mean.real mean.predval
   <fct>               <fct>         <fct> <fct>              <dbl>        <dbl>
 1 <b><span style='co… Uncolonized.… HNO   <b><span styl…     0.834        0.834
 2 <b><span style='co… Uncolonized.… Calu… <b><span styl…     0.912        0.912
 3 <b><span style='co… Uncolonized.… RPMI  <b><span styl…     0.883        0.883
 4 <i><b><span style=… Sau.HNO       HNO   <i><b><span s…     0.812        0.812
 5 <i><b><span style=… Sau.Calu-3    Calu… <i><b><span s…     0.998        0.998
 6 <i><b><span style=… Sau.RPMI      RPMI  <i><b><span s…     0.960        0.960
 7 <i><b><span style=… Dpi.HNO       HNO   <i><b><span s…     0.836        0.836
 8 <i><b><span style=… Dpi.Calu-3    Calu… <i><b><span s…     1.09         1.09 
 9 <i><b><span style=… Dpi.RPMI      RPMI  <i><b><span s…     1.07         1.07 
10 <i><b><span style=… Spn.HNO       HNO   <i><b><span s…     0.847        0.847
11 <i><b><span style=… Spn.Calu-3    Calu… <i><b><span s…     0.958        0.958
12 <i><b><span style=… Spn.RPMI      RPMI  <i><b><span s…     1.07         1.07 
# ℹ 3 more variables: mean.predval.se <dbl>, max <dbl>, min <dbl>

$plot_2