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")) # # TODO for debug # df_duplicates <- df %>% # group_by(OrfRep, Gene, num) %>% # filter(n() > 1) # interactions_duplicates <- interactions %>% # group_by(OrfRep, Gene, num) %>% # filter(n() > 1) # print(df_duplicates) # print(interactions_duplicates) 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(interactive) { # Use appropriate helper function based on plot type plot <- switch(config$plot_type, "scatter" = generate_scatter_plot(plot_base, config, interactive = 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(interactive = FALSE) # Generate the interactive plot interactive_plot <- generate_plot(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, 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 `interactive` is TRUE if (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()