fig_large_scale.Rmd

---
title: "Data munging and Figure 2 preparation"
author: "Shyam Saladi (saladi@caltech.edu)"
date: "09/20/2016"
output: html_document
---

## Load libraries
Load all libraries necessary for subsequent data munging and plotting.
Some defaults are set here to make things easier throughout.
```{r libraries}
library(stringr)
library(tidyverse)
library(dplyrExtras)
library(magrittr)

library(gridExtra)
# `uconv` for unit conversion
library(datamart)
library(RColorBrewer)
library(cowplot)

# devtools::install("myutils")
library(myutils)
```

### Load model data for prediction function
```{r}
model_env <- new.env()
load("model.RData", model_env)
svmpredict <- function(feature_df)
    prediction_fn(feature_df, model_env$daley_preprocess_params,
                  model_env$svm_model)
```


# Figure 2 - Large Scale Expression Prediction
## NYCOMPS Data Preparation
This block analyzes the model's predictive performance within the NYCOMPS
dataset and prepares the data for plotting.
```{r nycomps_data_prep, cache = TRUE}
# Load features as well as well as scores
feat_col_spec = cols(.default = col_double(), title = col_character())

# these are actually not the designations
# they were initially just grouped this way for convenience
nycomps_features <- c("only_pos", "only_neg", "mixed") %>%
    paste0("large-scale/nycomps_", .) %>%
    lapply(function(x) {
        read_csv(paste0(x, ".fna.allstats.csv.gz"),
                 col_types = feat_col_spec) %>%
            bind_cols(read_csv(paste0(x, "_phobius.phobiuscterm"),
                               col_names = "cterm"))
        }) %>%
    bind_rows %>%
    separate(title, c("gi", "id", "name"), convert = TRUE) %>%
    select(-gi)

nycomps_features$score <- svmpredict(nycomps_features)
nycomps_features$score_nornass <-
    svmpredict(nycomps_features %>% select(-avgRNAss, -minRNAss,
                                           -q25RNAss, -q75RNAss, -maxRNAss))

nycomps_scores <- read_csv("~/ml-datasets/nycomps_targettrack/nycomps_061116.fna.allstats.csv") %>%
    separate(title, c("id1", "id2", "orfno")) %>%
    select(-id1, -id2) %>%
    mutate(orfno = as.integer(orfno))

nycomps_scores$score <- svmpredict(nycomps_scores)

nycomps_actions <- read_csv("large-scale/nycomps-outcomes-all-actions-raw-dump.csv") %>%
    select(-itemno, -wgwellcoord, -previouswgwellcoord, -orflabid, -labprotocolno,
           -timestamp, -srcorgno, -labfundno) %>%
    filter(str_detect(type, "Expressed") |
               str_detect(type, "0.5 liter prep") |
               str_detect(type, "LargeScale Production")) %>%
    left_join(nycomps_scores %>% select(orfno, score), by = "orfno")
# what is the difference between 0.5 L prep and "Passed to Large-Scale"
# there are still some that dont have scores (797 genes)

posneg_all <- nycomps_actions %>%
    mutate(outcome = -1) %>%
    mutate_rows(action == "Coomassie Negative", outcome = 0) %>%
    mutate_rows(action == "Coomassie Positive", outcome = 1) %>%
    filter(outcome != -1) %>%
    group_by(type) %>%
    do(calc_auc_roc(.))

genes <- nycomps_actions %>%
    mutate(outcome = -1) %>%
    mutate_rows(action == "Coomassie Negative", outcome = 0) %>%
    mutate_rows(action == "Coomassie Positive", outcome = 1) %>%
    filter(outcome != -1) %>%
    group_by(score) %>%
    summarize(pos = sum(outcome),
              total = length(outcome),
              outcome = ifelse(pos == 0, 0, 1)) %>%
    mutate_rows(pos == total, outcome = 2) %>%
    ungroup() %>%
    do(calc_auc_roc(.))

nycomps_all_outcomes <- read_csv("large-scale/nycomps_outcomes.csv") %>%
    rename(id = X1,
           sum_outcomes = sum) %>%
    gather(plasmid_name, plasmid_outcome, N:MSGC.9_LDAO) %>%
    filter(Analyze == 1L) %>%
    select(id, plasmid_name, plasmid_outcome) %>%
    left_join(nycomps_features %>%
                  select(id, cterm, score, score_nornass),
              by = "id") %>%
    # 0L here means that condition wasn't tested
    filter(!is.na(score), plasmid_outcome != 0L)

gene_outcome_func <- function(x) {
    if (all(x == -1L)) {
        "neg"
    } else if (all(x == 1L)) {
        "pos"
    } else {
        "mix"
    }
}

nycomps_gene_outcomes <- nycomps_all_outcomes %>%
    group_by(id, cterm, score, score_nornass) %>%
    summarise(
        attempts = n(),
        gene_outcome = gene_outcome_func(plasmid_outcome)) %>%
    mutate(gene_outcome = factor(gene_outcome, levels = c("neg", "mix", "pos"),
                                 ordered = TRUE),
           posmix_outcome = gene_outcome != "neg") %>%
    ungroup

# for a consistent scale throughout
nycomps_xscale <- nycomps_gene_outcomes$score %>%
    quantile(probs = c(.0025, .9975), names = FALSE)
nycomps_xscale[[1]] <- nycomps_xscale[[1]] - 1
nycomps_xscale[[2]] <- nycomps_xscale[[2]] + 1
nycomps_breaks <- seq(-6, 6, 2)
```

## nycomps plotting
```{r nycomps_trial_counts}
gene_summary <- nycomps_gene_outcomes %>%
    group_by(gene_outcome) %>%
    summarize(count = n()) %>%
    spread(gene_outcome, count) %>%
    unlist

p2_trial_counts <- nycomps_gene_outcomes %>%
    mutate(count = factor(attempts, levels = 1:max(attempts),
                          ordered = TRUE),
           gene_outcome = factor(gene_outcome,
                                 levels(gene_outcome) %>% rev)) %>%
    ggplot(aes(x = count, fill = gene_outcome)) +
    geom_bar(position = "stack", width = 0.5) +
    annotate_text(x = 5, y = 4000, label = "Genes:") +
    annotate_text(x = 5, y = 3600,
                  label = paste(gene_summary['pos'], "Expressed"),
                  color = "#A10F1C") +
    annotate_text(x = 5, y = 3200,
                  label = paste(gene_summary['mix'], "Mixed"),
                  color = "#636EBC") +
    annotate_text(x = 5, y = 2800,
                  label = paste(gene_summary['neg'], "Not Expressed"),
                  color = "#A4A5A5") +
    xlab("Number of Trials") +
    ylab("Number of Different Genes") +
    scale_fill_manual(
        values = c(neg = "#A4A5A5", mix = "#636EBC", pos = "#A10F1C")) +
    scale_y_continuous() +
    scale_x_discrete(drop = FALSE)

p2_trial_counts
```

## NYCOMPS without RNA Secondary Structure predictions
This block analyzes the effect of not including full sequence-length RNA
secondary structure calculations on predictive performance on the NYCOMPS
dataset.
```{r nycomps_auc, cache=TRUE}
# make combos
nycomps_auc <- nycomps_all_outcomes %>%
    mutate(cut = "Trials") %>%
    bind_rows(mutate(., plasmid_name = "all")) %>%
    bind_rows(nycomps_gene_outcomes %>%
                  mutate(cut = "Genes")) %>%
    select(-attempts) %>%
    bind_rows(mutate(., cterm = "all")) %>%
    mutate(gene_outcome = as.numeric(gene_outcome)) %>%
    select(-posmix_outcome) %>%
    gather(outcome_type, outcome, gene_outcome, plasmid_outcome) %>%
    filter(!is.na(outcome)) %>%
    gather(score_type, score, score, score_nornass) %>%
    mutate(cut = factor(cut),
           cterm = factor(cterm),
           plasmid_name = factor(plasmid_name),
           outcome_type = factor(outcome_type),
           score_type = factor(score_type)) %>%
    group_by(cut, cterm, plasmid_name, outcome_type, score_type) %>% 
    do(calc_auc_roc(., method = "auc")) %>%
    ungroup()

nycomps_gene_outcomes %>%
    rename(outcome = gene_outcome) %>%
    group_by(outcome != "mix") %>%
    select(-score_nornass) %>%
    do(calc_auc_roc(., method = "auc")) %>%
    ungroup() %>% View

nycomps_auc %>%
    write_csv("nycomps_all_auc.csv")

# Remove No RNAss AUC
nycomps_auc %<>%
    filter(score_type == "score") %>%
    select(-score_type)
```

```{r nycomps_posmix_ppv}
p2_posmix_ppv_prep <- nycomps_gene_outcomes %>%
    rename(outcome = gene_outcome) %>%
    calc_auc_roc(method = "roc") %>%
    group_by(outcome) %>%
    arrange(desc(thresholds)) %>%
    do(prep_for_polygon(.)) %>%
    as_tibble %>%
    ungroup() %>%
    mutate(color = ifelse(y_greater_min, outcome_percentile, 0),
           outcome = factor(outcome, levels = c("pos", "mix")))

p2_posmix_ppv_prep_best <- p2_posmix_ppv_prep %>%
   filter(type == "ppv") %>%
   distinct(ppv, min_ppv, thresholds, threshold_percentiles, outcome) %>%
   inner_join(nycomps_auc %>%
                  filter(cut == "Genes", cterm == "all") %>%
                  mutate(outcome = factor(outcome, levels = c(2, 3),
                                          labels = c("mix", "pos"))),
              by = c("outcome", "thresholds" = "best_thresh")) %>%
    mutate(outcome = factor(as.character(outcome), levels = c("pos", "mix")),
           percent_imp = (ppv - min_ppv)/min_ppv)

p2_posmix_ppv <- p2_posmix_ppv_prep %>%
    ggplot(aes(x = threshold_percentiles, y = ppv)) +
    geom_vline(xintercept = 0.75, color = "#bcbf11") +
    geom_hline(aes(yintercept = ppv), color = "#bcbf11",
               data = p2_posmix_ppv_prep %>%
                   filter(type != "min_y",
                          find_closest(threshold_percentiles, 0.75)) %>%
                   distinct(outcome, threshold_percentiles, ppv)) +
    geom_polygon(aes(group = id, fill = color), size = 0.01) +
    geom_path(size = 0.25,
              data = p2_posmix_ppv_prep %>%
                  filter(type != "min_y") %>%
                  distinct(outcome, threshold_percentiles, ppv)) +
    # geom_point(size = 0.5, data = p2_posmix_ppv_prep_best) +
    geom_hline(aes(yintercept = min_ppv),
               data = p2_posmix_ppv_prep %>%
                   filter(type != "min_y") %>%
                   distinct(outcome, min_ppv)) +
    geom_text(aes(x = 0.25, y = 0.40, label = text,
                  color = as.numeric(outcome)),
              data = data_frame(
                  outcome = factor(c("pos", "mix")),
                  text = c("Only positive", "At least one positive"))) +
    scale_x_continuous(limits = c(0, 1),
                       labels = function(x) x*100) +
    scale_y_continuous(limits = c(0, .6001), breaks = seq(0, 1, .2),
                       labels = function(x) x*100) +
    xlab(expression("Percentile IMProve score")) +
    ylab("Positive Predictive Value") +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C",
                           limits = c(0, 2), guide = FALSE) +
    scale_fill_continuous(low = "#A4A5A5", high = "#A10F1C",
                          limits = c(0, 100), guide = FALSE) +
    theme(strip.text = element_blank()) +
    facet_wrap(~outcome, ncol = 1)

p2_posmix_ppv
```

```{r}
p2_posmix_ppv_chg <- p2_posmix_ppv_prep %>% 
    ungroup() %>%
    filter(type == "ppv") %>%
    mutate(percent_imp = ppv/min_ppv) %>%
    select(threshold_percentiles, percent_imp, outcome_percentile, id) %>% 
    group_by(outcome_percentile, id) %>%
    filter(threshold_percentiles == max(threshold_percentiles)) %>%
    # spread(outcome, ppv) %>%
    # mutate(ppv_imp = pos/mix) %>%
    ggplot(aes(x = threshold_percentiles, y = percent_imp)) +
    geom_segment(x = 0.75, xend = 0.75, y = 0, yend = 3, color = "#bcbf11",
                 linetype = "dashed", size = 0.25) +
    geom_hline(aes(yintercept = ppv/min_ppv), color = "#bcbf11",
               data = p2_posmix_ppv_prep %>%
                   filter(type != "min_y", outcome == "pos",
                          find_closest(threshold_percentiles, 0.75)) %>%
                   distinct(outcome, threshold_percentiles, min_ppv, ppv)) +
    geom_line(aes(color = outcome_percentile,
                  group = outcome_percentile)) +
    # geom_point(size = 0.5, data = p2_posmix_ppv_prep_best) +
    xlab(expression("Percentile IMProve score")) +
    ylab("Fold Change in PPV") +
    scale_x_continuous(limits = c(0, 1), labels = function(x) x*100) +
    scale_y_continuous() +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C",
                           limits = c(0, 100), guide = FALSE) +
    coord_cartesian(ylim = c(1, 6))

p2_posmix_ppv_chg
```

```{r}
p2_posmix_ppv_prep %>%
    filter(find_closest(thresholds, 0.5) |
               find_closest(threshold_percentiles, c(0.75, 0.97, 0.90)) |
               find_closest(ppv, min_ppv * 2 %>% unique) |
               thresholds %in% p2_posmix_ppv_prep_best$thresholds,
           type == "ppv") %>%
    distinct(outcome, threshold_percentiles, thresholds, ppv, min_ppv) %>%
    arrange(desc(threshold_percentiles), outcome) %>%
    transmute(outcome,
              percentile_threshold = round(threshold_percentiles*100, 0),
              thresholds = specify_decimal(thresholds, 3),
              ppv, min_ppv,
              ppv_boost = round((ppv - min_ppv)*100, 1),
              ppv_fold_chg = round(ppv / min_ppv, 3),
              additional_expressed = (ppv - min_ppv) %>%
                  multiply_by(nrow(nycomps_gene_outcomes)) %>%
                  round(1)) %>%
    # do this last, so the calculations above retain precision
    mutate(ppv = round(ppv*100, 1),
           min_ppv = round(min_ppv*100, 1))
```

```{r nycomps_posmix_histogram}
percent_rank_sep <- function(x, all) {
    lapply(x, function(x) myutils:::perc_rank_single(all, x)) %>% unlist
}

p2_posmix_sec_xaxis <-
    dup_axis(~ .,
             breaks = quantile(nycomps_gene_outcomes$score,
                               probs = c(.005, .02, .15, .50, .85, .98, .995),
                                   #seq(0, 1, .25),
                               names = FALSE),
             labels = function(x) {
                 percent_rank_sep(x, nycomps_gene_outcomes$score) %>%
                     round(1)
                 })

p2_posmix_histogram_meds <- nycomps_gene_outcomes %>%
    group_by(posmix_outcome) %>%
    summarize(med = median(score))

p2_posmix_histogram_bw <- fd_bw(nycomps_gene_outcomes$score)

p2_posmix_histogram <- nycomps_gene_outcomes %>%
    mutate(orig = TRUE) %>%
    bind_rows(. , mutate(., orig = FALSE) %>% filter(gene_outcome == "pos")) %>%
    mutate_rows(gene_outcome == "pos" & orig != TRUE, gene_outcome = "mix") %>%
    ggplot(aes(x = score)) +
    geom_vline(xintercept = 0.5, color = "#bcbf11") +
    geom_histogram(aes(fill = as.numeric(gene_outcome), group = gene_outcome),
                   binwidth = p2_posmix_histogram_bw,
                   alpha = 1, position = "identity") +
    # geom_step_hist(color = "#A10F1C", fill = "#A10F1C",
    #                binwidth = p2_posmix_histogram_bw,
    #                data = nycomps_gene_outcomes %>%
    #                    filter(gene_outcome == "pos")) +
    # geom_point(size = 0.5, aes(x = thresholds, y = 0),
    #            data = p2_posmix_ppv_prep_best) +
    scale_x_continuous(breaks = nycomps_breaks, limits = nycomps_xscale,
                       sec.axis = p2_posmix_sec_xaxis) +
    scale_y_continuous() +
    scale_fill_continuous(limits = c(1, 3), low = "#A4A5A5", high = "#A10F1C",
                          guide = FALSE) +
    ylab("Count") +
    xlab(expression("IMProve score")) +
    theme(legend.key.size = unit(0.25, "cm"),
          legend.position = c(0.08, 0.85),
          axis.title.x.top = element_blank())

# # plot medians
# p2_posmix_histogram <- p2_posmix_histogram +
#     geom_segment(aes(x = med, y = 0, xend = med, yend = 10),
#                  data = p2_posmix_histogram_meds) +
#     geom_segment(aes(x = `FALSE`, y = 5, xend = `TRUE`, yend = 5),
#                  data = spread(p2_posmix_histogram_meds, posmix_outcome, med))

# for ratios
p2_posmix_ratio_data <- ggplot_build(p2_posmix_histogram)$data[[2]]  %>%
    select(fill, x, count) %>%
    spread(fill, count) %>%
    rename(pos = `#A10F1C`,
           mix = `#AD665D`,
           neg = `#A4A5A5`) %>%
    mutate(total = pos + neg,
           ratio = pos / neg,
           pos_pct = pos / total) %>%
    as_tibble
# p2_posmix_ratio_data <- ggplot_build(p2_posmix_histogram)$data[[1]] %>%
#     filter(fill == "#A4A5A5") %>%
#     transmute(x, neg = count) %>%
#     inner_join(ggplot_build(p2_posmix_histogram)$data[[2]] %>%
#                    select(x, count), by = "x") %>%
#     as_tibble() %>%
#     rename(pos = count) %>%
#     mutate(total = pos + neg,
#            ratio = pos / neg,
#            pos_pct = pos / total,
#            pos_pct_hist_scale = max(c(pos, neg)) * pos_pct)
ggplot(p2_posmix_ratio_data, aes(x = x)) +
    geom_point(aes(color = factor(total > 20)), y = -.01) +
    geom_line(aes(y = pos_pct), color = "#784F36") +
    xlab("IMProve Score") +
    ylab("Only Positive/Total") +
    scale_y_continuous(labels = function(x) x*100) +
    theme(legend.position = c(.25, .75))

p2_posmix_histogram <- p2_posmix_histogram +
    geom_line(aes(x = x, y = pos_pct * 550), color = "#784F36",
              data = filter(p2_posmix_ratio_data, total > 20)) +
    scale_y_continuous(sec.axis = sec_axis(
        breaks = seq(0, 100, 10),
        trans = ~ . / 5.5,
        name = 'Only Positive/Total')) +
    theme(legend.position = c(.25, .75),
          axis.text.y.right = element_text(color = "#784F36"),
          axis.title.y.right = element_text(color = "#784F36"))

p2_posmix_histogram %<>% ggplot_build
p2_posmix_histogram$layout$panel_ranges[[1]]$y.sec.labels[seq(2,6,2)] <- ""

p2_posmix_histogram %<>% ggplot_gtable

ggdraw(p2_posmix_histogram)
```


```{r nycomps_posmix_roc}
format_pos_count <- function(df) {
    df %>%
        arrange(pos_count) %>%
        summarize(pct = paste0(round(pos_count/count*100), '%', collapse = ", "),
                  count = count[[1]]) %$%
        paste0(pct, " / ", count)
}

p2_posmix_roc_label <- nycomps_auc %>%
    filter(cut == "Genes", cterm == "all") %>%
    arrange(pos_count) %>%
    mutate(auc = specify_decimal(auc*100),
           lower95 = specify_decimal(lower95*100),
           upper95 = specify_decimal(upper95*100),
           ci = paste0(auc, " (", lower95, "-", upper95, ")"),
           y = c(.13, .04))

p2_posmix_roc <- nycomps_gene_outcomes %>%
    rename(outcome = gene_outcome) %>%
    calc_auc_roc(method = "roc") %>%
    group_by(outcome) %>%
    split_groups() %>%
    ggplot() +
    geom_abline(slope = 1, intercept = 0) +
    geom_path(aes(x = 1 - specificities, y = sensitivities,
                  color = as.character(outcome))) +
    geom_path(aes(x = 1 - specificities, y = sensitivities),
              linetype = "dashed", color = "#636ebc",
              data = nycomps_gene_outcomes %>%
                  rename(outcome = gene_outcome) %>%
                  filter(outcome != "pos") %>%
                  calc_auc_roc(method = "roc")) +
    geom_path(aes(x = 1 - specificities, y = sensitivities),
              linetype = "dashed", color = "#A10F1C",
              data = nycomps_gene_outcomes %>%
                  rename(outcome = gene_outcome) %>%
                  filter(outcome != "mix") %>%
                  calc_auc_roc(method = "roc")) +
    annotate_text(x = .70, y = .22,
                  label = format_pos_count(p2_posmix_roc_label)) +
    geom_text(aes(x = .70, y = y, label = ci, color = as.character(outcome)),
              data = p2_posmix_roc_label) +
    scale_y_continuous(labels = function(x) x*100) +
    scale_x_continuous(labels = function(x) x*100) +
    xlab("False Positive Rate") +
    ylab("True Positive Rate") +
    scale_color_manual(values = c('2' = "#636ebc", '3' = "#A10F1C",
                                  'mix' = "#636ebc", 'pos' = "#A10F1C")) +
    theme(aspect.ratio = 1)

p2_posmix_roc

nycomps_gene_outcomes %>%
    rename(outcome = gene_outcome) %>%
    filter(outcome != "mix") %>%
    calc_auc_roc()

nycomps_gene_outcomes %>%
    rename(outcome = gene_outcome) %>%
    filter(outcome != "pos") %>%
    calc_auc_roc()

ggsave(p2_posmix_roc, filename = "plots/p2_roc.pdf", height = 4, width = 4)
```

```{r}
p_nycomps_cterm <- nycomps_auc %>%
    filter(cut == "Trials", plasmid_name != "all") %>%
    mutate(
        plasmid_name = factor(plasmid_name,
               levels = c("all",
                          "N", "C",
                          "C_LDAO", "MSGC.24",
                          "MSGC.28", "MSGC.7",
                          "MSGC.9", "MSGC.9_LDAO"),
               labels = c("All\nConditions",
                          "His-FLAG-\nTEV-(N)", "-TEV-His\n(C)",
                          "-TEV-His\n(C,LDAO)", "His-GST-\nTEV-(24)",
                          "-TEV-His\n(28)", "His-TEV-\n(7)",
                          "His-MBP-\nTEV-(9)", "His-MBP-\nTEV-(9, LDAO)")))

p_nycomps_cterm <- p_nycomps_cterm %>%
    # calculate label positions
    mutate(lab_y = .35) %>%
    mutate_rows(cterm != "all", lab_y = lab_y - .032) %>%
    mutate_rows(cterm == "NON-CYTOPLASMIC.", lab_y = lab_y - .032) %>%
    ggplot() +
    geom_hline(yintercept = 0.50) +
    ggplot2::geom_text(aes(x = plasmid_name, y = lab_y,
                  color = cterm, label = fmt_cnt(pos_count, count)),
                  size = 5*5/14) +
    geom_pointrange(aes(x = plasmid_name, ymin = upper95, y = auc, ymax = lower95,
                        group = cterm, color = cterm),
                    fatten = 0.5,
                    position = position_dodge(width = 0.5)) +
    ylab("AUC") +
    scale_color_manual(
        values = setNames(c("black", brewer.pal(5L, "Dark2")[1:2]),
                          c("all", "CYTOPLASMIC.", "NON-CYTOPLASMIC.")),
                       labels = c(`all` = "All C-termini",
                                  `CYTOPLASMIC.` = "Cytoplasmic",
                                  `NON-CYTOPLASMIC.` = "Periplasmic")) +
    scale_x_discrete(limits = p_nycomps_cterm %>%
                         filter(cterm == "all") %>%
                         arrange(desc(count)) %$% unlist(plasmid_name)) +
    scale_y_continuous(labels = function(x) x*100) +
    coord_cartesian(ylim = c(0.45, .885)) +
    theme(legend.position = c(.06, .8),
          legend.key.size = unit(0.25, "cm"),
          axis.title.x = element_blank(),
          axis.text.x = element_text(size = 5),
          plot.margin = unit(c(0, 0, 21, 0), "pt"))

p_nycomps_cterm %<>%
    ggplot_build %>% ggplot_gtable

p_nycomps_cterm$layout %<>%
    mutate_rows(name == "panel", clip = "off")

# hide key
for (idx in 4:11)
    p_nycomps_cterm$grobs[[15L]]$grobs[[1L]]$grobs[[idx]] <- zeroGrob()

# fix colors of C-terminal labels in legend
p_nycomps_cterm$grobs[[15L]]$grobs[[1L]]$grobs[[13L]]$gp$col <-
    brewer.pal(3, "Dark2")[[1]]
p_nycomps_cterm$grobs[[15L]]$grobs[[1L]]$grobs[[14L]]$gp$col <-
    brewer.pal(3, "Dark2")[[2]]

ggdraw(p_nycomps_cterm)
```

## Compile all panels into full figure 2
A diagram with arrows showing a summary of the outcomes (with the three arrows)
was added in Adobe Illustrator after the fact.
```{r nycomps_fig_compilation}
plasmid_summary <- nycomps_all_outcomes %>%
    group_by(plasmid_outcome) %>%
    summarise(count = n()) %>%
    spread(plasmid_outcome, count) %>%
    unlist

p2_main <- ggdraw() +
    annotate_text(x = 0.15, y = 0.95,
                  label = paste(nrow(nycomps_gene_outcomes), "genes")) +
    annotate_text(x = 0.15, y = 0.875,
                  label = paste(nrow(nycomps_all_outcomes), "expression trials")) +
    annotate_text(x = 0.1, y = 0.8,
                  label = paste(plasmid_summary["1"], "\nsuccesses")) +
    annotate_text(x = 0.2, y = 0.8,
                  label = paste(plasmid_summary["-1"], "\nfailures")) +
    draw_plot(plot_grid(p2_posmix_ppv, p2_posmix_ppv_chg, nrow = 1,
                        labels = LETTERS[3:4], label_size = 8),
              x = .28, y = 0.6775, width = 1 - .28, height = .3225) +
    draw_plot(p2_posmix_histogram,
              x = 1/3, y = .355, width = 2/3, height = .3225) +
    draw_plot(p2_trial_counts, x = 0, y = .37, width = .33, height = .35) +
    draw_plot(p2_posmix_roc, x = -.12, y = 0, width = .5, height = .35) +
    draw_plot(p_nycomps_cterm, x = .27, y = 0, width = .715, height = .38) +
    draw_plot_label(label = c("A", "B", "E", "F", "G"),
                    x = c(0, 0, .3, 0, .272),
                    y = c(1, .75, .68, .385, .385))

save_plot(filename = "plots/fig2.pdf",
          plot = p2_main,
          base_width = uconv(136, "mm", "in", "Length"))

ggdraw(p2_main)
```

```{r session}
save(nycomps_features, nycomps_all_outcomes, nycomps_gene_outcomes,
     nycomps_xscale, nycomps_breaks,
     file = "large-scale.RData", compression_level = 9)
sessionInfo()
```