ソースを参照

Add cross-interaction linear model calculations

Bryan Roessler 7 ヶ月 前
コミット
9d23540157
1 ファイル変更106 行追加63 行削除
  1. 106 63
      qhtcp-workflow/apps/r/calculate_interaction_zscores.R

+ 106 - 63
qhtcp-workflow/apps/r/calculate_interaction_zscores.R

@@ -458,93 +458,126 @@ generate_and_save_plots <- function(output_dir, file_name, plot_configs, grid_la
 
 generate_scatter_plot <- function(plot, config) {
   
-  if (!is.null(config$delta_bg_point) && config$delta_bg_point) {
-    plot <- plot + geom_point(
-      shape = config$shape %||% 3,
-      size = config$size %||% 0.2
-    )
+  # 1. Determine Shape, Size, and Position for geom_point
+  shape <- if (!is.null(config$shape)) config$shape else 3
+  size <- if (!is.null(config$size)) {
+    config$size
+  } else {
+    if (!is.null(config$delta_bg_point) && config$delta_bg_point) 0.2
+    else if (!is.null(config$gene_point) && config$gene_point) 0.2
+    else 0.1
+  }
+  
+  position <- if (!is.null(config$delta_bg_point) && config$delta_bg_point) {
+    "identity"
   } else if (!is.null(config$gene_point) && config$gene_point) {
-    plot <- plot + geom_point(
-      shape = config$shape %||% 3,
-      size = config$size %||% 0.2,
-      position = "jitter"
-    )
+    "jitter"
   } else {
-    plot <- plot + geom_point(
-      shape = config$shape %||% 3,
-      position = if (!is.null(config$position) && config$position == "jitter") "jitter" else "identity",
-      size = config$size %||% 0.1
-    )
+    if (!is.null(config$position) && config$position == "jitter") "jitter" else "identity"
   }
   
-  # Add smooth line if specified
+  # 2. Add geom_point with determined parameters
+  plot <- plot + geom_point(shape = shape, size = size, position = position)
+  
+  # 3. Add Smooth Line if specified
   if (!is.null(config$add_smooth) && config$add_smooth) {
     if (!is.null(config$lm_line)) {
-      plot <- plot + geom_abline(intercept = config$lm_line$intercept, slope = config$lm_line$slope, color = "blue")
+      plot <- plot +
+        geom_abline(
+          intercept = config$lm_line$intercept,
+          slope = config$lm_line$slope,
+          color = "blue"
+        )
     } else {
-      plot <- plot + geom_smooth(method = "lm", se = FALSE, color = "blue")
+      plot <- plot +
+        geom_smooth(
+          method = "lm",
+          se = FALSE,
+          color = "blue"
+        )
     }
   }
   
-  # Add SD bands if specified
+  # 4. Add SD Bands if specified
   if (!is.null(config$sd_band_values)) {
     for (sd_band in config$sd_band_values) {
       plot <- plot +
-        annotate("rect", xmin = -Inf, xmax = Inf, ymin = sd_band, ymax = Inf, fill = "#542788", alpha = 0.3) +
-        annotate("rect", xmin = -Inf, xmax = Inf, ymin = -sd_band, ymax = -Inf, fill = "orange", alpha = 0.3) +
-        geom_hline(yintercept = c(-sd_band, sd_band), color = "gray")
+        annotate(
+          "rect",
+          xmin = -Inf, xmax = Inf,
+          ymin = sd_band, ymax = Inf,
+          fill = "#542788",
+          alpha = 0.3
+        ) +
+        annotate(
+          "rect",
+          xmin = -Inf, xmax = Inf,
+          ymin = -sd_band, ymax = -Inf,
+          fill = "orange",
+          alpha = 0.3
+        ) +
+        geom_hline(
+          yintercept = c(-sd_band, sd_band),
+          color = "gray"
+        )
     }
   }
   
-  # Add error bars if specified
+  # 5. 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
-    )
+    
+    plot <- plot +
+      geom_errorbar(
+        aes(
+          ymin = !!sym(y_mean_col) - !!sym(y_sd_col),
+          ymax = !!sym(y_mean_col) + !!sym(y_sd_col)
+        ),
+        alpha = 0.3
+      )
   }
   
-  # Add x-axis customization if specified
+  # 6. 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
-    )
+    plot <- plot +
+      scale_x_discrete(
+        name = config$x_label,
+        breaks = config$x_breaks,
+        labels = config$x_labels
+      )
   }
   
-  # Use coord_cartesian for zooming in without removing data outside the range
+  # 7. Apply coord_cartesian if specified
   if (!is.null(config$coord_cartesian)) {
     plot <- plot + coord_cartesian(ylim = config$coord_cartesian)
   }
   
-  # Use scale_y_continuous for setting the y-axis limits
+  # 8. 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
+  # 9. 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,
-        na.rm = TRUE
-      )
+      plot <- plot +
+        annotate(
+          "text",
+          x = annotation$x,
+          y = annotation$y,
+          label = annotation$label,
+          na.rm = TRUE
+        )
     }
   }
   
-  # Add titles and themes if specified
+  # 10. Add Title if specified
   if (!is.null(config$title)) {
     plot <- plot + ggtitle(config$title)
   }
   
+  # 11. Adjust Legend Position if specified
   if (!is.null(config$legend_position)) {
     plot <- plot + theme(legend.position = config$legend_position)
   }
@@ -552,7 +585,6 @@ generate_scatter_plot <- function(plot, config) {
   return(plot)
 }
 
-
 generate_box_plot <- function(plot, config) {
   plot <- plot + geom_boxplot()
   
@@ -794,30 +826,41 @@ generate_rank_plot_configs <- function(df_filtered, is_lm = FALSE, adjust = FALS
 }
 
 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()
 
-  variables <- c("r", "L", "K", "AUC")
+  for (rel in relationships) {
+    # Fit linear model
+    lm_model <- lm(as.formula(paste(rel$y, "~", rel$x)), data = df)
+    lm_summary <- summary(lm_model)
 
-  for (variable in variables) {
-    z_lm_var <- paste0("Z_lm_", variable)
-    avg_zscore_var <- paste0("Avg_Zscore_", variable)
-    lm_r_squared_col <- paste0("lm_R_squared_", variable)
-
-    configs[[length(configs) + 1]] <- list(
+    # Construct plot configuration
+    config <- list(
       df = df,
-      x_var = avg_zscore_var,
-      y_var = z_lm_var,
+      x_var = rel$x,
+      y_var = rel$y,
       plot_type = "scatter",
-      title = paste("Avg Zscore vs lm", variable),
-      color_var = "Overlap",
-      correlation_text = paste("R-squared =", round(df[[lm_r_squared_col]][1], 2)),
-      shape = 3,
-      geom_smooth = TRUE,
-      rect = list(xmin = -2, xmax = 2, ymin = -2, ymax = 2),  # To add the geom_rect layer
-      annotate_position = list(x = 0, y = 0),  # Position for the R-squared text
+      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 = 0, y = 0, label = paste("R-squared =", round(lm_summary$r.squared, 3)))
+      ),
+      add_smooth = TRUE,  # This flags that a geom_smooth layer should be added
+      lm_line = list(intercept = coef(lm_model)[1], slope = coef(lm_model)[2]),  # For direct geom_abline if needed
       legend_position = "right"
     )
+
+    configs[[length(configs) + 1]] <- config
   }
 
   return(configs)