diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index 9318a042..72268120 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -234,7 +234,7 @@ calculate_interaction_scores <- function(df, df_bg, type, overlap_threshold = 2)
 
   # Join WT statistics to df
   df <- df %>%
-    left_join(wt_stats, by = bg_group_vars)
+    left_join(wt_stats, by = c(bg_group_vars))
 
   # Compute mean values at zero concentration
   mean_zeroes <- df %>%
@@ -370,7 +370,7 @@ calculate_interaction_scores <- function(df, df_bg, type, overlap_threshold = 2)
       lm_slope_AUC = map_dbl(lm_AUC, "slope"),
       R_Squared_AUC = map_dbl(lm_AUC, "r_squared"),
       lm_Score_AUC = map_dbl(lm_AUC, "score")
-) %>%
+    ) %>%
     select(-lm_L, -lm_K, -lm_r, -lm_AUC) %>%
     ungroup()
 
@@ -730,8 +730,7 @@ calculate_interaction_scores <- function(df, df_bg, type, overlap_threshold = 2)
       "Raw_Shift_L", "Raw_Shift_K", "Raw_Shift_r", "Raw_Shift_AUC",
       "R_Squared_L", "R_Squared_K", "R_Squared_r", "R_Squared_AUC",
       # we need these for the interactions but the original code has the same names in both datasets
-      "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC",
-      "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC"
+      "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC"
     )))
 
   full_data <- calculations_no_overlap %>%
@@ -1435,6 +1434,12 @@ generate_correlation_plot_configs <- function(df, df_reference) {
       x_var <- paste0("Z_lm_", rel$x) # predictor
       y_var <- paste0("Z_lm_", rel$y) # response
 
+      # Skip this plot if there's no valid data
+      # if (all(is.na(df[[x_var]])) || all(is.na(df_reference[[x_var]])) ||
+      #   all(is.na(df[[y_var]])) || all(is.na(df_reference[[y_var]]))) {
+      #   next
+      # }
+
       # Find the max and min of both dataframes for printing linear regression line
       xmin <- min(c(min(df[[x_var]], na.rm = TRUE), min(df_reference[[x_var]], na.rm = TRUE)), na.rm = TRUE)
       xmax <- max(c(max(df[[x_var]], na.rm = TRUE), max(df_reference[[x_var]], na.rm = TRUE)), na.rm = TRUE)
@@ -1802,9 +1807,10 @@ main <- function() {
       
       message("Calculating reference strain interaction scores")
       reference_results <- calculate_interaction_scores(df_reference_interaction_stats, df_bg_stats, "reference")
+      df_reference_calculations <- reference_results$calculations
       df_reference_interactions_joined <- reference_results$full_data
       df_reference_interactions <- reference_results$interactions
-      write.csv(reference_results$calculations, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE)
+      write.csv(df_reference_calculations, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE)
       write.csv(df_reference_interactions, file = file.path(out_dir, "zscore_interactions_reference.csv"), row.names = FALSE)
 
       # message("Generating reference interaction plots")
@@ -1832,9 +1838,10 @@ main <- function() {
 
       message("Calculating deletion strain(s) interactions scores")
       deletion_results <- calculate_interaction_scores(df_deletion_stats, df_bg_stats, "deletion")
+      df_calculations <- deletion_results$calculations
       df_interactions <- deletion_results$interactions
       df_interactions_joined <- deletion_results$full_data
-      write.csv(deletion_results$calculations, file = file.path(out_dir, "zscore_calculations.csv"), row.names = FALSE)
+      write.csv(df_calculations, file = file.path(out_dir, "zscore_calculations.csv"), row.names = FALSE)
       write.csv(df_interactions, file = file.path(out_dir, "zscore_interactions.csv"), row.names = FALSE)
 
       # message("Generating deletion interaction plots")