diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index 0cc4252e..18436a59 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -304,6 +304,14 @@ calculate_interaction_scores <- function(df, max_conc, variables, group_vars = c
       R_Squared_K = summary(lm_K)$r.squared,
       R_Squared_r = summary(lm_r)$r.squared,
       R_Squared_AUC = summary(lm_AUC)$r.squared,
+      lm_intercept_L = coef(lm_L)[1],
+      lm_slope_L = coef(lm_L)[2],
+      lm_intercept_K = coef(lm_K)[1],
+      lm_slope_K = coef(lm_K)[2],
+      lm_intercept_r = coef(lm_r)[1],
+      lm_slope_r = coef(lm_r)[2],
+      lm_intercept_AUC = coef(lm_AUC)[1],
+      lm_slope_AUC = coef(lm_AUC)[2],
       NG = sum(NG, na.rm = TRUE),
       DB = sum(DB, na.rm = TRUE),
       SM = sum(SM, na.rm = TRUE)
@@ -342,12 +350,14 @@ calculate_interaction_scores <- function(df, max_conc, variables, group_vars = c
       "Zscore_L", "Zscore_K", "Zscore_r", "Zscore_AUC",
       "NG", "SM", "DB")
     
-  df <- df %>% select(-any_of(setdiff(names(calculations), c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))))
-  df <- left_join(df, calculations, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))
-  # df <- df %>% select(-any_of(setdiff(names(interactions), group_vars)))
-  # df <- left_join(df, interactions, by = group_vars)
+  calculations_joined <- df %>% select(-any_of(setdiff(names(calculations), c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))))
+  calculations_joined <- left_join(calculations_joined, calculations, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))
+  
+  interactions_joined <- df %>% select(-any_of(setdiff(names(interactions), c("OrfRep", "Gene", "num"))))
+  interactions_joined <- left_join(interactions_joined, interactions, by = c("OrfRep", "Gene", "num"))
 
-  return(list(calculations = calculations, interactions = interactions, joined = df))
+  return(list(calculations = calculations, interactions = interactions, interactions_joined = interactions_joined,
+    calculations_joined = calculations_joined))
 }
 
 generate_and_save_plots <- function(output_dir, file_name, plot_configs, grid_layout = NULL) {
@@ -554,7 +564,6 @@ generate_interaction_plot_configs <- function(df, variables) {
       ylim = limits_map[[variable]],
       lm_model = df[[paste0("lm_", variable)]][[1]],
       sd_col = paste0("WT_sd_", variable),
-      delta_var = paste0("Delta_", variable)
     )
 
     # Extract the precomputed linear model coefficients
@@ -563,8 +572,6 @@ generate_interaction_plot_configs <- function(df, variables) {
       slope = coef(var_info$lm_model)[2]
     )
 
-    # Dynamically create annotations based on variable
-    # Dynamically create annotations based on variable
     annotations <- lapply(names(annotation_positions[[variable]]), function(annotation_name) {
       y_pos <- annotation_positions[[variable]][[annotation_name]]
       
@@ -1030,23 +1037,23 @@ main <- function() {
           L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>%
         ungroup()
 
-      # Calculate interactions
-      interaction_vars <- c("L", "K", "r", "AUC")
       message("Calculating interaction scores")
-      # print("Reference strain:")
-      # print(head(reference_strain))
+      interaction_vars <- c("L", "K", "r", "AUC")
+
+      message("Calculating reference strain(s)")
       reference_results <- calculate_interaction_scores(reference_strain, max_conc, interaction_vars, group_vars = orf_group_vars)
-      # print("Deletion strains:")
-      # print(head(deletion_strains))
-      deletion_results <- calculate_interaction_scores(deletion_strains, max_conc, interaction_vars, group_vars = orf_group_vars)
-      
       zscores_calculations_reference <- reference_results$calculations
       zscores_interactions_reference <- reference_results$interactions
-      zscores_joined_reference <- reference_results$joined
+      zscores_calculations_reference_joined <- reference_results$calculations_joined
+      zscores_interactions_reference_joined <- reference_results$interactions_joined
+
+      message("Calculating deletion strain(s)")
+      deletion_results <- calculate_interaction_scores(deletion_strains, max_conc, interaction_vars, group_vars = orf_group_vars)
       zscores_calculations <- deletion_results$calculations
       zscores_interactions <- deletion_results$interactions
-      zscores_joined <- deletion_results$joined
-      
+      zscores_calculations_joined <- deletion_results$calculations_joined
+      zscores_interactions_joined <- deletion_results$interactions_joined
+
       # Writing Z-Scores to file
       write.csv(zscores_calculations_reference, file = file.path(out_dir, "RF_ZScores_Calculations.csv"), row.names = FALSE)
       write.csv(zscores_calculations, file = file.path(out_dir, "ZScores_Calculations.csv"), row.names = FALSE)
@@ -1055,8 +1062,8 @@ main <- function() {
 
       # Create interaction plots
       message("Generating interaction plot configurations")
-      reference_plot_configs <- generate_interaction_plot_configs(zscores_joined_reference, interaction_vars)
-      deletion_plot_configs <- generate_interaction_plot_configs(zscores_joined, interaction_vars)
+      reference_plot_configs <- generate_interaction_plot_configs(zscores_interactions_reference_joined, interaction_vars)
+      deletion_plot_configs <- generate_interaction_plot_configs(zscores_interactions_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))