all_genes <- read_tsv("data/refseq_genes.hg19.tsv",
                     col_names = c("transcript_id", "chromosome", "strand",
                                   "txStart", "txEnd", "gene_name"), skip = 1)

head(all_genes)

benchmark <- read_tsv("input/benchmark.tsv") %>% 
    left_join(all_genes, by = c("transcript_id", "gene_name")) %>% 
    # Define window around the TSS
  mutate(start = case_when(
    strand == "+" ~ as.numeric(txStart) - 2500,
    strand == "-" ~ as.numeric(txEnd) - 2500)) %>%
  mutate(end = case_when(
    strand == "+" ~ as.numeric(txStart) + 2500,
    strand == "-" ~ as.numeric(txEnd) + 2500)) %T>% 
    write_tsv("input/benchmark_coords.5kbp.tsv") %>% 
    makeGRangesFromDataFrame(keep.extra.columns = TRUE) %>% 
    keepStandardChromosomes()

benchmark %>%
  as.data.frame %>%
  dplyr::select(gene_name, transcript_id, chr = seqnames, start, end, switch) %>% 
  arrange(desc(switch), gene_name) %T>% 
  write_tsv("tables/tableS2.tsv") %>% 
  xtable()

refseq_unique <- all_genes %>% 
      # Keep the longest transcript per gene
      mutate(length = txEnd - txStart) %>%
      group_by(gene_name) %>%
      dplyr::slice(which.max(length)) %>% 
      ungroup() %>% 
      # Filter out genes where the length of the gene is less than half the
      # window size
      filter(length >= 2500)

getRefSeq <- function(window) {
    
  d <- window / 2
    
  refseq_unique %>%
      # Define window around the TSS
      mutate(start = case_when(
          strand == "+" ~ as.numeric(txStart) - d,
          strand == "-" ~ as.numeric(txEnd) - d)) %>%
      mutate(end = case_when(
          strand == "+" ~ as.numeric(txStart) + d,
          strand == "-" ~ as.numeric(txEnd) + d)) %T>% 
      write_tsv(paste0("input/benchmark_coords.", window, ".tsv")) %>% 
      makeGRangesFromDataFrame(keep.extra.columns = TRUE) %>% 
      keepStandardChromosomes(pruning.mode = "coarse")
}

refseq_5kbp <- getRefSeq(5000)

bench_s_k_raw <- mclapply(benchmark, safeCallSummary,
                       peaks = pk_k4me3,
                       metadata = meta,
                       mark = "H3K4me3",
                       filter = TRUE,
                       filter_columns = c("qValue", "signalValue"),
                       filter_thresholds = c(10, 4),
                       normalize = TRUE,
                       normalize_columns = c("qValue", "signalValue"),
                       summarize_columns = c("qValue", "signalValue"),
                       length = FALSE,
                       fraction = TRUE,
                       n = FALSE,
                       optimal_clusters = TRUE,
                       heatmap = FALSE,
                       BPPARAM = SerialParam(), mc.cores = n_cores) %>% 
  purrr::transpose()

bench_s_k <- bench_s_k_raw$result[map_lgl(bench_s_k_raw$error, is_null)] %>%
    bind_rows()

save(bench_s_k, file = "computations/bench_s_k.RData")

bench_s_d_raw <- mclapply(benchmark, safeCallSummary,
                               peaks = pk_dnase,
                               metadata = meta,
                               mark = "DNase",
                               filter = FALSE,
                               normalize = FALSE,
                               summarize_columns = NULL,
                               length = FALSE,
                               fraction = TRUE,
                               n = FALSE,
                               optimal_clusters = TRUE,
                               heatmap = FALSE,
                               BPPARAM = SerialParam(),
                               mc.cores = n_cores) %>% 
    purrr::transpose()

bench_s_d <- bench_s_d_raw$result[map_lgl(bench_s_d_raw$error, is_null)] %>%
    bind_rows()

save(bench_s_d, file = "computations/bench_s_d.RData")

bench_b_k_raw <- mclapply(benchmark,
                          safeCallBinary,
                          peaks = pk_k4me3,
                          metadata = meta,
                          filter = TRUE,
                          filter_columns = c("qValue","signalValue"),
                          filter_thresholds = c(10, 4),
                          reduce = TRUE,
                          gap = 300,
                          p = 0.4,
                          optimal_clusters = TRUE,
                          n_features = TRUE, mc.cores = n_cores) %>% 
  purrr::transpose()

bench_b_k <- bench_b_k_raw$result[map_lgl(bench_b_k_raw$error,
                                                is_null)] %>%
  bind_rows()

save(bench_b_k, file = "computations/bench_b_k.RData")

bench_b_d_raw <- mclapply(benchmark,
                          safeCallBinary,
                          peaks = pk_dnase,
                          metadata = meta,
                          filter = FALSE,
                          reduce = TRUE,
                          gap = 300,
                          p = 0.4,
                          optimal_clusters = TRUE,
                          n_features = TRUE, mc.cores = n_cores) %>% 
  purrr::transpose()

bench_b_d <- bench_b_d_raw$result[map_lgl(bench_b_d_raw$error, is_null)] %>%
  bind_rows()

save(bench_b_d, file = "computations/bench_b_d.RData")

bench_s_c_raw <- mclapply(benchmark,
                          safeCallSummary,
                          peaks = chmm_states,
                          metadata = meta,
                          mark = "State",
                          filter = FALSE,
                          normalize = FALSE,
                          summarize_columns = NULL,
                          length = FALSE,
                          fraction = TRUE,
                          n = FALSE,
                          heatmap = FALSE,
                          optimal_clusters = TRUE,
                          BPPARAM = SerialParam(),
                          mc.cores = n_cores) %>% 
  purrr::transpose()

bench_s_c <- bench_s_c_raw$result[map_lgl(bench_s_c_raw$error, is_null)] %>%
    bind_rows()

save(bench_s_c, file = "computations/bench_s_c.RData")

bench_b_c_raw <- mclapply(benchmark,
                          safeCallBinary,
                          peaks = chmm_states,
                          metadata = meta,
                          filter = FALSE,
                          reduce = FALSE,
                          p = 0.4,
                          optimal_clusters = TRUE,
                          n_features = TRUE,
                          mc.cores = n_cores) %>% 
  purrr::transpose()

bench_b_c <- bench_b_c_raw$result[map_lgl(bench_b_c_raw$error,
                                          is_null)] %>%
    bind_rows()

save(bench_b_c, file = "computations/bench_b_c.RData")

load("computations/bench_s_k.RData")
load("computations/bench_s_d.RData")
load("computations/bench_s_c.RData")
load("computations/bench_b_k.RData")
load("computations/bench_b_d.RData")
load("computations/bench_b_c.RData")

#' @param k, the number of query regions to make a prediction for
#' @param N, the number of times to draw
predictRandom <- function(k, N) {
    
    rep(k, N) %>%
        mclapply(runif, mc.cores = 2) %>% 
        data.frame() %>% 
        rowMeans()
}

# Repeat twice, to get one set of random predictions per method
draws <- replicate(2, predictRandom(120, 10000), simplify = FALSE)
save(draws, file = "computations/draws.RData")

load("computations/draws.RData")

bench_s_k <- bench_s_k %>%
    add_column(Random = draws[[1]], .after = "Consensus")

bench_b_d <- bench_b_d %>%
    add_column(Random = draws[[2]], .after = "Consensus")

predictShuffle <- function(clusters, N) {
  
  shuffle <- function(labels, clusters) {
    
    shuf <- sample(labels)
    contingency <- table(shuf, clusters)
    
    data.frame(
      ARI_shuffle          = mclust::adjustedRandIndex(clusters, shuf),
      NMI_shuffle          = chromswitch:::NMI(clusters, shuf),
      V_measure_shuffle    = chromswitch:::vMeasure(contingency)
    )
  }
  
  # Repeat N times
  replicate(N, shuffle(labels = meta$Condition, clust = unlist(clusters)),
            simplify = FALSE) %>%
    bind_rows() %>%
    colMeans()
}

shuffle_s_k <- bench_s_k %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

shuffle_s_d <- bench_s_d %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

shuffle_s_c <- bench_s_c %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

save(shuffle_s_k, shuffle_s_d, shuffle_s_c,
     file = "computations/shuffle_s.RData")

shuffle_b_k <- bench_b_k %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

shuffle_b_d <- bench_b_d %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

shuffle_b_c <- bench_b_c %>%
    dplyr::select(matches("^E[0-9]{3}")) %>%
    t() %>%
    as.data.frame %>%
    mclapply(predictShuffle, N = 10000, mc.cores = n_cores) %>% 
    as.data.frame() %>% 
    t() %>% 
    as.data.frame()

save(shuffle_b_c, shuffle_b_d, shuffle_b_k,
     file = "computations/shuffle_b.RData")

load("computations/shuffle_s.RData")
load("computations/shuffle_b.RData")

addShuf <- function(df) {
  df %>% 
    mutate(Consensus_shuffle =
               rowMeans(dplyr::select(., ARI_shuffle, NMI_shuffle, V_measure_shuffle)))
}

bench_s_k <- bind_cols(bench_s_k, shuffle_s_k) %>% addShuf()
bench_s_d <- bind_cols(bench_s_d, shuffle_s_d) %>% addShuf()
bench_s_c <- bind_cols(bench_s_c, shuffle_s_c) %>% addShuf()

bench_b_k <- bind_cols(bench_b_k, shuffle_b_k) %>% addShuf()
bench_b_d <- bind_cols(bench_b_d, shuffle_b_d) %>% addShuf()
bench_b_c <- bind_cols(bench_b_c, shuffle_b_c) %>% addShuf()

calcROC <- function(scores, method, input, smooth) {
    
    # Compute ROC info using each score as a predictor
    roc_out <- scores %>% 
        dplyr::select(
        matches("Consensus|Random")) %>% 
        map(~ pROC::roc(response = scores$switch,
                       predictor = as.numeric(.),
                       smooth = smooth))
    
    # Extract values to plot ROC curves
    roc_val <- roc_out %>%
        map(`[`, c("sensitivities", "specificities")) %>%
        map(as.data.frame) %>%
        map2_df(names(roc_out), function(df, index) mutate(df, index = index)) %>%
        mutate(TPR = sensitivities,
               FPR = 1 - specificities,
               Method = method,
               Input = input)

    return(roc_val)
}

benchmark_roc <- bind_rows(
  calcROC(bench_s_k, "Summary", "H3K4me3",  smooth = TRUE),
  calcROC(bench_s_d, "Summary", "DNase I",  smooth = TRUE),
  calcROC(bench_s_c, "Summary", "ChromHMM", smooth = TRUE),
  calcROC(bench_b_k, "Binary",  "H3K4me3",  smooth = TRUE),
  calcROC(bench_b_d, "Binary",  "DNase I",  smooth = TRUE),
  calcROC(bench_b_c, "Binary",  "ChromHMM", smooth = TRUE)) %>% 
  filter(grepl("Consensus|Random", index)) %>%
  filter(index != "Random" |
          (index == "Random") & (Method == "Binary") & (Input == "DNase I")) %>% 
  mutate(Input = ifelse(index == "Random", "NA", Input)) %T>% 
  write_tsv("tables/benchmark_roc.tsv")

roc_gg <- benchmark_roc %>% 
  makeLabel()

roc_gg %>%
    filter(Input != "DNase I") %>%
    filter(index != "Consensus_shuffle") %>%
    ggplot(aes(x = FPR, y = TPR, colour = Label, linetype = Label)) +
    geom_line(size = 1) +
    scale_colour_manual(name = "Strategy", values = pal_strategy) +
    scale_linetype_manual(name = "Strategy", values = linetype_strategy) +
    theme_min() +
    xlab("False Positive Rate") + ylab("True Positive Rate") +
    theme(legend.position = "bottom") +
    labs(colour = "") + labs(linetype = "") +
    guides(colour = guide_legend(ncol = 1))

roc_gg %>% 
    filter(index != "Random") %>% 
    ggplot(aes(x = FPR, y = TPR, colour = Label, linetype = Label)) +
    geom_line(size = 1) +
    scale_colour_manual(name = "Strategy", values = pal_strategy) +
    scale_linetype_manual(name = "Strategy", values = linetype_strategy2) +
    theme_min() +
    # theme(legend.position = "bottom") +
    xlab("False Positive Rate") + ylab("True Positive Rate") +
    theme(legend.position = "bottom") +
    facet_wrap(~ Input) +
    labs(colour = "") + labs(linetype = "") +
    guides(colour = guide_legend(ncol = 1))

calcAUC <- function(scores, method, input) {

    scores %>%
        dplyr::select(
        matches("Consensus|Random")) %>%
        map(~pROC::roc(response = scores$switch, predictor = as.numeric(.))) %>%
        map_df("auc") %>%
        mutate(Method = method,
               Input = input)
}

benchmark_auc <- bind_rows(
  calcAUC(bench_s_k, "Summary", "H3K4me3"),
  calcAUC(bench_s_d, "Summary", "DNase I"),
  calcAUC(bench_s_c, "Summary", "ChromHMM"),
  calcAUC(bench_b_k, "Binary", "H3K4me3"),
  calcAUC(bench_b_d, "Binary", "DNase I"),
  calcAUC(bench_b_c, "Binary", "ChromHMM")) %>% 
  select(Method, Input, Consensus, Consensus_shuffle, Random)

benchmark_auc

# A little manipulation to make the table for the supplementary

benchmark_auc %>%
  gather(Stat, Score, Consensus, Consensus_shuffle) %>%
  rowwise() %>%
  mutate(Method = ifelse(grepl("shuffle", Stat),
                         paste0(Method, " (shuffled)"), Method)) %>%
  ungroup() %>% 
  dplyr::select(Method, Input, Score) %>%
  spread(Input, Score) %T>%
  write_tsv("tables/benchmark_auc.tsv") %>% 
  xtable() # Generate latex code

getBest <- function(bench_out) {
    
    best_consensus <- pROC::roc(predictor = bench_out$Consensus,
                            response = bench_out$switch) %>%
    pROC::coords("best",
                 ret = c("threshold", "sensitivity",
                         "specificity", "accuracy", "ppv",
                 as.list = TRUE, drop = TRUE))
    
    best_threshold <- best_consensus["threshold"]
    return(best_threshold)
}

best_k <- getBest(bench_b_k)
best_k

grid <- expand.grid(gap = seq(200, 500, by = 50),
                    p_overlap = seq(0.2, 0.8, by = 0.1))

evaluateParams <- function(gap, p, best_threshold, ...) {
    
    out <- benchmark %>% 
        mclapply(safeCallBinary, gap = gap, p = p, ...) %>% 
        purrr::transpose()
    
    ok <- map_lgl(out$error, is_null)
    out <- out$result[ok] %>%
        bind_rows() %>% 
        mutate(pred = ifelse(Consensus >= best_threshold, 1, 0))
    
    # Now we can tally the results by comparing the predictions
    # to the known classes (switch or not)
    stats <- table(benchmark$switch[ok], out$pred)
    
    tn = fn = tp = fp = 0
    try(tn <- stats["0", "0"], silent = TRUE)
    try(fn <- stats["1", "0"], silent = TRUE)
    try(tp <- stats["1", "1"], silent = TRUE)
    try(fp <- stats["0", "1"], silent = TRUE)
    
    df <- data.frame(gap, p, tp, tn, fp, fn)
    return(list(predictions = out$pred, stats = df))
} 

safeEvalParams <- purrr::safely(evaluateParams)
grid_out <- Map(function(gap, p) safeEvalParams(gap, p,
                                                   best_threshold = best_k,
                                                   peaks = pk_k4me3,
                                                   metadata = meta,
                                                   filter = TRUE,
                                                   filter_columns = c("qValue", "signalValue"),
                                                   filter_thresholds = c(10, 4),
                                                   reduce = TRUE,
                                                   heatmap = FALSE,
                                                   BPPARAM = SerialParam(),
                                                   optimal_clusters = TRUE,
                                                   mc.cores = n_cores),
                   grid$gap, grid$p_overlap) %>%
    purrr::transpose()

save(grid_out, file = "computations/grid.RData")

load("computations/grid.RData")
makeGrid(grid_out, "H3K4me3")

load("tables/grid_scores.H3K4me3.RData")
plotGrid(grid_scores)

getSummaryFeatures <- function(query) {
  
  cols <- c("qValue", "signalValue")
  
  ft_mat <- pk_k4me3 %>% 
    filterPeaks(columns = cols,
                thresholds = c(10, 4)) %>% 
    normalizePeaks(columns = cols) %>% 
    retrievePeaks(metadata = meta, region = query) %>% 
    summarizePeaks(mark = "H3K4me3", cols = cols,
                   length = FALSE, fraction = TRUE, n = FALSE)
  
  return(ft_mat)
  
}

safeGetFeatures <- purrr::safely(getSummaryFeatures)

ft_s_k_raw <- benchmark %>% 
  mclapply(safeGetFeatures, mc.cores = n_cores) %>% 
  purrr::transpose()

ft_s_k <- ft_s_k_raw$result
save(ft_s_k, ft_s_k_raw,
     file = "computations/benchmark_features_s_k.RData")

getBinaryPeaks <- function(query) {
  
  pks <- pk_k4me3 %>% 
    retrievePeaks(metadata = meta, region = query) %>% 
    reducePeaks(gap = 300)
  
  return(pks)
  
}

safeGetBinaryPeaks <- purrr::safely(getBinaryPeaks)

pk_b_k_raw <- benchmark %>% 
  mclapply(safeGetBinaryPeaks, mc.cores = n_cores) %>% 
  purrr::transpose()
  
pk_b_k  <- pk_b_k_raw$result

# Extract the peaks from the localPeaks object
lpk_b_k <- lapply(pk_b_k, chromswitch:::lpkPeaks)

save(pk_b_k_raw, pk_b_k, lpk_b_k,
     file = "computations/benchmark_peaks_b_k.RData")

load("computations/benchmark_features_s_k.RData")
load("computations/benchmark_peaks_b_k.RData")

combos <- data.frame(
    n_other = c(2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 16, 16, 10, 12, 15, 16, 16),
    n_brain = c(2, 3, 4, 5, 6, 7, 7,  7,  7,  7,  7,  2,  3,  5,  4,  3,  4,  5)
)

combos

lvls_combos <- c(unlist(Map(function(x, y)
    paste0(x ," brain vs. ", y, " other"),
    combos$n_brain, combos$n_other)))

idx_other <- which(meta$Condition == "Other")
idx_brain <- which(meta$Condition == "Brain")

n_pos <- function(n_other, n_brain, idx_other) choose(length(idx_other), n_other) * choose(length(idx_brain), n_brain)

combos$n_pos <- unlist(Map(function(n_other, n_brain) n_pos(n_other, n_brain, idx_other),
                            combos$n_other, combos$n_brain))
combos

safeSubsampleFt <- purrr::safely(subsampleFromFeatures)
safeSubsamplePk <- purrr::safely(subsampleFromPeaks)

subsamp_s_raw <- mcMap(function(n_other, n_brain)
    safeSubsampleFt(ft_s_k, n_other, n_brain, idx_other),
    combos$n_other, combos$n_brain,
    mc.cores = n_cores) %>%
    purrr::transpose()

subsamp_s <- subsamp_s_raw$result
save(subsamp_s, subsamp_s_raw, file = "computations/subsamp_s.RData")

subsamp_b_raw <- mcMap(function(n_other, n_brain)
  safeSubsamplePk(pk_b_k, n_other, n_brain, idx_other),
    combos$n_other, combos$n_brain,
    mc.cores = n_cores) %>%
    purrr::transpose()

subsamp_b <- subsamp_b_raw$result
save(subsamp_b, subsamp_b_raw, file = "computations/subsamp_b.RData")

load("computations/subsamp_s.RData")
load("computations/subsamp_b.RData")

auc_s <- subsamp_s %>% map_df("auc") %>% 
  mutate(Method = "Summary", Input = "H3K4me3")

auc_b <- subsamp_b %>% map("scores") %>% 
  map_df("auc") %>% 
  mutate(Method = "Binary", Input = "H3K4me3")

# Simplify the facet labels
lvls_combos_new <- c(unlist(Map(function(x, y)
    paste0(x ," vs. ", y),
    combos$n_brain, combos$n_other)))
lvls_combos_new[1] <- "2 brain vs. 2 other"

# Calculate the mean AUC for each combination of n_brain, n_other
mean_auc <- bind_rows(auc_s, auc_b) %>% 
    rowwise() %>% 
    mutate(dataset = ifelse(n_brain == 2 && n_other == 2,
                            "2 brain vs. 2 other",
                            paste0(n_brain, " vs. ", n_other))) %>%
    ungroup() %>% 
    mutate(dataset = factor(dataset, levels = lvls_combos_new)) %>%
    dplyr::select(dataset, iteration, AUC = Consensus, Method, Input) %>%
    group_by(dataset, Method, Input) %>%
    summarise(mean_AUC = paste0("AUC: ", round(mean(AUC), 2))) %T>%
    write_tsv("tables/subsampling_auc.tsv")

roc_s <- subsamp_s %>% map_df("roc") %>%
  mutate(Method = "Summary", Input = "H3K4me3")

roc_b <- subsamp_b %>% map("scores") %>% 
  map_df("roc") %>% 
  mutate(Method = "Binary", Input = "H3K4me3")

roc_df <- bind_rows(roc_s, roc_b) %>% 
  rowwise() %>% 
  mutate(dataset = ifelse(n_brain == 2 && n_other == 2,
                          "2 brain vs. 2 other",
                          paste0(n_brain, " vs. ", n_other))) %>%
  ungroup() %>%
  mutate(dataset = factor(dataset, levels = lvls_combos_new)) %>%
  filter(index == "Consensus") %>% 
  dplyr::select(-index)

# Ensures that the iterations are plot in an order that looks nice
# with the colours of the lines
roc_df$iteration <- factor(roc_df$iteration,
                           levels = sort(unique(roc_df$iteration),
                                         decreasing = TRUE))

roc_df %T>%
    write_tsv("tables/subsampling_roc.tsv") %>% 
    filter(Method == "Summary" & Input == "H3K4me3") %>%
    plotSubsamp("red", "Summary", "H3K4me3")

roc_df %>%
    filter(Method == "Binary" & Input == "H3K4me3") %>%
    plotSubsamp("blue", "Binary", "H3K4me3")

window_genes <- c("GRIK2", "CHRM1", "TTYH1")

windows <- refseq_5kbp %>%
    as.data.frame %>%
    filter(gene_name %in% window_genes) %>%
    select(-start, -end)

get_window <- function(window_size, genes) {

    d <- window_size / 2

    genes %>%
        # Define window around the TSS
        mutate(start = case_when(
            strand == "+" ~ as.numeric(txStart) - d,
            strand == "-" ~ as.numeric(txEnd) - d)) %>%
        mutate(end = case_when(
            strand == "+" ~ as.numeric(txStart) + d,
            strand == "-" ~ as.numeric(txEnd) + d)) %>%
         mutate(gene_name = paste0(gene_name, "_", window_size/1000, "kbp"))
}

# Window sizes to test, in base pairs
sizes <- c(seq(1000, 10000, by = 1000),
           seq(10000, 50000, by = 10000),
           100000, 250000, 500000)

windows_gr <- lapply(sizes, get_window, windows) %>% bind_rows() %>%
  makeGRangesFromDataFrame(keep.extra.columns = TRUE)

win_s <- callSummary(windows_gr,
         peaks = pk_k4me3,
         metadata = meta,
         mark = "H3K4me3",
         filter = TRUE,
         filter_columns = c("qValue", "signalValue"),
         filter_thresholds = c(10, 4),
         normalize = TRUE,
         normalize_columns = c("qValue", "signalValue"),
         summarize_columns = c("qValue", "signalValue"),
         length = FALSE,
         fraction = TRUE,
         n = FALSE,
         optimal_clusters = TRUE,
         heatmap = TRUE,
         titles = paste0(windows_gr$gene_name, "_summary"),
         outdir = "figures/heatmaps")

win_s

save(win_s, file = "computations/windowsize_s.RData")

win_b <- callBinary(query = win_to_test,
                          peaks = pk_k4me3,
                          metadata = meta,
                          filter = TRUE,
                          filter_columns = c("qValue","signalValue"),
                          filter_thresholds = c(10, 4),
                          reduce = TRUE,
                          gap = 300,
                          p = 0.4,
                          optimal_clusters = TRUE,
                          heatmap = TRUE,
                          titles = paste0(windows_gr$gene_name, "_binary"),
                          outdir = "figures/heatmaps",
                          n_features = TRUE)

win_b

save(win_b, file = "computations/windowsize_b.RData")

load("computations/windowsize_s.RData")
load("computations/windowsize_b.RData")

windows <- bind_rows(mutate(win_s, method = "Summary"), mutate(win_b, method = "Binary")) %>% 
  select(gene_name, Consensus, method) %>% 
  separate(gene_name, into = c("gene_name", "window_size"), sep = "_") %>% 
  separate(window_size, into = c("window_size", "drop"), sep = "k") %>% 
  select(-drop) %>% 
  mutate(window_size = as.numeric(window_size)) %T>%
  write_tsv("tables/window_size.tsv")

# Plot TTYH1, "zoomed" in in the 0-50kbp range and in the 0-500kbp range
ttyh1_short <- windows %>% 
  filter(gene_name == "TTYH1") %>% 
  filter(window_size <= 50) %>% 
  mutate(range = "50kbp")

ttyh1_long <- windows %>% 
  filter(gene_name == "TTYH1") %>% 
  filter(window_size <= 500) %>% 
  mutate(range = "500kbp")

bind_rows(ttyh1_short, ttyh1_long) %>% 
  mutate(range = factor(range, levels = c("50kbp", "500kbp"))) %>% 
  ggplot(aes(x = window_size, y = Consensus)) +
  geom_point(aes(colour = method)) +
  geom_line(aes(colour = method), alpha = 0.8, size = 0.5) +
  scale_colour_manual(name = "Method", values = pal_methods) +
  xlab("Window Size (kbp)") + ylab("Consensus score") +
  facet_wrap(~ range, ncol = 1, scales = "free_x") +
  ggmin::theme_min() +
  theme(legend.position = c(0.9, 0.33)) +
  theme(strip.background = element_blank(), strip.text.x = element_blank())

keep_tissue <- c("Amygdala", "Anterior cingulate cortex", "Basal ganglia", 
                 "Cerebella hemisphere", "Cerebellum", "Cortex",
                 "Frontal cortex",
                 "Hippocampus", "Hypothalamus", "Substantia nigra", "Aorta",
                 "Sigmoid colon", "Transverse colon", "Esophagus",
                 "Heart (atrial appendage)", "Heart (left ventricle)",
                 "Kidney", "Liver", "Lung", "Skeletal muscle",
                 "Pancreas", "Small intestine", "Stomach")

brain <-  c("Amygdala", "Anterior cingulate cortex", "Basal ganglia", 
            "Cerebella hemisphere", "Cerebellum", "Cortex", "Frontal cortex",
            "Hippocampus", "Hypothalamus", "Substantia nigra")

other <- keep_tissue[!(keep_tissue %in% brain)]

gtex_rpkm <- read_tsv(
      "data/GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_median_rpkm.gct",
      skip = 2) %>%
    dplyr::rename(gene_name = Description) %>%
    unite(col = "unite", sep = "=", Name, gene_name) %>%
    gather(key = "Tissue", value = "mean_RPKM", 2:length(.)) %>% 
    mutate(Tissue = recode(.x = Tissue,
        `Brain - Amygdala` = "Amygdala",
        `Brain - Anterior cingulate cortex (BA24)` = "Anterior cingulate cortex",
        `Brain - Caudate (basal ganglia)` = "Basal ganglia",
        `Brain - Cerebellar Hemisphere` = "Cerebellar hemisphere",
        `Brain - Cerebellum` = "Cerebellum",
        `Brain - Cortex` = "Cortex",
        `Brain - Frontal Cortex (BA9)` = "Frontal cortex",
        `Brain - Hippocampus` = "Hippocampus",
        `Brain - Hypothalamus` = "Hypothalamus",
        `Brain - Nucleus accumbens (basal ganglia)` = "Nucleus accumbens",
        `Brain - Putamen (basal ganglia)` = "Putamen",
        `Brain - Substantia nigra` = "Substantia nigra",
        `Adipose - Subcutaneous` = "Subcutaneous adipose",
        `Adipose - Visceral (Omentum)` = "Omentum",
        `Artery - Aorta` = "Aorta",
        `Artery - Coronary` = "Coronary artery",
        `Artery - Tibial` = "Tibial artery",
        `Bladder` = "Bladder",
        `Breast - Mammary Tissue` = "Breast",
        `Cervix - Ectocervix` = "Ectocervix",
        `Cervix - Endocervix` = "Endocervix", 
        `Colon - Sigmoid` = "Sigmoid colon",
        `Colon - Transverse` = "Transverse colon",
        `Esophagus - Gastroesophageal Junction` = "Esophagus (GEJ)",
        `Esophagus - Mucosa` = "Esophageal mucosa",
        `Esophagus - Muscularis` = "Esophagus (Muscularis)",
        `Fallopian Tube` = "Fallopian tube",
        `Heart - Atrial Appendage` = "Heart (atrial appendage)",
        `Heart - Left Ventricle` = "Heart (left ventricle)",
        `Kidney - Cortex` = "Kidney",
        `Liver` = "Liver",
        `Lung` = "Lung",
        `Muscle - Skeletal` = "Skeletal muscle",
        `Pancreas` = "Pancreas",
        `Prostate` = "Prostate",
        `Skin - Nut Sun Exposed (Suprapubic)` = "Skin",
        `Small Intestine - Terminal Ileum` = "Small intestine",
        `Spleen` = "Spleen",
        `Stomach` = "Stomach",
        `Thyroid` = "Thyroid",
        `Uterus` = "Uterus",
        `Vagina` = "Vagina", .default = "drop")) %>% 
    filter(Tissue %in% keep_tissue) %>% 
    spread(Tissue, mean_RPKM) %>% 
    separate(unite, into = c("Name", "gene_name"), sep = "=")

gtex_fc <- gtex_rpkm

# For each gene, calculate the median RPKM across the tissues in each group
gtex_fc$median_brain_RPKM <- apply(select_(gtex_rpkm, brain), 1, median)
gtex_fc$median_other_RPKM <- apply(select_(gtex_rpkm, other), 1, median)

# Calculate log2 FC with a pseudocount
gtex_fc <- gtex_fc %>%
    mutate(FC     = (median_brain_RPKM + 1) / (median_other_RPKM + 1),
           log_FC = log(FC, base = 2))

joinGtex <- function(val_df) {
  
  val_df %>% 
    dplyr::select(gene_name, Consensus, state, k) %>% 
    inner_join(gtex_fc, by = "gene_name")
}

thresholdScores <- function(threshold, scores) {
    
    scores %>%
        dplyr::filter(Consensus >= threshold) %>%
        mutate(Threshold = threshold)
}

getCandidates <- function(df, method, input) {
  
  thresholds <- c(0.25, 0.5, 0.75)
  
  df <- df %>% filter(state == "ON")
  df_thresholded <- map_df(thresholds, thresholdScores, df) %>% 
    mutate(Threshold = paste0("Consensus >= ", Threshold),
           Method = method,
           Input = input)
}

val_out <- list(val_s_k, val_s_d, val_s_c,
                val_b_k, val_b_d, val_b_c)

eval_methods <- c(rep("Summary", 3), rep("Binary", 3))
eval_input <- rep(c("H3K4me3", "DNase I", "ChromHMM"), 2)

val_gtex <- val_out %>%
  mclapply(joinGtex, mc.cores = 2)

val_candidates <- val_gtex %>% 
  {mcMap(getCandidates, df = ., method = eval_methods, input = eval_input,
         mc.cores = 2)} %>% 
  bind_rows()

val_candidates_long <- val_candidates %>% 
  gather(Tissue, mean_RPKM, 6:26) %>% 
  filter(Threshold == "Consensus >= 0.75")

# Count candidates
counts <- ddply(.data = val_candidates,
                .(Method, Input, Threshold),
                summarize,
                # Define a string with the number of candidates which we can plot
                n = paste0(unique(Input), "\n(n=", length(gene_name), ")")) %>% 
    filter(Threshold == "Consensus >= 0.75")

val_candidates_long_clean <- val_candidates_long %>% 
    mutate(Label = ifelse(Tissue %in% brain, "Brain", "Other")) %>% 
    mutate(Tissue = factor(Tissue, levels = keep_tissue)) %>% 
    makeLabelSimple() %T>%
    write_tsv("tables/validation_boxplot.tsv")

val_candidates_long_clean %>%
    filter(Tissue %in% keep_tissue) %>%
    filter(Input == "H3K4me3") %>%
    mutate(Method = ifelse(Method == "Summary", "Summary (n=220)", 
                           "Binary (n=125)")) %>% 
    ggplot(aes(x = Tissue, y = mean_RPKM)) +
    geom_boxplot(width = 0.65, outlier.shape = NA, aes(fill = Method),
                 position = position_dodge(width = 0.8)) +
    scale_fill_manual(values = c("Summary (n=220)" = "red",
                                 "Binary (n=125)" = "#2f60af")) +
    facet_wrap(~ Input, ncol = 1) +
    ggmin::theme_min() +
    coord_cartesian(ylim = c(0, 55)) +
    ylab("Mean RPKM") +
    theme(text = element_text(size = 10),
          axis.text.x = element_text(angle = 90, hjust = 1),
      axis.ticks.x = element_blank(),
          legend.position = c(0.9, 0.7))

val_candidates_long_clean %>% 
    ggplot(aes(x = Tissue, y = mean_RPKM)) +
    geom_boxplot(outlier.shape = NA, aes(fill = Label), width = 0.65) +
    scale_fill_manual(values = pal_strategy) +
    ggmin::theme_min() +
    facet_grid(Input ~ Method) +
    coord_cartesian(ylim = c(0, 55)) +
    ylab("Mean RPKM") +
    theme(text = element_text(size = 10),
          axis.text.x = element_text(angle = 90, hjust = 1, size = 7),
          strip.text.y = element_blank(),
          legend.position = "none") +
    geom_text(data = counts, aes(x = 19, y = 52, label = n), size = 3,
                  colour = "black", inherit.aes = FALSE, parse = FALSE)

thresholds <- seq(0, 1, by = 0.1)

getFC <- function(threshold, scores) {

    scores %>%
        dplyr::filter(Consensus >= threshold) %>%
        `[[`("log_FC") %>%
        mean()
}

val_fc <- data.frame(thresholds)

val_fc$Summary_H3K4me3  <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[1]])

val_fc$Summary_DNase    <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[2]])

val_fc$Summary_ChromHMM <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[3]])

val_fc$Binary_H3K4me3   <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[4]])

val_fc$Binary_DNase     <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[5]])

val_fc$Binary_ChromHMM  <- data.frame(thresholds) %>% 
   apply(FUN = getFC, MARGIN = 1, val_gtex[[6]])

save(val_fc, file = "computations/validation_fc.RData")

val_fc_tidy <- val_fc %>% 
  gather(strategy, median_logFC, 2:length(.)) %>% 
  separate(strategy, into = c("Method", "Input"), sep = "_") %>% 
  mutate(Input = recode(.x = Input, DNase = "DNase I")) %>% 
  makeLabelSimple() %T>%
  write_tsv("tables/validation_foldchange.tsv")

val_fc_tidy %>% 
  filter(Input != "DNase I") %>% 
  ggplot(aes(x = thresholds, y = median_logFC)) +
  stat_smooth(aes(colour = Label, linetype = Label), se = FALSE) +
  scale_colour_manual(name = "Strategy", values = pal_strategy) +
  scale_linetype_manual(name = "Strategy", values = linetype_strategy) +
  ggmin::theme_min() +
  ylab("Mean log2 fold change of expression") + xlab("Threshold") +
  theme(legend.position = "bottom") +
  labs(colour = "") + labs(linetype = "") +
  guides(colour = guide_legend(ncol = 1))

val_fc_tidy %>% 
  ggplot(aes(x = thresholds, y = median_logFC)) +
  stat_smooth(aes(colour = Label, linetype = Label), se = FALSE) +
  scale_colour_manual(name = "Strategy", values = pal_strategy) +
  scale_linetype_manual(name = "Strategy", values = linetype_strategy) +
  ggmin::theme_min() +
  # ggtitle("Median Log FC of expression across candidate genes") +
  ylab("Mean log2 fold change of expression") + xlab("Threshold") +
  theme(legend.position = "bottom") +
  labs(colour = "") + labs(linetype = "") +
  guides(colour = guide_legend(ncol = 1))