diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index 65e2ffb1..93e25e03 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -377,7 +377,12 @@ calculate_interaction_scores <- function(df, max_conc, variables, group_vars = c
     arrange(desc(NG)) %>%
     ungroup()
 
-  return(list(calculations = calculations, interactions = interactions))
+  df <- df %>% select(-any_of(setdiff(names(calculations), group_vars)))
+  df <- left_join(df, calculations, by = group_vars)
+  df <- df %>% select(-any_of(setdiff(names(interactions), group_vars)))
+  df <- left_join(df, interactions, by = group_vars)
+
+  return(list(calculations = calculations, interactions = interactions, joined = df))
 }
 
 generate_and_save_plots <- function(output_dir, file_name, plot_configs, grid_layout = NULL) {
@@ -391,11 +396,20 @@ generate_and_save_plots <- function(output_dir, file_name, plot_configs, grid_la
     #   "delta_bg_tolerance", "delta_bg", "Gene", "L", "K", "r", "AUC", "NG", "DB"))), n = 5)
 
     # Plots are testy about missing aesthetics, so handle them here
-    aes_mapping <- aes(
-      x = !!sym(config$x_var),
-      y = if (!is.null(config$y_var)) !!sym(config$y_var) else NULL,
-      color = if (!is.null(config$color_var)) as.factor(!!sym(config$color_var)) else NULL
-    )
+    aes_mapping <-
+      if (is.null(config$color_var)) {
+        if (is.null(config$y_var)) {
+          aes(x = !!sym(config$x_var))
+        } else {
+          aes(x = !!sym(config$x_var), y = !!sym(config$y_var))
+        }
+      } else {
+        if (is.null(config$y_var)) {
+          aes(x = !!sym(config$x_var), color = as.factor(!!sym(config$color_var)))
+        } else {
+          aes(x = !!sym(config$x_var), y = !!sym(config$y_var), color = as.factor(!!sym(config$color_var)))
+        }
+      }
 
     # Start building the plot
     plot <- ggplot(df, aes_mapping)
@@ -553,10 +567,10 @@ generate_interaction_plot_configs <- function(df, variables) {
 
   # Define annotation positions based on the variable being plotted
   annotation_positions <- list(
-    L = list(ZShift = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
-    K = list(ZShift = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
-    r = list(ZShift = 0.45, lm_ZScore = 0.25, NG = -0.25, DB = -0.35, SM = -0.45),
-    AUC = list(ZShift = 4500, lm_ZScore = 2500, NG = -2500, DB = -3500, SM = -4500)
+    L = list(Z_Shift_L = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
+    K = list(Z_Shift_K = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
+    r = list(Z_Shift_r = 0.45, lm_ZScore = 0.25, NG = -0.25, DB = -0.35, SM = -0.45),
+    AUC = list(Z_Shift_AUC = 4500, lm_ZScore = 2500, NG = -2500, DB = -3500, SM = -4500)
   )
 
   # Define which annotations to include for each plot
@@ -586,7 +600,7 @@ generate_interaction_plot_configs <- function(df, variables) {
     # Dynamically generate the names of the columns
     var_info <- list(
       ylim = limits_map[[variable]],
-      lm_model = df[[paste0("lm_", variable)]][[1]],    # Access the precomputed linear model
+      lm_model = df[[paste0("lm_", variable)]][[1]],
       sd_col = paste0("WT_sd_", variable),
       delta_var = paste0("Delta_", variable)
     )
@@ -1086,8 +1100,10 @@ main <- function() {
       
       zscores_calculations_reference <- reference_results$calculations
       zscores_interactions_reference <- reference_results$interactions
+      zscores_joined_reference <- reference_results$joined
       zscores_calculations <- deletion_results$calculations
       zscores_interactions <- deletion_results$interactions
+      zscores_joined <- deletion_results$joined
       
       # Writing Z-Scores to file
       write.csv(zscores_calculations_reference, file = file.path(out_dir, "RF_ZScores_Calculations.csv"), row.names = FALSE)
@@ -1097,8 +1113,8 @@ main <- function() {
 
       # Create interaction plots
       message("Generating interaction plot configurations")
-      reference_plot_configs <- generate_interaction_plot_configs(df_reference, interaction_vars)
-      deletion_plot_configs <- generate_interaction_plot_configs(df_deletion, interaction_vars)
+      reference_plot_configs <- generate_interaction_plot_configs(zscores_joined_reference, interaction_vars)
+      deletion_plot_configs <- generate_interaction_plot_configs(zscores_joined, interaction_vars)
       message("Generating interaction plots")
       generate_and_save_plots(out_dir, "RF_interactionPlots", reference_plot_configs, grid_layout = list(ncol = 4, nrow = 3))
       generate_and_save_plots(out_dir, "InteractionPlots", deletion_plot_configs, grid_layout = list(ncol = 4, nrow = 3))