server-scripts/qhtcp-workflow/apps/r/calculate_interaction_zscores.R

suppressMessages({
  library(ggplot2)
  library(plotly)
  library(htmlwidgets)
  library(dplyr)
  library(ggthemes)
  library(data.table)
  library(unix)
})

options(warn = 2)
options(width = 10000)

# Set the memory limit to 30GB (30 * 1024 * 1024 * 1024 bytes)
soft_limit <- 30 * 1024 * 1024 * 1024
hard_limit <- 30 * 1024 * 1024 * 1024
rlimit_as(soft_limit, hard_limit)

# Constants for configuration
plot_width <- 14
plot_height <- 9
base_size <- 14

parse_arguments <- function() {
  args <- if (interactive()) {
    c(
      "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240116_jhartman2_DoxoHLD",
      "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/apps/r/SGD_features.tab",
      "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/easy/20240116_jhartman2_DoxoHLD/results_std.txt",
      "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240822_jhartman2_DoxoHLD/exp1",
      "Experiment 1: Doxo versus HLD",
      3,
      "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240822_jhartman2_DoxoHLD/exp2",
      "Experiment 2: HLD versus Doxo",
      3
    )
  } else {
    commandArgs(trailingOnly = TRUE)
  }
  
  # Extract paths, names, and standard deviations
  paths <- args[seq(4, length(args), by = 3)]
  names <- args[seq(5, length(args), by = 3)]
  sds <- as.numeric(args[seq(6, length(args), by = 3)])
  
  # Normalize paths
  normalized_paths <- normalizePath(paths, mustWork = FALSE)
  
  # Create named list of experiments
  experiments <- list()
  for (i in seq_along(paths)) {
    experiments[[names[i]]] <- list(
      path = normalized_paths[i],
      sd = sds[i]
    )
  }
  
  list(
    out_dir = normalizePath(args[1], mustWork = FALSE),
    sgd_gene_list = normalizePath(args[2], mustWork = FALSE),
    easy_results_file = normalizePath(args[3], mustWork = FALSE),
    experiments = experiments
  )
}

args <- parse_arguments()

# Should we keep output in exp dirs or combine in the study output dir?
# dir.create(file.path(args$out_dir, "zscores"), showWarnings = FALSE)
# dir.create(file.path(args$out_dir, "zscores", "qc"), showWarnings = FALSE)

# Define themes and scales
theme_publication <- function(base_size = 14, base_family = "sans", legend_position = "bottom") {
  theme_foundation <- ggplot2::theme_grey(base_size = base_size, base_family = base_family)
  
  theme_foundation %+replace%
    theme(
      plot.title = element_text(face = "bold", size = rel(1.2), hjust = 0.5),
      text = element_text(),
      panel.background = element_rect(colour = NA),
      plot.background = element_rect(colour = NA),
      panel.border = element_rect(colour = NA),
      axis.title = element_text(face = "bold", size = rel(1)),
      axis.title.y = element_text(angle = 90, vjust = 2),
      axis.title.x = element_text(vjust = -0.2),
      axis.line = element_line(colour = "black"),
      panel.grid.major = element_line(colour = "#f0f0f0"),
      panel.grid.minor = element_blank(),
      legend.key = element_rect(colour = NA),
      legend.position = legend_position,
      legend.direction = ifelse(legend_position == "right", "vertical", "horizontal"),
      plot.margin = unit(c(10, 5, 5, 5), "mm"),
      strip.background = element_rect(colour = "#f0f0f0", fill = "#f0f0f0"),
      strip.text = element_text(face = "bold")
    )
}

scale_fill_publication <- function(...) {
  discrete_scale("fill", "Publication", manual_pal(values = c(
    "#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506",
    "#a6cee3", "#fb9a99", "#984ea3", "#ffff33"
  )), ...)
}

scale_colour_publication <- function(...) {
  discrete_scale("colour", "Publication", manual_pal(values = c(
    "#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506",
    "#a6cee3", "#fb9a99", "#984ea3", "#ffff33"
  )), ...)
}

# Load the initial dataframe from the easy_results_file
load_and_process_data <- function(easy_results_file, sd = 3) {
  df <- read.delim(easy_results_file, skip = 2, as.is = TRUE, row.names = 1, strip.white = TRUE)
  
  df <- df %>%
    filter(!(.[[1]] %in% c("", "Scan"))) %>%
    filter(!is.na(ORF) & ORF != "" & !Gene %in% c("BLANK", "Blank", "blank") & Drug != "BMH21") %>%
    # Rename columns
    rename(L = l, num = Num., AUC = AUC96, scan = Scan, last_bg = LstBackgrd, first_bg = X1stBackgrd) %>%
    mutate(
      across(c(Col, Row, num, L, K, r, scan, AUC, last_bg, first_bg), as.numeric),
      delta_bg = last_bg - first_bg,
      delta_bg_tolerance = mean(delta_bg, na.rm = TRUE) + (sd * sd(delta_bg, na.rm = TRUE)),
      NG = if_else(L == 0 & !is.na(L), 1, 0),
      DB = if_else(delta_bg >= delta_bg_tolerance, 1, 0),
      SM = 0,
      OrfRep = if_else(ORF == "YDL227C", "YDL227C", OrfRep), # should these be hardcoded?
      conc_num = as.numeric(gsub("[^0-9\\.]", "", Conc)),
      conc_num_factor = as.numeric(as.factor(conc_num)) - 1
    )
  
  return(df)
}

# Update Gene names using the SGD gene list
update_gene_names <- function(df, sgd_gene_list) {
  # Load SGD gene list
  genes <- read.delim(file = sgd_gene_list,
    quote = "", header = FALSE,
    colClasses = c(rep("NULL", 3), rep("character", 2), rep("NULL", 11)))
  
  # Create a named vector for mapping ORF to GeneName
  gene_map <- setNames(genes$V5, genes$V4)
  # Vectorized match to find the GeneName from gene_map
  mapped_genes <- gene_map[df$ORF]
  # Replace NAs in mapped_genes with original Gene names (preserves existing Gene names if ORF is not found)
  updated_genes <- ifelse(is.na(mapped_genes) | df$OrfRep == "YDL227C", df$Gene, mapped_genes)
  # Ensure Gene is not left blank or incorrectly updated to "OCT1"
  df <- df %>%
    mutate(Gene = ifelse(updated_genes == "" | updated_genes == "OCT1", OrfRep, updated_genes))
  
  return(df)
}

# Calculate summary statistics for all variables
calculate_summary_stats <- function(df, variables, group_vars = c("OrfRep", "conc_num", "conc_num_factor")) {

  # Summarize the variables within the grouped data
  summary_stats <- df %>%
    group_by(across(all_of(group_vars))) %>%
    summarise(
      N = sum(!is.na(L)),
      across(all_of(variables), list(
        mean = ~mean(., na.rm = TRUE),
        median = ~median(., na.rm = TRUE),
        max = ~ ifelse(all(is.na(.)), NA, max(., na.rm = TRUE)),
        min = ~ ifelse(all(is.na(.)), NA, min(., na.rm = TRUE)),
        sd = ~sd(., na.rm = TRUE),
        se = ~ ifelse(all(is.na(.)), NA, sd(., na.rm = TRUE) / sqrt(sum(!is.na(.)) - 1))
      ), .names = "{.fn}_{.col}")
    )

  # print(summary_stats)

  # Prevent .x and .y suffix issues by renaming columns
  df_cleaned <- df %>%
    select(-any_of(setdiff(names(summary_stats), group_vars))) # Avoid duplicate columns in the final join
  
  # Join the stats back to the original dataframe
  df_with_stats <- left_join(df_cleaned, summary_stats, by = group_vars)

  return(list(summary_stats = summary_stats, df_with_stats = df_with_stats))
}

calculate_interaction_scores <- function(df, max_conc, variables, group_vars = c("OrfRep", "Gene", "num")) {

  # Calculate total concentration variables
  total_conc_num <- length(unique(df$conc_num))

  # Pull the background means and standard deviations from zero concentration
  bg_means <- list(
    L = df %>% filter(conc_num_factor == 0) %>% pull(mean_L) %>% first(),
    K = df %>% filter(conc_num_factor == 0) %>% pull(mean_K) %>% first(),
    r = df %>% filter(conc_num_factor == 0) %>% pull(mean_r) %>% first(),
    AUC = df %>% filter(conc_num_factor == 0) %>% pull(mean_AUC) %>% first()
  )

  bg_sd <- list(
    L = df %>% filter(conc_num_factor == 0) %>% pull(sd_L) %>% first(),
    K = df %>% filter(conc_num_factor == 0) %>% pull(sd_K) %>% first(),
    r = df %>% filter(conc_num_factor == 0) %>% pull(sd_r) %>% first(),
    AUC = df %>% filter(conc_num_factor == 0) %>% pull(sd_AUC) %>% first()
  )

  stats <- df %>%
    group_by(OrfRep, Gene, num, conc_num, conc_num_factor) %>%
    summarise(
      N = sum(!is.na(L)),
      NG = sum(NG, na.rm = TRUE),
      DB = sum(DB, na.rm = TRUE),
      SM = sum(SM, na.rm = TRUE),
      across(all_of(variables), list(
        mean = ~mean(., na.rm = TRUE),
        median = ~median(., na.rm = TRUE),
        max = ~ifelse(all(is.na(.)), NA, max(., na.rm = TRUE)),
        min = ~ifelse(all(is.na(.)), NA, min(., na.rm = TRUE)),
        sd = ~sd(., na.rm = TRUE),
        se = ~ifelse(sum(!is.na(.)) > 1, sd(., na.rm = TRUE) / sqrt(sum(!is.na(.)) - 1), NA)
      ), .names = "{.fn}_{.col}")
    )

  stats <- df %>%
    group_by(OrfRep, Gene, num) %>%
    mutate(
      WT_L = mean_L,
      WT_K = mean_K,
      WT_r = mean_r,
      WT_AUC = mean_AUC,
      WT_sd_L = sd_L,
      WT_sd_K = sd_K,
      WT_sd_r = sd_r,
      WT_sd_AUC = sd_AUC
    )

  stats <- stats %>%
    mutate(
      Raw_Shift_L = first(mean_L) - bg_means$L,
      Raw_Shift_K = first(mean_K) - bg_means$K,
      Raw_Shift_r = first(mean_r) - bg_means$r,
      Raw_Shift_AUC = first(mean_AUC) - bg_means$AUC,
      Z_Shift_L = first(Raw_Shift_L) / bg_sd$L,
      Z_Shift_K = first(Raw_Shift_K) / bg_sd$K,
      Z_Shift_r = first(Raw_Shift_r) / bg_sd$r,
      Z_Shift_AUC = first(Raw_Shift_AUC) / bg_sd$AUC
    )

  stats <- stats %>%
    mutate(
      Exp_L = WT_L + Raw_Shift_L,
      Exp_K = WT_K + Raw_Shift_K,
      Exp_r = WT_r + Raw_Shift_r,
      Exp_AUC = WT_AUC + Raw_Shift_AUC,
      Delta_L = mean_L - Exp_L,
      Delta_K = mean_K - Exp_K,
      Delta_r = mean_r - Exp_r,
      Delta_AUC = mean_AUC - Exp_AUC
    )

  stats <- stats %>%
    mutate(
      Delta_L = if_else(NG == 1, mean_L - WT_L, Delta_L),
      Delta_K = if_else(NG == 1, mean_K - WT_K, Delta_K),
      Delta_r = if_else(NG == 1, mean_r - WT_r, Delta_r),
      Delta_AUC = if_else(NG == 1, mean_AUC - WT_AUC, Delta_AUC),
      Delta_L = if_else(SM == 1, mean_L - WT_L, Delta_L)
    )

  stats <- stats %>%
    mutate(
      Zscore_L = Delta_L / WT_sd_L,
      Zscore_K = Delta_K / WT_sd_K,
      Zscore_r = Delta_r / WT_sd_r,
      Zscore_AUC = Delta_AUC / WT_sd_AUC
    )

  # Calculate linear models
  lm_L <- lm(Delta_L ~ conc_num_factor, data = stats)
  lm_K <- lm(Delta_K ~ conc_num_factor, data = stats)
  lm_r <- lm(Delta_r ~ conc_num_factor, data = stats)
  lm_AUC <- lm(Delta_AUC ~ conc_num_factor, data = stats)

  interactions <- stats %>%
    group_by(OrfRep, Gene, num) %>%
    summarise(
      OrfRep = first(OrfRep),
      Gene = first(Gene),
      num = first(num),
      conc_num = first(conc_num),
      conc_num_factor = first(conc_num_factor),
      Raw_Shift_L = first(Raw_Shift_L),
      Raw_Shift_K = first(Raw_Shift_K),
      Raw_Shift_r = first(Raw_Shift_r),
      Raw_Shift_AUC = first(Raw_Shift_AUC),
      Z_Shift_L = first(Z_Shift_L),
      Z_Shift_K = first(Z_Shift_K),
      Z_Shift_r = first(Z_Shift_r),
      Z_Shift_AUC = first(Z_Shift_AUC),
      Sum_Zscore_L = sum(Zscore_L, na.rm = TRUE),
      Sum_Zscore_K = sum(Zscore_K, na.rm = TRUE),
      Sum_Zscore_r = sum(Zscore_r, na.rm = TRUE),
      Sum_Zscore_AUC = sum(Zscore_AUC, na.rm = TRUE),
      lm_Score_L = max_conc * coef(lm_L)[2] + coef(lm_L)[1],
      lm_Score_K = max_conc * coef(lm_K)[2] + coef(lm_K)[1],
      lm_Score_r = max_conc * coef(lm_r)[2] + coef(lm_r)[1],
      lm_Score_AUC = max_conc * coef(lm_AUC)[2] + coef(lm_AUC)[1],
      R_Squared_L = summary(lm_L)$r.squared,
      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)
    )

  num_non_removed_concs <- total_conc_num - sum(stats$DB, na.rm = TRUE) - 1
  
  interactions <- interactions %>%
    mutate(
      Avg_Zscore_L = Sum_Zscore_L / num_non_removed_concs,
      Avg_Zscore_K = Sum_Zscore_K / num_non_removed_concs,
      Avg_Zscore_r = Sum_Zscore_r / (total_conc_num - 1),
      Avg_Zscore_AUC = Sum_Zscore_AUC / (total_conc_num - 1),
      Z_lm_L = (lm_Score_L - mean(lm_Score_L, na.rm = TRUE)) / sd(lm_Score_L, na.rm = TRUE),
      Z_lm_K = (lm_Score_K - mean(lm_Score_K, na.rm = TRUE)) / sd(lm_Score_K, na.rm = TRUE),
      Z_lm_r = (lm_Score_r - mean(lm_Score_r, na.rm = TRUE)) / sd(lm_Score_r, na.rm = TRUE),
      Z_lm_AUC = (lm_Score_AUC - mean(lm_Score_AUC, na.rm = TRUE)) / sd(lm_Score_AUC, na.rm = TRUE)
    ) %>%
    arrange(desc(Z_lm_L)) %>%
    arrange(desc(NG))

  # Declare column order for output
  calculations <- stats %>%
    select(
      "OrfRep", "Gene", "conc_num", "conc_num_factor", "N",
      "mean_L", "mean_K", "mean_r", "mean_AUC",
      "median_L", "median_K", "median_r", "median_AUC",
      "sd_L", "sd_K", "sd_r", "sd_AUC",
      "se_L", "se_K", "se_r", "se_AUC",
      "Raw_Shift_L", "Raw_Shift_K", "Raw_Shift_r", "Raw_Shift_AUC",
      "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC",
      "WT_L", "WT_K", "WT_r", "WT_AUC",
      "WT_sd_L", "WT_sd_K", "WT_sd_r", "WT_sd_AUC",
      "Exp_L", "Exp_K", "Exp_r", "Exp_AUC",
      "Delta_L", "Delta_K", "Delta_r", "Delta_AUC",
      "Zscore_L", "Zscore_K", "Zscore_r", "Zscore_AUC",
      "NG", "SM", "DB")
    
  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", "conc_num", "conc_num_factor"))))
  interactions_joined <- left_join(interactions_joined, interactions, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor"))

  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) {
  message("Generating html and pdf plots for: ", file_name)

  # Prepare lists to collect plots
  static_plots <- list()
  plotly_plots <- list()

  for (i in seq_along(plot_configs)) {
    config <- plot_configs[[i]]
    df <- config$df

    # Build the aes_mapping based on config
    aes_mapping <- if (is.null(config$color_var)) {
      if (is.null(config$y_var)) {
        aes(x = .data[[config$x_var]])
      } else {
        aes(x = .data[[config$x_var]], y = .data[[config$y_var]])
      }
    } else {
      if (is.null(config$y_var)) {
        aes(x = .data[[config$x_var]], color = as.factor(.data[[config$color_var]]))
      } else {
        aes(x = .data[[config$x_var]], y = .data[[config$y_var]], color = as.factor(.data[[config$color_var]]))
      }
    }

    # Start building the plot with aes_mapping
    plot_base <- ggplot(df, aes_mapping)

    # Function to generate the plot
    generate_plot <- function(is_interactive) {
      # Use appropriate helper function based on plot type
      plot <- switch(config$plot_type,
        "scatter" = generate_scatter_plot(plot_base, config, is_interactive = is_interactive),
        "box" = generate_box_plot(plot_base, config),
        "density" = plot_base + geom_density(),
        "bar" = plot_base + geom_bar(),
        plot_base  # default case if no type matches
      )

      # Apply additional settings if provided
      if (!is.null(config$legend_position)) {
        plot <- plot + theme(legend.position = config$legend_position)
      }

      # Add title and labels if provided
      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)
      }

      # Return the plot
      plot
    }

    # Generate the static plot
    static_plot <- generate_plot(is_interactive = FALSE)

    # Generate the interactive plot
    interactive_plot <- generate_plot(is_interactive = TRUE)

    # Convert to plotly object
    plotly_plot <- ggplotly(interactive_plot, tooltip = "text")
    if (!is.null(config$legend_position) && config$legend_position == "bottom") {
      plotly_plot <- plotly_plot %>% layout(legend = list(orientation = "h"))
    }

    # Add plots to lists
    static_plots[[i]] <- static_plot
    plotly_plots[[i]] <- plotly_plot
  }

  # PDF saving logic
  pdf(file.path(output_dir, paste0(file_name, ".pdf")), width = 14, height = 9)
  lapply(static_plots, print)
  dev.off()

  # Combine and save interactive plots
  combined_plot <- subplot(plotly_plots, nrows = grid_layout$nrow %||% length(plotly_plots), margin = 0.05)
  saveWidget(combined_plot, file = file.path(output_dir, paste0(file_name, ".html")), selfcontained = TRUE)
}

generate_scatter_plot <- function(plot, config, is_interactive = FALSE) {
  # Check for missing or out-of-range data
  missing_data <- config$df %>%
    filter(
      is.na(!!sym(config$x_var)) | is.na(!!sym(config$y_var)) |
        !!sym(config$y_var) < min(config$ylim_vals, na.rm = TRUE) |
        !!sym(config$y_var) > max(config$ylim_vals, na.rm = TRUE)
    )

  # Print the rows with missing or out-of-range data if any
  if (nrow(missing_data) > 0) {
    message("Missing or out-of-range data for ", config$title, ":")
    print(
      missing_data %>% select(any_of(
        c(
          "OrfRep",
          "Gene",
          "num",
          "conc_num",
          "conc_num_factor",
          config$x_var,
          config$y_var
        )
      )),
      n = 100
    )
  }

  # Add the interactive `text` aesthetic if `is_interactive` is TRUE
  if (is_interactive) {
    if (!is.null(config$delta_bg_point) && config$delta_bg_point) {
      plot <- plot + geom_point(
        aes(text = paste("ORF:", OrfRep, "Gene:", Gene, "delta_bg:", delta_bg)),
        shape = config$shape %||% 3,
        size = config$size %||% 0.2
      )
    } else if (!is.null(config$gene_point) && config$gene_point) {
      plot <- plot + geom_point(
        aes(text = paste("ORF:", OrfRep, "Gene:", Gene)),
        shape = config$shape %||% 3,
        size = config$size %||% 0.2,
        position = "jitter"
      )
    } else {
      plot <- plot + geom_point(
        aes(text = paste("ORF:", OrfRep, "Gene:", Gene)),
        shape = config$shape %||% 3,
        size = config$size %||% 0.2
      )
    }
  } else {
    # For non-interactive plots, just add `geom_point` without `text` aesthetic
    plot <- plot + geom_point(
      shape = config$shape %||% 3,
      size = config$size %||% 0.2,
      position = if (!is.null(config$position) && config$position == "jitter") "jitter" else "identity"
    )
  }

  # Add smooth line if specified
  if (!is.null(config$add_smooth) && config$add_smooth) {
    plot <- if (!is.null(config$lm_line)) {
      plot + geom_abline(intercept = config$lm_line$intercept, slope = config$lm_line$slope)
    } else {
      plot + geom_smooth(method = "lm", se = FALSE)
    }
  }

  # Add SD bands (iterate over sd_band only here)
  if (!is.null(config$sd_band)) {
    for (i in config$sd_band) {
      plot <- plot +
        annotate("rect", xmin = -Inf, xmax = Inf, ymin = i, ymax = Inf, fill = "#542788", alpha = 0.3) +
        annotate("rect", xmin = -Inf, xmax = Inf, ymin = -i, ymax = -Inf, fill = "orange", alpha = 0.3) +
        geom_hline(yintercept = c(-i, i), color = "gray")
    }
  }

  # Add error bars if specified
  if (!is.null(config$error_bar) && config$error_bar) {
    y_mean_col <- paste0("mean_", config$y_var)
    y_sd_col <- paste0("sd_", config$y_var)
    plot <- plot + geom_errorbar(
      aes(
        ymin = !!sym(y_mean_col) - !!sym(y_sd_col),
        ymax = !!sym(y_mean_col) + !!sym(y_sd_col)
      ),
      alpha = 0.3
    )
  }

  # Add x-axis customization if specified
  if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) {
    plot <- plot + scale_x_continuous(
      name = config$x_label,
      breaks = config$x_breaks,
      labels = config$x_labels
    )
  }

  # Use coord_cartesian for zooming in without removing data outside the range
  if (!is.null(config$coord_cartesian)) {
    plot <- plot + coord_cartesian(ylim = config$coord_cartesian)
  }

  # Use scale_y_continuous for setting the y-axis limits
  if (!is.null(config$ylim_vals)) {
    plot <- plot + scale_y_continuous(limits = config$ylim_vals)
  }

  # Add annotations if provided
  if (!is.null(config$annotations)) {
    for (annotation in config$annotations) {
      plot <- plot + annotate("text", x = annotation$x, y = annotation$y, label = annotation$label)
    }
  }

  # Add titles and themes if specified
  if (!is.null(config$title)) {
    plot <- plot + ggtitle(config$title)
  }

  if (!is.null(config$legend_position)) {
    plot <- plot + theme(legend.position = config$legend_position)
  }

  return(plot)
}

generate_box_plot <- function(plot, config) {
  plot <- plot + geom_boxplot()
  
  if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) {
    plot <- plot + scale_x_discrete(
      name = config$x_label,
      breaks = config$x_breaks,
      labels = config$x_labels
    )
  }

  if (!is.null(config$coord_cartesian)) {
    plot <- plot + coord_cartesian(ylim = config$coord_cartesian)
  }

  return(plot)
}

# Adjust missing values and calculate ranks
adjust_missing_and_rank <- function(df, variables) {
  
  # Loop over each variable
  for (var in variables) {
    # Construct column names
    avg_zscore_col <- paste0("Avg_Zscore_", var)
    z_lm_col <- paste0("Z_lm_", var)
    rank_col <- paste0("Rank_", var)
    rank_lm_col <- paste0("Rank_lm_", var)
    
    # Check if the columns exist in the data frame
    if (all(c(avg_zscore_col, z_lm_col) %in% names(df))) {
      # Adjust missing values by replacing NA with 0.001
      df[[avg_zscore_col]] <- if_else(is.na(df[[avg_zscore_col]]), 0.001, df[[avg_zscore_col]])
      df[[z_lm_col]] <- if_else(is.na(df[[z_lm_col]]), 0.001, df[[z_lm_col]])
      
      # Compute ranks and create new columns
      df[[rank_col]] <- rank(df[[avg_zscore_col]], na.last = "keep")
      df[[rank_lm_col]] <- rank(df[[z_lm_col]], na.last = "keep")
    } else {
      warning(paste("Columns", avg_zscore_col, "or", z_lm_col, "not found in data frame"))
    }
  }
  
  return(df)
}

generate_interaction_plot_configs <- function(df, variables) {
  configs <- list()

  # Define common y-limits and other attributes for each variable dynamically
  limits_map <- list(L = c(-65, 65), K = c(-65, 65), r = c(-0.65, 0.65), AUC = c(-6500, 6500))

  # Define annotation positions based on the variable being plotted
  annotation_positions <- list(
    L = list(Z_Shift_L = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
    K = list(Z_Shift_K = 45, lm_ZScore = 25, NG = -25, DB = -35, SM = -45),
    r = list(Z_Shift_r = 0.45, lm_ZScore = 0.25, NG = -0.25, DB = -0.35, SM = -0.45),
    AUC = list(Z_Shift_AUC = 4500, lm_ZScore = 2500, NG = -2500, DB = -3500, SM = -4500)
  )

  # Define which annotations to include for each plot
  annotation_labels <- list(
    ZShift = function(df, var) {
      val <- df[[paste0("Z_Shift_", var)]]
      if (is.numeric(val)) {
        paste("ZShift =", round(val, 2))
      } else {
        paste("ZShift =", val)
      }
    },
    lm_ZScore = function(df, var) {
      val <- df[[paste0("Z_lm_", var)]]
      if (is.numeric(val)) {
        paste("lm ZScore =", round(val, 2))
      } else {
        paste("lm ZScore =", val)
      }
    },
    NG = function(df, var) paste("NG =", df$NG),
    DB = function(df, var) paste("DB =", df$DB),
    SM = function(df, var) paste("SM =", df$SM)
  )

  for (variable in variables) {
    # Dynamically generate the names of the columns
    var_info <- list(
      ylim = limits_map[[variable]],
      sd_col = paste0("WT_sd_", variable)
    )

    # Extract the precomputed linear model coefficients
    lm_line <- list(
      intercept = df[[paste0("lm_intercept_", variable)]],
      slope = df[[paste0("lm_slope_", variable)]]
    )

    annotations <- lapply(names(annotation_positions[[variable]]), function(annotation_name) {
      message("Processing annotation: ", annotation_name, " for variable: ", variable)
      y_pos <- annotation_positions[[variable]][[annotation_name]]
      
      # Check if the annotation_name exists in annotation_labels
      if (!is.null(annotation_labels[[annotation_name]])) {
        label <- annotation_labels[[annotation_name]](df, variable)
        list(x = 1, y = y_pos, label = label)
      } else {
        message(paste("Warning: No annotation function found for", annotation_name))
        NULL
      }
    })

    # Filter out any NULL annotations
    annotations <- Filter(Negate(is.null), annotations)

    # Add scatter plot configuration for this variable
    configs[[length(configs) + 1]] <- list(
      df = df,
      x_var = "conc_num_factor",
      y_var = variable,
      plot_type = "scatter",
      title = sprintf("%s      %s", df$OrfRep[1], df$Gene[1]),
      ylim_vals = var_info$ylim,
      annotations = annotations,
      lm_line = lm_line,  # Precomputed linear model
      error_bar = TRUE,
      x_breaks = unique(df$conc_num_factor),
      x_labels = unique(as.character(df$conc_num)),
      x_label = unique(df$Drug[1]),
      position = "jitter",
      coord_cartesian = c(0, max(var_info$ylim))  # You can customize this per plot as needed
    )

    # Add box plot configuration for this variable
    configs[[length(configs) + 1]] <- list(
      df = df,
      x_var = "conc_num_factor",
      y_var = variable,
      plot_type = "box",
      title = sprintf("%s      %s (Boxplot)", df$OrfRep[1], df$Gene[1]),
      ylim_vals = var_info$ylim,
      annotations = annotations,
      error_bar = FALSE,
      x_breaks = unique(df$conc_num_factor),
      x_labels = unique(as.character(df$conc_num)),
      x_label = unique(df$Drug[1]),
      coord_cartesian = c(0, max(var_info$ylim))  # Customize this as needed
    )
  }

  return(configs)
}

generate_rank_plot_configs <- function(df, rank_var, zscore_var, var, is_lm = FALSE) {
  configs <- list()
  
  plot_title_prefix <- if (is_lm) "Interaction Z score vs. Rank for" else "Average Z score vs. Rank for"

  # Single config with all sd bands
  configs[[length(configs) + 1]] <- list(
    df = df,
    x_var = rank_var,
    y_var = zscore_var,
    plot_type = "scatter",
    title = paste(plot_title_prefix, var, "Rank Plot"),
    sd_band = c(1, 2, 3),  # Pass all sd bands at once
    enhancer_label = list(
      x = nrow(df) / 2, y = 10,
      label = paste("Deletion Enhancers =", nrow(df[df[[zscore_var]] >= 1, ]))  # Example for the first SD band
    ),
    suppressor_label = list(
      x = nrow(df) / 2, y = -10,
      label = paste("Deletion Suppressors =", nrow(df[df[[zscore_var]] <= -1, ]))
    ),
    shape = 3,
    size = 0.1
  )

  # Non-annotated version
  configs[[length(configs) + 1]] <- list(
    df = df,
    x_var = rank_var,
    y_var = zscore_var,
    plot_type = "scatter",
    title = paste(plot_title_prefix, var, "Rank Plot No Annotations"),
    sd_band = c(1, 2, 3),
    enhancer_label = NULL,
    suppressor_label = NULL,
    shape = 3,
    size = 0.1,
    position = "jitter"
  )

  return(configs)
}

generate_correlation_plot_configs <- function(df, variables) {
  configs <- list()

  for (variable in variables) {
    z_lm_var <- paste0("Z_lm_", variable)
    avg_zscore_var <- paste0("Avg_Zscore_", variable)
    lm_r_squared_col <- paste0("lm_R_squared_", variable)

    configs[[length(configs) + 1]] <- list(
      df = df,
      x_var = avg_zscore_var,
      y_var = z_lm_var,
      plot_type = "scatter",
      title = paste("Avg Zscore vs lm", variable),
      color_var = "Overlap",
      correlation_text = paste("R-squared =", round(df[[lm_r_squared_col]][1], 2)),
      shape = 3,
      geom_smooth = TRUE,
      rect = list(xmin = -2, xmax = 2, ymin = -2, ymax = 2),  # To add the geom_rect layer
      annotate_position = list(x = 0, y = 0),  # Position for the R-squared text
      legend_position = "right"
    )
  }

  return(configs)
}

main <- function() {
  lapply(names(args$experiments), function(exp_name) {
    exp <- args$experiments[[exp_name]]
    exp_path <- exp$path
    exp_sd <- exp$sd
    out_dir <- file.path(exp_path, "zscores")
    out_dir_qc <- file.path(exp_path, "zscores", "qc")
    dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
    dir.create(out_dir_qc, recursive = TRUE, showWarnings = FALSE)

    summary_vars <- c("L", "K", "r", "AUC", "delta_bg") # fields to filter and calculate summary stats across
    group_vars <- c("OrfRep", "conc_num", "conc_num_factor") # default fields to group by
    orf_group_vars <- c("OrfRep", "Gene", "num")
    print_vars <- c("OrfRep", "Plate", "scan", "Col", "Row", "num", "OrfRep", "conc_num", "conc_num_factor",
      "delta_bg_tolerance", "delta_bg", "Gene", "L", "K", "r", "AUC", "NG", "DB")
    
    message("Loading and filtering data")
    df <- load_and_process_data(args$easy_results_file, sd = exp_sd)
    df <- update_gene_names(df, args$sgd_gene_list)
    df <- as_tibble(df)

    # Filter rows that are above tolerance for quality control plots
    df_above_tolerance <- df %>% filter(DB == 1)

    # Set L, r, K, AUC (and delta_bg?) to NA for rows that are above tolerance
    df_na <- df %>% mutate(across(all_of(summary_vars), ~ ifelse(DB == 1, NA, .)))

    # Remove rows with 0 values in L
    df_no_zeros <- df_na %>% filter(L > 0)
    
    # Save some constants
    max_conc <- max(df$conc_num_factor)
    l_half_median <- (median(df_above_tolerance$L, na.rm = TRUE)) / 2
    k_half_median <- (median(df_above_tolerance$K, na.rm = TRUE)) / 2

    message("Calculating summary statistics before quality control")
    ss <- calculate_summary_stats(df, summary_vars, group_vars = group_vars)
    # df_ss <- ss$summary_stats
    df_stats <- ss$df_with_stats
    df_filtered_stats <- df_stats %>%
      {
        non_finite_rows <- filter(., if_any(c(L), ~ !is.finite(.)))
        if (nrow(non_finite_rows) > 0) {
          message("Filtering out the following non-finite rows:")
          print(non_finite_rows %>% select(any_of(print_vars)), n = 200)
        }
        filter(., if_all(c(L), is.finite))
      }

    message("Calculating summary statistics after quality control")
    ss <- calculate_summary_stats(df_na, summary_vars, group_vars = group_vars)
    df_na_ss <- ss$summary_stats
    df_na_stats <- ss$df_with_stats
    write.csv(df_na_ss, file = file.path(out_dir, "summary_stats_all_strains.csv"), row.names = FALSE)
    # Filter out non-finite rows for plotting
    df_na_filtered_stats <- df_na_stats %>%
      {
        non_finite_rows <- filter(., if_any(c(L), ~ !is.finite(.)))
        if (nrow(non_finite_rows) > 0) {
          message("Removed the following non-finite rows:")
          print(non_finite_rows %>% select(any_of(print_vars)), n = 200)
        }
        filter(., if_all(c(L), is.finite))
      }

    message("Calculating summary statistics after quality control excluding zero values")
    ss <- calculate_summary_stats(df_no_zeros, summary_vars, group_vars = group_vars)
    df_no_zeros_stats <- ss$df_with_stats
    df_no_zeros_filtered_stats <- df_no_zeros_stats %>%
      {
        non_finite_rows <- filter(., if_any(c(L), ~ !is.finite(.)))
        if (nrow(non_finite_rows) > 0) {
          message("Removed the following non-finite rows:")
          print(non_finite_rows %>% select(any_of(print_vars)), n = 200)
        }
        filter(., if_all(c(L), is.finite))
      }

    message("Filtering by 2SD of K")
    df_na_within_2sd_k <- df_na_stats %>%
      filter(K >= (mean_K - 2 * sd_K) & K <= (mean_K + 2 * sd_K))
    df_na_outside_2sd_k <- df_na_stats %>%
      filter(K < (mean_K - 2 * sd_K) | K > (mean_K + 2 * sd_K))

    message("Calculating summary statistics for L within 2SD of K")
    # TODO We're omitting the original z_max calculation, not sure if needed?
    ss <- calculate_summary_stats(df_na_within_2sd_k, "L", group_vars = c("conc_num", "conc_num_factor"))
    l_within_2sd_k_ss <- ss$summary_stats
    df_na_l_within_2sd_k_stats <- ss$df_with_stats
    write.csv(l_within_2sd_k_ss,
      file = file.path(out_dir_qc, "max_observed_L_vals_for_spots_within_2sd_K.csv"), row.names = FALSE)
    
    message("Calculating summary statistics for L outside 2SD of K")
    ss <- calculate_summary_stats(df_na_outside_2sd_k, "L", group_vars = c("conc_num", "conc_num_factor"))
    l_outside_2sd_k_ss <- ss$summary_stats
    df_na_l_outside_2sd_k_stats <- ss$df_with_stats
    write.csv(l_outside_2sd_k_ss,
      file = file.path(out_dir, "max_observed_L_vals_for_spots_outside_2sd_K.csv"), row.names = FALSE)

    # Each plots list corresponds to a file
    message("Generating quality control plot configurations")
    l_vs_k_plots <- list(
      list(
        df = df,
        x_var = "L",
        y_var = "K",
        plot_type = "scatter",
        delta_bg_point = TRUE,
        title = "Raw L vs K before quality control",
        color_var = "conc_num",
        error_bar = FALSE,
        legend_position = "right"
      )
    )

    frequency_delta_bg_plots <- list(
      list(
        df = df_filtered_stats,
        x_var = "delta_bg",
        y_var = NULL,
        plot_type = "density",
        title = "Plate analysis by Drug Conc for Delta Background before quality control",
        color_var = "conc_num",
        x_label = "Delta Background",
        y_label = "Density",
        error_bar = FALSE,
        legend_position = "right"),
      list(
        df = df_filtered_stats,
        x_var = "delta_bg",
        y_var = NULL,
        plot_type = "bar",
        title = "Plate analysis by Drug Conc for Delta Background before quality control",
        color_var = "conc_num",
        x_label = "Delta Background",
        y_label = "Count",
        error_bar = FALSE,
        legend_position = "right")
    )

    above_threshold_plots <- list(
      list(
        df = df_above_tolerance,
        x_var = "L",
        y_var = "K",
        plot_type = "scatter",
        delta_bg_point = TRUE,
        title = paste("Raw L vs K for strains above Delta Background threshold of",
          df_above_tolerance$delta_bg_tolerance[[1]], "or above"),
        color_var = "conc_num",
        position = "jitter",
        annotations = list(
          x = l_half_median,
          y = k_half_median,
          label = paste("# strains above Delta Background tolerance =", nrow(df_above_tolerance))
        ),
        error_bar = FALSE,
        legend_position = "right"
      )
    )

    plate_analysis_plots <- list()
    for (var in summary_vars) {
      for (stage in c("before", "after")) {
        if (stage == "before") {
          df_plot <- df_filtered_stats
        } else {
          df_plot <- df_na_filtered_stats
        }
        
        config <- list(
          df = df_plot,
          x_var = "scan",
          y_var = var,
          plot_type = "scatter",
          title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"),
          error_bar = TRUE,
          color_var = "conc_num",
          position = "jitter")

        plate_analysis_plots <- append(plate_analysis_plots, list(config))
      }
    }

    plate_analysis_boxplots <- list()
    for (var in summary_vars) {
      for (stage in c("before", "after")) {
        if (stage == "before") {
          df_plot <- df_filtered_stats
        } else {
          df_plot <- df_na_filtered_stats
        }
        
        config <- list(
          df = df_plot,
          x_var = "scan",
          y_var = var,
          plot_type = "box",
          title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"),
          error_bar = FALSE,
          color_var = "conc_num")

        plate_analysis_boxplots <- append(plate_analysis_boxplots, list(config))
      }
    }

    plate_analysis_no_zeros_plots <- list()
    for (var in summary_vars) {
      config <- list(
        df = df_no_zeros_filtered_stats,
        x_var = "scan",
        y_var = var,
        plot_type = "scatter",
        title = paste("Plate analysis by Drug Conc for", var, "after quality control"),
        error_bar = TRUE,
        color_var = "conc_num",
        position = "jitter")

      plate_analysis_no_zeros_plots <- append(plate_analysis_no_zeros_plots, list(config))
    }

    plate_analysis_no_zeros_boxplots <- list()
    for (var in summary_vars) {
      config <- list(
        df = df_no_zeros_filtered_stats,
        x_var = "scan",
        y_var = var,
        plot_type = "box",
        title = paste("Plate analysis by Drug Conc for", var, "after quality control"),
        error_bar = FALSE,
        color_var = "conc_num"
      )
      plate_analysis_no_zeros_boxplots <- append(plate_analysis_no_zeros_boxplots, list(config))
    }

    l_outside_2sd_k_plots <- list(
      list(
        df = df_na_l_outside_2sd_k_stats,
        x_var = "L",
        y_var = "K",
        plot_type = "scatter",
        delta_bg_point = TRUE,
        title = "Raw L vs K for strains falling outside 2SD of the K mean at each Conc",
        color_var = "conc_num",
        position = "jitter",
        legend_position = "right"
      )
    )

    delta_bg_outside_2sd_k_plots <- list(
      list(
        df = df_na_l_outside_2sd_k_stats,
        x_var = "delta_bg",
        y_var = "K",
        plot_type = "scatter",
        gene_point = TRUE,
        title = "Delta Background vs K for strains falling outside 2SD of the K mean at each Conc",
        color_var = "conc_num",
        position = "jitter",
        legend_position = "right"
      )
    )

    # message("Generating quality control plots")
    # generate_and_save_plots(out_dir_qc, "L_vs_K_before_quality_control", l_vs_k_plots)
    # generate_and_save_plots(out_dir_qc, "frequency_delta_background", frequency_delta_bg_plots)
    # generate_and_save_plots(out_dir_qc, "L_vs_K_above_threshold", above_threshold_plots)
    # generate_and_save_plots(out_dir_qc, "plate_analysis", plate_analysis_plots)
    # generate_and_save_plots(out_dir_qc, "plate_analysis_boxplots", plate_analysis_boxplots)
    # generate_and_save_plots(out_dir_qc, "plate_analysis_no_zeros", plate_analysis_no_zeros_plots)
    # generate_and_save_plots(out_dir_qc, "plate_analysis_no_zeros_boxplots", plate_analysis_no_zeros_boxplots)
    # generate_and_save_plots(out_dir_qc, "L_vs_K_for_strains_2SD_outside_mean_K", l_outside_2sd_k_plots)
    # generate_and_save_plots(out_dir_qc, "delta_background_vs_K_for_strains_2sd_outside_mean_K", delta_bg_outside_2sd_k_plots)

    # Clean up
    rm(df, df_above_tolerance, df_no_zeros, df_no_zeros_stats, df_no_zeros_filtered_stats, ss)
    gc()

    # TODO: Originally this filtered L NA's
    # Let's try to avoid for now since stats have already been calculated

    # Process background strains
    bg_strains <- c("YDL227C")
    lapply(bg_strains, function(strain) {
      
      message("Processing background strain: ", strain)
      
      # Handle missing data by setting zero values to NA
      # and then removing any rows with NA in L col
      df_bg <- df_na %>%
        filter(OrfRep == strain) %>%
        mutate(
          L = if_else(L == 0, NA, L),
          K = if_else(K == 0, NA, K),
          r = if_else(r == 0, NA, r),
          AUC = if_else(AUC == 0, NA, AUC)
        ) %>%
        filter(!is.na(L))
      
      # Recalculate summary statistics for the background strain
      message("Calculating summary statistics for background strain")
      ss_bg <- calculate_summary_stats(df_bg, summary_vars, group_vars = group_vars)
      summary_stats_bg <- ss_bg$summary_stats
      # df_bg_stats <- ss_bg$df_with_stats
      write.csv(summary_stats_bg,
        file = file.path(out_dir, paste0("SummaryStats_BackgroundStrains_", strain, ".csv")),
        row.names = FALSE)
      
      # Filter reference and deletion strains
      # Formerly X2_RF (reference strains)
      df_reference <- df_na_stats %>%
        filter(OrfRep == strain) %>%
        mutate(SM = 0)
      
      # Formerly X2 (deletion strains)
      df_deletion <- df_na_stats %>%
        filter(OrfRep != strain) %>%
        mutate(SM = 0)

      # Set the missing values to the highest theoretical value at each drug conc for L
      # Leave other values as 0 for the max/min
      reference_strain <- df_reference %>%
        group_by(conc_num) %>%
        mutate(
          max_l_theoretical = max(max_L, na.rm = TRUE),
          L = ifelse(L == 0 & !is.na(L) & conc_num > 0, max_l_theoretical, L),
          SM = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, 1, SM),
          L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>%
        ungroup()

      # Ditto for deletion strains
      deletion_strains <- df_deletion %>%
        group_by(conc_num) %>%
        mutate(
          max_l_theoretical = max(max_L, na.rm = TRUE),
          L = ifelse(L == 0 & !is.na(L) & conc_num > 0, max_l_theoretical, L),
          SM = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, 1, SM),
          L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>%
        ungroup()

      message("Calculating interaction scores")
      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)
      zscores_calculations_reference <- reference_results$calculations
      zscores_interactions_reference <- reference_results$interactions
      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_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)
      write.csv(zscores_interactions_reference, file = file.path(out_dir, "RF_ZScores_Interaction.csv"), row.names = FALSE)
      write.csv(zscores_interactions, file = file.path(out_dir, "ZScores_Interaction.csv"), row.names = FALSE)

      # Create interaction plots
      message("Generating interaction plot configurations")
      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))

      # Define conditions for enhancers and suppressors
      # TODO Add to study config file?
      threshold <- 2
      enhancer_condition_L <- zscores_interactions$Avg_Zscore_L >= threshold
      suppressor_condition_L <- zscores_interactions$Avg_Zscore_L <= -threshold
      enhancer_condition_K <- zscores_interactions$Avg_Zscore_K >= threshold
      suppressor_condition_K <- zscores_interactions$Avg_Zscore_K <= -threshold
      
      # Subset data
      enhancers_L <- zscores_interactions[enhancer_condition_L, ]
      suppressors_L <- zscores_interactions[suppressor_condition_L, ]
      enhancers_K <- zscores_interactions[enhancer_condition_K, ]
      suppressors_K <- zscores_interactions[suppressor_condition_K, ]
      
      # Save enhancers and suppressors
      message("Writing enhancer/suppressor csv files")
      write.csv(enhancers_L, file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_L.csv"), row.names = FALSE)
      write.csv(suppressors_L, file = file.path(out_dir, "ZScores_Interaction_Deletion_Suppressors_L.csv"), row.names = FALSE)
      write.csv(enhancers_K, file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_K.csv"), row.names = FALSE)
      write.csv(suppressors_K, file = file.path(out_dir, "ZScores_Interaction_Deletion_Suppressors_K.csv"), row.names = FALSE)
      
      # Combine conditions for enhancers and suppressors
      enhancers_and_suppressors_L <- zscores_interactions[enhancer_condition_L | suppressor_condition_L, ]
      enhancers_and_suppressors_K <- zscores_interactions[enhancer_condition_K | suppressor_condition_K, ]
      
      # Save combined enhancers and suppressors
      write.csv(enhancers_and_suppressors_L,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_and_Suppressors_L.csv"), row.names = FALSE)
      write.csv(enhancers_and_suppressors_K,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_and_Suppressors_K.csv"), row.names = FALSE)
      
      # Handle linear model based enhancers and suppressors
      lm_threshold <- 2
      enhancers_lm_L <- zscores_interactions[zscores_interactions$Z_lm_L >= lm_threshold, ]
      suppressors_lm_L <- zscores_interactions[zscores_interactions$Z_lm_L <= -lm_threshold, ]
      enhancers_lm_K <- zscores_interactions[zscores_interactions$Z_lm_K >= lm_threshold, ]
      suppressors_lm_K <- zscores_interactions[zscores_interactions$Z_lm_K <= -lm_threshold, ]
      
      # Save linear model based enhancers and suppressors
      message("Writing linear model enhancer/suppressor csv files")
      write.csv(enhancers_lm_L,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_L_lm.csv"), row.names = FALSE)
      write.csv(suppressors_lm_L,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Suppressors_L_lm.csv"), row.names = FALSE)
      write.csv(enhancers_lm_K,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Enhancers_K_lm.csv"), row.names = FALSE)
      write.csv(suppressors_lm_K,
        file = file.path(out_dir, "ZScores_Interaction_Deletion_Suppressors_K_lm.csv"), row.names = FALSE)

      # TODO needs explanation
      zscores_interactions_adjusted <- adjust_missing_and_rank(zscores_interactions)
      
      rank_plot_configs <- c(
        generate_rank_plot_configs(zscores_interactions_adjusted, "Rank_L", "Avg_Zscore_L", "L"),
        generate_rank_plot_configs(zscores_interactions_adjusted, "Rank_K", "Avg_Zscore_K", "K")
      )
      generate_and_save_plots(output_dir = out_dir, file_name = "RankPlots",
        plot_configs = rank_plot_configs, grid_layout = list(ncol = 3, nrow = 2))

      rank_lm_plot_config <- c(
        generate_rank_plot_configs(zscores_interactions_adjusted, "Rank_lm_L", "Z_lm_L", "L", is_lm = TRUE),
        generate_rank_plot_configs(zscores_interactions_adjusted, "Rank_lm_K", "Z_lm_K", "K", is_lm = TRUE)
      )
      generate_and_save_plots(output_dir = out_dir, file_name = "RankPlots_lm",
        plot_configs = rank_lm_plot_config, grid_layout = list(ncol = 3, nrow = 2))
      
      # Formerly X_NArm
      zscores_interactions_filtered <- zscores_interactions %>%
        group_by(across(all_of(orf_group_vars))) %>%
          filter(!is.na(Z_lm_L) | !is.na(Avg_Zscore_L))
      
      # Final filtered correlation calculations and plots
      lm_results <- zscores_interactions_filtered %>%
        summarise(
          lm_R_squared_L = if (n() > 1) summary(lm(Z_lm_L ~ Avg_Zscore_L))$r.squared else NA,
          lm_R_squared_K = if (n() > 1) summary(lm(Z_lm_K ~ Avg_Zscore_K))$r.squared else NA,
          lm_R_squared_r = if (n() > 1) summary(lm(Z_lm_r ~ Avg_Zscore_r))$r.squared else NA,
          lm_R_squared_AUC = if (n() > 1) summary(lm(Z_lm_AUC ~ Avg_Zscore_AUC))$r.squared else NA
        )

      zscores_interactions_filtered <- zscores_interactions_filtered %>%
        left_join(lm_results, by = orf_group_vars) %>%
        mutate(
          Overlap = case_when(
            Z_lm_L >= 2 & Avg_Zscore_L >= 2 ~ "Deletion Enhancer Both",
            Z_lm_L <= -2 & Avg_Zscore_L <= -2 ~ "Deletion Suppressor Both",
            Z_lm_L >= 2 & Avg_Zscore_L <= 2 ~ "Deletion Enhancer lm only",
            Z_lm_L <= -2 & Avg_Zscore_L >= -2 ~ "Deletion Suppressor lm only",
            Z_lm_L >= 2 & Avg_Zscore_L <= -2 ~ "Deletion Enhancer lm, Deletion Suppressor Avg Z score",
            Z_lm_L <= -2 & Avg_Zscore_L >= 2 ~ "Deletion Suppressor lm, Deletion Enhancer Avg Z score",
            TRUE ~ "No Effect"
          )
        ) %>%
        ungroup()

      rank_plot_configs <- c(
        generate_rank_plot_configs(zscores_interactions_filtered, "Rank_L", "Avg_Zscore_L", "L"),
        generate_rank_plot_configs(zscores_interactions_filtered, "Rank_K", "Avg_Zscore_K", "K")
      )
      generate_and_save_plots(output_dir = out_dir, file_name = "RankPlots",
        plot_configs = rank_plot_configs, grid_layout = list(ncol = 3, nrow = 2))
          
      rank_lm_plot_configs <- c(
        generate_rank_plot_configs(zscores_interactions_filtered, "Rank_lm_L", "Z_lm_L", "L", is_lm = TRUE),
        generate_rank_plot_configs(zscores_interactions_filtered, "Rank_lm_K", "Z_lm_K", "K", is_lm = TRUE)
      )
      generate_and_save_plots(output_dir = out_dir, file_name = "RankPlots_lm",
        plot_configs = rank_lm_plot_configs, grid_layout = list(ncol = 3, nrow = 2))

      correlation_plot_configs <- generate_correlation_plot_configs(zscores_interactions_filtered, interaction_vars)
      generate_and_save_plots(output_dir = out_dir, file_name = "Avg_Zscore_vs_lm_NA_rm",
        plot_configs = correlation_plot_configs, grid_layout = list(ncol = 2, nrow = 2))
    })
  })
}
main()