diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index b217802a..a68fe8c6 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -717,7 +717,8 @@ generate_scatter_plot <- function(plot, config) {
   plot <- plot + geom_point(
     shape = ifelse(!is.null(config$shape), config$shape, 3),
     size = ifelse(!is.null(config$size), config$size, 1.5),
-    color = ifelse(!is.null(config$color), config$color, "black"),
+    # Do not uncomment, something is broken in ggplot aes() even if this makes logical sense
+    # color = ifelse(!is.null(config$color), config$color, "black"),
     position = position
   )
 
@@ -807,14 +808,6 @@ generate_scatter_plot <- function(plot, config) {
     # Extract necessary values
     intercept <- config$lm_line$intercept # required
     slope <- config$lm_line$slope # required
-
-    message("----- Generating lm_line -----")
-    message("Intercept: ", intercept)
-    message("Slope: ", slope)
-    message("Class: ", class(config$df[[config$x_var]]))
-    message("------------------------------")
-    
-
     xmin <- ifelse(!is.null(config$lm_line$xmin), config$lm_line$xmin, min(as.numeric(config$df[[config$x_var]])))
     xmax <- ifelse(!is.null(config$lm_line$xmax), config$lm_line$xmax, max(as.numeric(config$df[[config$x_var]])))
     color <- ifelse(!is.null(config$lm_line$color), config$lm_line$color, "blue")
@@ -1542,30 +1535,30 @@ main <- function() {
 
     plot_configs <- list(
       list(out_dir = out_dir_qc, filename = "L_vs_K_before_quality_control",
-        plot_configs = l_vs_k_plot_configs, page_width = 12, page_height = 8),
-      list(out_dir = out_dir_qc, filename = "frequency_delta_background",
-        plot_configs = frequency_delta_bg_plot_configs, page_width = 12, page_height = 8),
-      list(out_dir = out_dir_qc, filename = "L_vs_K_above_threshold",
-        plot_configs = above_threshold_plot_configs, page_width = 12, page_height = 8),
-      list(out_dir = out_dir_qc, filename = "plate_analysis",
-        plot_configs = plate_analysis_plot_configs, page_width = 14, page_height = 9),
-      list(out_dir = out_dir_qc, filename = "plate_analysis_boxplots",
-        plot_configs = plate_analysis_boxplot_configs, page_width = 18, page_height = 9),
-      list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros",
-        plot_configs = plate_analysis_no_zeros_plot_configs, page_width = 14, page_height = 9),
-      list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros_boxplots",
-        plot_configs = plate_analysis_no_zeros_boxplot_configs, page_width = 18, page_height = 9),
-      list(out_dir = out_dir_qc, filename = "L_vs_K_for_strains_2SD_outside_mean_K",
-        plot_configs = l_outside_2sd_k_plot_configs, page_width = 10, page_height = 8),
-      list(out_dir = out_dir_qc, filename = "delta_background_vs_K_for_strains_2SD_outside_mean_K",
-        plot_configs = delta_bg_outside_2sd_k_plot_configs, page_width = 10, page_height = 8)
+        plot_configs = l_vs_k_plot_configs, page_width = 12, page_height = 8)
+      # list(out_dir = out_dir_qc, filename = "frequency_delta_background",
+      #   plot_configs = frequency_delta_bg_plot_configs, page_width = 12, page_height = 8),
+      # list(out_dir = out_dir_qc, filename = "L_vs_K_above_threshold",
+      #   plot_configs = above_threshold_plot_configs, page_width = 12, page_height = 8),
+      # list(out_dir = out_dir_qc, filename = "plate_analysis",
+      #   plot_configs = plate_analysis_plot_configs, page_width = 14, page_height = 9),
+      # list(out_dir = out_dir_qc, filename = "plate_analysis_boxplots",
+      #   plot_configs = plate_analysis_boxplot_configs, page_width = 18, page_height = 9),
+      # list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros",
+      #   plot_configs = plate_analysis_no_zeros_plot_configs, page_width = 14, page_height = 9),
+      # list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros_boxplots",
+      #   plot_configs = plate_analysis_no_zeros_boxplot_configs, page_width = 18, page_height = 9),
+      # list(out_dir = out_dir_qc, filename = "L_vs_K_for_strains_2SD_outside_mean_K",
+      #   plot_configs = l_outside_2sd_k_plot_configs, page_width = 10, page_height = 8),
+      # list(out_dir = out_dir_qc, filename = "delta_background_vs_K_for_strains_2SD_outside_mean_K",
+      #   plot_configs = delta_bg_outside_2sd_k_plot_configs, page_width = 10, page_height = 8)
     )
 
     # Parallelize background and quality control plot generation
-    # furrr::future_map(plot_configs, function(config) {
-    #   generate_and_save_plots(config$out_dir, config$filename, config$plot_configs,
-    #     page_width = config$page_width, page_height = config$page_height)
-    # }, .options = furrr_options(seed = TRUE))
+    furrr::future_map(plot_configs, function(config) {
+      generate_and_save_plots(config$out_dir, config$filename, config$plot_configs,
+        page_width = config$page_width, page_height = config$page_height)
+    }, .options = furrr_options(seed = TRUE))
 
     # Loop over background strains
     # TODO currently only tested against one strain, if we want to do multiple strains we'll
@@ -1646,9 +1639,9 @@ main <- function() {
       write.csv(calculations_reference_csv, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE)
       write.csv(interactions_reference_csv, file = file.path(out_dir, "zscore_interactions_reference.csv"), row.names = FALSE)
 
-      # message("Generating reference interaction plots")
-      # reference_plot_configs <- generate_interaction_plot_configs(stats_reference, interactions_reference_joined, "reference")
-      # generate_and_save_plots(out_dir, "interaction_plots_reference", reference_plot_configs, page_width = 16, page_height = 16)
+      message("Generating reference interaction plots")
+      reference_plot_configs <- generate_interaction_plot_configs(stats_reference, interactions_reference_joined, "reference")
+      generate_and_save_plots(out_dir, "interaction_plots_reference", reference_plot_configs, page_width = 16, page_height = 16)
 
       message("Setting missing deletion values to the highest theoretical value at each drug conc for L")
       deletion_adjusted <- stats_na %>% # formerly X2
@@ -1679,9 +1672,9 @@ main <- function() {
       write.csv(calculations_csv, file = file.path(out_dir, "zscore_calculations.csv"), row.names = FALSE)
       write.csv(interactions_csv, file = file.path(out_dir, "zscore_interactions.csv"), row.names = FALSE)
 
-      # message("Generating deletion interaction plots")
-      # deletion_plot_configs <- generate_interaction_plot_configs(stats_reference, interactions_joined, "deletion")
-      # generate_and_save_plots(out_dir, "interaction_plots", deletion_plot_configs, page_width = 16, page_height = 16)
+      message("Generating deletion interaction plots")
+      deletion_plot_configs <- generate_interaction_plot_configs(stats_reference, interactions_joined, "deletion")
+      generate_and_save_plots(out_dir, "interaction_plots", deletion_plot_configs, page_width = 16, page_height = 16)
 
       message("Writing enhancer/suppressor csv files")
       interaction_threshold <- 2  # TODO add to study config?