plots.R

## See the distribution of % Overlap vs. wells

suppressPackageStartupMessages({
    library(ggplot2)
    library(tidyverse)
    library(readr)
    library(forcats)
    library(varhandle)
})

## Do not generate Rplots.pdf
pdf(NULL)

wells_all <- read_csv("../results/tables/rnai-p_values.csv") %>%
    mutate(significant = p_value < 0.05,
           plate = factor(plate))
binwidth <- 0.025

theme_set(theme_bw())

plot_save <- function(prefix, dir = "../results/plots",
                      extensions = c("png", "pdf"),
                      width = 8, height = 5) {
    invisible(sapply(file.path(dir, paste(prefix, extensions, sep = ".")),
                     ggplot2::ggsave, width = width, height = height,
                     units = "in"))
}

ggplot(filter(wells_all, is.na(control)),
       aes(coloc, fill = significant)) +
    geom_histogram(binwidth = binwidth) +
    scale_fill_manual(values = c("grey", "blue")) +
    ggtitle("Histogram of RNAi wells with significant overlap (p-value < 5%)") +
    xlab("% colocalization")
plot_save("rnai-hist-signif")

## 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") +
    facet_wrap(~ plate) +
    ggtitle("Control colocalization by plate") +
    ylab("colocalization fold change")
plot_save("controls-by_plate", width = 12, height = 8)

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("Control colocalization and p-value across plates")
plot_save("controls-all_plates", width = 12, height = 8)

symbols_significant <- filter(wells_all, significant == TRUE)$symbol %>%
  unique
rnai_replicates <- filter(wells_all,
                          is.na(control),
                          symbol %in% symbols_significant) %>%
  group_by(symbol) %>%
  filter(n() > 1)  %>%
  mutate(fold_change = coloc / coloc_lacz,
         replicates = factor(n())) %>%
  gather(key = yaxis, value = y, fold_change, p_value_control)
rnai_fold_change <- rnai_replicates %>% filter(yaxis == "fold_change")
rnai_p_value <- rnai_replicates %>% filter(yaxis == "p_value_control")
ggplot(rnai_replicates, aes(x = symbol, y = y)) +
  facet_grid(yaxis ~ ., scale = "free") +
  geom_boxplot(data = rnai_fold_change, aes(linetype = replicates)) +
  geom_boxplot(data = rnai_p_value, aes(linetype = replicates)) +
  geom_jitter(alpha = 0.5, aes(color = plate)) +
  scale_color_brewer(palette = "Set1") +
  ggtitle("RNAi replicates") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_save("rnai-replicates", width = 12, height = 8)

wells_all <- wells_all %>%
    mutate(control = ifelse(is.na(control), "RNAi", control))

ggplot(wells_all, aes(coloc, color = control)) +
    geom_freqpoly(binwidth = binwidth) +
    scale_y_log10() +
    facet_wrap(~ control) +
    ggtitle("Histogram of all well types") +
    xlab("% colocalization")
plot_save("all-hist")

ggplot(wells_all, aes(coloc, color = control)) +
    geom_density() +
    facet_wrap(~ control) +
    ggtitle("Density plot of wells to show distribution") +
    xlab("% colocalization")
plot_save("all-density")