fig_features.Rmd

---
title: "Data munging and figure 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 setup}
# if TRUE, then data will be written to file

library(datamart) # `uconv` for unit conversion

library(iterators)
library(tidyverse)
library(magrittr)
library(dplyrExtras)

library(caret)
library(pROC)

library(corrplot)
library(RColorBrewer)
library(viridis)

library(cowplot)

# devtools::install("myutils")
## ggcorr modified from GGally package
library(myutils)
```

### Load data from training process
```{r}
training_env <- new.env()
load("training.RData", training_env)
```

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

### Load NYCOMPS data
```{r}
nycomps_env <- new.env()
load("large-scale.RData", nycomps_env)
```

### Load small-scale data
```{r}
small_env <- new.env()
load("small-scale.RData", small_env)
```

### Define auc function for the subsequent calculations
```{r}
auc_w_ci <- partial(ci.auc, direction = "<",
                    method = "delong"#,
                    # boot.n = 1000,
                    # progress = "none",
                    # parallel = TRUE
                    )

transform_features <- function(ref_feat, weights) {
    # scale and center features
    feat_weights <- setNames(weights$weight, as.character(weights$old_name))
        
    trans_feat <- preprocess_vars(ref_feat, model_env$daley_preprocess_params)
    # reorder vector
    feat_weights <- feat_weights[colnames(trans_feat)]
    for (i in seq_along(trans_feat))
        set(trans_feat, j = i, value = trans_feat[[i]] * feat_weights[i])

    trans_feat[is.na(trans_feat)] <- 0

    trans_feat
}
```

```{r assessment_function}
# data is allstats with outcomes (grouped)
group_stat <- function(exp_data, feat_grp, SUMMARY_FUN) {
    feat_grp %>%
        do({
            feat_df <- .
            keep_feats <- feat_df[, 'old_name'] %>% unlist

            exp_data %>%
                do({
                    data_df <- .
                    data_df$score <- data_df %>%
                        select(one_of(keep_feats)) %>%
                        svmpredict()
                    
                    data_df %>%
                        do(SUMMARY_FUN(.))
                })
            
        })
}
```

## Load SVM weights and calculate correlations
This block loads a csv file of the model's weights and renames certain features
whose original identifiers could be more accurately written by the feature
calculation Python scripts.
```{r load_weights}
# name corrections (see smsaladi/ml-expression#16)
rename_feats <- function(df) {
    rename(df,
           avgCU_first5 = minCU_win5,
           avgCU_first10 = minCU_win10,
           nterm_hydOCT = nterm_hydGES,
           avgRONNTM1_2 = RONNlongestloop,
           RONNlongestloop = avgRONNTM1_2)
}

# load weights
svmweights <- model_env$svm_model %>%
    enframe(name = "old_name", value = "weight") %>%
    mutate(abs_weight = abs(weight),
           rel = weight / mean(weight)) %>%
    inner_join(read_csv("training/feature_key.csv"), by = "old_name") #%>%
    # mutate_rows(old_name == "X40deltaG", old_name = "40deltaG") %>%
    # mutate_rows(old_name == "X40freqens", old_name = "40freqens")
```

## Calculate feature correlations
Calculates the Pearson correlation coefficients between each feature within a
reference dataset, i.e. all forward predicted genomes.
```{r feature_correlation}
# Add feature groupings back to feature weights dataframe
# with manual category summary information (i.e. number of
# features that compose the grouping).
svmweights <- svmweights %>%
    group_by(category) %>%
    summarize(man_weight = sum(abs(weight)),
              man_count = n()) %>%
    right_join(svmweights, by = c("category")) %>%
    mutate(man_rank = dense_rank(-man_weight),
           man_alpha = LETTERS[man_rank])
```

# order of labels throughout figure
```{r}
category_color <- svmweights %>%
    group_by(prot_nuc, category) %>%
    summarize(weight = sum(weight)) %>%
    arrange(weight) %>%
    mutate_rows(prot_nuc == "Nucleotide",
                prot_nuc = brewer.pal(5, "Dark2")[[5]]) %>%
    mutate_rows(prot_nuc == "Protein",
                prot_nuc = brewer.pal(7, "Dark2")[[7]]) %$%
    setNames(prot_nuc, category)
```


# correlations between feature groups after adding by feature weights
```{r}
p_feat_category_corr <- nycomps_env$nycomps_features %>%
    # rename(`40deltaG` = X40deltaG,
    #        `40freqens` = X40freqens) %>%
    semi_join(nycomps_env$nycomps_gene_outcomes, "id")

p_feat_category_corr <- p_feat_category_corr %>%
    transform_features(svmweights) %>%
    bind_cols(p_feat_category_corr %>% select(id)) %>%
    gather(feature, value, -id) %>%
    left_join(svmweights %>% select(category, old_name),
              by = c("feature" = "old_name")) %>%
    group_by(id, category) %>%
    summarize(value = sum(value, na.omit = TRUE)) %>%
    spread(category, value) %>%
    ungroup

p_feat_category_corr <- p_feat_category_corr[, names(category_color)] %>%
    ggcorr(method = c("pairwise.complete.obs", "pearson"),
           hjust = 1, size = 2,
           diagColors = svmweights %>%
           distinct(prot_nuc, category) %>%
           rename(feature = category,
                  colorgrp = prot_nuc) %>%
           data.frame
           ) +
    labs(fill = "Pearson's r") +
    theme(legend.position = c(-.07, .83),
          legend.title = element_text(size = 6),
          legend.text = element_text(size = 5),
          legend.key.height = unit(5, "pt"),
          legend.key.width = unit(8, "pt"),
          plot.margin = unit(rep(0, 4), "pt")) +
    scale_color_manual(values = setNames(brewer.pal(7, "Dark2")[c(5, 7)],
                                         c("Nucleotide", "Protein")),
                       guide = FALSE) +
    scale_fill_gradient2(breaks = c(-1, 0, 1), limits = c(-1, 1),
        guide = guide_colourbar(title = "Pearson's r",
                                title.position = "top",
                                title.hjust = 0.5,
                                nbin = 100,
                                raster = FALSE,
                                direction = "horizontal",
                                ticks = FALSE))

p_feat_category_corr %<>% ggplot_build %>% ggplot_gtable
p_feat_category_corr$layout$clip <- "off"

ggdraw(p_feat_category_corr)
```


```{r}
geom_bar_def <-
    partial(geom_bar, position = "stack", stat = "identity",
            width = 0.5, size = 0.5, fill = "white", color = "black")

p_weight_man_sum <- svmweights %>%
    group_by(category) %>%
    summarize(count = n(),
              tot_weight = sum(weight)) %>%
    ggplot(aes(x = category, y = weight)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "lightgrey") +
    geom_bar_def(data = svmweights %>%
                     filter(weight > 0)) +
    geom_bar_def(data = svmweights %>%
                     filter(weight <= 0)) +
    geom_text(aes(y = 1.6, label = count), size = 5*5/14) +
    geom_point(aes(y = tot_weight), color = "red") +
    ylab("Weight") +
    scale_x_discrete(limits = names(category_color)) +
    coord_flip() +
    theme(legend.position = "none",
          axis.title.y = element_blank(),
          axis.text.y = element_text(color = category_color),
          plot.margin = unit(rep(0, 4), "pt"))

p_weight_man_sum %<>% ggplot_build %>% ggplot_gtable
p_weight_man_sum$layout$clip <- "off"

ggdraw(p_weight_man_sum)
```

## feature importances to e coli data
```{r training_imp}
training_data <- training_env$ecoli_training %>%
    left_join(training_env$daley_allstats, by = "id") %>%
    filter(normalized, !is.na(score)) %>%
    group_by(cterm)

boot_perm_partial <- partial(boot_perm, fn_prefix = "__",
                             cores = 4, count = 5050, force = TRUE)

ecoli_all <- training_data %>%
    group_stat(svmweights, boot_perm_partial)

ecoli_protnuc <- training_data %>%
    group_stat(svmweights %>% group_by(prot_nuc), boot_perm_partial)

ecoli_cate <- training_data %>%
    group_stat(svmweights %>% group_by(category), boot_perm_partial)

ecoli_all_auc <- training_data %>%
    group_stat(svmweights, calc_auc_roc)

ecoli_protnuc_auc <- training_data %>%
    group_stat(svmweights %>% group_by(prot_nuc), calc_auc_roc)

small_scale <- small_env$small_scale %>%
    group_by(dataset)

small_all <- group_stat(small_scale, svmweights, calc_auc_roc)

small_protnuc <- group_stat(small_scale,
                            svmweights %>% group_by(prot_nuc), calc_auc_roc)

small_cate <- group_stat(small_scale,
                         svmweights %>% group_by(category), calc_auc_roc)
```

```{r ec_cate}
p_ec_cate <- ecoli_all %>%
    mutate(category = "IMProve Score") %>%
    bind_rows(ecoli_cate) %>%
    ungroup %>%
    mutate(cterm = factor(cterm, levels = c("fluman", "out", "in"))) %>%
    ggplot(aes(x = category, y = t0, ymin = bca_low, ymax = bca_high)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
    geom_pointrange(aes(color = cterm, group = cterm),
                    position = position_dodge(width = 0.5),
                    size = 0.5, fatten = 0.5) +
    ylab("Spearman's rho") +
    scale_color_manual(values = setNames(brewer.pal(3, "Dark2"),
                                         c("in", "out", "fluman")),
                       guide = FALSE) +
    scale_y_continuous(expand = c(0, 0),
                       breaks = seq(-.6, .6, .2) %>% round(1),
                       limits = c(-.6, .6)) +
    scale_x_discrete(limits = names(category_color)) +
    coord_flip() +
    theme(axis.text.y = element_text(color = category_color),
          axis.title.y = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

ggdraw(p_ec_cate)
```

## plot AUC within NYCOMPS with subsets of features
```{r nycomps_subsets}
nycomps_outcomes <- nycomps_env$nycomps_all_outcomes %>%
    mutate(cut = "Trials") %>%
    bind_rows(mutate(., plasmid_name = "all")) %>%
    bind_rows(nycomps_env$nycomps_gene_outcomes %>%
                  mutate(cut = "Genes")) %>%
    select(-attempts) %>%
    mutate(gene_outcome = as.numeric(gene_outcome)) %>%
    select(-posmix_outcome, -score, -score_nornass ,-cterm) %>%
    gather(outcome_type, outcome, gene_outcome, plasmid_outcome) %>%
    filter(!is.na(outcome)) %>%
    mutate(cut = factor(cut),
           plasmid_name = factor(plasmid_name),
           outcome_type = factor(outcome_type)) %>%
    group_by(cut, plasmid_name, outcome_type) %>%
    left_join(nycomps_env$nycomps_features %>%
                  select(-name, -cterm, -score, -score_nornass), by = "id")

svmweights %>%
    group_by(category) %>%
    do({
        df <- .
        svmweights %>%
            anti_join(df, by = "category") %>%
            mutate(category = unique(df$category))
    }) 

nycomps_all_auc <- bind_rows(
    nycomps_outcomes %>%
    group_stat(svmweights, calc_auc_roc) %>%
    mutate(feat_type = "all"),
    nycomps_outcomes %>%
    group_stat(svmweights %>% group_by(prot_nuc), calc_auc_roc) %>%
    mutate(feat_type = "protnuc"),
    nycomps_outcomes %>%
    group_stat(svmweights %>% group_by(category), calc_auc_roc) %>%
    mutate(feat_type = "cate"),
    nycomps_outcomes %>%
    group_stat(svmweights %>%
                   group_by(category) %>%
                   do({
                       df <- .
                       svmweights %>%
                           anti_join(df, by = "category") %>%
                           mutate(category = unique(df$category))
                       }), calc_auc_roc) %>%
    mutate(feat_type = "cate_drop"),
    nycomps_outcomes %>%
    group_stat(svmweights %>% group_by(old_name), calc_auc_roc) %>%
    mutate(feat_type = "indiv"))
```

## indiv weights
```{r}
geom_pointrange_half <- partial(geom_pointrange, size = 0.5, fatten = 0.5)

feat_order <- nycomps_all_auc %>%
    filter(feat_type %in% c("all", "indiv"), outcome == 3) %>%
    mutate_rows(feat_type == "all", old_name = "IMProve Score") %>%
    arrange(old_name != "IMProve Score", desc(auc)) %$%
    factor(old_name, old_name)
    
p_indiv_feats <- nycomps_all_auc %>%
    filter(feat_type %in% c("all", "indiv")) %>%
    mutate_rows(feat_type == "all", old_name = "IMProve Score")

p_indiv_feats <- p_indiv_feats %>%
    ggplot(aes(x = old_name, y = auc, ymin = lower95, ymax = upper95,
               color = outcome)) +
    geom_hline(yintercept = 0.5, linetype = "dashed") +
    # Have to do this in separate steps to get layering correct
    geom_pointrange_half(data = filter(p_indiv_feats, outcome == 1)) +
    geom_pointrange_half(data = filter(p_indiv_feats, outcome == 2)) +
    geom_pointrange_half(data = filter(p_indiv_feats, outcome == 3)) +
    coord_flip() +
    ylab("AUC") +
    scale_x_discrete(limits = rev(feat_order)) +
    scale_y_continuous(expand = c(0, 0), labels = function(x) x*100) +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C", guide = FALSE) +
    theme(axis.title.y = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

p_indiv_feats
```

```{r}
category_color %<>%
    c(`IMProve Score` = "#000000")

make_incl_excl <- function(df)
    df %>%
    ggplot(aes(x = category, y = auc, ymin = lower95, ymax = upper95,
               color = outcome)) +
    geom_hline(yintercept = 0.5, linetype = "dashed") +
    # Have to do this in separate steps to get layering correct
    geom_pointrange_half(data = filter(df, outcome == 1)) +
    geom_pointrange_half(data = filter(df, outcome == 2)) +
    geom_pointrange_half(data = filter(df, outcome == 3)) +
    coord_flip() +
    ylab("AUC") +
    scale_x_discrete(limits = names(category_color)) +
    scale_y_continuous(expand = c(0, 0), labels = function(x) x*100) +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C", guide = FALSE) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_text(color = category_color),
          plot.margin = unit(rep(0, 4), "pt"))
    theme(legend.position = "none",
          axis.text.y = element_text(color = category_color),
          axis.title.y = element_blank(),
          strip.text = element_blank(),
          strip.background = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

p_category_incl <- nycomps_all_auc %>%
    filter(feat_type %in% c("cate", "all"),
           cut == "Genes" | plasmid_name != "all") %>%
    mutate_rows(feat_type == "all", category = "IMProve Score") %>%
    make_incl_excl() %>%
    ggdraw() +
    draw_plot_label("Only Single\nCategory", x = .23, y = .92,
                    fontface = "bold", size = 6, hjust = 0, vjust = 0)

p_category_excl <- nycomps_all_auc %>%
    filter(feat_type %in% c("cate_drop", "all"),
           cut == "Genes" | plasmid_name != "all") %>%
    mutate_rows(feat_type == "all", category = "IMProve Score") %>%
    make_incl_excl() %>%
    add(theme(axis.text.y = element_blank())) %>%
    ggdraw() +
    draw_plot_label("Excluding Single\nCategory", x = .02, y = .92,
                    fontface = "bold", size = 6, hjust = 0, vjust = 0)


p_category_spl <- nycomps_all_auc %>%
    filter(feat_type == "cate",
           cut == "Genes" | plasmid_name != "all",
           category %in% c("5' Codon Usage",
                           "SD-like Sites",
                           "Loop Length",
                           "Length/pI")) %>% 
    mutate_rows(feat_type == "all", category = "IMProve Score")

p_category_spl <- NULL %>%
    ggplot(aes(x = category, y = auc, ymin = lower95, ymax = upper95,
                     color = outcome, group = plasmid_name)) +
    geom_hline(yintercept = 0.5, linetype = "dashed") +
    # Have to do this in separate steps to get layering correct
    geom_pointrange_half(position = position_dodge(width = .7),
                         data = filter(p_category_spl, outcome == 1)) +
    geom_pointrange_half(data = filter(p_category_spl, outcome == 2)) +
    geom_pointrange_half(data = filter(p_category_spl, outcome == 3)) +
    coord_flip() +
    ylab("AUC") +
    scale_x_discrete(limits = p_category_spl$category %>% unique %>% rev) +
    scale_y_continuous(expand = c(0, 0), labels = function(x) x*100) +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C", guide = FALSE) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_text(color = category_color[p_category_spl$category %>% unique %>% rev]),
          plot.margin = unit(rep(0, 4), "pt"))


ggdraw(p_category_incl)
ggdraw(p_category_excl)

ggdraw(p_category_spl)
```


```{r}
p_weight_protnuc <- svmweights %>%
    group_by(prot_nuc) %>%
    summarize(
        count = n(),
        tot_weight = sum(weight)) %>%
    ggplot(aes(x = prot_nuc, y = weight)) +
    geom_bar_def(data = svmweights %>%
                     filter(weight > 0)) +
    geom_bar_def(data = svmweights %>%
                     filter(weight <= 0)) +
    geom_text(aes(y = 7, label = count),
              size = 5*5/14) +
    geom_point(aes(y = tot_weight), color = "red") +
    ylab("Weight") +
    coord_flip() +
    theme(axis.title.y = element_blank(),
          axis.text.y =
              element_text(color = brewer.pal(7, "Dark2")[c(5, 7)]),
          plot.margin = unit(rep(0, 4), "pt"))

p_weight_protnuc %<>% ggplot_build %>% ggplot_gtable
p_weight_protnuc$layout$clip <- "off"

ggdraw(p_weight_protnuc)
```

# protein or nucleotide auc plot
```{r}
p_auc_protnuc <- nycomps_all_auc %>%
    filter(feat_type %in% c("all", "protnuc"),
           cut == "Genes" | plasmid_name != "all") %>%
    mutate_rows(feat_type == "all", prot_nuc = "IMProve")

p_auc_protnuc <- p_auc_protnuc %>%
    ggplot(aes(x = prot_nuc, y = auc, ymin = lower95, ymax = upper95,
               color = outcome, group = plasmid_name)) +
    geom_hline(yintercept = 0.5, linetype = "dashed") +
    # Have to do this in separate steps to get layering correct
    geom_pointrange(size = 0.5, fatten = 0.5,
                    position = position_dodge(width = .7),
                    data = filter(p_auc_protnuc, outcome == 1)) +
    geom_pointrange(size = 0.5, fatten = 0.5,
                    data = filter(p_auc_protnuc, outcome == 2)) +
    geom_pointrange(size = 0.5, fatten = 0.5,
                    data = filter(p_auc_protnuc, outcome == 3)) +
    ylab("AUC") +
    scale_color_continuous(low = "#A4A5A5", high = "#A10F1C", guide = FALSE) +
    scale_y_continuous(expand = c(0, 0), limits = c(.16, .86),
                       labels = function(x) x*100) +
    scale_x_discrete(limits = c("Nucleotide", "Protein", "IMProve"),
                     labels = c(Nucleotide = "Nucleotide\nOnly",
                                Protein = "Protein\nOnly",
                                IMProve = "IMProve\nScore")) +
    coord_flip() +
    theme(legend.position = "none",
          axis.text.y = element_text(
              color = c(brewer.pal(7, "Dark2")[c(5, 7)], "#000000")),
          axis.title.y = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

p_auc_protnuc
```

### Load data from training process
```{r}
training_env <- new.env()
load("training.RData", training_env)
```

## Load Fluman, et al., aSD mutant data (Figure 6 in their paper)
```{r fluman_data_prep, cache=TRUE}
fluman_asd <- read_csv("forward-genomes/fluman_asd_mutant.fna.allstats.csv") %>%
    separate(title, c("gene")) %>%
    mutate(gene = tolower(gene),
           construct = "MUT")

fluman_asd <- training_env$daley_allstats %>%
    filter(id %in% fluman_asd$gene) %>%
    rename(gene = id) %>%
    mutate(construct = "WT") %>%
    bind_rows(fluman_asd)

fluman_asd$score <- fluman_asd %>%
    # rename(`40deltaG` = X40deltaG,
    #        `40freqens` = X40freqens) %>%
    svmpredict

fluman_asd <- read_csv("forward-genomes/fig6-analysis.csv") %>%
    filter(band == "FOLDED") %>%
    rename(outcome = `Area from Image`) %>%
    mutate(gene = tolower(gene)) %>%
    select(-`Percent (Just an identifier)`, -band) %>%
    left_join(fluman_asd, by = c("gene", "construct")) %>%
    gather(metric, value, -gene, -construct) %>%
    spread(construct, value) %>%
    mutate(rel_diff = rel_diff(MUT, WT))
```

# fluman asd plotting
```{r fluman_asd_plotting}
metrics_of_interest <- data.frame(
    row.names = c("SD-like Sites",  "Folded Protein", "IMProve Score"),
    colors = c("#e7298a", "#7570b3", "#e6ab02"),
    short = c("totalSDsites", "outcome", "score")) %>% t

p4_fluman_asd <- fluman_asd %>%
    filter(metric %in% metrics_of_interest['short',]) %>%
    mutate(gene = factor(gene,
                         levels = c("ygdd", "brnq", "ybjj"),
                         labels = c("YgdD", "BrnQ", "YbjJ")),
           metric = factor(metric,
                           levels = metrics_of_interest['short',],
                           labels = colnames(metrics_of_interest))) %>%
    ggplot(aes(x = gene, y = rel_diff, fill = metric)) +
    geom_bar(stat = "identity", position = "dodge") +
    ylab("Relative Difference") +
    scale_fill_manual(values = metrics_of_interest['colors',]) +
    scale_y_continuous(labels = c(-1, -.5, 0, .25)) +
    theme(axis.title.x = element_blank(),
          legend.title = element_blank(),
          legend.key = element_blank(),
          legend.key.size = unit(0.25, "cm"),
          legend.position = c(.60, .23),
          plot.margin = unit(rep(0, 4), "pt"))

p4_fluman_asd %<>% ggplotGrob()

# remove keystones and set legend text label colors
# for (idx in c(4L, 6L, 8L))
#     #p4_fluman_asd$grobs[[8L]]$grobs[[1L]]$grobs[[idx]] <- zeroGrob()
#     p4_fluman_asd[[1]][[15]][[1]][[1]][[idx]]
# 
# p4_fluman_asd$grobs[[8L]]$grobs[[1L]]$grobs[[9L]]$gp$col <-
#     metrics_of_interest['colors', 'SD-like Sites']
# 
# p4_fluman_asd$grobs[[8L]]$grobs[[1L]]$grobs[[10L]]$gp$col <-
#     metrics_of_interest['colors', 'Folded Protein']
# 
# p4_fluman_asd$grobs[[8L]]$grobs[[1L]]$grobs[[11L]]$gp$col <-
#     metrics_of_interest['colors', 'IMProve Score']

ggdraw(p4_fluman_asd)
```


## Figure importance plot
This block plots the entire figure relating to feature importances.
```{r plot_feature_importances}
plt_lab <- partial(draw_plot_label, size = 8)

p5_weight_importance <-
    ggdraw() +
    draw_plot(p_weight_man_sum, x = .34, y = .7, width = .66, height = .3) +
    draw_plot(p_feat_category_corr,
              x = 0.05, y = .71, width = .29, height = .29) +
    # draw_plot(
    #     plot_grid(p_category_incl, p_category_excl,
    #               nrow = 1L, rel_widths = c(2.6, 2)),
    #     x = 0, y = .21, width = 1, height = .49) +
    draw_plot(
        plot_grid(p_ec_cate, p_category_spl,
                  nrow = 1L, rel_widths = c(2.6, 2)),
        x = 0, y = .21, width = 1, height = .49) +
    draw_plot(
        plot_grid(p_weight_protnuc, p_auc_protnuc, p4_fluman_asd, nrow = 1L,
                  align = "h", rel_widths = c(2, 2, 1),
                  labels = LETTERS[5:7], label_size = 8, vjust = 0.5),
        x = 0, y = 0, width = 1, height = .2) +
    plt_lab(x = 0, y = 1, label = "A") +
    plt_lab(x = .31, y = 1, label = "B") +
    plt_lab(x = 0, y = .71, label = "C") +
    plt_lab(x = 0.53, y = .71, label = "D")
    # plt_lab(x = 0, y = .22, label = "e") +
    # plt_lab(x = 0.53, y = .22, label = "f")

p5_weight_importance %<>% ggplotGrob

# # dont clip e, f, g labels
# p5_weight_importance$grobs[[4L]]$children[[5L]]$children$layout$layout %<>%
#     mutate_rows(name == "panel", clip = "off")

save_plot(filename = "plots/fig6_fluman.pdf", plot = p5_weight_importance,
      base_width = uconv(183, "mm", "in", "Length"),
      base_height = 5.5)

ggdraw(p5_weight_importance)
```

```{r plot_big_corr}
p_feat_corr_labels <- svmweights %>%
    group_by(prot_nuc, category) %>%
    summarize(catweight = sum(weight)) %>%
    left_join(svmweights %>%
                  select(feature, category, weight),
              by = "category") %>%
    ungroup %>%
    arrange(catweight, weight) %>%
    mutate_rows(prot_nuc == "Nucleotide",
                prot_nuc_col = brewer.pal(5, "Dark2")[[5]]) %>%
    mutate_rows(prot_nuc == "Protein",
                prot_nuc_col = brewer.pal(7, "Dark2")[[7]]) %>%
    mutate_rows(feature == "X40deltaG", feature = "40deltaG") %>%
    mutate_rows(feature == "X40freqens", feature = "40freqens")

p_feat_corr_category <- p_feat_corr_labels %>%
    mutate(rank = row_number()) %>%
    group_by(prot_nuc, category) %>%
    (function(df)
        df %>%
         filter(weight == min(weight)) %>%
         rename(lower_rank = rank) %>%
         inner_join(df %>%
                        filter(weight == max(weight)) %>%
                        rename(upper_rank = rank),
                    by = "category")) %>%
    ungroup %>%
    transmute(prot_nuc = prot_nuc.x,
              category, lower_rank, upper_rank) %>%
    rowwise() %>%
    mutate(avg_rank = mean(c(lower_rank, upper_rank))) %>%
    ungroup

p_feat_corr <- nycomps_env$nycomps_features %>%
    rename(`40deltaG` = X40deltaG,
           `40freqens` = X40freqens) %>%
    select(-score, -score_nornass) %>%
    preprocess_vars(model_env$daley_preprocess_params) %>%
    rename_feats %>%
    cor(use = "pairwise.complete.obs") %>%
    as.data.frame %>%
    rownames_to_column(var = "x") %>%
    gather(y, cor, -x) %>%
    # adapted from ggcorr
    ggplot(aes()) +
    geom_tile(aes(x = x, y = y, fill = cor), color = "white") +
    geom_rect(aes(xmin = lower_rank - 0.5, ymin = lower_rank - 0.5,
                  xmax = upper_rank + 0.5, ymax = upper_rank + 0.5),
              color = "black", fill = NA,
              data = p_feat_corr_category) +
    geom_segment(aes(x = 89.9, y = lower_rank,
                     xend = 89.9, yend = upper_rank),
              color = "black",
              data = p_feat_corr_category) +
    geom_segment(aes(x = lower_rank, y = 89,
                     xend = upper_rank, yend = 89),
              color = "black",
              data = p_feat_corr_category) + 
    geom_text(aes(x = 89.9, y = avg_rank, label = category, color = prot_nuc),
              size = 5*5/14, nudge_x = 1, hjust = "left",
              data = p_feat_corr_category %>%
                  mutate(category = gsub(" ", "\n", category))) +
    geom_text(aes(x = avg_rank, y = 89, label = category, color = prot_nuc),
              size = 5*5/14,
              data = p_feat_corr_category) +
    labs(x = NULL, y = NULL, fill = "Pearson's r") +
    coord_equal() +
    scale_x_discrete(limits = p_feat_corr_labels$feature) +
    scale_y_discrete(limits = p_feat_corr_labels$feature) +
    scale_color_manual(values = p_feat_corr_labels %>%
                           distinct(prot_nuc_col, prot_nuc) %$%
                           setNames(prot_nuc_col, prot_nuc),
                       guide = FALSE) +
    scale_fill_gradient2(breaks = c(-1, 0, 1),
                         limits = c(-1, 1),
                         guide = guide_colourbar(title.position = "top",
                                                 title.hjust = 0.5,
                                                 nbin = 100,
                                                 raster = FALSE,
                                                 direction = "horizontal",
                                                 ticks = FALSE)) +
    theme(aspect.ratio = 1,
          axis.text.x = element_text(color = p_feat_corr_labels$prot_nuc_col,
                                     angle = 45, hjust = 1, vjust = 1),
          axis.text.y = element_blank(),
            #element_text(color = p_feat_corr_labels$prot_nuc_col),
          panel.background = element_blank(),
          legend.position = c(.5, 1.2),
          legend.key.height = unit(5, "pt"),
          legend.key.width = unit(10, "pt"),
          plot.margin = unit(rep(0, 4), "pt"))

p_feat_corr %<>% ggplot_build %>% ggplot_gtable
p_feat_corr$layout$clip <- "off"

ggdraw(p_feat_corr)
```

## Plot weights
This block plots the weights for the SVM model trained in the order of the
hclust clustering of the correlation matrix (plotted directly above).
```{r plot_weights}
p_svmweights <- ggplot(p_feat_corr_labels) +
    geom_bar(aes(x = feature, y = weight), stat = "identity") +
    scale_x_discrete(limits = p_feat_corr_labels$feature) +
    draw_label(x = -2, y = 0, label = "Weight", size = 6) +
    coord_flip() +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_text(color = p_feat_corr_labels$prot_nuc_col),
          axis.line.y = element_blank(),
          axis.ticks.y = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

p_svmweights %<>% ggplot_build %>% ggplot_gtable
p_svmweights$layout$clip <- "off"

ggdraw(p_svmweights)

#rm(p_feat_corr_labels, p_feat_corr_category)
```

### Plot all features
```{r}
# vertical alignment, based on tutorial from
# https://cran.r-project.org/web/packages/cowplot/vignettes/plot_grid.html
plot_grid(p_svmweights, p_feat_corr,
          rel_widths = c(1, 5), nrow = 1, align = "h") %T>%
    save_plot(filename = "plots/figS3_feature_correlations.pdf",
              plot = .,
              base_width = 9, base_height = 8) %>%
    ggdraw()
```


```{r ec_protnuc}
p_ec_kendall <- ecoli_all %>%
    mutate(prot_nuc = "IMProve") %>%
    bind_rows(ecoli_protnuc) %>%
    ungroup %>%
    mutate(cterm = factor(cterm, levels = c("fluman", "out", "in"))) %>%
    ggplot(aes(x = prot_nuc, y = t0, ymin = bca_low, ymax = bca_high)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
    geom_pointrange(aes(color = cterm, group = cterm),
                    position = position_dodge(width = 0.5),
                    size = 0.5, fatten = 0.5) +
    ylab("Spearman's rho") +
    scale_color_manual(values = setNames(brewer.pal(3, "Dark2"),
                                         c("in", "out", "fluman"))) +
    scale_y_continuous(expand = c(0, 0),
                       breaks = seq(-.6, .6, .2) %>% round(1),
                       limits = c(-.6, .6)) +
    scale_x_discrete(limits = c("Nucleotide", "Protein", "IMProve"),
                     labels = c(Nucleotide = "Nucleotide\nOnly",
                                Protein = "Protein\nOnly",
                                IMProve = "IMProve\nScore")) +
    coord_flip() +
    theme(legend.position = "none",
          axis.text.y = element_text(
              color = c(brewer.pal(7, "Dark2")[c(5, 7)], "#000000")),
          axis.title.y = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

ggdraw(p_ec_kendall)
```


```{r train_all_auc}
ecoli_protnuc_auc %<>%
    mutate(group = factor(cterm, levels = c("in", "out", "fluman"),
                          ordered = TRUE))

p_ec_auc_base <- ecoli_all_auc %>%
    ggplot(aes(x = outcome_percentile, y = auc, color = cterm)) +
    geom_hline(yintercept = .5, linetype = "dashed", color = "darkgrey") +
    #geom_line(linetype = "dotdash", alpha = 0.75) +
    scale_y_continuous(expand = c(0L, 0L), breaks = seq(0, 1, .20),
                       limits = c(0, 1),
                       labels = function(x) x*100) +
    scale_x_continuous(expand = c(0L, 0L), breaks = seq(0L, 100L, 25),
                       limits = c(0, 100)) +
    ylab("AUC") +
    xlab("Activity Percentile Rank") +
    scale_color_manual(values = setNames(brewer.pal(3, "Dark2"),
                                         c("in", "out", "fluman")),
                       guide = FALSE) +
    theme(strip.background = element_blank(),
          plot.margin = unit(rep(0, 4), "pt"))

p_ec_auc_prot <- p_ec_auc_base +
    geom_line(aes(x = outcome_percentile, y = auc, color = cterm),
              data = ecoli_protnuc_auc %>% filter(prot_nuc == "Protein")) +
    annotate("text", x = .1, y = .9, label = "Protein", vjust = "left",
              size = 7*5/14, fontface = "bold", color = "#A6761D")

p_ec_auc_nuc <- p_ec_auc_base +
    geom_line(aes(x = outcome_percentile, y = auc, color = cterm),
              data = ecoli_protnuc_auc %>% filter(prot_nuc == "Nucleotide")) +
    annotate("text", x = .1, y = .9, label = "Nucleotide", vjust = "left",
              size = 7*5/14, fontface = "bold", color = "#66A61E") +
    theme(axis.title.y = element_blank(),
         # axis.text.y = element_blank(),
          plot.margin = unit(c(0, 0, 0, 10), "pt"))

plot_grid(p_ec_auc_prot, p_ec_auc_nuc, nrow = 1)
```

### small scale feature plotting
```{r small_scale_plots}
p_small_protnuc <- small_all %>%
    group_by(dataset) %>%
    mutate(prot_nuc = "IMProve") %>%
    filter(auc == max(auc)) %>%
    # keep condition that corresponds to overall max AUC
    bind_rows(semi_join(small_protnuc, .,
                        by = c("dataset", "outcome"))) %>%
    ungroup() %>%
    mutate(dataset = factor(dataset,
                            levels = c("surade",
                                       "tmaritima",
                                       "hpylori",
                                       "lundstrom_ec",
                                       "tcross",
                                       "archaea"),
                            labels = c("Microbial\nSecondary\nTransporters",
                                       "T. maritima",
                                       "H. pylori",
                                       "Mammalian\nGPCRs\n(E. coli)",
                                       "M. tuberculosis",
                                       "Archaeal\nTransporters")),
           prot_nuc = factor(prot_nuc,
                             levels = c("Nucleotide", "Protein", "IMProve"))) %>%
    ggplot() +
    geom_hline(aes(yintercept = .5), linetype = "dashed", color = "darkgrey") +
    geom_pointrange(aes(x = dataset,
                        y = auc, ymin = lower95, ymax = upper95,
                        color = prot_nuc, group = prot_nuc),
                    position = position_dodge(width = 0.5),
                    size = 0.5, fatten = 0.5,) +
    ylab("AUC") +
    coord_flip() +
    scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, .25),
                       limits = c(0, 1), labels = function(x){x*100}) +
    scale_color_manual(
        values = setNames(c("#000000", brewer.pal(7, "Dark2")[c(5, 7)]),
                          c("IMProve", "Nucleotide", "Protein")),
        guide = FALSE) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_text(face = c("plain", "italic", "italic",
                                              "plain", "italic", "plain")),
          plot.margin = unit(rep(0, 4), "pt"))

ggdraw(p_small_protnuc)
```

```{r ec_small_compile}
p_s4_feats <-
    # ggdraw() +
    # plot_grid(p_ec_kendall, p_ec_auc_prot, p_ec_auc_nuc,
    #           align = "h", nrow = 1,
    #           labels = LETTERS[1:3], label_size = 8, vjust = 0.5) %>%
    # draw_plot(x = 0, y = .51, width = 1, height = 0.49) +
    # plot_grid(p_ec_cate, p_small_protnuc,
    #           nrow = 1, align = "h",
    #           labels = LETTERS[4:5], label_size = 8, vjust = 0.5) %>%
    # draw_plot(x = 0, y = 0, width = 1, height = 0.5)
    plot_grid(p_category_incl, p_category_excl,
              p_ec_kendall, p_small_protnuc,
              p_ec_auc_prot, p_ec_auc_nuc,
              ncol = 2, rel_heights = c(5, 4, 3),
              labels = LETTERS, label_size = 8)

# p_s4_feats %<>% ggplotGrob

save_plot(filename = "plots/figS4_feats.pdf", plot = p_s4_feats,
      base_width = uconv(183, "mm", "in", "Length"),
      base_height = 8)
    
ggdraw(p_s4_feats)
```

### Load data from training process
```{r, eval = FALSE}
training_env <- new.env()
load("training.RData", training_env)
```

# function to train a single model
```{r, eval = FALSE}
train_model(daley_to_train)

train_model <- function(df, tune = 4, ext = tune,
                        binary = "/ul/saladi/apps/svm_rank/svm_rank_learn") {
    
    data_fn <- df %>%
        select(-cterm, -id) %>%
        write_dataset(paste0("training/daley_gfp_phoa.", ext))
    model_fn <- paste0(data_fn, ".model")
    
    # train_model
    args <- paste("-c ", tune, " -t 0 -m 4000", data_fn, model_fn)
    system2(binary, args, stdout = FALSE)
    
    # read model
    model <- df %>%
        ungroup %>%
        select(-cterm, -id, -outcome, -grouping) %>%
        colnames %>%
        read_svmlight_model(model_fn, feat_names = .)
    
    # delete files
    file.remove(data_fn, model_fn)
    
    model
}

test_model <- function(model, data, xtrans = NULL) {
    data$score <- prediction_fn(data, xtrans = xtrans, model)

    data %>%
        group_by(cterm, grouping) %>%
        partition(cterm, grouping) %>%
        do(ConDis_fast(.$score, .$outcome)) %>%
        collect() %>%
        group_by(cterm) %>%
        summarize(total = sum(valid_count),
                  kendall = sum(con_dis)/total) %>%
        ungroup %>%
        bind_cols(model %>% t %>% as_data_frame %>% bind_rows(., ., .))
}
```

### Test each model
against the Training (Daley, Rapp, et al.) as well as Fluman, et al. datasets
```{r, test_models, eval = FALSE}
n_models <- 10

models <- seq.int(n_models) %>%
    foreach(idx = ., .combine = bind_rows, .multicombine = TRUE) %dopar% {
        df <- daley_to_train %>%
        #training_env$daley_to_train %>%
            sample_frac(0.95) %>%
            train_model() %>%
            test_model(data = daley_allstats,
                       xtrans = daley_preprocess_params)
}
```


```{r session}
sessionInfo()
```