diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index e09c199a..ef960e45 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -294,12 +294,16 @@ calculate_interaction_scores <- function(df, max_conc, bg_stats, group_vars, ove
   lm_means_sds <- calculations %>%
     group_by(across(all_of(group_vars))) %>%
     summarise(
+      mean_mean_L = mean(mean_L, na.rm = TRUE),
       mean_lm_L = mean(lm_Score_L, na.rm = TRUE),
       sd_lm_L = sd(lm_Score_L, na.rm = TRUE),
+      mean_mean_K = mean(mean_K, na.rm = TRUE),
       mean_lm_K = mean(lm_Score_K, na.rm = TRUE),
       sd_lm_K = sd(lm_Score_K, na.rm = TRUE),
+      mean_mean_r = mean(mean_r, na.rm = TRUE),
       mean_lm_r = mean(lm_Score_r, na.rm = TRUE),
       sd_lm_r = sd(lm_Score_r, na.rm = TRUE),
+      mean_mean_AUC = mean(mean_AUC, na.rm = TRUE),
       mean_lm_AUC = mean(lm_Score_AUC, na.rm = TRUE),
       sd_lm_AUC = sd(lm_Score_AUC, na.rm = TRUE)
     )
@@ -408,8 +412,10 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
     plot_configs  # Multiple groups
   }
 
-  for (group in plot_groups) {
+  # Open the PDF device once for all plots
+  pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 16, height = 9)
 
+  for (group in plot_groups) {
     static_plots <- list()
     plotly_plots <- list()
 
@@ -420,31 +426,26 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
       config <- plots[[i]]
       df <- config$df
 
-      if (config$plot_type == "bar") {
+      # Set up aes mapping based on plot type
+      aes_mapping <- if (config$plot_type == "bar" || config$plot_type == "density") {
         if (!is.null(config$color_var)) {
-          aes_mapping <- aes(x = .data[[config$x_var]], fill = .data[[config$color_var]], color = .data[[config$color_var]])
+          aes(x = .data[[config$x_var]], fill = .data[[config$color_var]], color = .data[[config$color_var]])
         } else {
-          aes_mapping <- aes(x = .data[[config$x_var]])
-        }
-      } else if (config$plot_type == "density") {
-        if (!is.null(config$color_var)) {
-          aes_mapping <- aes(x = .data[[config$x_var]], color = .data[[config$color_var]])
-        } else {
-          aes_mapping <- aes(x = .data[[config$x_var]])
+          aes(x = .data[[config$x_var]])
         }
       } else {
-        # For scatter and other plot types
         if (!is.null(config$y_var) && !is.null(config$color_var)) {
-          aes_mapping <- aes(x = .data[[config$x_var]], y = .data[[config$y_var]], color = .data[[config$color_var]])
+          aes(x = .data[[config$x_var]], y = .data[[config$y_var]], color = .data[[config$color_var]])
         } else if (!is.null(config$y_var)) {
-          aes_mapping <- aes(x = .data[[config$x_var]], y = .data[[config$y_var]])
+          aes(x = .data[[config$x_var]], y = .data[[config$y_var]])
         } else {
-          aes_mapping <- aes(x = .data[[config$x_var]])
+          aes(x = .data[[config$x_var]])
         }
       }
 
       plot <- ggplot(df, aes_mapping) + theme_publication(legend_position = config$legend_position)
 
+      # Add appropriate plot layer based on plot type
       plot <- switch(config$plot_type,
         "scatter" = generate_scatter_plot(plot, config),
         "box" = generate_boxplot(plot, config),
@@ -453,6 +454,7 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
         plot  # default (unused)
       )
 
+      # Add labels and title
       if (!is.null(config$title)) plot <- plot + ggtitle(config$title)
       if (!is.null(config$x_label)) plot <- plot + xlab(config$x_label)
       if (!is.null(config$y_label)) plot <- plot + ylab(config$y_label)
@@ -460,56 +462,66 @@ generate_and_save_plots <- function(out_dir, filename, plot_configs) {
 
       # Add error bars if specified
       if (!is.null(config$error_bar) && config$error_bar) {
-        error_bar_color <- if (!is.null(config$error_bar_params$color)) {
-          config$error_bar_params$color
-        } else {
-          "red"
+        error_bar_color <- config$error_bar_params$color %||% "red"
+        y_mean_col <- paste0("mean_", config$y_var)
+        y_sd_col <- paste0("sd_", config$y_var)
+
+        if (!is.null(config$error_bar_params$center_point)) {
+          plot <- plot + geom_point(aes(
+            x = .data[[config$x_var]],
+            y = .data[[y_mean_col]]),
+            color = error_bar_color,
+            shape = 16)
         }
 
         # Use error_bar_params if provided, otherwise calculate from mean and sd
         if (!is.null(config$error_bar_params$ymin) && !is.null(config$error_bar_params$ymax)) {
           plot <- plot + geom_errorbar(aes(
             ymin = config$error_bar_params$ymin,
-            ymax = config$error_bar_params$ymax,
-            color = error_bar_color))
+            ymax = config$error_bar_params$ymax),
+            color = error_bar_color)
         } else {
-          y_mean_col <- paste0("mean_", config$y_var)
-          y_sd_col <- paste0("sd_", config$y_var)
           plot <- plot + geom_errorbar(aes(
             ymin = .data[[y_mean_col]] - .data[[y_sd_col]],
-            ymax = .data[[y_mean_col]] + .data[[y_sd_col]],
-            color = error_bar_color))
+            ymax = .data[[y_mean_col]] + .data[[y_sd_col]]),
+            color = error_bar_color)
         }
       }
 
+      # Convert ggplot to plotly for interactive version
       plotly_plot <- suppressWarnings(plotly::ggplotly(plot))
+
+      # Store both static and interactive versions
       static_plots[[i]] <- plot
       plotly_plots[[i]] <- plotly_plot
     }
 
-    pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 16, height = 9)
-
+    # Print the plots to the PDF (one page per plot or in a grid)
     if (is.null(grid_layout)) {
+      # Print each plot individually on separate pages if no grid layout is specified
       for (plot in static_plots) {
         print(plot)
       }
     } else {
+      # Arrange plots in grid layout on a single page
       grid.arrange(
         grobs = static_plots,
         ncol = grid_layout$ncol,
         nrow = grid_layout$nrow
       )
     }
-
-    dev.off()
-
-    # out_html_file <- file.path(out_dir, paste0(filename, ".html"))
-    # message("Saving combined HTML file: ", out_html_file)
-    # htmltools::save_html(
-    #   htmltools::tagList(plotly_plots),
-    #   file = out_html_file
-    # )
   }
+
+  # Close the PDF device after all plots are done
+  dev.off()
+
+  # Optional: Uncomment and save the interactive HTML version if needed
+  # out_html_file <- file.path(out_dir, paste0(filename, ".html"))
+  # message("Saving combined HTML file: ", out_html_file)
+  # htmltools::save_html(
+  #   htmltools::tagList(plotly_plots),
+  #   file = out_html_file
+  # )
 }
 
 generate_scatter_plot <- function(plot, config) {
@@ -519,7 +531,7 @@ generate_scatter_plot <- function(plot, config) {
   size <- if (!is.null(config$size)) config$size else 1.5
   position <-
     if (!is.null(config$position) && config$position == "jitter") {
-      position_jitter(width = 0.2, height = 0)
+      position_jitter(width = 0.3, height = 0.1)
     } else {
       "identity"
     }
@@ -727,47 +739,66 @@ generate_interaction_plot_configs <- function(df, type) {
     r = c(-0.6, 0.6),
     AUC = c(-6000, 6000)
   )
-
-  overall_plot_configs <- list()
+  
+  stats_plot_configs <- list()
+  stats_boxplot_configs <- list()
   delta_plot_configs <- list()
 
   # Overall statistics plots
   OrfRep <- first(df$OrfRep) # this should correspond to the reference strain
 
-  for (var in names(limits_map)) {
-    y_limits <- limits_map[[var]]
+  for (plot_type in c("scatter", "box")) {
 
-    # Use the pre-calculated lm intercept and slope from the dataframe
-    lm_intercept_col <- paste0("lm_intercept_", var)
-    lm_slope_col <- paste0("lm_slope_", var)
+    for (var in names(limits_map)) {
+      y_limits <- limits_map[[var]]
 
-    # Ensure no NA or invalid values in lm_line calculations
-    intercept_value <- mean(df[[lm_intercept_col]], na.rm = TRUE)
-    slope_value <- mean(df[[lm_slope_col]], na.rm = TRUE)
+      # Use the pre-calculated lm intercept and slope from the dataframe
+      lm_intercept_col <- paste0("lm_intercept_", var)
+      lm_slope_col <- paste0("lm_slope_", var)
 
-    plot_config <- list(
-      df = df,
-      plot_type = "scatter",
-      x_var = "conc_num_factor_factor",
-      y_var = var,
-      shape = 16,
-      x_label = unique(df$Drug)[1],
-      title = sprintf("%s Scatter RF for %s with SD", OrfRep, var),
-      coord_cartesian = y_limits,
-      error_bar = TRUE,
-      error_bar_params = list(
-        color = "red"
-      ),
-      x_breaks = unique(df$conc_num_factor_factor),
-      x_labels = as.character(unique(df$conc_num)),
-      position = "jitter",
-      smooth = TRUE,
-      lm_line = list(
-        intercept = intercept_value,
-        slope = slope_value
+      # Ensure no NA or invalid values in lm_line calculations
+      intercept_value <- mean(df[[lm_intercept_col]], na.rm = TRUE)
+      slope_value <- mean(df[[lm_slope_col]], na.rm = TRUE)
+
+      # Common plot configuration
+      plot_config <- list(
+        df = df,
+        x_var = "conc_num_factor_factor",
+        y_var = var,
+        shape = 16,
+        x_label = unique(df$Drug)[1],
+        coord_cartesian = y_limits,
+        x_breaks = unique(df$conc_num_factor_factor),
+        x_labels = as.character(unique(df$conc_num)),
+        lm_line = list(
+          intercept = intercept_value,
+          slope = slope_value
+        )
       )
-    )
-    overall_plot_configs <- append(overall_plot_configs, list(plot_config))
+
+      # Add specific configurations for scatter and box plots
+      if (plot_type == "scatter") {
+        plot_config$plot_type <- "scatter"
+        plot_config$title <- sprintf("%s Scatter RF for %s with SD", OrfRep, var)
+        plot_config$error_bar = TRUE
+        plot_config$error_bar_params <- list(
+          color = "red",
+          center_point = TRUE
+        )
+        plot_config$position <- "jitter"
+        
+        # Append to scatter plot configurations
+        stats_plot_configs <- append(stats_plot_configs, list(plot_config))
+
+      } else if (plot_type == "box") {
+        plot_config$plot_type <- "box"
+        plot_config$title <- sprintf("%s Boxplot RF for %s with SD", OrfRep, var)
+        plot_config$position <- "dodge"  # Boxplots don't need jitter, use dodge instead
+
+        # Append to boxplot configurations
+        stats_boxplot_configs <- append(stats_boxplot_configs, list(plot_config))
+      }
+    }
   }
 
   # Delta interaction plots
@@ -850,7 +881,8 @@ generate_interaction_plot_configs <- function(df, type) {
   grid_nrow <- ceiling(num_plots / grid_ncol)
 
   return(list(
-    list(grid_layout = list(ncol = 2, nrow = 2), plots = overall_plot_configs),
+    list(grid_layout = list(ncol = 2, nrow = 2), plots = stats_plot_configs),
+    list(grid_layout = list(ncol = 2, nrow = 2), plots = stats_boxplot_configs),
     list(grid_layout = list(ncol = 4, nrow = grid_nrow), plots = delta_plot_configs)
   ))
 }