Refactor interaction plots to handle both reference and deletion scoring

qhtcp-workflow/apps/r/calculate_interaction_zscores.R

@@ -150,8 +150,9 @@ load_and_filter_data <- function(easy_results_file, sd = 3) {
       SM = 0,
       OrfRep = if_else(ORF == "YDL227C", "YDL227C", OrfRep), # should these be hardcoded?
       conc_num = as.numeric(gsub("[^0-9\\.]", "", Conc)),
-      conc_num_factor = factor(as.numeric(factor(conc_num)) - 1),
-      conc_num_factor_num = as.numeric(conc_num_factor)
+      conc_num_factor_new = factor(conc_num),
+      conc_num_factor_zeroed = factor(as.numeric(conc_num_factor2) - 1),
+      conc_num_factor = as.numeric(conc_num_factor_zeroed) # for legacy purposes, neither conc_num nor a factor
     )
     
     return(df)
@@ -250,10 +251,10 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats,
       Zscore_AUC = Delta_AUC / WT_sd_AUC,
 
       # Fit linear models and store in list-columns
-      gene_lm_L = list(lm(Delta_L ~ conc_num_factor_num, data = pick(everything()))),
-      gene_lm_K = list(lm(Delta_K ~ conc_num_factor_num, data = pick(everything()))),
-      gene_lm_r = list(lm(Delta_r ~ conc_num_factor_num, data = pick(everything()))),
-      gene_lm_AUC = list(lm(Delta_AUC ~ conc_num_factor_num, data = pick(everything()))),
+      gene_lm_L = list(lm(Delta_L ~ conc_num_factor_zeroed_num, data = pick(everything()))),
+      gene_lm_K = list(lm(Delta_K ~ conc_num_factor_zeroed_num, data = pick(everything()))),
+      gene_lm_r = list(lm(Delta_r ~ conc_num_factor_zeroed_num, data = pick(everything()))),
+      gene_lm_AUC = list(lm(Delta_AUC ~ conc_num_factor_zeroed_num, data = pick(everything()))),
 
       # Extract coefficients using purrr::map_dbl
       lm_intercept_L = map_dbl(gene_lm_L, ~ coef(.x)[1]),
@@ -293,12 +294,12 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats,
     )
 
   calculations <- calculations %>%
-  mutate(
-    Z_lm_L = (lm_Score_L - lm_means_sds$lm_mean_L) / lm_means_sds$lm_sd_L,
-    Z_lm_K = (lm_Score_K - lm_means_sds$lm_mean_K) / lm_means_sds$lm_sd_K,
-    Z_lm_r = (lm_Score_r - lm_means_sds$lm_mean_r) / lm_means_sds$lm_sd_r,
-    Z_lm_AUC = (lm_Score_AUC - lm_means_sds$lm_mean_AUC) / lm_means_sds$lm_sd_AUC
-  )
+    mutate(
+      Z_lm_L = (lm_Score_L - lm_means_sds$lm_mean_L) / lm_means_sds$lm_sd_L,
+      Z_lm_K = (lm_Score_K - lm_means_sds$lm_mean_K) / lm_means_sds$lm_sd_K,
+      Z_lm_r = (lm_Score_r - lm_means_sds$lm_mean_r) / lm_means_sds$lm_sd_r,
+      Z_lm_AUC = (lm_Score_AUC - lm_means_sds$lm_mean_AUC) / lm_means_sds$lm_sd_AUC
+    )
 
   # Summarize some of the stats
   interactions <- calculations %>%
@@ -321,7 +322,7 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats,
       Sum_Z_Score_K = sum(Zscore_K),
       Sum_Z_Score_r = sum(Zscore_r),
       Sum_Z_Score_AUC = sum(Zscore_AUC),
-      
+
       # Calculate Average Z-scores
       Avg_Zscore_L = Sum_Z_Score_L / num_non_removed_concs,
       Avg_Zscore_K = Sum_Z_Score_K / num_non_removed_concs,
@@ -346,7 +347,8 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats,
       "Exp_L", "Exp_K", "Exp_r", "Exp_AUC",
       "Delta_L", "Delta_K", "Delta_r", "Delta_AUC",
       "Zscore_L", "Zscore_K", "Zscore_r", "Zscore_AUC",
-      "NG", "SM", "DB")
+      "NG", "SM", "DB"
+    )
 
   interactions <- interactions %>%
     select(
@@ -357,30 +359,49 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats,
       "Avg_Zscore_L", "Avg_Zscore_K", "Avg_Zscore_r", "Avg_Zscore_AUC",
       "lm_Score_L", "lm_Score_K", "lm_Score_r", "lm_Score_AUC",
       "R_Squared_L", "R_Squared_K", "R_Squared_r", "R_Squared_AUC",
-      "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC")
-    
+      "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC"
+    )
+
+  # Clean the original dataframe by removing overlapping columns
   cleaned_df <- df %>%
     select(-any_of(
-      setdiff(intersect(names(df), names(interactions)),
+      setdiff(intersect(names(df), names(calculations)),
+      c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))))
+
+  # Join the original dataframe with calculations
+  df_with_calculations <- left_join(cleaned_df, calculations, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))
+
+  # Remove overlapping columns between df_with_calculations and interactions
+  df_with_calculations_clean <- df_with_calculations %>%
+    select(-any_of(
+      setdiff(intersect(names(df_with_calculations), names(interactions)),
       c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))))
 
-  interactions_joined <- left_join(cleaned_df, interactions, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))
+  # Join with interactions to create the full dataset
+  full_data <- left_join(df_with_calculations_clean, interactions, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))
 
   return(list(
     calculations = calculations,
     interactions = interactions,
-    interactions_joined = interactions_joined))
+    full_data = full_data
+  ))
 }
 
 generate_and_save_plots <- function(out_dir, filename, plot_configs) {
   message("Generating ", filename, ".pdf and ", filename, ".html")
 
-  # Iterate through the plot_configs (which contain both plots and grid_layout)
   for (config_group in plot_configs) {
     plot_list <- config_group$plots
     grid_nrow <- config_group$grid_layout$nrow
     grid_ncol <- config_group$grid_layout$ncol
 
+    # Set defaults if nrow or ncol are not provided
+    if (is.null(grid_nrow) || is.null(grid_ncol)) {
+      num_plots <- length(plot_list)
+      grid_nrow <- ifelse(is.null(grid_nrow), 1, grid_nrow)
+      grid_ncol <- ifelse(is.null(grid_ncol), num_plots, grid_ncol)
+    }
+
     # Prepare lists to collect static and interactive plots
     static_plots <- list()
     plotly_plots <- list()
@@ -419,11 +440,11 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
 
       # Use appropriate helper function based on plot type
       plot <- switch(config$plot_type,
-                     "scatter" = generate_scatter_plot(plot, config),
-                     "box" = generate_box_plot(plot, config),
-                     "density" = plot + geom_density(),
-                     "bar" = plot + geom_bar(),
-                     plot  # default case if no type matches
+        "scatter" = generate_scatter_plot(plot, config),
+        "box" = generate_box_plot(plot, config),
+        "density" = plot + geom_density(),
+        "bar" = plot + geom_bar(),
+        plot  # default case if no type matches
       )
 
       # Add title and labels
@@ -462,9 +483,9 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
       # Convert to plotly object and suppress warnings here
       plotly_plot <- suppressWarnings({
         if (length(tooltip_vars) > 0) {
-          ggplotly(plot, tooltip = tooltip_vars)
+          plotly::ggplotly(plot, tooltip = tooltip_vars)
         } else {
-          ggplotly(plot, tooltip = "none")
+          plotly::ggplotly(plot, tooltip = "none")
         }
       })
 
@@ -483,8 +504,7 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
     grid.arrange(grobs = static_plots, ncol = grid_ncol, nrow = grid_nrow)
     dev.off()
 
-    # Combine and save interactive HTML plot(s)
-    combined_plot <- subplot(
+    combined_plot <- plotly::subplot(
       plotly_plots,
       nrows = grid_nrow,
       ncols = grid_ncol,
@@ -492,7 +512,8 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
     )
 
     # Save combined HTML plot(s)
-    saveWidget(combined_plot, file = file.path(out_dir, paste0(filename, ".html")), selfcontained = TRUE)
+    html_file <- file.path(out_dir, paste0(filename, ".html"))
+    saveWidget(combined_plot, file = html_file, selfcontained = TRUE)
   }
 }
 
@@ -635,8 +656,10 @@ generate_scatter_plot <- function(plot, config) {
 }
 
 generate_box_plot <- function(plot, config) {
-  plot <- plot + geom_boxplot()
-  
+  # Convert x_var to a factor within aes mapping
+  plot <- plot + geom_boxplot(aes(x = factor(.data[[config$x_var]])))
+
+  # Apply scale_x_discrete for breaks, labels, and axis label if provided
   if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) {
     plot <- plot + scale_x_discrete(
       name = config$x_label,
@@ -681,7 +704,7 @@ generate_plate_analysis_plot_configs <- function(variables, stages = c("before",
         plot_type = plot_type,
         title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"),
         error_bar = error_bar,
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         position = position,
         size = 0.2
       )
@@ -691,75 +714,127 @@ generate_plate_analysis_plot_configs <- function(variables, stages = c("before",
   return(plots)
 }
 
-generate_interaction_plot_configs <- function(df, limits_map = NULL, stats_df = NULL) {
+generate_interaction_plot_configs <- function(df, limits_map = NULL, plot_type = "reference") {
+  # Define limits if not provided
   if (is.null(limits_map)) {
     limits_map <- list(
       L = c(0, 130),
       K = c(-20, 160),
       r = c(0, 1),
-      AUC = c(0, 12500)
+      AUC = c(0, 12500),
+      Delta_L = c(-60, 60),
+      Delta_K = c(-60, 60),
+      Delta_r = c(-0.6, 0.6),
+      Delta_AUC = c(-6000, 6000)
     )
   }
 
-  # Ensure proper grouping by OrfRep, Gene, and num
-  df_filtered <- df %>%
-    filter(
-      !is.na(L) & L >= limits_map$L[1] & L <= limits_map$L[2],
-      !is.na(K) & K >= limits_map$K[1] & K <= limits_map$K[2],
-      !is.na(r) & r >= limits_map$r[1] & r <= limits_map$r[2],
-      !is.na(AUC) & AUC >= limits_map$AUC[1] & AUC <= limits_map$AUC[2]
-    ) %>%
-    group_by(OrfRep, Gene, num)  # Group by OrfRep, Gene, and num
+  # Define grouping variables and filter data based on plot type
+  if (plot_type == "reference") {
+    group_vars <- c("OrfRep", "Gene", "num")
+    df_filtered <- df %>%
+      mutate(
+        OrfRepCombined = paste(OrfRep, Gene, num, sep = "_")
+      )
+  } else if (plot_type == "deletion") {
+    group_vars <- c("OrfRep", "Gene")
+    df_filtered <- df %>%
+      mutate(
+        OrfRepCombined = paste(OrfRep, Gene, sep = "_") # Compare by OrfRep and Gene for deletion
+      )
+  }
 
-  scatter_configs <- list()
-  box_configs <- list()
+  # Create a list to store all configs
+  configs <- list()
 
-  # Generate scatter and box plots for each variable (L, K, r, AUC)
-  for (var in names(limits_map)) {
-    scatter_configs[[length(scatter_configs) + 1]] <- list(
+  # Generate the first 8 scatter/box plots for L, K, r, AUC (shared between reference and deletion)
+  overall_vars <- c("L", "K", "r", "AUC")
+  for (var in overall_vars) {
+    y_limits <- limits_map[[var]]
+
+    config <- list(
       df = df_filtered,
-      x_var = "conc_num",  # X-axis variable
-      y_var = var,         # Y-axis variable (Delta_L, Delta_K, Delta_r, Delta_AUC)
       plot_type = "scatter",
+      x_var = "conc_num_factor_new",
+      y_var = var,
+      x_label = unique(df_filtered$Drug)[1],
       title = sprintf("Scatter RF for %s with SD", var),
-      coord_cartesian = limits_map[[var]],  # Set limits for Y-axis
-      annotations = list(
-        list(x = -0.25, y = 10, label = "NG"),
-        list(x = -0.25, y = 5, label = "DB"),
-        list(x = -0.25, y = 0, label = "SM")
-      ),
-      grid_layout = list(ncol = 4, nrow = 3)
-    )
-    box_configs[[length(box_configs) + 1]] <- list(
-      df = df_filtered,
-      x_var = "conc_num",  # X-axis variable
-      y_var = var,         # Y-axis variable (Delta_L, Delta_K, Delta_r, Delta_AUC)
-      plot_type = "box",
-      title = sprintf("Boxplot RF for %s with SD", var),
-      coord_cartesian = limits_map[[var]],
-      grid_layout = list(ncol = 4, nrow = 3)
+      coord_cartesian = y_limits,
+      error_bar = TRUE,
+      x_breaks = unique(df_filtered$conc_num_factor_new),
+      x_labels = as.character(unique(df_filtered$conc_num)),
+      grid_layout = list(ncol = 2, nrow = 2)
     )
+    configs <- append(configs, list(config))
   }
 
-  # Combine scatter and box plots into grids
-  configs <- list(
-    list(
-      grid_layout = list(nrow = 2, ncol = 2),  # Scatter plots in a 2x2 grid (for the 8 plots)
-      plots = scatter_configs[1:4]
-    ),
-    list(
-      grid_layout = list(nrow = 2, ncol = 2),  # Box plots in a 2x2 grid (for the 8 plots)
-      plots = box_configs
-    ),
-    list(
-      grid_layout = list(nrow = 3, ncol = 4),  # Delta_ plots in a 3x4 grid
-      plots = scatter_configs
-    ),
-    list(
-      grid_layout = list(nrow = 3, ncol = 4),  # Delta_ box plots in a 3x4 grid
-      plots = box_configs
-    )
-  )
+  # Generate Delta comparison plots (4x3 grid for deletion and reference)
+  unique_groups <- df_filtered %>% select(all_of(group_vars)) %>% distinct()
+
+  for (i in seq_len(nrow(unique_groups))) {
+    group <- unique_groups[i, ]
+    group_data <- df_filtered %>% filter(across(all_of(group_vars), ~ . == group[[cur_column()]]))
+
+    OrfRep <- as.character(group$OrfRep)
+    Gene <- if ("Gene" %in% names(group)) as.character(group$Gene) else ""
+    num <- if ("num" %in% names(group)) as.character(group$num) else ""
+    OrfRepCombined <- paste(OrfRep, Gene, num, sep = "_")
+
+    # Generate plots for Delta variables
+    delta_vars <- c("Delta_L", "Delta_K", "Delta_r", "Delta_AUC")
+    for (var in delta_vars) {
+      y_limits <- limits_map[[var]]
+      upper_y <- y_limits[2]
+      lower_y <- y_limits[1]
+      y_span <- upper_y - lower_y
+
+      # Get WT_sd_var for error bar calculations
+      WT_sd_var <- paste0("WT_sd_", sub("Delta_", "", var))
+      WT_sd_value <- group_data[[WT_sd_var]][1]
+      error_bar_ymin <- 0 - (2 * WT_sd_value)
+      error_bar_ymax <- 0 + (2 * WT_sd_value)
+
+      # Set annotations (Z_Shifts, lm Z-Scores, NG, DB, SM values)
+      Z_Shift_var <- paste0("Z_Shift_", sub("Delta_", "", var))
+      Z_lm_var <- paste0("Z_lm_", sub("Delta_", "", var))
+      Z_Shift_value <- round(group_data[[Z_Shift_var]][1], 2)
+      Z_lm_value <- round(group_data[[Z_lm_var]][1], 2)
+      NG_value <- group_data$NG[1]
+      DB_value <- group_data$DB[1]
+      SM_value <- group_data$SM[1]
+
+      annotations <- list(
+        list(x = 1, y = upper_y - 0.2 * y_span, label = paste("ZShift =", Z_Shift_value)),
+        list(x = 1, y = upper_y - 0.3 * y_span, label = paste("lm ZScore =", Z_lm_value)),
+        list(x = 1, y = lower_y + 0.2 * y_span, label = paste("NG =", NG_value)),
+        list(x = 1, y = lower_y + 0.1 * y_span, label = paste("DB =", DB_value)),
+        list(x = 1, y = lower_y, label = paste("SM =", SM_value))
+      )
+
+      # Create configuration for each Delta plot
+      config <- list(
+        df = group_data,
+        plot_type = "scatter",
+        x_var = "conc_num",
+        y_var = var,
+        x_label = unique(group_data$Drug)[1],
+        title = paste(OrfRep, Gene, sep = "      "),
+        coord_cartesian = y_limits,
+        annotations = annotations,
+        error_bar = TRUE,
+        error_bar_params = list(
+          ymin = error_bar_ymin,
+          ymax = error_bar_ymax
+        ),
+        lm_smooth = TRUE,
+        x_breaks = unique(group_data$conc_num_factor_new),
+        x_labels = as.character(unique(group_data$conc_num)),
+        ylim_vals = y_limits,
+        grid_layout = list(ncol = 4, nrow = 3) # Adjust grid layout for gene-gene comparisons
+      )
+      configs <- append(configs, list(config))
+    }
+  }
 
   return(configs)
 }
@@ -826,14 +901,14 @@ generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FA
         annotations = list(
           list(
             x = median(df_ranked[[rank_var]], na.rm = TRUE),
-            y = 10,
+            y = max(df_ranked[[zscore_var]], na.rm = TRUE) * 0.9,
             label = paste("Deletion Enhancers =", num_enhancers),
             hjust = 0.5,
             vjust = 1
           ),
           list(
             x = median(df_ranked[[rank_var]], na.rm = TRUE),
-            y = -10,
+            y = min(df_ranked[[zscore_var]], na.rm = TRUE) * 0.9,
             label = paste("Deletion Suppressors =", num_suppressors),
             hjust = 0.5,
             vjust = 0
@@ -870,7 +945,7 @@ generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FA
     }
   }
     
-  # Avg ZScore and Rank Avg ZScore Plots for r, L, K, and AUC
+  # Avg ZScore and Rank Avg ZScore Plots for variables
   for (variable in variables) {
     for (plot_type in c("Avg Zscore vs lm", "Rank Avg Zscore vs lm")) {
 
@@ -894,32 +969,30 @@ generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FA
       # Fit the linear model
       lm_model <- lm(as.formula(paste(y_var, "~", x_var)), data = df_ranked)
       
-      # Extract intercept and slope from the model coefficients
+      # Extract intercept, slope, and R-squared from the model coefficients
       intercept <- coef(lm_model)[1]
       slope <- coef(lm_model)[2]
-        
+      r_squared <- summary(lm_model)$r.squared
+
+      # Annotations: include R-squared in the plot
+      annotations <- list(
+        list(
+          x = mean(range(df_ranked[[x_var]], na.rm = TRUE)),
+          y = mean(range(df_ranked[[y_var]], na.rm = TRUE)),
+          label = paste("R-squared =", round(r_squared, 2)),
+          hjust = 0.5,
+          vjust = 1,
+          size = 5
+        )
+      )
+
       configs[[length(configs) + 1]] <- list(
         df = df_ranked,
         x_var = x_var,
         y_var = y_var,
         plot_type = "scatter",
         title = title,
-        annotations = list(
-          list(
-            x = median(df_ranked[[rank_var]], na.rm = TRUE),
-            y = 10,
-            label = paste("Deletion Enhancers =", num_enhancers),
-            hjust = 0.5,
-            vjust = 1
-          ),
-          list(
-            x = median(df_ranked[[rank_var]], na.rm = TRUE),
-            y = -10,
-            label = paste("Deletion Suppressors =", num_suppressors),
-            hjust = 0.5,
-            vjust = 0
-          )
-        ),
+        annotations = annotations,
         shape = 3,
         size = 0.25,
         smooth = TRUE,
@@ -936,7 +1009,7 @@ generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FA
   return(configs)
 }
 
-generate_correlation_plot_configs <- function(df) {
+generate_correlation_plot_configs <- function(df, highlight_cyan = FALSE) {
   # Define relationships for plotting
   relationships <- list(
     list(x = "Z_lm_L", y = "Z_lm_K", label = "Interaction L vs. Interaction K"),
@@ -953,6 +1026,9 @@ generate_correlation_plot_configs <- function(df) {
     # Fit linear model
     lm_model <- lm(as.formula(paste(rel$y, "~", rel$x)), data = df)
     lm_summary <- summary(lm_model)
+    intercept <- coef(lm_model)[1]
+    slope <- coef(lm_model)[2]
+    r_squared <- lm_summary$r.squared
 
     # Construct plot configuration
     config <- list(
@@ -965,18 +1041,18 @@ generate_correlation_plot_configs <- function(df) {
       y_label = paste("z-score", gsub("Z_lm_", "", rel$y)),
       annotations = list(
         list(
-          x = Inf,
-          y = Inf,
-          label = paste("R-squared =", round(lm_summary$r.squared, 3)),
-          hjust = 1.1,
-          vjust = 2,
-          size = 4,
+          x = mean(range(df[[rel$x]], na.rm = TRUE)),
+          y = mean(range(df[[rel$y]], na.rm = TRUE)),
+          label = paste("R-squared =", round(r_squared, 3)),
+          hjust = 0.5,
+          vjust = 1,
+          size = 5,
           color = "black"
         )
       ),
       smooth = TRUE,
       smooth_color = "tomato3",
-      lm_line = list(intercept = coef(lm_model)[1], slope = coef(lm_model)[2]),
+      lm_line = list(intercept = intercept, slope = slope),
       legend_position = "right",
       shape = 3,
       size = 0.5,
@@ -987,8 +1063,7 @@ generate_correlation_plot_configs <- function(df) {
           fill = NA, color = "grey20", alpha = 0.1
         )
       ),
-      cyan_points = TRUE,
-      grid_layout = list(ncol = 2, nrow = 2)
+      cyan_points = highlight_cyan,  # Toggle cyan point highlighting
     )
 
     configs[[length(configs) + 1]] <- config
@@ -1023,7 +1098,7 @@ main <- function() {
     df_no_zeros <- df_na %>% filter(L > 0) # formerly X_noZero
     
     # Save some constants
-    max_conc <- max(df$conc_num_factor_num)
+    max_conc <- max(df$conc_num_factor_zeroed_num)
     l_half_median <- (median(df_above_tolerance$L, na.rm = TRUE)) / 2
     k_half_median <- (median(df_above_tolerance$K, na.rm = TRUE)) / 2
 
@@ -1106,7 +1181,7 @@ main <- function() {
         plot_type = "scatter",
         delta_bg_point = TRUE,
         title = "Raw L vs K before quality control",
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         error_bar = FALSE,
         legend_position = "right"
       )
@@ -1119,7 +1194,7 @@ main <- function() {
         y_var = NULL,
         plot_type = "density",
         title = "Density plot for Delta Background by [Drug] (All Data)",
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         x_label = "Delta Background",
         y_label = "Density",
         error_bar = FALSE,
@@ -1130,7 +1205,7 @@ main <- function() {
         y_var = NULL,
         plot_type = "bar",
         title = "Bar plot for Delta Background by [Drug] (All Data)",
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         x_label = "Delta Background",
         y_label = "Count",
         error_bar = FALSE,
@@ -1146,7 +1221,7 @@ main <- function() {
         delta_bg_point = TRUE,
         title = paste("Raw L vs K for strains above Delta Background threshold of",
           round(df_above_tolerance$delta_bg_tolerance[[1]], 3), "or above"),
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         position = "jitter",
         annotations = list(
           list(
@@ -1194,7 +1269,7 @@ main <- function() {
         plot_type = "scatter",
         delta_bg_point = TRUE,
         title = "Raw L vs K for strains falling outside 2SD of the K mean at each Conc",
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         position = "jitter",
         legend_position = "right"
       )
@@ -1208,7 +1283,7 @@ main <- function() {
         plot_type = "scatter",
         gene_point = TRUE,
         title = "Delta Background vs K for strains falling outside 2SD of the K mean at each Conc",
-        color_var = "conc_num",
+        color_var = "conc_num_factor_zeroed",
         position = "jitter",
         legend_position = "right"
       )
@@ -1305,7 +1380,7 @@ main <- function() {
       reference_results <- calculate_interaction_scores(df_reference_stats, max_conc, bg_stats, group_vars = c("OrfRep", "Gene", "num"))
       zscore_calculations_reference <- reference_results$calculations
       zscore_interactions_reference <- reference_results$interactions
-      zscore_interactions_reference_joined <- reference_results$interactions_joined
+      zscore_interactions_reference_joined <- reference_results$full_data
 
       message("Calculating deletion strain(s) interactions scores")
       df_deletion_stats <- calculate_summary_stats(
@@ -1316,7 +1391,7 @@ main <- function() {
       deletion_results <- calculate_interaction_scores(df_deletion_stats, max_conc, bg_stats, group_vars = c("OrfRep"))
       zscore_calculations <- deletion_results$calculations
       zscore_interactions <- deletion_results$interactions
-      zscore_interactions_joined <- deletion_results$interactions_joined
+      zscore_interactions_joined <- deletion_results$full_data
 
       # Writing Z-Scores to file
       write.csv(zscore_calculations_reference, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE)
@@ -1326,11 +1401,11 @@ main <- function() {
 
       # Create interaction plots
       message("Generating reference interaction plots")
-      reference_plot_configs <- generate_interaction_plot_configs(zscore_interactions_reference_joined)
+      reference_plot_configs <- generate_interaction_plot_configs(zscore_interactions_reference_joined, plot_type = "reference")
       generate_and_save_plots(out_dir, "interaction_plots_reference", reference_plot_configs)
 
       message("Generating deletion interaction plots")
-      deletion_plot_configs <- generate_interaction_plot_configs(zscore_interactions_joined)
+      deletion_plot_configs <- generate_interaction_plot_configs(zscore_interactions_joined, plot_type = "deletion")
       generate_and_save_plots(out_dir, "interaction_plots", deletion_plot_configs)
 
       # Define conditions for enhancers and suppressors