Analysis for joxm

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE,
                      fig.height = 10, fig.width = 14)
dir.create("figures/donor_specific", showWarnings = FALSE, recursive = TRUE)
library(tidyverse)
library(scater)
library(ggridges)
library(GenomicRanges)
library(RColorBrewer)
library(edgeR)
library(ggrepel)
library(rlang)
library(limma)
library(org.Hs.eg.db)
library(ggforce)
library(cardelino)
library(cowplot)
library(IHW)
library(viridis)
library(ggthemes)
library(superheat)
options(stringsAsFactors = FALSE)

load("data/human_c6_v5p2.rdata")
load("data/human_H_v5p2.rdata")
load("data/human_c2_v5p2.rdata")

vep_best <- read_tsv("data/exome-point-mutations/high-vs-low-exomes.v62.ft.alldonors-filt_lenient.all_filt_sites.vep_most_severe_csq.txt")
colnames(vep_best)[1] <- "Uploaded_variation"
## deduplicate dataframe
vep_best <- vep_best[!duplicated(vep_best[["Uploaded_variation"]]),]

exome_sites <- read_tsv("data/exome-point-mutations/high-vs-low-exomes.v62.ft.filt_lenient-alldonors.txt.gz",
    col_types = "ciccdcciiiiccccccccddcdcll", comment = "#",
    col_names = TRUE)
exome_sites <- dplyr::mutate(
    exome_sites, 
    chrom = paste0("chr", gsub("chr", "", chrom)),
    var_id = paste0(chrom, ":", pos, "_", ref, "_", alt))
## deduplicate sites list
exome_sites <- exome_sites[!duplicated(exome_sites[["var_id"]]),]

vep_best[["var_id"]] <- paste0("chr", vep_best[["Uploaded_variation"]])
exome_sites <- inner_join(exome_sites, 
                           vep_best[, c("var_id", "Location", "Consequence")], 
                           by = "var_id")

cell_assign_joxm <- readRDS(file.path("data/cell_assignment", 
        paste0("cardelino_results.joxm.filt_lenient.cell_coverage_sites.rds")))

params <- list()
params$callset <- "filt_lenient.cell_coverage_sites"
fls <- list.files("data/sces")
fls <- fls[grepl(params$callset, fls)]
donors <- gsub(".*ce_([a-z]+)_.*", "\\1", fls)

sce_unst_list <- list()
for (don in donors) {
    sce_unst_list[[don]] <- readRDS(file.path("data/sces", 
        paste0("sce_", don, "_with_clone_assignments.", params$callset, ".rds")))
    cat(paste("reading", don, ":   ", ncol(sce_unst_list[[don]]), "cells.\n"))
}

assignments_lst <- list()
for (don in donors) {
    assignments_lst[[don]] <- as_data_frame(
        colData(sce_unst_list[[don]])[, 
                                      c("donor_short_id", "highest_prob", 
                                        "assigned", "total_features",
                                        "total_counts_endogenous", "num_processed")])
}
assignments <- do.call("rbind", assignments_lst)

sce_joxm <- readRDS("data/sces/sce_joxm_with_clone_assignments.filt_lenient.cell_coverage_sites.rds")
sce_joxm

table(sce_joxm$assigned)

de_res <- readRDS("data/de_analysis_FTv62/filt_lenient.cell_coverage_sites.de_results_unstimulated_cells.rds")

plot_tree <- function(tree, orient="h") {
  node_total <- max(tree$edge)
  node_shown <- length(tree$P[, 1])
  node_hidden <- node_total - node_shown
  
  prevalence <- c(tree$P[, 1]*100, rep(0, node_hidden))
  # node_size <- c(rep(20, node_shown), rep(0, node_hidden))
  
  mut_ids <- 0
  mut_id_all <- tree$Z %*% (2**seq(ncol(tree$Z),1))
  mut_id_all <- seq(length(unique(mut_id_all)),1)[as.factor(mut_id_all)]
  
  branch_ids <- NULL
  for (i in seq_len(node_total)) {
    if (i <= node_shown) {
      tree$tip.label[i] = paste0("C", i, ": ", round(prevalence[i], digits = 0),
                                 "%")
    }
    mut_num = sum(tree$sna[,3] == i)
    if (mut_num == 0) {
      if (i == node_shown + 1) {branch_ids = c(branch_ids, "Root")}
      else{branch_ids = c(branch_ids, "")} #NA
    }
    else {
      vaf <- mean(tree$VAF[tree$sna[,3] == i])
      mut_ids <- mut_ids + 1
      mut_ids <- mean(mut_id_all[tree$sna[,3] == i])
      branch_ids <- c(branch_ids, paste0(mut_num, " mutations"))
    }
  }
  pt <- ggtree::ggtree(tree)
  pt <- pt + ggplot2::geom_label(ggplot2::aes_string(x = "branch"), 
                                 label = branch_ids, color = "firebrick", size = 6)
  pt <- pt + ggplot2::xlim(-0, node_hidden + 0.5) + ggplot2::ylim(0.8, node_shown + 0.5) #the degree may not be 3
  if (orient == "v") {
    pt <- pt + ggtree::geom_tiplab(hjust = 0.39, vjust = 1.0, size = 6) + 
        ggplot2::scale_x_reverse() + ggplot2::coord_flip() 
  } else {
    pt <- pt + ggtree::geom_tiplab(hjust = 0.0, vjust = 0.5, size = 6)
  }
  pt
}

fig_tree <- plot_tree(cell_assign_joxm$full_tree, orient = "v") + 
    xlab("Clonal tree") +
    cardelino:::heatmap.theme(size = 16) +
    theme(axis.text.x = element_blank(), axis.title.y = element_text(size = 20))

prob_to_plot <- cell_assign_joxm$prob_mat[
    colnames(sce_joxm)[sce_joxm$well_condition == "unstimulated"], ]
hc <- hclust(dist(prob_to_plot))

clone_ids <- colnames(prob_to_plot)
clone_frac <- colMeans(prob_to_plot[matrixStats::rowMaxs(prob_to_plot) > 0.5,])
clone_perc <- paste0(clone_ids, ": ", 
                          round(clone_frac*100, digits = 1), "%")

colnames(prob_to_plot) <- clone_perc
nba.m <- as_data_frame(prob_to_plot[hc$order,]) %>%
    dplyr::mutate(cell = rownames(prob_to_plot[hc$order,])) %>%
    gather(key = "clone", value = "prob", -cell)
nba.m <- dplyr::mutate(nba.m, cell = factor(
    cell, levels = rownames(prob_to_plot[hc$order,])))
fig_assign <- ggplot(nba.m, aes(clone, cell, fill = prob)) + 
    geom_tile(show.legend = TRUE) +
    # scale_fill_gradient(low = "white", high = "firebrick4",
    #                     name = "posterior probability of assignment") +
    scico::scale_fill_scico(palette = "oleron", direction = 1) +
    ylab(paste("Single cells")) + 
    cardelino:::heatmap.theme(size = 16) + #cardelino:::pub.theme() +
    theme(axis.title.y = element_text(size = 20), legend.position = "bottom",
          legend.text = element_text(size = 12), legend.key.size = unit(0.05, "npc"))

plot_grid(fig_tree, fig_assign, nrow = 2, rel_heights = c(0.46, 0.52))

ggsave("figures/donor_specific/joxm_tree_probmat.png", height = 10, width = 7.5)
ggsave("figures/donor_specific/joxm_tree_probmat.pdf", height = 10, width = 7.5)
ggsave("figures/donor_specific/joxm_tree_probmat.svg", height = 10, width = 7.5)

ggsave("figures/donor_specific/joxm_tree_probmat_wide.png", height = 9, width = 10)
ggsave("figures/donor_specific/joxm_tree_probmat_wide.pdf", height = 9, width = 10)
ggsave("figures/donor_specific/joxm_tree_probmat_wide.svg", height = 9, width = 10)

joxm_config <- as_data_frame(cell_assign_joxm$full_tree$Z)
joxm_config[["var_id"]] <- rownames(cell_assign_joxm$full_tree$Z)
exome_sites_joxm <- inner_join(exome_sites, joxm_config)
exome_sites_joxm[["clone_presence"]] <- ""
for (cln in colnames(cell_assign_joxm$full_tree$Z)[-1]) {
    exome_sites_joxm[["clone_presence"]][
        as.logical(exome_sites_joxm[[cln]])] <- paste(
            exome_sites_joxm[["clone_presence"]][
                as.logical(exome_sites_joxm[[cln]])], cln, sep = "&")
}
exome_sites_joxm[["clone_presence"]] <- gsub("^&", "",
                                        exome_sites_joxm[["clone_presence"]])

exome_sites_joxm %>% group_by(Consequence, clone_presence) %>%
    summarise(n_vars = n()) %>%
ggplot(aes(x = n_vars, y = reorder(Consequence, n_vars, max), 
       colour = reorder(Consequence, n_vars, max))) +
    geom_point(size = 5) +
    geom_segment(aes(x = 0, y = Consequence, xend = n_vars, yend = Consequence)) +
    facet_wrap(~clone_presence) +
#    scale_color_brewer(palette = "Set2") +
    scale_color_manual(values = colorRampPalette(brewer.pal(8, "Accent"))(12)) +
    guides(colour = FALSE) +
    ggtitle("joxm clone tagging variants by consequence class") +
    xlab("number of variants") + ylab("consequence") +
    theme_bw(16)

## filter out any remaining ERCC genes
sce_joxm <- sce_joxm[!rowData(sce_joxm)$is_feature_control,]
sce_joxm_gr <- makeGRangesFromDataFrame(rowData(sce_joxm),
                                   start.field = "start_position",
                                   end.field = "end_position",
                                   keep.extra.columns = TRUE)
exome_sites_joxm[["chrom"]] <- gsub("chr", "", exome_sites_joxm[["chrom"]])
exome_sites_joxm_gr <- makeGRangesFromDataFrame(exome_sites_joxm,
                                   start.field = "pos",
                                   end.field = "pos",
                                   keep.extra.columns = TRUE)
# find overlaps
ov_joxm <- findOverlaps(sce_joxm_gr, exome_sites_joxm_gr)
tmp_cols <- colnames(mcols(exome_sites_joxm_gr))
tmp_cols <- tmp_cols[grepl("clone", tmp_cols)]
tmp_cols <- c("Consequence", tmp_cols, "var_id")
mut_genes_exprs_joxm <- logcounts(sce_joxm)[queryHits(ov_joxm),]
mut_genes_df_joxm <- as_data_frame(mut_genes_exprs_joxm)
mut_genes_df_joxm[["gene"]] <- rownames(mut_genes_exprs_joxm)
mut_genes_df_joxm <- bind_cols(mut_genes_df_joxm,
                          as_data_frame(
                              exome_sites_joxm_gr[
                                  subjectHits(ov_joxm)])[, tmp_cols]
)

knitr::opts_chunk$set(fig.height = 6, fig.width = 8)
summary(decideTests(de_res$qlf_list$joxm, p.value = 0.1))

topTags(de_res$qlf_list$joxm)

plotBCV(de_res$dge_list$joxm)

plotQLDisp(de_res$qlf_list$joxm)

as.data.frame(colData(de_res[["sce_list_unst"]][["joxm"]])) %>%
    dplyr::mutate(Cell_Cycle = factor(cyclone_phase, levels = c("G2M", "S", "G1")),
                  assigned = factor(assigned, levels = c("clone3", "clone2", "clone1"))) %>%
  ggplot(aes(x = assigned, fill = Cell_Cycle)) +
  geom_bar() +
  scale_fill_manual(values = c("#ff6a5c", "#ccdfcb", "#414141")) +
  coord_flip() + 
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y = element_blank())

as.data.frame(colData(de_res[["sce_list_unst"]][["joxm"]])) %>%
    dplyr::mutate(Cell_Cycle = factor(cyclone_phase, levels = c("G2M", "S", "G1")),
                  assigned = factor(assigned, levels = c("clone3", "clone2", "clone1"))) %>%
  ggplot(aes(x = assigned, fill = Cell_Cycle)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("#ff6a5c", "#ccdfcb", "#414141")) +
  coord_flip() + 
  ylab("proportion") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y = element_blank())

freqs <- as.matrix(table(
    de_res[["sce_list_unst"]][["joxm"]]$assigned,  
    de_res[["sce_list_unst"]][["joxm"]]$cyclone_phase))
fisher.test(freqs)

fisher.test(freqs[-3,])

choose_joxm_cells <- (sce_joxm$well_condition == "unstimulated" &
                          sce_joxm$assigned != "unassigned")
pca_unst <- reducedDim(runPCA(sce_joxm[, choose_joxm_cells], 
                              ntop = 500, ncomponents = 10), "PCA")
pca_unst <- data.frame(
    PC1 = pca_unst[, 1], PC2 = pca_unst[, 2], 
    PC3 = pca_unst[, 3], PC4 = pca_unst[, 4],
    PC5 = pca_unst[, 5], PC6 = pca_unst[, 6],
    clone = sce_joxm[, choose_joxm_cells]$assigned,
    nvars_cloneid = sce_joxm[, choose_joxm_cells]$nvars_cloneid,
    cyclone_phase = sce_joxm[, choose_joxm_cells]$cyclone_phase,
    G1 = sce_joxm[, choose_joxm_cells]$G1,
    G2M = sce_joxm[, choose_joxm_cells]$G2M,
    S = sce_joxm[, choose_joxm_cells]$S,
    clone1_prob = sce_joxm[, choose_joxm_cells]$clone1_prob,
    clone2_prob = sce_joxm[, choose_joxm_cells]$clone2_prob,
    clone3_prob = sce_joxm[, choose_joxm_cells]$clone3_prob,
    RPS6KA2 = as.vector(logcounts(sce_joxm[grep("RPS6KA2", rownames(sce_joxm)), choose_joxm_cells]))
    )

ggplot(pca_unst, aes(x = PC1, y = PC2, fill = clone)) +
    geom_point(pch = 21, size = 4, colour = "gray30") +
    scale_fill_brewer(palette = "Accent", name = "assigned\nclone") +
    theme_classic(14)

ggplot(pca_unst, aes(x = PC2, y = PC3, fill = clone)) +
    geom_point(pch = 21, size = 4, colour = "gray30") +
    scale_fill_brewer(palette = "Accent", name = "assigned\nclone") +
    theme_classic(14)

ggplot(pca_unst, aes(x = PC2, y = PC4, fill = clone)) +
    geom_point(pch = 21, size = 4, colour = "gray30") +
    scale_fill_brewer(palette = "Accent", name = "assigned\nclone") +
    theme_classic(14)

ggplot(pca_unst, aes(x = PC3, y = PC4, fill = clone)) +
    geom_point(pch = 21, size = 4, colour = "gray30") +
    scale_fill_brewer(palette = "Accent", name = "assigned\nclone") +
    theme_classic(14)

ppc1 <- ggplot(pca_unst, aes(x = PC1, y = PC2, fill = clone1_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone1\nposterior\nprobability") +
  theme_classic(14)

ppc2 <- ggplot(pca_unst, aes(x = PC1, y = PC2, fill = clone2_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone2\nposterior\nprobability") +
  theme_classic(14)

ppc3 <- ggplot(pca_unst, aes(x = PC1, y = PC2, fill = clone3_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone3\nposterior\nprobability") +
  theme_classic(14)

plot_grid(ppc1, ppc2, ppc3, labels = "auto", ncol = 1)

ggsave("figures/donor_specific/joxm_pca_pc1_pc2_clone_probs.png", height = 11, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_pc1_pc2_clone_probs.pdf", height = 11, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_pc1_pc2_clone_probs.svg", height = 11, width = 9.5)

ppc1 <- ggplot(pca_unst, aes(x = PC2, y = PC3, fill = clone1_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone1\nposterior\nprobability") +
  theme_classic(14)

ppc2 <- ggplot(pca_unst, aes(x = PC2, y = PC3, fill = clone2_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone2\nposterior\nprobability") +
  theme_classic(14)

ppc3 <- ggplot(pca_unst, aes(x = PC2, y = PC3, fill = clone3_prob)) +
  geom_point(pch = 21, size = 4, colour = "gray30") +
  scale_fill_viridis(option = "A", name = "clone3\nposterior\nprobability") +
  theme_classic(14)

plot_grid(ppc1, ppc2, ppc3, labels = "auto", ncol = 1)

ggsave("figures/donor_specific/joxm_pca_pc2_pc3_clone_probs.png", height = 11, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_pc2_pc3_clone_probs.pdf", height = 11, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_pc2_pc3_clone_probs.svg", height = 11, width = 9.5)

pca_unst$cyclone_phase <- factor(pca_unst$cyclone_phase, levels = c("G1", "S", "G2M"))
ggplot(pca_unst, aes(x = PC1, y = PC2, colour = cyclone_phase,
                     shape = clone)) +
    geom_point(size = 6) +
    scale_color_manual(values = magma(6)[c(1, 3, 5)], name = "cell cycle\nphase") +
    xlab("PC1 (10% variance)") +
    ylab("PC2 (5% variance)") +
    theme_classic(18)

ggsave("figures/donor_specific/joxm_pca.png", height = 6, width = 9.5)
ggsave("figures/donor_specific/joxm_pca.pdf", height = 6, width = 9.5)
ggsave("figures/donor_specific/joxm_pca.svg", height = 6, width = 9.5)


pca_unst$cyclone_phase <- factor(pca_unst$cyclone_phase, levels = c("G1", "S", "G2M"))
ggplot(pca_unst, aes(x = PC2, y = PC3, colour = cyclone_phase,
                     shape = clone)) +
    geom_point(size = 6) +
    scale_color_manual(values = magma(6)[c(1, 3, 5)], name = "cell cycle\nphase") +
    xlab("PC2 (5% variance)") +
    ylab("PC3 (3% variance)") +
    theme_classic(18)

ggplot(pca_unst, aes(x = PC1, y = PC2, fill = G2M,
                     shape = clone)) +
    geom_point(colour = "gray50", size = 5) +
    scale_shape_manual(values = c(21, 23, 25), name = "clone") +
    scico::scale_fill_scico(palette = "bilbao", name  = "G2/M score") +
    scale_size_continuous(range = c(4, 6)) +
    xlab("PC1 (10% variance)") +
    ylab("PC2 (5% variance)") +
    theme_classic(18)

ggsave("figures/donor_specific/joxm_pca_g2m_score.png", height = 6, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_g2m_score.pdf", height = 6, width = 9.5)
ggsave("figures/donor_specific/joxm_pca_g2m_score.svg", height = 6, width = 9.5)

ggplot(pca_unst, aes(x = PC1, y = PC2, colour = S,
                     shape = clone)) +
    geom_point(size = 5) +
    scale_color_viridis(option = "B") +
    xlab("PC1 (10% variance)") +
    ylab("PC2 (5% variance") +
    theme_classic(18)

ggplot(pca_unst, aes(x = PC1, y = PC2, colour = G1,
                     shape = clone)) +
    geom_point(size = 5) +
    scale_color_viridis(option = "B") +
    xlab("PC1 (10% variance)") +
    ylab("PC2 (5% variance") +
    theme_classic(18)

ggplot(pca_unst, aes(x = PC1, y = PC2, fill = clone2_prob, size = nvars_cloneid)) +
    geom_point(pch = 21, colour = "gray30") +
    scale_fill_viridis(option = "B", name = "clone2\nprobability") +
    scale_size_continuous(name = "# variants\nfor clone ID") +
    theme_classic(14)

de_res_cc <- readRDS("data/de_analysis_FTv62/cellcycle_analyses/filt_lenient.cell_coverage_sites.de_results_unstimulated_cells.cc.rds")
de_joxm_cellcycle_only <- de_res_cc$cellcycle_only$qlf_list$joxm
topTags(de_res_cc$cellcycle_only$qlf_list$joxm)

summary(decideTests(de_res_cc$cellcycle_only$qlf_list$joxm, p.value = 0.1))

summary(decideTests(de_res_cc$cellcycle_clone$qlf_list$joxm, p.value = 0.1))

topTags(de_res_cc$cellcycle_clone$qlf_list$joxm)

## accounting for cell cycle in model
head(de_res_cc$camera$H$joxm$clone2_clone1$logFC)

df <- data_frame(
  pval_clone = de_res$qlf_list$joxm$table$PValue,
  fdr_clone = p.adjust(de_res$qlf_list$joxm$table$PValue, method = "BH"),
  pval_cellcycle_only = de_res_cc$cellcycle_only$qlf_list$joxm$table$PValue,
  pval_cellcycle_clone = de_res_cc$cellcycle_clone$qlf_list$joxm$table$PValue)

ggplot(df, aes(-log10(pval_clone), -log10(pval_cellcycle_clone),
               colour = fdr_clone < 0.05)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) +
  xlab("P-value for clone DE, not accounting for cell cycle (-log10 scale)") +
  ylab("P-value for clone DE, accounting for cell cycle (-log10 scale)") +
  theme_classic()

## tree and cell assignment
fig_tree <- plot_tree(cell_assign_joxm$full_tree, orient = "v") + 
    xlab("Clonal tree") +
    cardelino:::heatmap.theme(size = 16) +
    theme(axis.text.x = element_blank(), axis.title.y = element_text(size = 20))

prob_to_plot <- cell_assign_joxm$prob_mat[
    colnames(sce_joxm)[sce_joxm$well_condition == "unstimulated"], ]
hc <- hclust(dist(prob_to_plot))

clone_ids <- colnames(prob_to_plot)
clone_frac <- colMeans(prob_to_plot[matrixStats::rowMaxs(prob_to_plot) > 0.5,])
clone_perc <- paste0(clone_ids, ": ", 
                          round(clone_frac*100, digits = 1), "%")

colnames(prob_to_plot) <- clone_perc
nba.m <- as_data_frame(prob_to_plot[hc$order,]) %>%
    dplyr::mutate(cell = rownames(prob_to_plot[hc$order,])) %>%
    gather(key = "clone", value = "prob", -cell)
nba.m <- dplyr::mutate(nba.m, cell = factor(
    cell, levels = rownames(prob_to_plot[hc$order,])))
fig_assign <- ggplot(nba.m, aes(clone, cell, fill = prob)) + 
    geom_tile(show.legend = TRUE) +
    # scale_fill_gradient(low = "white", high = "firebrick4",
    #                     name = "posterior probability of assignment") +
    scico::scale_fill_scico(palette = "oleron", direction = 1) +
    ylab(paste("Single cells")) + 
    cardelino:::heatmap.theme(size = 16) + #cardelino:::pub.theme() +
    theme(axis.title.y = element_text(size = 20), legend.position = "bottom",
          legend.text = element_text(size = 12), legend.key.size = unit(0.05, "npc"))

p_tree <- plot_grid(fig_tree, fig_assign, nrow = 2, rel_heights = c(0.46, 0.52))


## cell cycle barplot
p_bar <- as.data.frame(colData(de_res[["sce_list_unst"]][["joxm"]])) %>%
  dplyr::mutate(Cell_Cycle = factor(cyclone_phase, levels = c("G2M", "S", "G1")),
                assigned = factor(assigned, levels = c("clone3", "clone2", "clone1"))) %>%
  ggplot(aes(x = assigned, fill = Cell_Cycle, group = Cell_Cycle)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("#ff6a5c", "#ccdfcb", "#414141")) +
  coord_flip() + 
  ylab("proportion of cells") +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(axis.title.y = element_blank())

## effects on mutated clone
df_to_plot <- mut_genes_df_allcells %>%
  dplyr::filter(!is.na(logFC), donor == "joxm") %>%
  dplyr::filter(!duplicated(gene)) %>%
  dplyr::mutate(
    FDR = p.adjust(PValue, method = "BH"),
    consequence_simplified = factor(
      consequence_simplified, 
      levels(as.factor(consequence_simplified))[order(tmp4[["med"]])]),
    de  = ifelse(FDR < 0.2, "FDR < 0.2", "FDR > 0.2"))


p_mutated_clone <- ggplot(df_to_plot, aes(y = logFC, x = consequence_simplified)) +
    geom_hline(yintercept = 0, linetype = 1, colour = "black") +
    geom_boxplot(outlier.size = 0, outlier.alpha = 0, fill = "gray90",
                 colour = "firebrick4", width = 0.2, size = 1) +
    ggbeeswarm::geom_quasirandom(aes(fill = -log10(PValue)), 
                                 colour = "gray40", pch = 21, size = 4) +
    geom_segment(aes(y = -0.25, x = 0, yend = -1, xend = 0), 
                 colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
    annotate("text", y = -3, x = 0, size = 6, label = "lower in mutated clone") +
    geom_segment(aes(y = 0.25, x = 0, yend = 1, xend = 0), 
                 colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
    annotate("text", y = 3, x = 0, size = 6, label = "higher in mutated clone") +
    scale_x_discrete(expand = c(0.1, .05), name = "consequence") +
    scale_y_continuous(expand = c(0.1, 0.1), name = "logFC") +
    expand_limits(x = c(-0.75, 8)) +
    theme_ridges(22) +
    coord_flip() +
    scale_fill_viridis(option = "B", name = "-log10(P)") +
    theme(strip.background = element_rect(fill = "gray90"),
          legend.position = "right") +
    guides(color = FALSE)

## PCA
p_pca <- ggplot(pca_unst, aes(x = PC2, y = PC3, fill = clone,
                     shape = clone)) +
    geom_point(colour = "gray50", size = 5) +
    scale_shape_manual(values = c(21, 23, 25, 22, 24, 26), name = "clone") +
    # scico::scale_fill_scico(palette = "bilbao", name  = "G2/M score") +
    ggthemes::scale_fill_canva(palette = "Surf and turf") +
    scale_size_continuous(range = c(4, 6)) +
    xlab("PC2 (5% variance)") +
    ylab("PC3 (3% variance)") +
    theme_classic(18)

 # ggplot(pca_unst, aes(x = PC2, y = PC3, colour = clone,
 #                     shape = cyclone_phase)) +
 #    geom_point(alpha = 0.9, size = 5) +
 #    scale_shape_manual(values = c(15, 17, 19), name = "clone") +
 #    # scico::scale_fill_scico(palette = "bilbao", name  = "G2/M score") +
 #    ggthemes::scale_color_canva(palette = "Surf and turf") +
 #    scale_size_continuous(range = c(4, 6)) +
 #    xlab("PC2 (5% variance)") +
 #    ylab("PC3 (3% variance)") +
 #    theme_classic(18)

## volcano
de_joxm_cl2_vs_cl1 <- de_res$qlf_pairwise$joxm$clone2_clone1$table
de_joxm_cl2_vs_cl1[["gene"]] <- rownames(de_joxm_cl2_vs_cl1)
de_joxm_cl2_vs_cl1 <- de_joxm_cl2_vs_cl1 %>% 
  dplyr::mutate(FDR = adj_pvalues(ihw(PValue ~ logCPM, alpha = 0.1)), 
                sig = FDR < 0.1,
                signed_F = sign(logFC) * F) 
de_joxm_cl2_vs_cl1[["lab"]] <- ""
int_genes_entrezid <- c(Hs.H$HALLMARK_G2M_CHECKPOINT, Hs.H$HALLMARK_E2F_TARGETS,
                        Hs.H$HALLMARK_MITOTIC_SPINDLE)
mm <- match(int_genes_entrezid, de_joxm_cl2_vs_cl1$entrezid)
mm <- mm[!is.na(mm)]
int_genes_hgnc <- de_joxm_cl2_vs_cl1$hgnc_symbol[mm]
genes_to_label <- (de_joxm_cl2_vs_cl1[["hgnc_symbol"]] %in% int_genes_hgnc &
                     de_joxm_cl2_vs_cl1[["FDR"]] < 0.01)
de_joxm_cl2_vs_cl1[["lab"]][genes_to_label] <-
  de_joxm_cl2_vs_cl1[["hgnc_symbol"]][genes_to_label]
de_joxm_cl2_vs_cl1[["cell_cycle_gene"]] <- (de_joxm_cl2_vs_cl1$entrezid %in% 
                                              int_genes_entrezid)

p_volcano <- ggplot(de_joxm_cl2_vs_cl1, aes(x = logFC, y = -log10(PValue), 
                                            fill = sig, label = lab)) +
  geom_point(aes(size = sig), pch = 21, colour = "gray40") +
  geom_label_repel(show.legend = FALSE, 
                   arrow = arrow(type = "closed", length = unit(0.25, "cm")), 
                   nudge_x = 0.2, nudge_y = 0.3, fill = "gray95") +
  geom_segment(aes(x = -1, y = 0, xend = -4, yend = 0), 
               colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
  annotate("text", x = -4, y = -0.5, label = "higher in clone1", size = 6) +
  geom_segment(aes(x = 1, y = 0, xend = 4, yend = 0), 
               colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
  annotate("text", x = 4, y = -0.5, label = "higher in clone2", size = 6) +
  scale_fill_manual(values = c("gray60", "firebrick"), 
                    label = c("N.S.", "FDR < 10%"), name = "") +
  scale_size_manual(values = c(1, 3), guide = FALSE) +
  ylab(expression(-"log"[10](P))) +
  xlab(expression("log"[2](FC))) +
  guides(alpha = FALSE,
         fill = guide_legend(override.aes = list(size = 5))) +
  theme_classic(20) + theme(legend.position = "right")

# ggplot(de_joxm_cl2_vs_cl1, aes(x = logFC, y = -log10(PValue), 
#                                fill = cell_cycle_gene, label = lab)) +
#   geom_point(aes(size = sig), pch = 21, colour = "gray40") +
#   geom_point(aes(size = sig), pch = 21, colour = "gray40", 
#              data = dplyr::filter(de_joxm_cl2_vs_cl1, cell_cycle_gene)) +
#   geom_label_repel(show.legend = FALSE, 
#                    arrow = arrow(type = "closed", length = unit(0.25, "cm")), 
#                    nudge_x = 0.2, nudge_y = 0.3, fill = "gray95") +
#   geom_segment(aes(x = -1, y = 0, xend = -4, yend = 0), 
#                colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
#   annotate("text", x = -4, y = -0.5, label = "higher in clone1", size = 6) +
#   geom_segment(aes(x = 1, y = 0, xend = 4, yend = 0), 
#                colour = "black", size = 1, arrow = arrow(length = unit(0.5, "cm"))) +
#   annotate("text", x = 4, y = -0.5, label = "higher in clone2", size = 6) +
#   scale_fill_manual(values = c("gray60", "firebrick"), 
#                     label = c("N.S.", "FDR < 10%"), name = "") +
#   scale_size_manual(values = c(1, 3), guide = FALSE) +
#   guides(alpha = FALSE) +
#   theme_classic(20) + theme(legend.position = "right")

## genesets
cam_H_pw <- de_res$camera$H$joxm$clone2_clone1$logFC
cam_H_pw[["geneset"]] <- rownames(cam_H_pw)
cam_H_pw <- cam_H_pw %>% 
  dplyr::mutate(sig = FDR < 0.05) 
cam_H_pw[["lab"]] <- ""
cam_H_pw[["lab"]][1:3] <-
    cam_H_pw[["geneset"]][1:3]

cam_H_pw <- dplyr::mutate(
    cam_H_pw,
    Direction = gsub("clone4", "clone2", Direction)
)
p_genesets <- cam_H_pw %>% 
  ggplot(aes(x = Direction, y = -log10(PValue), colour = sig, 
             label = lab)) +
  ggbeeswarm::geom_quasirandom(aes(size = NGenes)) +
  geom_label_repel(show.legend = FALSE,
                   nudge_y = 0.3, nudge_x = 0.3, fill = "gray95") +
  scale_colour_manual(values = c("gray50", "firebrick"), 
                      label = c("N.S.", "FDR < 5%"), name = "") +
  guides(alpha = FALSE,
         colour = guide_legend(override.aes = list(size = 5))) +
  xlab("Gene set enrichment direction") +
  ylab(expression(-"log"[10](P))) +
  theme_classic(20) + theme(legend.position = "right")


## produce combined fig
## combine pca and barplot
p_bar_pca <- ggdraw() +
  draw_plot(p_pca + theme(legend.justification = "bottom"), 0, 0, 1, 1) +
  draw_plot(p_bar, x = 0.48, 0.65, height = 0.35, width = 0.52, scale = 1)

ggdraw() +
    draw_plot(p_tree, x = 0,  y = 0.45, width = 0.48, height = 0.55, scale = 1) +
    draw_plot(p_bar_pca, x = 0.52, y = 0.45, width = 0.48, height = 0.55, scale = 1) +
    draw_plot(p_volcano, x = 0,  y = 0, width = 0.48, height = 0.45, scale = 1) +
    draw_plot(p_genesets, x = 0.52,  y = 0, width = 0.48, height = 0.45, scale = 1) +
    draw_plot_label(letters[1:4], x = c(0, 0.5, 0, 0.5), 
                    y = c(1, 1, 0.45, 0.45), size = 36)

ggsave("figures/donor_specific/joxm_combined_fig.png", 
       height = 16, width = 18)
ggsave("figures/donor_specific/joxm_combined_fig.pdf", 
       height = 16, width = 18)


## plots for talk
ggsave("figures/donor_specific/joxm_bar_pca.png", plot = p_bar_pca,
       height = 7, width = 10)
ggsave("figures/donor_specific/joxm_volcano.png", plot = p_volcano,
       height = 6, width = 10)
ggsave("figures/donor_specific/joxm_genesets.png", plot = p_genesets,
       height = 6, width = 10)
ggsave("figures/donor_specific/joxm_mutated_clone.png", plot = p_mutated_clone,
       height = 6, width = 14)


# ggdraw() +
#     draw_plot(p_tree, x = 0,  y = 0.57, width = 0.48, height = 0.43, scale = 1) +
#     draw_plot(p_bar_pca, x = 0.52, y = 0.57, width = 0.48, height = 0.43, scale = 1) +
#     draw_plot(p_volcano, x = 0,  y = 0.3, width = 0.48, height = 0.27, scale = 1) +
#     draw_plot(p_genesets, x = 0.52,  y = 0.3, width = 0.48, height = 0.27, scale = 1) +
#     #draw_plot(p_table, x = 0,  y = 0.2, width = 1, height = 0.15, scale = 1) +
#     draw_plot(p_mutated_clone, x = 0.05,  y = 0, width = 0.9, height = 0.3, scale = 1) +
#     draw_plot_label(letters[1:5], x = c(0, 0.5, 0, 0.5, 0), 
#                     y = c(1, 1, 0.57, 0.57, 0.3), size = 36)
# ggsave("figures/donor_specific/joxm_combined_fig.png", 
#        height = 20, width = 19)
# ggsave("figures/donor_specific/joxm_combined_fig.pdf", 
#        height = 20, width = 19)

devtools::session_info()