diff --git a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
index 2cd527af..950ceb30 100644
--- a/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
+++ b/qhtcp-workflow/apps/r/calculate_interaction_zscores.R
@@ -806,7 +806,7 @@ generate_scatter_plot <- function(plot, config) {
       ) +
       geom_hline(
         yintercept = c(-config$sd_band, config$sd_band),
-        color = ifelse(!is.null(config$hl_color), config$hl_color, "gray")
+        color = ifelse(!is.null(config$hl_color), config$hl_color, "black")
       )
   }
 
@@ -1081,11 +1081,10 @@ generate_interaction_plot_configs <- function(df_summary, df_interactions, type)
     }
   }
 
-  # Group delta plots in chunks of 12
+  # Group delta plots in chunks of 12 per page
   chunk_size <- 12
   delta_plot_chunks <- split(delta_plot_configs, ceiling(seq_along(delta_plot_configs) / chunk_size))
 
-  # TODO, only return first page of plots for testing, remove this later
   return(c(
     list(list(grid_layout = list(ncol = 2), plots = stats_plot_configs)),
     list(list(grid_layout = list(ncol = 2), plots = stats_boxplot_configs)),
@@ -1109,11 +1108,11 @@ generate_rank_plot_configs <- function(df, is_lm = FALSE, adjust = FALSE, overla
     df[[paste0("Rank_lm_", variable)]] <- rank(df[[paste0("Z_lm_", variable)]], na.last = "keep")
   }
 
-  # Helper function to create a plot configuration
+  # Helper function to create a rank plot configuration
   create_plot_config <- function(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = TRUE) {
     num_enhancers <- sum(df[[zscore_var]] >= sd_band, na.rm = TRUE)
     num_suppressors <- sum(df[[zscore_var]] <= -sd_band, na.rm = TRUE)
-    
+
     # Default plot config
     plot_config <- list(
       df = df,
@@ -1121,7 +1120,7 @@ generate_rank_plot_configs <- function(df, is_lm = FALSE, adjust = FALSE, overla
       y_var = zscore_var,
       x_label = "Rank",
       plot_type = "scatter",
-      title = paste(y_label, "vs. Rank for", variable, "above", sd_band),
+      title = paste(y_label, "vs. Rank for", variable, "above", sd_band, "SD"),
       sd_band = sd_band,
       fill_positive = "#542788",
       fill_negative = "orange",
@@ -1134,8 +1133,8 @@ generate_rank_plot_configs <- function(df, is_lm = FALSE, adjust = FALSE, overla
       legend_position = "none"
     )
     
+    # Selectively add annotations
     if (with_annotations) {
-      # Add specific annotations for plots with annotations
       plot_config$annotations <- list(
         list(
           x = nrow(df) / 2,
@@ -1169,7 +1168,13 @@ generate_rank_plot_configs <- function(df, is_lm = FALSE, adjust = FALSE, overla
     }
   }
 
-  return(list(grid_layout = list(ncol = 3), plots = plot_configs))
+  # Group delta plots in chunks of 6 per page
+  chunk_size <- 6
+  plot_chunks <- split(plot_configs, ceiling(seq_along(plot_configs) / chunk_size))
+
+  return(c(
+    lapply(plot_chunks, function(chunk) list(grid_layout = list(ncol = 3), plots = chunk))
+  ))
 }
 
 generate_correlation_plot_configs <- function(df) {
@@ -1591,11 +1596,9 @@ main <- function() {
       write.csv(deletion_results$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)
 
-      print(df_interactions_joined, n = 20, width = 100000)
-
-      message("Generating deletion interaction plots")
-      deletion_plot_configs <- generate_interaction_plot_configs(df_reference_summary_stats, df_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(df_reference_summary_stats, df_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?
@@ -1631,7 +1634,7 @@ main <- function() {
 
       message("Generating rank plots")
       rank_plot_configs <- generate_rank_plot_configs(
-        df_interactions_joined,
+        df_interactions,
         is_lm = FALSE,
         adjust = TRUE
       )
@@ -1640,7 +1643,7 @@ main <- function() {
 
       message("Generating ranked linear model plots")
       rank_lm_plot_configs <- generate_rank_plot_configs(
-        df_interactions_joined,
+        df_interactions,
         is_lm = TRUE,
         adjust = TRUE
       )
@@ -1649,17 +1652,17 @@ main <- function() {
 
       message("Generating filtered ranked plots")
       rank_plot_filtered_configs <- generate_rank_plot_configs(
-        df_interactions_joined,
+        df_interactions,
         is_lm = FALSE,
         adjust = FALSE,
         overlap_color = TRUE
       )
-      generate_and_save_plots(out_dir, "RankPlots_na_rm", rank_plot_filtered_configs,
+      generate_and_save_plots(out_dir, "rank_plots_na_rm", rank_plot_filtered_configs,
         page_width = 18, page_height = 12)
 
       message("Generating filtered ranked linear model plots")
       rank_plot_lm_filtered_configs <- generate_rank_plot_configs(
-        df_interactions_joined,
+        df_interactions,
         is_lm = TRUE,
         adjust = FALSE,
         overlap_color = TRUE