p50 Binding Sites

genes_p50_path <- file.path(PROJHOME, "reference", 
                            "genes.protein_coding.1kb_flank.p50_overlap.bed")

genes_p50 <- 
  readr::read_tsv(genes_p50_path, col_names = c("chrom", "start", "end", "name")) %>% 
  tidyr::separate(name, c("gene_id", "gene"), sep = "_")

is_defined <- col_data$coo_ns %in% c("ABC", "GCB")
texpr_p50 <- texpr[ , is_defined]
gois <- rownames(texpr)[gene_ids_filt %in% genes_p50$gene_id]

nfkbiz_mutants <- muts %>% dplyr::filter(gene == "NFKBIZ") %$% patient %>% unique()

data <- 
  texpr_p50 %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column("gene") %>% 
  tidyr::gather(patient, expr, -gene) %>% 
  dplyr::left_join(col_data, by = "patient") %>% 
  dplyr::mutate(nfkbiz = ifelse(patient %in% nfkbiz_mutants, "mutated", "not_mutated"))

p50_abc_vs_gcb <- 
  data %>% 
  dplyr::select(-patient) %>% 
  dplyr::group_by(coo_ns, gene, nfkbiz) %>% 
  tidyr::nest() %>% 
  tidyr::spread(nfkbiz, data) %>% 
  dplyr::mutate(
    wilcox = purrr::map2(mutated, not_mutated, ~wilcox.test(.x$expr, .y$expr)), 
    pval = purrr::map_dbl(wilcox, "p.value")) %>% 
  dplyr::group_by(coo_ns) %>% 
  dplyr::mutate(qval = p.adjust(pval, "BH")) %>% 
  dplyr::select(coo_ns, gene, pval, qval) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(pval)

knitr::kable(head(p50_abc_vs_gcb, 30))

p50_abc_vs_gcb %>% 
  dplyr::mutate(is_goi = gene %in% gois) %>% 
  ggplot(aes(x = pval, fill = is_goi)) +
  geom_histogram(bins = 25) +
  scale_fill_manual(values = c(`FALSE` = "grey", `TRUE` = "red")) +
  facet_grid(~coo_ns)