ENH: Facet plots and normalize to fold change.

plots.R

@@ -2,23 +2,27 @@
 
 suppressPackageStartupMessages({
     library(ggplot2)
-    library(dplyr)
+    library(tidyverse)
     library(readr)
     library(forcats)
     library(varhandle)
 })
 
+## Do not generate Rplots.pdf
+pdf(NULL)
+
 wells_all <- read_csv("rnai-p_values.csv") %>%
     mutate(significant = p_value < 0.05,
-           plate = factor(plate) %>% fct_relevel("504"))
+           plate = factor(plate))
 binwidth <- 0.025
 
 theme_set(theme_bw())
 
 plot_save <- function(prefix, dir = "plots",
-                      extensions = c("png", "pdf", "svg")) {
+                      extensions = c("png", "pdf", "svg"),
+                      width = 8, height = 5) {
     invisible(sapply(file.path(dir, paste(prefix, extensions, sep = ".")),
-                     ggplot2::ggsave, width = 8, height = 5,
+                     ggplot2::ggsave, width = width, height = height,
                      units = "in"))
 }
 
@@ -30,52 +34,57 @@ ggplot(filter(wells_all, is.na(control)),
     xlab("% colocalization")
 plot_save("rnai-hist-signif")
 
-ggplot(filter(wells_all,
-              !is.na(control),
-              plate == "504"),
-       aes(coloc, fill = symbol)) +
-    geom_histogram(binwidth = binwidth) +
-    facet_wrap(~ control) +
-    ggtitle(paste("Histogram of control wells in plate 504")) +
-    xlab("% colocalization")
-
-ggplot(filter(wells_all, !is.na(control)) %>%
-       group_by(well) %>% mutate(ymin = min(coloc), ymax = max(coloc)),
-       aes(well, coloc, color = plate)) +
-    geom_errorbar(aes(ymin = ymin, ymax = ymax), color = "gray") +
-    geom_jitter(aes(shape = plate)) +
+## Normalize fold change to average "LacZ" colocalization value.
+##
+## FIXME: I'm sure there's a more mathemmatically sound way to
+## calculate fold change.
+coloc_lacz <- filter(wells_all, symbol == "LacZ") %>%
+    select(n_coloc, n) %>%
+    summarise(coloc = sum(n_coloc) / sum(n)) %>%
+    unlist
+data <- filter(wells_all, !is.na(control)) %>%
+    group_by(well) %>%
+    mutate(fold_change = coloc / coloc_lacz,
+           well_symbol = table(symbol) %>% sort() %>% head(1) %>%
+               names %>% paste(well, ., sep = "\n"))
+ggplot(data,
+       aes(x = symbol, y = fold_change)) +
+    geom_boxplot() +
+    geom_jitter(alpha = 0.5) +
     scale_color_brewer(palette = "Set1") +
-    ggtitle(paste("Scatter plot of control well",
-                  "colocalization across plates")) +
-    xlab("well") + ylab("% colocalization")
+    facet_wrap(~ plate) +
+    ggtitle("Control colocalization by plate") +
+    ylab("colocalization fold change")
+plot_save("controls-by_plate", width = 12, height = 8)
 
-ggplot(filter(wells_all, !is.na(control)) %>%
-       group_by(well) %>% mutate(ymin = min(p_value_control),
-                                 ymax = max(p_value_control),
-                                 ytext = plate %>%
-                                     lvls_revalue(seq(1.5, 1.1, -0.1) %>%
-                                                  as.character) %>%
-                                     unfactor),
-       aes(well, p_value_control, color = plate)) +
-    geom_errorbar(aes(ymin = ymin, ymax = ymax, linetype = control),
-                  color = "grey") +
-    geom_jitter(aes(color = plate), size = 3) +
-    geom_text(aes(label = symbol, y = ytext, color = plate)) +
+data_wells <- data %>%
+    select(symbol, plate, fold_change, p_value_control) %>%
+    gather(key = yaxis, value = y, fold_change, p_value_control)
+data_fold_change <- data_wells %>% filter(yaxis == "fold_change")
+data_p_value <- data_wells %>% filter(yaxis == "p_value_control")
+ggplot(data_wells, aes(symbol, y, color = plate)) +
+    facet_grid(yaxis ~ ., scale = "free") +
+    geom_boxplot(data = select(data_fold_change, -plate), color = "gray") +
+    geom_boxplot(data = select(data_p_value, -plate), color = "gray") +
+    geom_jitter(data = data_wells, alpha = 0.5) +
     scale_color_brewer(palette = "Set1") +
-    ggtitle(paste("Scatter plot of control well p-values across plates")) +
-    xlab("well") + ylab("p-value against LacZ")
+    ggtitle("Control colocalization and p-value across plates")
+plot_save("controls-all_plates", width = 12, height = 8)
+
+wells_all <- wells_all %>%
+    mutate(control = ifelse(is.na(control), "RNAi", control))
 
-ggplot(wells_all, aes(coloc, color = type)) +
+ggplot(wells_all, aes(coloc, color = control)) +
     geom_freqpoly(binwidth = binwidth) +
     scale_y_log10() +
-    facet_wrap(~ type) +
+    facet_wrap(~ control) +
     ggtitle("Histogram of all well types") +
     xlab("% colocalization")
 plot_save("all-hist")
 
-ggplot(wells_all, aes(coloc, color = type)) +
+ggplot(wells_all, aes(coloc, color = control)) +
     geom_density() +
-    facet_wrap(~ type) +
+    facet_wrap(~ control) +
     ggtitle("Density plot of wells to show distribution") +
     xlab("% colocalization")
 plot_save("all-density")