suppressMessages({ library("ggplot2") library("plotly") library("htmlwidgets") library("dplyr") library("rlang") library("ggthemes") library("data.table") library("future") library("furrr") library("purrr") }) # These parallelization libraries are very noisy suppressPackageStartupMessages({ library("future") library("furrr") library("purrr") }) options(warn = 2) # 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) } out_dir <- normalizePath(args[1], mustWork = FALSE) sgd_gene_list <- normalizePath(args[2], mustWork = FALSE) easy_results_file <- normalizePath(args[3], mustWork = FALSE) # The remaining arguments should be in groups of 3 exp_args <- args[-(1:3)] if (length(exp_args) %% 3 != 0) { stop("Experiment arguments should be in groups of 3: path, name, sd.") } experiments <- list() for (i in seq(1, length(exp_args), by = 3)) { exp_name <- exp_args[i + 1] experiments[[exp_name]] <- list( path = normalizePath(exp_args[i], mustWork = FALSE), sd = as.numeric(exp_args[i + 2]) ) } list( out_dir = out_dir, sgd_gene_list = sgd_gene_list, easy_results_file = easy_results_file, 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, size = 18), axis.title.x = element_text(vjust = -0.2, size = 18), axis.line = element_line(colour = "black"), axis.text.x = element_text(size = 16), axis.text.y = element_text(size = 16), 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_filter_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 = factor(as.numeric(factor(conc_num)) - 1), conc_num_factor_num = as.numeric(conc_num_factor) ) 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_stats <- function(df, variables, group_vars) { summary_stats <- df %>% group_by(across(all_of(group_vars))) %>% summarise( N = n(), 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 = ~sd(., na.rm = TRUE) / sqrt(N - 1) # TODO non-standard SE, needs explanation ), .names = "{.fn}_{.col}" ), .groups = "drop" ) # Create a cleaned version of df that doesn't overlap with summary_stats cleaned_df <- df %>% select(-any_of(setdiff(intersect(names(df), names(summary_stats)), group_vars))) df_joined <- left_join(cleaned_df, summary_stats, by = group_vars) return(list(summary_stats = summary_stats, df_with_stats = df_joined)) } calculate_interaction_scores <- function(df, max_conc, bg_stats, group_vars = c("OrfRep", "Gene", "num")) { # Calculate total concentration variables total_conc_num <- length(unique(df$conc_num)) calculations <- df %>% group_by(across(all_of(group_vars))) %>% mutate( 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(DB, na.rm = TRUE) - 1, # Calculate raw data Raw_Shift_L = first(mean_L) - bg_stats$mean_L, Raw_Shift_K = first(mean_K) - bg_stats$mean_K, Raw_Shift_r = first(mean_r) - bg_stats$mean_r, Raw_Shift_AUC = first(mean_AUC) - bg_stats$mean_AUC, Z_Shift_L = first(Raw_Shift_L) / bg_stats$sd_L, Z_Shift_K = first(Raw_Shift_K) / bg_stats$sd_K, Z_Shift_r = first(Raw_Shift_r) / bg_stats$sd_r, Z_Shift_AUC = first(Raw_Shift_AUC) / bg_stats$sd_AUC, 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, 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), # Calculate Z-scores 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, # Fit linear models and store in list-columns gene_lm_L = list(lm(Delta_L ~ conc_num_factor_num, data = pick(everything()))), gene_lm_K = list(lm(Delta_K ~ conc_num_factor_num, data = pick(everything()))), gene_lm_r = list(lm(Delta_r ~ conc_num_factor_num, data = pick(everything()))), gene_lm_AUC = list(lm(Delta_AUC ~ conc_num_factor_num, data = pick(everything()))), # Extract coefficients using purrr::map_dbl lm_intercept_L = map_dbl(gene_lm_L, ~ coef(.x)[1]), lm_slope_L = map_dbl(gene_lm_L, ~ coef(.x)[2]), lm_intercept_K = map_dbl(gene_lm_K, ~ coef(.x)[1]), lm_slope_K = map_dbl(gene_lm_K, ~ coef(.x)[2]), lm_intercept_r = map_dbl(gene_lm_r, ~ coef(.x)[1]), lm_slope_r = map_dbl(gene_lm_r, ~ coef(.x)[2]), lm_intercept_AUC = map_dbl(gene_lm_AUC, ~ coef(.x)[1]), lm_slope_AUC = map_dbl(gene_lm_AUC, ~ coef(.x)[2]), # Calculate lm_Score_* based on coefficients lm_Score_L = max_conc * lm_slope_L + lm_intercept_L, lm_Score_K = max_conc * lm_slope_K + lm_intercept_K, lm_Score_r = max_conc * lm_slope_r + lm_intercept_r, lm_Score_AUC = max_conc * lm_slope_AUC + lm_intercept_AUC, # Calculate R-squared values R_Squared_L = map_dbl(gene_lm_L, ~ summary(.x)$r.squared), R_Squared_K = map_dbl(gene_lm_K, ~ summary(.x)$r.squared), R_Squared_r = map_dbl(gene_lm_r, ~ summary(.x)$r.squared), R_Squared_AUC = map_dbl(gene_lm_AUC, ~ summary(.x)$r.squared) ) %>% ungroup() # Calculate overall mean and SD for lm_Score_* variables lm_means_sds <- calculations %>% summarise( lm_mean_L = mean(lm_Score_L, na.rm = TRUE), lm_sd_L = sd(lm_Score_L, na.rm = TRUE), lm_mean_K = mean(lm_Score_K, na.rm = TRUE), lm_sd_K = sd(lm_Score_K, na.rm = TRUE), lm_mean_r = mean(lm_Score_r, na.rm = TRUE), lm_sd_r = sd(lm_Score_r, na.rm = TRUE), lm_mean_AUC = mean(lm_Score_AUC, na.rm = TRUE), lm_sd_AUC = sd(lm_Score_AUC, na.rm = TRUE) ) calculations <- calculations %>% mutate( Z_lm_L = (lm_Score_L - lm_means_sds$lm_mean_L) / lm_means_sds$lm_sd_L, Z_lm_K = (lm_Score_K - lm_means_sds$lm_mean_K) / lm_means_sds$lm_sd_K, Z_lm_r = (lm_Score_r - lm_means_sds$lm_mean_r) / lm_means_sds$lm_sd_r, Z_lm_AUC = (lm_Score_AUC - lm_means_sds$lm_mean_AUC) / lm_means_sds$lm_sd_AUC ) # Summarize some of the stats interactions <- calculations %>% group_by(across(all_of(group_vars))) %>% mutate( # Calculate raw shifts 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), # Calculate Z-shifts 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 Z-scores Sum_Z_Score_L = sum(Zscore_L), Sum_Z_Score_K = sum(Zscore_K), Sum_Z_Score_r = sum(Zscore_r), Sum_Z_Score_AUC = sum(Zscore_AUC), # Calculate Average Z-scores Avg_Zscore_L = Sum_Z_Score_L / num_non_removed_concs, Avg_Zscore_K = Sum_Z_Score_K / num_non_removed_concs, Avg_Zscore_r = Sum_Z_Score_r / (total_conc_num - 1), Avg_Zscore_AUC = Sum_Z_Score_AUC / (total_conc_num - 1) ) %>% arrange(desc(Z_lm_L), desc(NG)) %>% ungroup() # Declare column order for output calculations <- calculations %>% select( "OrfRep", "Gene", "num", "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") interactions <- interactions %>% select( "OrfRep", "Gene", "conc_num", "conc_num_factor", "num", "NG", "DB", "SM", "Raw_Shift_L", "Raw_Shift_K", "Raw_Shift_r", "Raw_Shift_AUC", "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC", "Sum_Z_Score_L", "Sum_Z_Score_K", "Sum_Z_Score_r", "Sum_Z_Score_AUC", "Avg_Zscore_L", "Avg_Zscore_K", "Avg_Zscore_r", "Avg_Zscore_AUC", "lm_Score_L", "lm_Score_K", "lm_Score_r", "lm_Score_AUC", "R_Squared_L", "R_Squared_K", "R_Squared_r", "R_Squared_AUC", "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC") cleaned_df <- df %>% select(-any_of( setdiff(intersect(names(df), names(interactions)), c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor")))) interactions_joined <- left_join(cleaned_df, interactions, by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor")) return(list( calculations = calculations, interactions = interactions, interactions_joined = interactions_joined)) } generate_and_save_plots <- function(out_dir, filename, plot_configs, grid_layout = NULL) { message("Generating ", filename, ".pdf and ", filename, ".html") # 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 aes_mapping <- if (config$plot_type == "bar" || config$plot_type == "density") { if (is.null(config$color_var)) { aes(x = .data[[config$x_var]]) } else { aes(x = .data[[config$x_var]], color = as.factor(.data[[config$color_var]])) } } else if (is.null(config$color_var)) { aes(x = .data[[config$x_var]], y = .data[[config$y_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) # Use appropriate helper function based on plot type plot <- switch(config$plot_type, "scatter" = generate_scatter_plot(plot_base, config), "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 (!is.null(config$legend_position)) { plot <- plot + theme(legend.position = config$legend_position) } # Add title and labels 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) } # Apply scale_color_discrete(guide = FALSE) when color_var is NULL if (is.null(config$color_var)) { plot <- plot + scale_color_discrete(guide = "none") } # Add interactive tooltips for plotly tooltip_vars <- c() if (config$plot_type == "scatter") { if (!is.null(config$delta_bg_point) && config$delta_bg_point) { tooltip_vars <- c(tooltip_vars, "OrfRep", "Gene", "delta_bg") } else if (!is.null(config$gene_point) && config$gene_point) { tooltip_vars <- c(tooltip_vars, "OrfRep", "Gene") } else if (!is.null(config$y_var) && !is.null(config$x_var)) { tooltip_vars <- c(config$x_var, config$y_var) } } # Convert to plotly object if (length(tooltip_vars) > 0) { plotly_plot <- ggplotly(plot, tooltip = tooltip_vars) } else { plotly_plot <- ggplotly(plot, tooltip = "none") } # Adjust legend position if specified 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]] <- plot plotly_plots[[i]] <- plotly_plot } # Save static PDF plots pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 14, height = 9) lapply(static_plots, print) dev.off() # Combine and save interactive HTML plots combined_plot <- subplot( plotly_plots, nrows = if (!is.null(grid_layout) && !is.null(grid_layout$nrow)) { grid_layout$nrow } else { # Calculate nrow based on the length of plotly_plots (default 1 row if only one plot) ceiling(length(plotly_plots) / ifelse(!is.null(grid_layout) && !is.null(grid_layout$ncol), grid_layout$ncol, 1)) }, margin = 0.05 ) saveWidget(combined_plot, file = file.path(out_dir, paste0(filename, ".html")), selfcontained = TRUE) } generate_scatter_plot <- function(plot, config) { shape <- if (!is.null(config$shape)) config$shape else 3 size <- if (!is.null(config$size)) config$size else 0.1 position <- if (!is.null(config$position) && config$position == "jitter") { position_jitter(width = 0.1, height = 0) } else { "identity" } plot <- plot + geom_point( shape = shape, size = size, position = position ) if (!is.null(config$cyan_points) && config$cyan_points) { plot <- plot + geom_point( aes(x = .data[[config$x_var]], y = .data[[config$y_var]]), color = "cyan", shape = 3, size = 0.5 ) } # Add Smooth Line if specified if (!is.null(config$smooth) && config$smooth) { smooth_color <- if (!is.null(config$smooth_color)) config$smooth_color else "blue" if (!is.null(config$lm_line)) { plot <- plot + geom_abline( intercept = config$lm_line$intercept, slope = config$lm_line$slope, color = smooth_color ) } else { plot <- plot + geom_smooth( method = "lm", se = FALSE, color = smooth_color ) } } # Add SD Bands if specified if (!is.null(config$sd_band)) { plot <- plot + annotate( "rect", xmin = -Inf, xmax = Inf, ymin = config$sd_band, ymax = Inf, fill = ifelse(!is.null(config$fill_positive), config$fill_positive, "#542788"), alpha = ifelse(!is.null(config$alpha_positive), config$alpha_positive, 0.3) ) + annotate( "rect", xmin = -Inf, xmax = Inf, ymin = -config$sd_band, ymax = -Inf, fill = ifelse(!is.null(config$fill_negative), config$fill_negative, "orange"), alpha = ifelse(!is.null(config$alpha_negative), config$alpha_negative, 0.3) ) + geom_hline( yintercept = c(-config$sd_band, config$sd_band), color = ifelse(!is.null(config$hl_color), config$hl_color, "gray") ) } # Add Rectangles if specified if (!is.null(config$rectangles)) { for (rect in config$rectangles) { plot <- plot + annotate( "rect", xmin = rect$xmin, xmax = rect$xmax, ymin = rect$ymin, ymax = rect$ymax, fill = ifelse(is.null(rect$fill), NA, rect$fill), color = ifelse(is.null(rect$color), "black", rect$color), alpha = ifelse(is.null(rect$alpha), 0.1, rect$alpha) ) } } # Add Error Bars if specified if (!is.null(config$error_bar) && config$error_bar && !is.null(config$y_var)) { 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 ) } # Customize X-axis if specified 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 ) } # Apply coord_cartesian if specified if (!is.null(config$coord_cartesian)) { plot <- plot + coord_cartesian(ylim = config$coord_cartesian) } # Set Y-axis limits if specified if (!is.null(config$ylim_vals)) { plot <- plot + scale_y_continuous(limits = config$ylim_vals) } # Add Annotations if specified if (!is.null(config$annotations)) { for (annotation in config$annotations) { plot <- plot + annotate( "text", x = annotation$x, y = annotation$y, label = annotation$label, hjust = ifelse(is.null(annotation$hjust), 0.5, annotation$hjust), vjust = ifelse(is.null(annotation$vjust), 0.5, annotation$vjust), size = ifelse(is.null(annotation$size), 4, annotation$size), color = ifelse(is.null(annotation$color), "black", annotation$color) ) } } # Add Title if specified if (!is.null(config$title)) { plot <- plot + ggtitle(config$title) } # Adjust Legend Position if specified 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) } generate_plate_analysis_plot_configs <- function(variables, stages = c("before", "after"), df_before = NULL, df_after = NULL, plot_type = "scatter") { plots <- list() for (var in variables) { for (stage in stages) { df_plot <- if (stage == "before") df_before else df_after # Adjust settings based on plot_type if (plot_type == "scatter") { error_bar <- TRUE position <- "jitter" } else if (plot_type == "box") { error_bar <- FALSE position <- NULL } config <- list( df = df_plot, x_var = "scan", y_var = var, plot_type = plot_type, title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"), error_bar = error_bar, color_var = "conc_num_factor", position = position ) plots <- append(plots, list(config)) } } return(plots) } generate_interaction_plot_configs <- function(df, variables) { configs <- list() # Set the y-limits for each variable df_filtered <- df %>% filter( !is.na(L) & between(L, -65, 65), !is.na(K) & between(K, -65, 65), !is.na(r) & between(r, -0.65, 0.65), !is.na(AUC) & between(AUC, -6500, 6500) ) # Define annotation label functions generate_annotation_labels <- function(df, var, annotation_name) { switch(annotation_name, ZShift = paste("ZShift =", round(df[[paste0("Z_Shift_", var)]], 2)), lm_ZScore = paste("lm ZScore =", round(df[[paste0("Z_lm_", var)]], 2)), NG = paste("NG =", df$NG), DB = paste("DB =", df$DB), SM = paste("SM =", df$SM), NULL # Default case for unrecognized annotation names ) } # Define annotation positions relative to the y-axis range calculate_annotation_positions <- function(y_range) { y_min <- min(y_range) y_max <- max(y_range) y_span <- y_max - y_min list( ZShift = y_max - 0.1 * y_span, lm_ZScore = y_max - 0.2 * y_span, NG = y_min + 0.2 * y_span, DB = y_min + 0.1 * y_span, SM = y_min + 0.05 * y_span ) } # Create configurations for each variable for (variable in variables) { y_range <- limits_map[[variable]] annotation_positions <- calculate_annotation_positions(y_range) lm_line <- list( intercept = df_filtered[[paste0("lm_intercept_", variable)]], slope = df_filtered[[paste0("lm_slope_", variable)]] ) # Determine x-axis midpoint num_levels <- length(levels(df_filtered$conc_num_factor)) x_pos <- (1 + num_levels) / 2 # Midpoint of x-axis # Generate annotations annotations <- lapply(names(annotation_positions), function(annotation_name) { label <- generate_annotation_labels(df_filtered, variable, annotation_name) y_pos <- annotation_positions[[annotation_name]] if (!is.null(label)) { list(x = x_pos, y = y_pos, label = label) } else { message(paste("Warning: No annotation found for", annotation_name)) NULL } }) # Remove NULL annotations annotations <- Filter(Negate(is.null), annotations) # Shared plot settings plot_settings <- list( df = df_filtered, x_var = "conc_num_factor", y_var = variable, ylim_vals = y_range, annotations = annotations, lm_line = lm_line, x_breaks = levels(df_filtered$conc_num_factor), x_labels = levels(df_filtered$conc_num_factor), x_label = unique(df_filtered$Drug[1]), coord_cartesian = y_range # Use the actual y-limits ) # Scatter plot config configs[[length(configs) + 1]] <- modifyList(plot_settings, list( plot_type = "scatter", title = sprintf("%s %s", df_filtered$OrfRep[1], df_filtered$Gene[1]), error_bar = TRUE, position = "jitter" )) # Box plot config configs[[length(configs) + 1]] <- modifyList(plot_settings, list( plot_type = "box", title = sprintf("%s %s (box plot)", df_filtered$OrfRep[1], df_filtered$Gene[1]), error_bar = FALSE )) } return(configs) } generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FALSE, overlap_color = FALSE) { sd_bands <- c(1, 2, 3) avg_zscore_cols <- paste0("Avg_Zscore_", variables) z_lm_cols <- paste0("Z_lm_", variables) rank_avg_zscore_cols <- paste0("Rank_", variables) rank_z_lm_cols <- paste0("Rank_lm_", variables) configs <- list() if (adjust) { message("Replacing NA with 0.001 for Avg_Zscore_ and Z_lm_ columns for ranks") df <- df %>% mutate( across(all_of(avg_zscore_cols), ~ifelse(is.na(.), 0.001, .)), across(all_of(z_lm_cols), ~ifelse(is.na(.), 0.001, .)) ) } message("Calculating ranks for Avg_Zscore and Z_lm columns") rank_col_mapping <- setNames(rank_avg_zscore_cols, avg_zscore_cols) df_ranked <- df %>% mutate(across(all_of(avg_zscore_cols), ~rank(., na.last = "keep"), .names = "{rank_col_mapping[.col]}")) rank_lm_col_mapping <- setNames(rank_z_lm_cols, z_lm_cols) df_ranked <- df_ranked %>% mutate(across(all_of(z_lm_cols), ~rank(., na.last = "keep"), .names = "{rank_lm_col_mapping[.col]}")) # SD-based plots for L and K for (variable in c("L", "K")) { if (is_lm) { rank_var <- paste0("Rank_lm_", variable) zscore_var <- paste0("Z_lm_", variable) y_label <- paste("Int Z score", variable) } else { rank_var <- paste0("Rank_", variable) zscore_var <- paste0("Avg_Zscore_", variable) y_label <- paste("Avg Z score", variable) } for (sd_band in sd_bands) { num_enhancers <- sum(df_ranked[[zscore_var]] >= sd_band, na.rm = TRUE) num_suppressors <- sum(df_ranked[[zscore_var]] <= -sd_band, na.rm = TRUE) # Annotated plot configuration configs[[length(configs) + 1]] <- list( df = df_ranked, x_var = rank_var, y_var = zscore_var, plot_type = "scatter", title = paste(y_label, "vs. Rank for", variable, "above", sd_band, "SD"), sd_band = sd_band, fill_positive = "#542788", fill_negative = "orange", alpha_positive = 0.3, alpha_negative = 0.3, annotations = list( list( x = median(df_ranked[[rank_var]], na.rm = TRUE), y = 10, label = paste("Deletion Enhancers =", num_enhancers), hjust = 0.5, vjust = 1 ), list( x = median(df_ranked[[rank_var]], na.rm = TRUE), y = -10, label = paste("Deletion Suppressors =", num_suppressors), hjust = 0.5, vjust = 0 ) ), shape = 3, size = 0.1, y_label = y_label, x_label = "Rank", legend_position = "none" ) # Non-Annotated Plot Configuration configs[[length(configs) + 1]] <- list( df = df_ranked, x_var = rank_var, y_var = zscore_var, plot_type = "scatter", title = paste(y_label, "vs. Rank for", variable, "above", sd_band, "SD No Annotations"), sd_band = sd_band, fill_positive = "#542788", fill_negative = "orange", alpha_positive = 0.3, alpha_negative = 0.3, annotations = NULL, shape = 3, size = 0.1, y_label = y_label, x_label = "Rank", legend_position = "none" ) } } # Avg ZScore and Rank Avg ZScore Plots for r, L, K, and AUC for (variable in variables) { for (plot_type in c("Avg Zscore vs lm", "Rank Avg Zscore vs lm")) { title <- paste(plot_type, variable) # Define specific variables based on plot type if (plot_type == "Avg Zscore vs lm") { x_var <- paste0("Avg_Zscore_", variable) y_var <- paste0("Z_lm_", variable) rectangles <- list( list(xmin = -2, xmax = 2, ymin = -2, ymax = 2, fill = NA, color = "grey20", alpha = 0.1 ) ) } else if (plot_type == "Rank Avg Zscore vs lm") { x_var <- paste0("Rank_", variable) y_var <- paste0("Rank_lm_", variable) rectangles <- NULL } # Fit the linear model lm_model <- lm(as.formula(paste(y_var, "~", x_var)), data = df_ranked) # Extract intercept and slope from the model coefficients intercept <- coef(lm_model)[1] slope <- coef(lm_model)[2] configs[[length(configs) + 1]] <- list( df = df_ranked, x_var = x_var, y_var = y_var, plot_type = "scatter", title = title, annotations = list( list( x = median(df_ranked[[rank_var]], na.rm = TRUE), y = 10, label = paste("Deletion Enhancers =", num_enhancers), hjust = 0.5, vjust = 1 ), list( x = median(df_ranked[[rank_var]], na.rm = TRUE), y = -10, label = paste("Deletion Suppressors =", num_suppressors), hjust = 0.5, vjust = 0 ) ), shape = 3, size = 0.1, smooth = TRUE, smooth_color = "black", lm_line = list(intercept = intercept, slope = slope), legend_position = "right", color_var = if (overlap_color) "Overlap" else NULL, x_label = x_var, y_label = y_var, rectangles = rectangles ) } } return(configs) } generate_correlation_plot_configs <- function(df) { # Define relationships for plotting relationships <- list( list(x = "Z_lm_L", y = "Z_lm_K", label = "Interaction L vs. Interaction K"), list(x = "Z_lm_L", y = "Z_lm_r", label = "Interaction L vs. Interaction r"), list(x = "Z_lm_L", y = "Z_lm_AUC", label = "Interaction L vs. Interaction AUC"), list(x = "Z_lm_K", y = "Z_lm_r", label = "Interaction K vs. Interaction r"), list(x = "Z_lm_K", y = "Z_lm_AUC", label = "Interaction K vs. Interaction AUC"), list(x = "Z_lm_r", y = "Z_lm_AUC", label = "Interaction r vs. Interaction AUC") ) configs <- list() for (rel in relationships) { # Fit linear model lm_model <- lm(as.formula(paste(rel$y, "~", rel$x)), data = df) lm_summary <- summary(lm_model) # Construct plot configuration config <- list( df = df, x_var = rel$x, y_var = rel$y, plot_type = "scatter", title = rel$label, x_label = paste("z-score", gsub("Z_lm_", "", rel$x)), y_label = paste("z-score", gsub("Z_lm_", "", rel$y)), annotations = list( list( x = Inf, y = Inf, label = paste("R-squared =", round(lm_summary$r.squared, 3)), hjust = 1.1, vjust = 2, size = 4, color = "black" ) ), smooth = TRUE, smooth_color = "tomato3", lm_line = list(intercept = coef(lm_model)[1], slope = coef(lm_model)[2]), legend_position = "right", shape = 3, size = 0.5, color_var = "Overlap", rectangles = list( list( xmin = -2, xmax = 2, ymin = -2, ymax = 2, fill = NA, color = "grey20", alpha = 0.1 ) ), cyan_points = TRUE ) configs[[length(configs) + 1]] <- config } 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 interaction_vars <- c("L", "K", "r", "AUC") # fields to calculate interaction z-scores 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 for experiment: ", exp_name) df <- load_and_filter_data(args$easy_results_file, sd = exp_sd) %>% update_gene_names(args$sgd_gene_list) %>% as_tibble() # Filter rows above delta background tolerance df_above_tolerance <- df %>% filter(DB == 1) df_na <- df %>% mutate(across(all_of(summary_vars), ~ ifelse(DB == 1, NA, .))) # formerly X df_no_zeros <- df_na %>% filter(L > 0) # formerly X_noZero # Save some constants max_conc <- max(df$conc_num_factor_num) 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") df_stats <- calculate_summary_stats( df = df, variables = summary_vars, group_vars = c("conc_num", "conc_num_factor"))$df_with_stats message("Filtering non-finite data") message("Calculating summary statistics after quality control") ss <- calculate_summary_stats( df = df_na, variables = summary_vars, group_vars = c("conc_num", "conc_num_factor")) 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) df_na_stats <- df_na_stats %>% 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 ) # Pull the background means and standard deviations from zero concentration for interactions bg_stats <- df_na_stats %>% filter(conc_num == 0) %>% summarise( mean_L = first(mean_L), mean_K = first(mean_K), mean_r = first(mean_r), mean_AUC = first(mean_AUC), sd_L = first(sd_L), sd_K = first(sd_K), sd_r = first(sd_r), sd_AUC = first(sd_AUC) ) message("Calculating summary statistics after quality control excluding zero values") df_no_zeros_stats <- calculate_summary_stats( df = df_no_zeros, variables = summary_vars, group_vars = c("conc_num", "conc_num_factor") )$df_with_stats 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"))$summary_stats write.csv(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")) df_na_l_outside_2sd_k_stats <- ss$df_with_stats write.csv(ss$summary_stats, 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 l_vs_k_plot_configs <- 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_factor", error_bar = FALSE, legend_position = "right" ) ) frequency_delta_bg_plot_configs <- list( list( df = df_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_factor", x_label = "Delta Background", y_label = "Density", error_bar = FALSE, legend_position = "right"), list( df = df_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_factor", x_label = "Delta Background", y_label = "Count", error_bar = FALSE, legend_position = "right") ) above_threshold_plot_configs <- 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_factor", position = "jitter", annotations = list( 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_plot_configs <- generate_plate_analysis_plot_configs( variables = summary_vars, df_before = df_stats, df_after = df_na_stats, ) plate_analysis_boxplot_configs <- generate_plate_analysis_plot_configs( variables = summary_vars, df_before = df_stats, df_after = df_na_stats, plot_type = "box" ) plate_analysis_no_zeros_plot_configs <- generate_plate_analysis_plot_configs( variables = summary_vars, stages = c("after"), # Only after QC df_after = df_no_zeros_stats, ) plate_analysis_no_zeros_boxplot_configs <- generate_plate_analysis_plot_configs( variables = summary_vars, stages = c("after"), # Only after QC df_after = df_no_zeros_stats, plot_type = "box" ) l_outside_2sd_k_plot_configs <- 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_factor", position = "jitter", legend_position = "right" ) ) delta_bg_outside_2sd_k_plot_configs <- 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_factor", position = "jitter", legend_position = "right" ) ) message("Generating quality control plots") # TODO trying out some parallelization # future::plan(future::multicore, workers = parallel::detectCores()) future::plan(future::multisession, workers = 3) plot_configs <- list( list(out_dir = out_dir_qc, filename = "L_vs_K_before_quality_control", plot_configs = l_vs_k_plot_configs), list(out_dir = out_dir_qc, filename = "frequency_delta_background", plot_configs = frequency_delta_bg_plot_configs), list(out_dir = out_dir_qc, filename = "L_vs_K_above_threshold", plot_configs = above_threshold_plot_configs), list(out_dir = out_dir_qc, filename = "plate_analysis", plot_configs = plate_analysis_plot_configs), list(out_dir = out_dir_qc, filename = "plate_analysis_boxplots", plot_configs = plate_analysis_boxplot_configs), list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros", plot_configs = plate_analysis_no_zeros_plot_configs), list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros_boxplots", plot_configs = plate_analysis_no_zeros_boxplot_configs), list(out_dir = out_dir_qc, name = "L_vs_K_for_strains_2SD_outside_mean_K", plot_configs = l_outside_2sd_k_plot_configs), list(out_dir = out_dir_qc, name = "delta_background_vs_K_for_strains_2sd_outside_mean_K", plot_configs = delta_bg_outside_2sd_k_plot_configs) ) # Generating quality control plots in parallel # furrr::future_map(plot_configs, function(config) { # generate_and_save_plots(config$out_dir, config$filename, config$plot_configs) # }, .options = furrr_options(seed = TRUE)) # 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 = c("OrfRep", "conc_num", "conc_num_factor")) summary_stats_bg <- ss_bg$summary_stats write.csv(summary_stats_bg, file = file.path(out_dir, paste0("SummaryStats_BackgroundStrains_", strain, ".csv")), row.names = FALSE) # Set the missing values to the highest theoretical value at each drug conc for L # Leave other values as 0 for the max/min df_reference <- df_na_stats %>% # formerly X2_RF filter(OrfRep == strain) %>% filter(!is.na(L)) %>% group_by(conc_num, conc_num_factor) %>% 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, 0), L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>% ungroup() # Ditto for deletion strains df_deletion <- df_na_stats %>% # formerly X2 filter(OrfRep != strain) %>% filter(!is.na(L)) %>% mutate(SM = 0) %>% group_by(conc_num, conc_num_factor) %>% 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 reference strain interaction scores") df_reference_stats <- calculate_summary_stats( df = df_reference, variables = interaction_vars, group_vars = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor") )$df_with_stats reference_results <- calculate_interaction_scores(df_reference_stats, max_conc, bg_stats, group_vars = c("OrfRep", "Gene", "num")) zscores_calculations_reference <- reference_results$calculations zscores_interactions_reference <- reference_results$interactions zscores_interactions_reference_joined <- reference_results$interactions_joined message("Calculating deletion strain(s) interactions scores") df_deletion_stats <- calculate_summary_stats( df = df_deletion, variables = interaction_vars, group_vars = c("OrfRep", "Gene", "conc_num", "conc_num_factor") )$df_with_stats deletion_results <- calculate_interaction_scores(df_deletion_stats, max_conc, bg_stats, group_vars = c("OrfRep")) zscores_calculations <- deletion_results$calculations zscores_interactions <- deletion_results$interactions 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 reference interaction plots") reference_plot_configs <- generate_interaction_plot_configs(zscores_interactions_reference_joined, interaction_vars) generate_and_save_plots(out_dir, "RF_interactionPlots", reference_plot_configs, grid_layout = list(ncol = 4, nrow = 3)) message("Generating deletion interaction plots") deletion_plot_configs <- generate_interaction_plot_configs(zscores_interactions_joined, interaction_vars) 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? 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) message("Generating rank plots") rank_plot_configs <- generate_rank_plot_configs( df = zscores_interactions_joined, variables = interaction_vars, is_lm = FALSE, adjust = TRUE ) generate_and_save_plots(out_dir = out_dir, filename = "RankPlots", plot_configs = rank_plot_configs, grid_layout = list(ncol = 3, nrow = 2)) message("Generating ranked linear model plots") rank_lm_plot_configs <- generate_rank_plot_configs( df = zscores_interactions_joined, variables = interaction_vars, is_lm = TRUE, adjust = TRUE ) generate_and_save_plots(out_dir = out_dir, filename = "RankPlots_lm", plot_configs = rank_lm_plot_configs, grid_layout = list(ncol = 3, nrow = 2)) message("Filtering and reranking plots") # Formerly X_NArm zscores_interactions_filtered <- zscores_interactions_joined %>% filter(!is.na(Z_lm_L) & !is.na(Avg_Zscore_L)) %>% 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 Enhancer Avg Zscore only", Z_lm_L <= -2 & Avg_Zscore_L >= -2 ~ "Deletion Suppressor lm only", Z_lm_L >= -2 & Avg_Zscore_L <= -2 ~ "Deletion Suppressor Avg Zscore 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" ), lm_R_squared_L = summary(lm(Z_lm_L ~ Avg_Zscore_L))$r.squared, lm_R_squared_K = summary(lm(Z_lm_K ~ Avg_Zscore_K))$r.squared, lm_R_squared_r = summary(lm(Z_lm_r ~ Avg_Zscore_r))$r.squared, lm_R_squared_AUC = summary(lm(Z_lm_AUC ~ Avg_Zscore_AUC))$r.squared ) message("Generating filtered ranked plots") rank_plot_filtered_configs <- generate_rank_plot_configs( df = zscores_interactions_filtered_ranked, variables = interaction_vars, is_lm = FALSE, adjust = FALSE, overlap_color = TRUE ) generate_and_save_plots( out_dir = out_dir, filename = "RankPlots_na_rm", plot_configs = rank_plot_filtered_configs, grid_layout = list(ncol = 3, nrow = 2)) message("Generating filtered ranked linear model plots") rank_plot_lm_filtered_configs <- generate_rank_plot_configs( df = zscores_interactions_filtered_ranked, variables = interaction_vars, is_lm = TRUE, adjust = FALSE, overlap_color = TRUE ) generate_and_save_plots( out_dir = out_dir, filename = "RankPlots_lm_na_rm", plot_configs = rank_plot_lm_filtered_configs, grid_layout = list(ncol = 3, nrow = 2)) message("Generating correlation plots") correlation_plot_configs <- generate_correlation_plot_configs(zscores_interactions_filtered) generate_and_save_plots( out_dir = out_dir, filename = "Avg_Zscore_vs_lm_NA_rm", plot_configs = correlation_plot_configs, grid_layout = list(ncol = 2, nrow = 2)) }) }) } main()