suppressMessages({ library("ggplot2") library("plotly") library("htmlwidgets") library("htmltools") library("dplyr") library("rlang") library("ggthemes") library("data.table") library("gridExtra") 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.") } # Extract the experiments into a list 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]) ) } # Extract the trailing number from each path trailing_numbers <- sapply(experiments, function(x) { path <- x$path nums <- gsub("[^0-9]", "", basename(path)) as.integer(nums) }) # Sort the experiments based on the trailing numbers sorted_experiments <- experiments[order(trailing_numbers)] list( out_dir = out_dir, sgd_gene_list = sgd_gene_list, easy_results_file = easy_results_file, experiments = sorted_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) theme_publication <- function(base_size = 14, base_family = "sans", legend_position = NULL) { # Ensure that legend_position has a valid value or default to "none" legend_position <- if (is.null(legend_position) || length(legend_position) == 0) "none" else legend_position theme_foundation <- ggthemes::theme_foundation(base_size = base_size, base_family = base_family) theme_foundation %+replace% theme( plot.title = element_text(face = "bold", size = rel(1.6), hjust = 0.5), text = element_text(), panel.background = element_blank(), plot.background = element_blank(), panel.border = element_blank(), axis.title = element_text(face = "bold", size = rel(1.4)), axis.title.y = element_text(angle = 90, vjust = 2), axis.text = element_text(size = rel(1.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 = if (legend_position == "right") { "vertical" } else if (legend_position == "bottom") { "horizontal" } else { NULL # No legend direction if position is "none" or other values }, legend.spacing = unit(0, "cm"), legend.title = element_text(face = "italic", size = rel(1.3)), legend.text = element_text(size = rel(1.2)), plot.margin = unit(c(10, 5, 5, 5), "mm") ) } 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 = as.numeric(as.factor(conc_num)) - 1, # for legacy purposes conc_num_factor_factor = as.factor(conc_num) ) 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(.x, na.rm = TRUE), median = ~ median(.x, na.rm = TRUE), max = ~ ifelse(all(is.na(.x)), NA, max(.x, na.rm = TRUE)), min = ~ ifelse(all(is.na(.x)), NA, min(.x, na.rm = TRUE)), sd = ~ sd(.x, na.rm = TRUE), se = ~ sd(.x, na.rm = TRUE) / sqrt(n() - 1) ), .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, overlap_threshold = 2) { # Calculate total concentration variables total_conc_num <- length(unique(df$conc_num)) # Initial calculations 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 = Raw_Shift_L / bg_stats$sd_L, Z_Shift_K = Raw_Shift_K / bg_stats$sd_K, Z_Shift_r = Raw_Shift_r / bg_stats$sd_r, Z_Shift_AUC = Raw_Shift_AUC / bg_stats$sd_AUC, # Expected values 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, # Deltas 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 ) %>% group_modify(~ { # Perform linear models lm_L <- lm(Delta_L ~ conc_num_factor, data = .x) lm_K <- lm(Delta_K ~ conc_num_factor, data = .x) lm_r <- lm(Delta_r ~ conc_num_factor, data = .x) lm_AUC <- lm(Delta_AUC ~ conc_num_factor, data = .x) .x %>% mutate( lm_intercept_L = coef(lm_L)[1], lm_slope_L = coef(lm_L)[2], R_Squared_L = summary(lm_L)$r.squared, lm_Score_L = max_conc * lm_slope_L + lm_intercept_L, lm_intercept_K = coef(lm_K)[1], lm_slope_K = coef(lm_K)[2], R_Squared_K = summary(lm_K)$r.squared, lm_Score_K = max_conc * lm_slope_K + lm_intercept_K, lm_intercept_r = coef(lm_r)[1], lm_slope_r = coef(lm_r)[2], R_Squared_r = summary(lm_r)$r.squared, lm_Score_r = max_conc * lm_slope_r + lm_intercept_r, lm_intercept_AUC = coef(lm_AUC)[1], lm_slope_AUC = coef(lm_AUC)[2], R_Squared_AUC = summary(lm_AUC)$r.squared, lm_Score_AUC = max_conc * lm_slope_AUC + lm_intercept_AUC ) }) %>% ungroup() # Summary statistics for lm scores lm_means_sds <- calculations %>% group_by(across(all_of(group_vars))) %>% summarise( mean_mean_L = mean(mean_L, na.rm = TRUE), mean_lm_L = mean(lm_Score_L, na.rm = TRUE), sd_lm_L = sd(lm_Score_L, na.rm = TRUE), mean_mean_K = mean(mean_K, na.rm = TRUE), mean_lm_K = mean(lm_Score_K, na.rm = TRUE), sd_lm_K = sd(lm_Score_K, na.rm = TRUE), mean_mean_r = mean(mean_r, na.rm = TRUE), mean_lm_r = mean(lm_Score_r, na.rm = TRUE), sd_lm_r = sd(lm_Score_r, na.rm = TRUE), mean_mean_AUC = mean(mean_AUC, na.rm = TRUE), mean_lm_AUC = mean(lm_Score_AUC, na.rm = TRUE), sd_lm_AUC = sd(lm_Score_AUC, na.rm = TRUE) ) # Continue with gene Z-scores and interactions calculations <- calculations %>% left_join(lm_means_sds, by = group_vars) %>% group_by(across(all_of(group_vars))) %>% mutate( Z_lm_L = (lm_Score_L - mean_lm_L) / sd_lm_L, Z_lm_K = (lm_Score_K - mean_lm_K) / sd_lm_K, Z_lm_r = (lm_Score_r - mean_lm_r) / sd_lm_r, Z_lm_AUC = (lm_Score_AUC - mean_lm_AUC) / sd_lm_AUC ) # Build summary stats (interactions) interactions <- calculations %>% group_by(across(all_of(group_vars))) %>% summarise( Avg_Zscore_L = sum(Zscore_L, na.rm = TRUE) / first(num_non_removed_concs), Avg_Zscore_K = sum(Zscore_K, na.rm = TRUE) / first(num_non_removed_concs), Avg_Zscore_r = sum(Zscore_r, na.rm = TRUE) / first(num_non_removed_concs), Avg_Zscore_AUC = sum(Zscore_AUC, na.rm = TRUE) / first(num_non_removed_concs), # Interaction Z-scores Z_lm_L = first(Z_lm_L), Z_lm_K = first(Z_lm_K), Z_lm_r = first(Z_lm_r), Z_lm_AUC = first(Z_lm_AUC), # 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), # 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), # NG, DB, SM values NG = first(NG), DB = first(DB), SM = first(SM) ) # Creating the final calculations and interactions dataframes with only required columns for csv output calculations_df <- calculations %>% select( all_of(group_vars), conc_num, conc_num_factor, conc_num_factor_factor, N, NG, DB, SM, mean_L, median_L, sd_L, se_L, mean_K, median_K, sd_K, se_K, mean_r, median_r, sd_r, se_r, mean_AUC, median_AUC, sd_AUC, 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 ) interactions_df <- interactions %>% select( all_of(group_vars), NG, DB, SM, Avg_Zscore_L, Avg_Zscore_K, Avg_Zscore_r, Avg_Zscore_AUC, Z_lm_L, Z_lm_K, Z_lm_r, Z_lm_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 ) calculations_no_overlap <- calculations %>% # DB, NG, SM are same as in interactions, the rest may be different and need to be checked select(-any_of(c( "DB", "NG", "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", "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC" ))) # Use left_join to avoid dimension mismatch issues full_data <- calculations_no_overlap %>% left_join(interactions, by = group_vars) # Return full_data and the two required dataframes (calculations and interactions) return(list( calculations = calculations_df, interactions = interactions_df, full_data = full_data )) } generate_and_save_plots <- function(out_dir, filename, plot_configs) { message("Generating ", filename, ".pdf and ", filename, ".html") # Check if we're dealing with multiple plot groups plot_groups <- if ("plots" %in% names(plot_configs)) { list(plot_configs) # Single group } else { plot_configs # Multiple groups } # Open the PDF device once for all plots pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 16, height = 9) for (group in plot_groups) { static_plots <- list() plotly_plots <- list() grid_layout <- group$grid_layout plots <- group$plots for (i in seq_along(plots)) { config <- plots[[i]] df <- config$df # Set up aes mapping based on plot type aes_mapping <- if (config$plot_type == "bar" || config$plot_type == "density") { if (!is.null(config$color_var)) { aes(x = .data[[config$x_var]], fill = .data[[config$color_var]], color = .data[[config$color_var]]) } else { aes(x = .data[[config$x_var]]) } } else { if (!is.null(config$y_var) && !is.null(config$color_var)) { aes(x = .data[[config$x_var]], y = .data[[config$y_var]], color = .data[[config$color_var]]) } else if (!is.null(config$y_var)) { aes(x = .data[[config$x_var]], y = .data[[config$y_var]]) } else { aes(x = .data[[config$x_var]]) } } plot <- ggplot(df, aes_mapping) + theme_publication(legend_position = config$legend_position) # Add appropriate plot layer based on plot type plot <- switch(config$plot_type, "scatter" = generate_scatter_plot(plot, config), "box" = generate_boxplot(plot, config), "density" = plot + geom_density(), "bar" = plot + geom_bar(), plot # default (unused) ) # Add labels and title if (!is.null(config$title)) plot <- plot + ggtitle(config$title) if (!is.null(config$x_label)) plot <- plot + xlab(config$x_label) if (!is.null(config$y_label)) plot <- plot + ylab(config$y_label) if (!is.null(config$coord_cartesian)) plot <- plot + coord_cartesian(ylim = config$coord_cartesian) # Add error bars if specified if (!is.null(config$error_bar) && config$error_bar) { error_bar_color <- config$error_bar_params$color %||% "red" y_mean_col <- paste0("mean_", config$y_var) y_sd_col <- paste0("sd_", config$y_var) if (!is.null(config$error_bar_params$center_point)) { plot <- plot + geom_point(aes( x = .data[[config$x_var]], y = .data[[y_mean_col]]), color = error_bar_color, shape = 16) } # Use error_bar_params if provided, otherwise calculate from mean and sd if (!is.null(config$error_bar_params$ymin) && !is.null(config$error_bar_params$ymax)) { plot <- plot + geom_errorbar(aes( ymin = config$error_bar_params$ymin, ymax = config$error_bar_params$ymax), color = error_bar_color) } else { plot <- plot + geom_errorbar(aes( ymin = .data[[y_mean_col]] - .data[[y_sd_col]], ymax = .data[[y_mean_col]] + .data[[y_sd_col]]), color = error_bar_color) } } # Convert ggplot to plotly for interactive version plotly_plot <- suppressWarnings(plotly::ggplotly(plot)) # Store both static and interactive versions static_plots[[i]] <- plot plotly_plots[[i]] <- plotly_plot } # Print the plots to the PDF (one page per plot or in a grid) if (is.null(grid_layout)) { # Print each plot individually on separate pages if no grid layout is specified for (plot in static_plots) { print(plot) } } else { # Arrange plots in grid layout on a single page grid.arrange( grobs = static_plots, ncol = grid_layout$ncol, nrow = grid_layout$nrow ) } } # Close the PDF device after all plots are done dev.off() # Optional: Uncomment and save the interactive HTML version if needed # out_html_file <- file.path(out_dir, paste0(filename, ".html")) # message("Saving combined HTML file: ", out_html_file) # htmltools::save_html( # htmltools::tagList(plotly_plots), # file = out_html_file # ) } generate_scatter_plot <- function(plot, config) { # Define the points shape <- if (!is.null(config$shape)) config$shape else 3 size <- if (!is.null(config$size)) config$size else 1.5 position <- if (!is.null(config$position) && config$position == "jitter") { position_jitter(width = 0.3, height = 0.1) } 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) ) } } # Customize X-axis if specified if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) { # Check if x_var is factor or character (for discrete x-axis) if (is.factor(plot$data[[config$x_var]]) || is.character(plot$data[[config$x_var]])) { plot <- plot + scale_x_discrete( name = config$x_label, breaks = config$x_breaks, labels = config$x_labels ) } else { plot <- plot + scale_x_continuous( name = config$x_label, breaks = config$x_breaks, labels = config$x_labels ) } } # 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), 6, annotation$size), color = ifelse(is.null(annotation$color), "black", annotation$color) ) } } return(plot) } generate_boxplot <- function(plot, config) { # Convert x_var to a factor within aes mapping plot <- plot + geom_boxplot(aes(x = factor(.data[[config$x_var]]))) # Customize X-axis if specified if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) { # Check if x_var is factor or character (for discrete x-axis) if (is.factor(plot$data[[config$x_var]]) || is.character(plot$data[[config$x_var]])) { plot <- plot + scale_x_discrete( name = config$x_label, breaks = config$x_breaks, labels = config$x_labels ) } else { plot <- plot + scale_x_continuous( name = config$x_label, breaks = config$x_breaks, labels = config$x_labels ) } } return(plot) } generate_plate_analysis_plot_configs <- function(variables, df_before = NULL, df_after = NULL, plot_type = "scatter", stages = c("before", "after")) { plot_configs <- list() for (var in variables) { for (stage in stages) { df_plot <- if (stage == "before") df_before else df_after # Check for non-finite values in the y-variable df_plot_filtered <- df_plot %>% filter(is.finite(!!sym(var))) # Adjust settings based on plot_type plot_config <- list( df = df_plot_filtered, x_var = "scan", y_var = var, plot_type = plot_type, title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"), color_var = "conc_num_factor_factor", position = if (plot_type == "scatter") "jitter" else NULL, size = 0.2, error_bar = (plot_type == "scatter") ) # Add config to plots list plot_configs <- append(plot_configs, list(plot_config)) } } return(list(plots = plot_configs)) } generate_interaction_plot_configs <- function(df, type) { # Set group_vars based on the type (reference or deletion) if (type == "reference") { group_vars <- c("OrfRep", "Gene", "num") } else if (type == "deletion") { group_vars <- c("OrfRep", "Gene") } # Define the limits for the plots limits_map <- list( L = c(0, 130), K = c(-20, 160), r = c(0, 1), AUC = c(0, 12500) ) delta_limits_map <- list( L = c(-60, 60), K = c(-60, 60), r = c(-0.6, 0.6), AUC = c(-6000, 6000) ) stats_plot_configs <- list() stats_boxplot_configs <- list() delta_plot_configs <- list() # Overall statistics plots OrfRep <- first(df$OrfRep) # this should correspond to the reference strain for (plot_type in c("scatter", "box")) { for (var in names(limits_map)) { y_limits <- limits_map[[var]] # Use the pre-calculated lm intercept and slope from the dataframe lm_intercept_col <- paste0("lm_intercept_", var) lm_slope_col <- paste0("lm_slope_", var) # Ensure no NA or invalid values in lm_line calculations intercept_value <- mean(df[[lm_intercept_col]], na.rm = TRUE) slope_value <- mean(df[[lm_slope_col]], na.rm = TRUE) # Common plot configuration plot_config <- list( df = df, x_var = "conc_num_factor_factor", y_var = var, shape = 16, x_label = unique(df$Drug)[1], coord_cartesian = y_limits, x_breaks = unique(df$conc_num_factor_factor), x_labels = as.character(unique(df$conc_num)), lm_line = list( intercept = intercept_value, slope = slope_value ) ) # Add specific configurations for scatter and box plots if (plot_type == "scatter") { plot_config$plot_type <- "scatter" plot_config$title <- sprintf("%s Scatter RF for %s with SD", OrfRep, var) plot_config$error_bar = TRUE plot_config$error_bar_params <- list( color = "red", center_point = TRUE ) plot_config$position <- "jitter" # Append to scatter plot configurations stats_plot_configs <- append(stats_plot_configs, list(plot_config)) } else if (plot_type == "box") { plot_config$plot_type <- "box" plot_config$title <- sprintf("%s Boxplot RF for %s with SD", OrfRep, var) plot_config$position <- "dodge" # Boxplots don't need jitter, use dodge instead # Append to boxplot configurations stats_boxplot_configs <- append(stats_boxplot_configs, list(plot_config)) } } } # Delta interaction plots grouped_data <- df %>% group_by(across(all_of(group_vars))) %>% group_split() for (group_data in grouped_data) { OrfRep <- first(group_data$OrfRep) Gene <- first(group_data$Gene) num <- if ("num" %in% names(group_data)) first(group_data$num) else "" if (type == "reference") { OrfRepTitle <- paste(OrfRep, Gene, num, sep = "_") } else if (type == "deletion") { OrfRepTitle <- OrfRep } for (var in names(delta_limits_map)) { y_limits <- delta_limits_map[[var]] y_span <- y_limits[2] - y_limits[1] # Error bars WT_sd_value <- first(group_data[[paste0("WT_sd_", var)]], default = 0) # Z_Shift and lm values Z_Shift_value <- round(first(group_data[[paste0("Z_Shift_", var)]], default = 0), 2) Z_lm_value <- round(first(group_data[[paste0("Z_lm_", var)]], default = 0), 2) R_squared_value <- round(first(group_data[[paste0("R_Squared_", var)]], default = 0), 2) # NG, DB, SM values NG_value <- first(group_data$NG, default = 0) DB_value <- first(group_data$DB, default = 0) SM_value <- first(group_data$SM, default = 0) annotations <- list( list(x = 1, y = y_limits[2] - 0.2 * y_span, label = paste("ZShift =", Z_Shift_value)), list(x = 1, y = y_limits[2] - 0.3 * y_span, label = paste("lm ZScore =", Z_lm_value)), list(x = 1, y = y_limits[2] - 0.4 * y_span, label = paste("R-squared =", R_squared_value)), list(x = 1, y = y_limits[1] + 0.2 * y_span, label = paste("NG =", NG_value)), list(x = 1, y = y_limits[1] + 0.1 * y_span, label = paste("DB =", DB_value)), list(x = 1, y = y_limits[1], label = paste("SM =", SM_value)) ) # lm_line for delta plots lm_intercept_value <- first(group_data[[lm_intercept_col]], default = 0) lm_slope_value <- first(group_data[[lm_slope_col]], default = 0) plot_config <- list( df = group_data, plot_type = "scatter", x_var = "conc_num_factor_factor", y_var = var, x_label = unique(group_data$Drug)[1], title = paste(OrfRepTitle, Gene, num, sep = " "), coord_cartesian = y_limits, annotations = annotations, error_bar = TRUE, error_bar_params = list( ymin = 0 - (2 * WT_sd_value), ymax = 0 + (2 * WT_sd_value), color = "black" ), smooth = TRUE, x_breaks = unique(group_data$conc_num_factor_factor), x_labels = as.character(unique(group_data$conc_num)), ylim_vals = y_limits, lm_line = list( intercept = lm_intercept_value, slope = lm_slope_value ) ) delta_plot_configs <- append(delta_plot_configs, list(plot_config)) } } # Calculate dynamic grid layout grid_ncol <- 4 num_plots <- length(delta_plot_configs) grid_nrow <- ceiling(num_plots / grid_ncol) return(list( list(grid_layout = list(ncol = 2, nrow = 2), plots = stats_plot_configs), list(grid_layout = list(ncol = 2, nrow = 2), plots = stats_boxplot_configs), list(grid_layout = list(ncol = 4, nrow = grid_nrow), plots = delta_plot_configs) )) } generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FALSE, overlap_color = FALSE) { sd_bands <- c(1, 2, 3) plot_configs <- list() variables <- c("L", "K") # Adjust (if necessary) and rank columns for (variable in variables) { if (adjust) { df[[paste0("Avg_Zscore_", variable)]] <- ifelse(is.na(df[[paste0("Avg_Zscore_", variable)]]), 0.001, df[[paste0("Avg_Zscore_", variable)]]) df[[paste0("Z_lm_", variable)]] <- ifelse(is.na(df[[paste0("Z_lm_", variable)]]), 0.001, df[[paste0("Z_lm_", variable)]]) } df[[paste0("Rank_", variable)]] <- rank(df[[paste0("Avg_Zscore_", variable)]], na.last = "keep") df[[paste0("Rank_lm_", variable)]] <- rank(df[[paste0("Z_lm_", variable)]], na.last = "keep") } # Helper function to create a plot configuration create_plot_config <- function(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = TRUE) { num_enhancers <- sum(df[[zscore_var]] >= sd_band, na.rm = TRUE) num_suppressors <- sum(df[[zscore_var]] <= -sd_band, na.rm = TRUE) # Default plot config plot_config <- list( df = df, x_var = rank_var, y_var = zscore_var, plot_type = "scatter", title = paste(y_label, "vs. Rank for", variable, "above", sd_band), 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" ) if (with_annotations) { # Add specific annotations for plots with annotations plot_config$annotations <- list( list( x = median(df[[rank_var]], na.rm = TRUE), y = max(df[[zscore_var]], na.rm = TRUE) * 0.9, label = paste("Deletion Enhancers =", num_enhancers) ), list( x = median(df[[rank_var]], na.rm = TRUE), y = min(df[[zscore_var]], na.rm = TRUE) * 0.9, label = paste("Deletion Suppressors =", num_suppressors) ) ) } return(plot_config) } # Generate plots for each variable for (variable in variables) { rank_var <- if (is_lm) paste0("Rank_lm_", variable) else paste0("Rank_", variable) zscore_var <- if (is_lm) paste0("Z_lm_", variable) else paste0("Avg_Zscore_", variable) y_label <- if (is_lm) paste("Int Z score", variable) else paste("Avg Z score", variable) # Loop through SD bands for (sd_band in sd_bands) { # Create plot with annotations plot_configs[[length(plot_configs) + 1]] <- create_plot_config(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = TRUE) # Create plot without annotations plot_configs[[length(plot_configs) + 1]] <- create_plot_config(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = FALSE) } } # Calculate dynamic grid layout based on the number of plots grid_ncol <- 3 num_plots <- length(plot_configs) grid_nrow <- ceiling(num_plots / grid_ncol) # Automatically calculate the number of rows return(list(grid_layout = list(ncol = grid_ncol, nrow = grid_nrow), plots = plot_configs)) } generate_correlation_plot_configs <- function(df, correlation_stats) { # Define relationships for different-variable correlations relationships <- list( list(x = "L", y = "K"), list(x = "L", y = "r"), list(x = "L", y = "AUC"), list(x = "K", y = "r"), list(x = "K", y = "AUC"), list(x = "r", y = "AUC") ) plot_configs <- list() # Iterate over the option to highlight cyan points (TRUE/FALSE) highlight_cyan_options <- c(FALSE, TRUE) for (highlight_cyan in highlight_cyan_options) { for (rel in relationships) { # Extract relevant variable names for Z_lm values x_var <- paste0("Z_lm_", rel$x) y_var <- paste0("Z_lm_", rel$y) # Access the correlation statistics from the correlation_stats list relationship_name <- paste0(rel$x, "_vs_", rel$y) # Example: L_vs_K stats <- correlation_stats[[relationship_name]] intercept <- stats$intercept slope <- stats$slope r_squared <- stats$r_squared # Generate the label for the plot plot_label <- paste("Interaction", rel$x, "vs.", rel$y) # Construct plot config plot_config <- list( df = df, x_var = x_var, y_var = y_var, plot_type = "scatter", title = plot_label, annotations = list( list( x = mean(df[[x_var]], na.rm = TRUE), y = mean(df[[y_var]], na.rm = TRUE), label = paste("R-squared =", round(r_squared, 3)) ) ), smooth = TRUE, smooth_color = "tomato3", lm_line = list( intercept = intercept, slope = slope ), shape = 3, size = 0.5, color_var = "Overlap", cyan_points = highlight_cyan # Include cyan points or not based on the loop ) plot_configs <- append(plot_configs, list(plot_config)) } } return(list(plots = plot_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 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) 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", "conc_num_factor_factor"))$df_with_stats 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", "conc_num_factor_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) # For plotting (ggplot warns on NAs) df_na_stats_filtered <- df_na_stats %>% filter(if_all(all_of(summary_vars), is.finite)) 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", "conc_num_factor_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", "conc_num_factor_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", "conc_num_factor_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 list of plots corresponds to a file l_vs_k_plot_configs <- list( plots = list( list( df = df, x_var = "L", y_var = "K", plot_type = "scatter", tooltip_vars = c("OrfRep", "Gene", "delta_bg"), title = "Raw L vs K before quality control", color_var = "conc_num_factor_factor", error_bar = FALSE, legend_position = "right" ) ) ) frequency_delta_bg_plot_configs <- list( plots = list( list( df = df_stats, x_var = "delta_bg", y_var = NULL, plot_type = "density", title = "Density plot for Delta Background by [Drug] (All Data)", color_var = "conc_num_factor_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 = "Bar plot for Delta Background by [Drug] (All Data)", color_var = "conc_num_factor_factor", x_label = "Delta Background", y_label = "Count", error_bar = FALSE, legend_position = "right" ) ) ) above_threshold_plot_configs <- list( plots = list( list( df = df_above_tolerance, x_var = "L", y_var = "K", plot_type = "scatter", tooltip_vars = c("OrfRep", "Gene", "delta_bg"), title = paste("Raw L vs K for strains above Delta Background threshold of", round(df_above_tolerance$delta_bg_tolerance[[1]], 3), "or above"), color_var = "conc_num_factor_factor", position = "jitter", annotations = list( list( x = median(df_above_tolerance$L, na.rm = TRUE) / 2, y = median(df_above_tolerance$K, na.rm = TRUE) / 2, 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_filtered ) plate_analysis_boxplot_configs <- generate_plate_analysis_plot_configs( variables = summary_vars, df_before = df_stats, df_after = df_na_stats_filtered, 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( plots = list( list( df = df_na_l_outside_2sd_k_stats, x_var = "L", y_var = "K", plot_type = "scatter", title = "Raw L vs K for strains falling outside 2SD of the K mean at each Conc", color_var = "conc_num_factor_factor", position = "jitter", tooltip_vars = c("OrfRep", "Gene", "delta_bg"), annotations = list( list( x = median(df_na_l_outside_2sd_k_stats$L, na.rm = TRUE) / 2, y = median(df_na_l_outside_2sd_k_stats$K, na.rm = TRUE) / 2, label = paste("Total strains:", nrow(df_na_l_outside_2sd_k_stats)) ) ), error_bar = FALSE, legend_position = "right" ) ) ) delta_bg_outside_2sd_k_plot_configs <- list( plots = list( list( df = df_na_l_outside_2sd_k_stats, x_var = "delta_bg", y_var = "K", plot_type = "scatter", title = "Delta Background vs K for strains falling outside 2SD of the K mean at each Conc", color_var = "conc_num_factor_factor", position = "jitter", tooltip_vars = c("OrfRep", "Gene", "delta_bg"), annotations = list( list( x = median(df_na_l_outside_2sd_k_stats$delta_bg, na.rm = TRUE) / 2, y = median(df_na_l_outside_2sd_k_stats$K, na.rm = TRUE) / 2, label = paste("Total strains:", nrow(df_na_l_outside_2sd_k_stats)) ) ), error_bar = FALSE, legend_position = "right" ) ) ) message("Generating quality control plots in parallel") # future::plan(future::multicore, workers = parallel::detectCores()) future::plan(future::multisession, workers = 3) # generate 3 plots in parallel 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, filename = "L_vs_K_for_strains_2SD_outside_mean_K", plot_configs = l_outside_2sd_k_plot_configs), list(out_dir = out_dir_qc, filename = "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", "conc_num_factor_factor")) summary_stats_bg <- ss_bg$summary_stats write.csv(summary_stats_bg, file = file.path(out_dir, paste0("summary_stats_background_strain_", 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, conc_num_factor_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, conc_num_factor_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")) zscore_calculations_reference <- reference_results$calculations zscore_interactions_reference <- reference_results$interactions zscore_interactions_reference_joined <- reference_results$full_data 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", "Gene")) zscore_calculations <- deletion_results$calculations zscore_interactions <- deletion_results$interactions zscore_interactions_joined <- deletion_results$full_data # Writing Z-Scores to file write.csv(zscore_calculations_reference, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE) write.csv(zscore_calculations, file = file.path(out_dir, "zscore_calculations.csv"), row.names = FALSE) write.csv(zscore_interactions_reference, file = file.path(out_dir, "zscore_interactions_reference.csv"), row.names = FALSE) write.csv(zscore_interactions, file = file.path(out_dir, "zscore_interactions.csv"), row.names = FALSE) # Create interaction plots message("Generating reference interaction plots") reference_plot_configs <- generate_interaction_plot_configs(zscore_interactions_reference_joined, "reference") generate_and_save_plots(out_dir, "interaction_plots_reference", reference_plot_configs) message("Generating deletion interaction plots") deletion_plot_configs <- generate_interaction_plot_configs(zscore_interactions_joined, "deletion") generate_and_save_plots(out_dir, "interaction_plots", deletion_plot_configs) # Define conditions for enhancers and suppressors # TODO Add to study config? threshold <- 2 enhancer_condition_L <- zscore_interactions$Avg_Zscore_L >= threshold suppressor_condition_L <- zscore_interactions$Avg_Zscore_L <= -threshold enhancer_condition_K <- zscore_interactions$Avg_Zscore_K >= threshold suppressor_condition_K <- zscore_interactions$Avg_Zscore_K <= -threshold # Subset data enhancers_L <- zscore_interactions[enhancer_condition_L, ] suppressors_L <- zscore_interactions[suppressor_condition_L, ] enhancers_K <- zscore_interactions[enhancer_condition_K, ] suppressors_K <- zscore_interactions[suppressor_condition_K, ] # Save enhancers and suppressors message("Writing enhancer/suppressor csv files") write.csv(enhancers_L, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_L.csv"), row.names = FALSE) write.csv(suppressors_L, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_L.csv"), row.names = FALSE) write.csv(enhancers_K, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_K.csv"), row.names = FALSE) write.csv(suppressors_K, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_K.csv"), row.names = FALSE) # Combine conditions for enhancers and suppressors enhancers_and_suppressors_L <- zscore_interactions[enhancer_condition_L | suppressor_condition_L, ] enhancers_and_suppressors_K <- zscore_interactions[enhancer_condition_K | suppressor_condition_K, ] # Save combined enhancers and suppressors write.csv(enhancers_and_suppressors_L, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_and_suppressors_L.csv"), row.names = FALSE) write.csv(enhancers_and_suppressors_K, file = file.path(out_dir, "zscore_interaction_deletion_enhancers_and_suppressors_K.csv"), row.names = FALSE) # Handle linear model based enhancers and suppressors lm_threshold <- 2 # TODO add to study config? enhancers_lm_L <- zscore_interactions[zscore_interactions$Z_lm_L >= lm_threshold, ] suppressors_lm_L <- zscore_interactions[zscore_interactions$Z_lm_L <= -lm_threshold, ] enhancers_lm_K <- zscore_interactions[zscore_interactions$Z_lm_K >= lm_threshold, ] suppressors_lm_K <- zscore_interactions[zscore_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, "zscore_interactions_deletion_enhancers_lm_L.csv"), row.names = FALSE) write.csv(suppressors_lm_L, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_lm_L.csv"), row.names = FALSE) write.csv(enhancers_lm_K, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_lm_K.csv"), row.names = FALSE) write.csv(suppressors_lm_K, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_lm_K.csv"), row.names = FALSE) message("Generating rank plots") rank_plot_configs <- generate_rank_plot_configs( df = zscore_interactions_joined, is_lm = FALSE, adjust = TRUE ) generate_and_save_plots(out_dir = out_dir, filename = "rank_plots", plot_configs = rank_plot_configs) message("Generating ranked linear model plots") rank_lm_plot_configs <- generate_rank_plot_configs( df = zscore_interactions_joined, is_lm = TRUE, adjust = TRUE ) generate_and_save_plots(out_dir = out_dir, filename = "rank_plots_lm", plot_configs = rank_lm_plot_configs) message("Generating filtered ranked plots") rank_plot_filtered_configs <- generate_rank_plot_configs( df = zscore_interactions_filtered, 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) message("Generating filtered ranked linear model plots") rank_plot_lm_filtered_configs <- generate_rank_plot_configs( df = zscore_interactions_filtered, is_lm = TRUE, adjust = FALSE, overlap_color = TRUE ) generate_and_save_plots( out_dir = out_dir, filename = "rank_plots_lm_na_rm", plot_configs = rank_plot_lm_filtered_configs) message("Generating correlation curve parameter pair plots") correlation_plot_configs <- generate_correlation_plot_configs(zscore_interactions_filtered) generate_and_save_plots( out_dir = out_dir, filename = "correlation_cpps", plot_configs = correlation_plot_configs, ) }) }) } main() # For future simplification of joined dataframes # df_joined <- left_join(cleaned_df, summary_stats, by = group_vars, suffix = c("_original", "_stats"))