Mutated gene analysis across all lines

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"]]),]

table(vep_best[["Consequence"]])

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::select(exome_sites,
                             -c(clinsing, SIFT_Pdel, PolyPhen_Pdel, gene_name,
                                bcftools_csq))
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"]])
vep_best[, c("var_id", "Location", "Consequence")]

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

# com_vars <- intersect(exome_sites$var_id, csq_annos$var_id)
# length(com_vars)
# mean(exome_sites$var_id %in% csq_annos$var_id)
# mean(exome_sites$pos %in% csq_annos$Pos)

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

sce_unst_list <- list()
for (don in lines) {
    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 lines) {
    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)

de_res <- readRDS(paste0("data/de_analysis_FTv62/",
                         params$callset, 
                         ".de_results_unstimulated_cells.rds"))

cell_assign_list <- list()
for (don in lines) {
    cell_assign_list[[don]] <- readRDS(file.path("data/cell_assignment", 
        paste0("cardelino_results.", don, ".", params$callset, ".rds")))
    cat(paste("reading", don, "\n"))
} 

get_sites_by_line <- function(sites_df, sce_list, assign_list) {
    if (!identical(sort(names(sce_list)), sort(names(assign_list))))
        stop("lines do not match between sce_list and assign_list.")
    sites_by_line <- list()
    for (don in names(sce_list)) {
        config <- as_data_frame(assign_list[[don]]$tree$Z)
        config[["var_id"]] <- rownames(assign_list[[don]]$tree$Z)
        sites_line <- inner_join(sites_df, config)
        sites_line[["clone_presence"]] <- ""
        for (cln in colnames(assign_list[[don]]$tree$Z)[-1]) {
            sites_line[["clone_presence"]][
                as.logical(sites_line[[cln]])] <- paste(
                    sites_line[["clone_presence"]][
                        as.logical(sites_line[[cln]])], cln, sep = "&")
        }
        sites_line[["clone_presence"]] <- gsub("^&", "",
                                                sites_line[["clone_presence"]])
        ## drop config columns as these won't match up between lines
        keep_cols <- grep("^clone[0-9]$", colnames(sites_line), invert = TRUE)
        sites_by_line[[don]] <- sites_line[, keep_cols]
    }
    do.call("bind_rows", sites_by_line)
}

sites_by_line <- get_sites_by_line(exome_sites, sce_unst_list, cell_assign_list)
sites_by_line <- dplyr::mutate(sites_by_line, 
                                Consequence = gsub("_variant", "", Consequence))

## run DE for mutated cells vs unmutated cells using existing DE results
## filter out any remaining ERCC genes
for (don in names(de_res[["sce_list_unst"]]))
    de_res[["sce_list_unst"]][[don]] <- de_res[["sce_list_unst"]][[don]][
        !rowData(de_res[["sce_list_unst"]][[don]])$is_feature_control,]  
sce_de_list_gr <- list()
for (don in names(de_res[["sce_list_unst"]])) {
    sce_de_list_gr[[don]] <- makeGRangesFromDataFrame(
        rowData(de_res[["sce_list_unst"]][[don]]),
        start.field = "start_position",
        end.field = "end_position",
        keep.extra.columns = TRUE)
    seqlevelsStyle(sce_de_list_gr[[don]]) <- "Ensembl"
}
sites_by_line <- dplyr::mutate(sites_by_line, end_pos = pos)
sites_by_line_gr <- makeGRangesFromDataFrame(
        sites_by_line,
        start.field = "pos",
        end.field = "end_pos",
        keep.extra.columns = TRUE)
seqlevelsStyle(sites_by_line_gr) <- "Ensembl"

mut_genes_df_allcells_list <- list()
for (don in names(de_res[["sce_list_unst"]])) {
    cat("working on ", don, "\n")
    sites_tmp <- sites_by_line_gr[sites_by_line_gr$donor_short_id == don]
    ov_tmp <- findOverlaps(sce_de_list_gr[[don]], sites_tmp)
    sce_tmp <- de_res[["sce_list_unst"]][[don]][queryHits(ov_tmp),]
    sites_tmp <- sites_tmp[subjectHits(ov_tmp)]
    sites_tmp$gene <- rownames(sce_tmp)
    dge_tmp <- de_res[["dge_list"]][[don]]
    dge_tmp <- dge_tmp[intersect(rownames(dge_tmp), sites_tmp$gene),]
    base_design <- dge_tmp$design[, !grepl("assigned", colnames(dge_tmp$design))]
    de_tbl_tmp <- data.frame(line = don,
                             gene = sites_tmp$gene, 
                             hgnc_symbol = gsub(".*_", "", sites_tmp$gene),
                             ensembl_gene_id = gsub("_.*", "", sites_tmp$gene),
                             var_id = sites_tmp$var_id,
                             location = sites_tmp$Location,
                             consequence = sites_tmp$Consequence,
                             clone_presence = sites_tmp$clone_presence,
                             logFC = NA, logCPM = NA, F = NA, PValue = NA,
                             comment = "", stringsAsFactors = FALSE)
    for (i in seq_len(length(sites_tmp))) {
        clones_tmp <- strsplit(sites_tmp$clone_presence[i], split = "&")[[1]]
        mutatedclone <- as.numeric(sce_tmp$assigned %in% clones_tmp)
        dsgn_tmp <- cbind(base_design, data.frame(mutatedclone))
        if (sites_tmp$gene[i] %in% rownames(dge_tmp) && is.fullrank(dsgn_tmp)) {
            qlfit_tmp <- glmQLFit(dge_tmp[sites_tmp$gene[i],], dsgn_tmp)
            de_tmp <- glmQLFTest(qlfit_tmp, coef = ncol(dsgn_tmp))
            de_tbl_tmp$logFC[i] <- de_tmp$table$logFC
            de_tbl_tmp$logCPM[i] <- de_tmp$table$logCPM
            de_tbl_tmp$F[i] <- de_tmp$table$F
            de_tbl_tmp$PValue[i] <- de_tmp$table$PValue
        }
        if (!(sites_tmp$gene[i] %in% rownames(dge_tmp)))
            de_tbl_tmp$comment[i] <- "gene did not pass DE filters"
        if (!is.fullrank(dsgn_tmp))
            de_tbl_tmp$comment[i] <- "insufficient cells assigned to clone"
    }
    mut_genes_df_allcells_list[[don]] <- de_tbl_tmp
}

mut_genes_df_allcells <- do.call("bind_rows", mut_genes_df_allcells_list)
dim(mut_genes_df_allcells)

head(mut_genes_df_allcells)

## add FDRs for genes tested here for DE
ihw_res_all <- ihw(PValue ~ logCPM, data = mut_genes_df_allcells, alpha = 0.2)
mut_genes_df_allcells$FDR <- adj_pvalues(ihw_res_all)
## add simplified consequence categories
mut_genes_df_allcells$consequence_simplified <- 
    mut_genes_df_allcells$consequence
mut_genes_df_allcells$consequence_simplified[
    mut_genes_df_allcells$consequence_simplified %in% 
        c("stop_retained", "start_lost", "stop_lost", "stop_gained")] <- "nonsense"
mut_genes_df_allcells$consequence_simplified[
    mut_genes_df_allcells$consequence_simplified %in% 
        c("splice_donor", "splice_acceptor", "splice_region")] <- "splicing"
table(mut_genes_df_allcells$consequence_simplified)

summary(mut_genes_df_allcells$FDR)

dplyr::arrange(mut_genes_df_allcells, FDR) %>% dplyr::select(-location) %>% 
  DT::datatable(., options = list(pageLength = 20))

dplyr::filter(mut_genes_df_allcells, !is.na(logFC)) %>%
ggplot(aes(x = logFC, y = -log10(PValue), 
                         fill = FDR < 0.1)) +
    geom_point(colour = "gray40", pch = 21) +
    geom_rug(sides = "b", alpha = 0.6) +
    geom_vline(xintercept = 0, linetype = 2) +
    # geom_text_repel(show.legend = FALSE) +
    scale_fill_manual(values = c("gray70", "firebrick"),
                        label = c("N.S.", "FDR < 10%"), 
                        name = "DE result") +
    facet_wrap(~consequence, ncol = 5) +
    guides(alpha = FALSE) +    
    theme_classic(16) +
    theme(strip.background = element_rect(fill = "gray90"),
          legend.position = c(0.9, 0.1))

ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_vep_anno_allcells.png", 
       height = 7, width = 16)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_vep_anno_allcells.pdf", 
       height = 7, width = 16)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_vep_anno_allcells.svg", 
       height = 7, width = 16)  

dplyr::filter(mut_genes_df_allcells, !is.na(logFC)) %>%
ggplot(aes(x = logFC, y = -log10(PValue), 
                         fill = FDR < 0.2)) +
    geom_point(colour = "gray40", pch = 21) +
    geom_rug(sides = "b", alpha = 0.6) +
    geom_vline(xintercept = 0, linetype = 2) +
    # geom_text_repel(show.legend = FALSE) +
    scale_fill_manual(values = c("gray70", "firebrick"),
                        label = c("N.S.", "FDR < 20%"), 
                        name = "DE result") +
    facet_wrap(~consequence_simplified, ncol = 3) +
    guides(alpha = FALSE) +    
    theme_classic(20) +
    theme(strip.background = element_rect(fill = "gray90"),
          strip.text = element_text(size = 13),
          legend.position = c(0.8, 0.15))

ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells.png", 
       height = 7, width = 10.5)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells.pdf", 
       height = 7, width = 10.5)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells.svg", 
       height = 7, width = 10.5)  

tmp2 <- mut_genes_df_allcells %>%
    group_by(consequence_simplified) %>%
    summarise(med = median(logFC, na.rm = TRUE),
              nvars = n())
dplyr::arrange(tmp2, med)

plist <- list()
simple_cons_vec <- c("nonsense", "splicing", "missense", 
                     "synonymous", "intron", 
                     "non_coding_transcript_exon",
                     "5_prime_UTR", "3_prime_UTR")
for (cons in simple_cons_vec) {
    p <- dplyr::filter(mut_genes_df_allcells, !is.na(logFC), 
                       consequence_simplified == cons) %>%
        ggplot(aes(x = logFC, y = -log10(PValue), 
                   fill = FDR < 0.2)) +
        geom_point(colour = "gray40", pch = 21, size = 2) +
        geom_rug(sides = "b", alpha = 0.6) +
        geom_vline(xintercept = 0, linetype = 2) +
        scale_fill_manual(values = c("gray70", "firebrick"),
                          label = c("N.S.", "FDR < 20%"), 
                          name = "DE result") +
        guides(alpha = FALSE) +    
        theme_classic(14) +
        ggtitle(cons) +
        xlim(1.05 * range(mut_genes_df_allcells$logFC, na.rm = TRUE)) +
        ylim(c(0, 1.05 * max(-log10(mut_genes_df_allcells$PValue), na.rm = TRUE))) +
        theme(strip.background = element_rect(fill = "gray90"),
              strip.text = element_text(size = 12),
              legend.position = "none",
              plot.title = element_text(face = "bold"))
    plist[[cons]] <- ggExtra::ggMarginal(
        p, type = "boxplot", margins = "x", fill = "wheat", 
        colour = "firebrick4", outlier.size = 0, outlier.alpha = 0)
}
cowplot::plot_grid(plotlist = plist, nrow = 2)

ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells_with_boxplot.png", 
       height = 7, width = 18)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells_with_boxplot.pdf", 
       height = 7, width = 18)
ggsave("figures/mutated_genes/alllines_mutgenes_volcano_by_simple_vep_anno_allcells_with_boxplot.svg", 
       height = 7, width = 18)  

tmp2 <- mut_genes_df_allcells %>%
    group_by(consequence) %>%
    summarise(med = median(logFC, na.rm = TRUE),
              nvars = n())

mut_genes_df_allcells %>%
    dplyr::filter(!is.na(logFC)) %>%
    dplyr::mutate(consequence = factor(
        consequence, levels(as.factor(consequence))[order(tmp2[["med"]])]),
        de  = ifelse(FDR < 0.2, "FDR < 0.2", "FDR > 0.2")) %>%
ggplot(aes(y = logFC, x = consequence)) +
    geom_hline(yintercept = 0, linetype = 1, colour = "black") +
    ggbeeswarm::geom_quasirandom(colour = "gray40", fill = "gray70", pch = 21) +
    geom_boxplot(outlier.size = 0, outlier.alpha = 0, fill = "wheat",
                 colour = "firebrick4", width = 0.2, size = 1) +
    scale_x_discrete(expand = c(.01, 0)) +
    scale_y_continuous(expand = c(0, 0), name = "logFC") +
    theme_ridges(16) +
    coord_flip() +
    theme(strip.background = element_rect(fill = "gray90")) +
    guides(fill = FALSE, color = FALSE, point_color = FALSE)

ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_vep_anno_allcells.png", 
       height = 9, width = 10)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_vep_anno_allcells.pdf", 
       height = 9, width = 10)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_vep_anno_allcells.svg", 
       height = 9, width = 10)

tmp3 <- mut_genes_df_allcells %>%
    group_by(consequence_simplified) %>%
    summarise(med = median(logFC, na.rm = TRUE),
              nvars = n(),
              t_coef = t.test(logFC, mu = 0)$estimate,
              t_stat = t.test(logFC, mu = 0)$statistic,
              t_pval = t.test(logFC, mu = 0)$p.value)
tmp3

mut_genes_df_allcells %>%
    dplyr::filter(!is.na(logFC)) %>%
    dplyr::mutate(consequence_simplified = factor(
        consequence_simplified, 
        levels(as.factor(consequence_simplified))[order(tmp3[["med"]])]),
        de  = ifelse(FDR < 0.2, "FDR < 0.2", "FDR > 0.2")) %>%
ggplot(aes(y = logFC, x = consequence_simplified)) +
    geom_hline(yintercept = 0, linetype = 1, colour = "black") +
    ggbeeswarm::geom_quasirandom(colour = "gray40", fill = "gray70", pch = 21) +
    geom_boxplot(outlier.size = 0, outlier.alpha = 0, fill = "wheat",
                 colour = "firebrick4", width = 0.2, size = 1) +
    scale_x_discrete(expand = c(.01, 0), name = "consequence") +
    scale_y_continuous(expand = c(0, 0), name = "logFC") +
    theme_ridges(20) +
    coord_flip() +
    theme(strip.background = element_rect(fill = "gray90")) +
    guides(fill = FALSE, color = FALSE, point_color = FALSE)

ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.png", 
       height = 7, width = 10)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.pdf", 
       height = 7, width = 10)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.svg", 
       height = 7, width = 10)

mut_genes_df_allcells_filt_list <- list()
for (don in names(de_res[["sce_list_unst"]])) {
  cat("working on ", don, "\n")
  sites_tmp <- sites_by_line_gr[sites_by_line_gr$donor_short_id == don]
  ov_tmp <- findOverlaps(sce_de_list_gr[[don]], sites_tmp)
  sce_tmp <- de_res[["sce_list_unst"]][[don]][queryHits(ov_tmp),]
  sites_tmp <- sites_tmp[subjectHits(ov_tmp)]
  sites_tmp$gene <- rownames(sce_tmp)
  dge_tmp <- de_res[["dge_list"]][[don]]
  dge_tmp <- dge_tmp[intersect(rownames(dge_tmp), sites_tmp$gene),]
  base_design <- dge_tmp$design[, !grepl("assigned", colnames(dge_tmp$design))]
  de_tbl_tmp <- data.frame(line = don,
                           gene = sites_tmp$gene, 
                           hgnc_symbol = gsub(".*_", "", sites_tmp$gene),
                           ensembl_gene_id = gsub("_.*", "", sites_tmp$gene),
                           var_id = sites_tmp$var_id,
                           location = sites_tmp$Location,
                           consequence = sites_tmp$Consequence,
                           clone_presence = sites_tmp$clone_presence,
                           logFC = NA, logCPM = NA, F = NA, PValue = NA,
                           comment = "", stringsAsFactors = FALSE)
  for (i in seq_len(length(sites_tmp))) {
    clones_tmp <- strsplit(sites_tmp$clone_presence[i], split = "&")[[1]]
    mutatedclone <- as.numeric(sce_tmp$assigned %in% clones_tmp)
    dsgn_tmp <- cbind(base_design, data.frame(mutatedclone))
    if (sites_tmp$gene[i] %in% rownames(dge_tmp) && is.fullrank(dsgn_tmp)
        && (sum(mutatedclone) > 9.5) && (sum(!mutatedclone) > 9.5)) {
      qlfit_tmp <- glmQLFit(dge_tmp[sites_tmp$gene[i],], dsgn_tmp)
      de_tmp <- glmQLFTest(qlfit_tmp, coef = ncol(dsgn_tmp))
      de_tbl_tmp$logFC[i] <- de_tmp$table$logFC
      de_tbl_tmp$logCPM[i] <- de_tmp$table$logCPM
      de_tbl_tmp$F[i] <- de_tmp$table$F
      de_tbl_tmp$PValue[i] <- de_tmp$table$PValue
    }
    if ((sum(mutatedclone) < 9.5) || (sum(!mutatedclone) < 9.5))
      de_tbl_tmp$comment[i] <- "minimum group size < 10 cells"
    if (!(sites_tmp$gene[i] %in% rownames(dge_tmp)))
      de_tbl_tmp$comment[i] <- "gene did not pass DE filters"
    if (!is.fullrank(dsgn_tmp))
      de_tbl_tmp$comment[i] <- "insufficient cells assigned to clone"
  }
  mut_genes_df_allcells_filt_list[[don]] <- de_tbl_tmp
}

mut_genes_df_allcells_filt <- do.call("bind_rows", mut_genes_df_allcells_filt_list)

## add FDRs for genes tested here for DE
ihw_res_all <- ihw(PValue ~ logCPM, data = mut_genes_df_allcells_filt, alpha = 0.2)
mut_genes_df_allcells_filt$FDR <- adj_pvalues(ihw_res_all)
## add simplified consequence categories
mut_genes_df_allcells_filt$consequence_simplified <- 
    mut_genes_df_allcells_filt$consequence
mut_genes_df_allcells_filt$consequence_simplified[
    mut_genes_df_allcells_filt$consequence_simplified %in% 
        c("stop_retained", "start_lost", "stop_lost", "stop_gained")] <- "nonsense"
mut_genes_df_allcells_filt$consequence_simplified[
    mut_genes_df_allcells_filt$consequence_simplified %in% 
        c("splice_donor", "splice_acceptor", "splice_region")] <- "splicing"
table(mut_genes_df_allcells_filt$consequence_simplified)

dplyr::arrange(mut_genes_df_allcells_filt, FDR) %>% 
  dplyr::select(-location, line, hgnc_symbol, var_id, consequence, FDR, PValue,
                everything()) %>% 
  DT::datatable(., options = list(pageLength = 20))

tmp5 <- mut_genes_df_allcells_filt %>%
    group_by(consequence_simplified) %>%
    summarise(med = median(logFC, na.rm = TRUE),
              nvars = n()
    )
dplyr::arrange(tmp5, med)

mut_genes_df_allcells_filt %>%
    dplyr::filter(!is.na(logFC)) %>%
    dplyr::mutate(consequence_simplified = factor(
        consequence_simplified, 
        levels(as.factor(consequence_simplified))[order(tmp5[["med"]])]),
        de  = ifelse(FDR < 0.2, "FDR < 0.2", "FDR > 0.2")) %>%
ggplot(aes(y = logFC, x = consequence_simplified, fill = consequence)) +
    geom_hline(yintercept = 0, linetype = 1, colour = "black") +
    ggbeeswarm::geom_quasirandom(colour = "gray40", pch = 21) +
    geom_boxplot(outlier.size = 0, outlier.alpha = 0, fill = "wheat", alpha = 0.5,
                 colour = "firebrick4", width = 0.2, size = 1) +
    scale_x_discrete(expand = c(0.1, 0.1), name = "consequence") +
    scale_y_continuous(expand = c(0.1, 0.1), name = "logFC") +
    theme_ridges(20) +
    coord_flip() +
    theme(strip.background = element_rect(fill = "gray90")) +
    guides(color = FALSE, point_color = FALSE) +
    scale_fill_brewer(palette = "Paired") +
    theme(legend.position = "bottom")

ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.png", 
       height = 9, width = 15)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.pdf", 
       height = 9, width = 15)
ggsave("figures/mutated_genes/alllines_mutgenes_logfc-box_by_simple_vep_anno_allcells.svg", 
       height = 9, width = 15)

## run DE for mutated cells vs unmutated cells using existing DE results
## fitting PCs as covariates in DE model
mut_genes_df_allcells_list_PCreg <- list()
for (don in names(de_res[["sce_list_unst"]])) {
    cat("working on ", don, "\n")
    sites_tmp <- sites_by_line_gr[sites_by_line_gr$donor_short_id == don]
    ov_tmp <- findOverlaps(sce_de_list_gr[[don]], sites_tmp)
    sce_tmp <- de_res[["sce_list_unst"]][[don]][queryHits(ov_tmp),]
    sites_tmp <- sites_tmp[subjectHits(ov_tmp)]
    sites_tmp$gene <- rownames(sce_tmp)
    dge_tmp <- de_res[["dge_list"]][[don]]
    dge_tmp <- dge_tmp[intersect(rownames(dge_tmp), sites_tmp$gene),]
    base_design <- dge_tmp$design[, !grepl("assigned", colnames(dge_tmp$design))]
    de_tbl_tmp <- data.frame(line = don,
                             gene = sites_tmp$gene, 
                             hgnc_symbol = gsub(".*_", "", sites_tmp$gene),
                             ensembl_gene_id = gsub("_.*", "", sites_tmp$gene),
                             var_id = sites_tmp$var_id,
                             location = sites_tmp$Location,
                             consequence = sites_tmp$Consequence,
                             clone_presence = sites_tmp$clone_presence,
                             logFC = NA, logCPM = NA, F = NA, PValue = NA,
                             comment = "", stringsAsFactors = FALSE)
    pca_tmp <- reducedDim(runPCA(
        de_res[["sce_list_unst"]][[don]], ntop = 500, ncomponents = 5), "PCA")
    for (i in seq_len(length(sites_tmp))) {
        clones_tmp <- strsplit(sites_tmp$clone_presence[i], split = "&")[[1]]
        mutatedclone <- as.numeric(sce_tmp$assigned %in% clones_tmp)
        dsgn_tmp <- cbind(base_design, pca_tmp[, 1], data.frame(mutatedclone))
        if (sites_tmp$gene[i] %in% rownames(dge_tmp) && is.fullrank(dsgn_tmp)) {
            qlfit_tmp <- glmQLFit(dge_tmp[sites_tmp$gene[i],], dsgn_tmp)
            de_tmp <- glmQLFTest(qlfit_tmp, coef = ncol(dsgn_tmp))
            de_tbl_tmp$logFC[i] <- de_tmp$table$logFC
            de_tbl_tmp$logCPM[i] <- de_tmp$table$logCPM
            de_tbl_tmp$F[i] <- de_tmp$table$F
            de_tbl_tmp$PValue[i] <- de_tmp$table$PValue
        }
        if (!(sites_tmp$gene[i] %in% rownames(dge_tmp)))
            de_tbl_tmp$comment[i] <- "gene did not pass DE filters"
        if (!is.fullrank(dsgn_tmp))
            de_tbl_tmp$comment[i] <- "insufficient cells assigned to clone"
    }
    mut_genes_df_allcells_list_PCreg[[don]] <- de_tbl_tmp
}

mut_genes_df_allcells_PCreg <- do.call("bind_rows", mut_genes_df_allcells_list_PCreg)
## add FDRs for genes tested here for DE
ihw_res_all <- ihw(PValue ~ logCPM, data = mut_genes_df_allcells_PCreg, alpha = 0.2)
mut_genes_df_allcells_PCreg$FDR <- adj_pvalues(ihw_res_all)
## add simplified consequence categories
mut_genes_df_allcells_PCreg$consequence_simplified <- 
    mut_genes_df_allcells_PCreg$consequence
mut_genes_df_allcells_PCreg$consequence_simplified[
    mut_genes_df_allcells_PCreg$consequence_simplified %in% 
        c("stop_retained", "start_lost", "stop_lost", "stop_gained")] <- "nonsense"
mut_genes_df_allcells_PCreg$consequence_simplified[
    mut_genes_df_allcells_PCreg$consequence_simplified %in% 
        c("splice_donor", "splice_acceptor", "splice_region")] <- "splicing"
table(mut_genes_df_allcells_PCreg$consequence_simplified)

dplyr::arrange(mut_genes_df_allcells_PCreg, FDR) %>% 
  dplyr::select(-location) %>% 
   DT::datatable(., options = list(pageLength = 20))

data_frame(PValue_default = mut_genes_df_allcells[["PValue"]],
           PValue_PCreg = mut_genes_df_allcells_PCreg[["PValue"]],
           sig_default = mut_genes_df_allcells[["FDR"]] < 0.1,
           sig_PCreg = mut_genes_df_allcells_PCreg[["FDR"]] < 0.1) %>%
    ggplot(aes(x = -log10(PValue_default), y = -log10(PValue_PCreg),
               shape = sig_default, colour = sig_PCreg)) +
    geom_point() +
    geom_smooth(aes(group = 1), colour = "firebrick") +
    geom_abline(slope = 1, intercept = 0) +
    scale_color_manual(values = c("gray60", "black"))

## run DE for mutated cells vs unmutated cells using existing DE results
## fitting PCs as covariates in DE model
mut_genes_df_allcells_list_CCreg <- list()
for (don in names(de_res[["sce_list_unst"]])) {
    cat("working on ", don, "\n")
    sites_tmp <- sites_by_line_gr[sites_by_line_gr$donor_short_id == don]
    ov_tmp <- findOverlaps(sce_de_list_gr[[don]], sites_tmp)
    sce_tmp <- de_res[["sce_list_unst"]][[don]][queryHits(ov_tmp),]
    sites_tmp <- sites_tmp[subjectHits(ov_tmp)]
    sites_tmp$gene <- rownames(sce_tmp)
    dge_tmp <- de_res[["dge_list"]][[don]]
    dge_tmp <- dge_tmp[intersect(rownames(dge_tmp), sites_tmp$gene),]
    base_design <- dge_tmp$design[, !grepl("assigned", colnames(dge_tmp$design))]
    de_tbl_tmp <- data.frame(line = don,
                             gene = sites_tmp$gene, 
                             hgnc_symbol = gsub(".*_", "", sites_tmp$gene),
                             ensembl_gene_id = gsub("_.*", "", sites_tmp$gene),
                             var_id = sites_tmp$var_id,
                             location = sites_tmp$Location,
                             consequence = sites_tmp$Consequence,
                             clone_presence = sites_tmp$clone_presence,
                             logFC = NA, logCPM = NA, F = NA, PValue = NA,
                             comment = "", stringsAsFactors = FALSE)
    G1 <- de_res[["sce_list_unst"]][[don]]$G1
    S <- de_res[["sce_list_unst"]][[don]]$S
    G2M <- de_res[["sce_list_unst"]][[don]]$G2M
    for (i in seq_len(length(sites_tmp))) {
        clones_tmp <- strsplit(sites_tmp$clone_presence[i], split = "&")[[1]]
        mutatedclone <- as.numeric(sce_tmp$assigned %in% clones_tmp)
        dsgn_tmp <- cbind(base_design, G1, S, G2M, data.frame(mutatedclone))
        if (sites_tmp$gene[i] %in% rownames(dge_tmp) && is.fullrank(dsgn_tmp)) {
            qlfit_tmp <- glmQLFit(dge_tmp[sites_tmp$gene[i],], dsgn_tmp)
            de_tmp <- glmQLFTest(qlfit_tmp, coef = ncol(dsgn_tmp))
            de_tbl_tmp$logFC[i] <- de_tmp$table$logFC
            de_tbl_tmp$logCPM[i] <- de_tmp$table$logCPM
            de_tbl_tmp$F[i] <- de_tmp$table$F
            de_tbl_tmp$PValue[i] <- de_tmp$table$PValue
        }
        if (!(sites_tmp$gene[i] %in% rownames(dge_tmp)))
            de_tbl_tmp$comment[i] <- "gene did not pass DE filters"
        if (!is.fullrank(dsgn_tmp))
            de_tbl_tmp$comment[i] <- "insufficient cells assigned to clone"
    }
    mut_genes_df_allcells_list_CCreg[[don]] <- de_tbl_tmp
}

mut_genes_df_allcells_CCreg <- do.call("bind_rows", mut_genes_df_allcells_list_CCreg)
## add FDRs for genes tested here for DE
ihw_res_all <- ihw(PValue ~ logCPM, data = mut_genes_df_allcells_CCreg, alpha = 0.2)
mut_genes_df_allcells_CCreg$FDR <- adj_pvalues(ihw_res_all)
## add simplified consequence categories
mut_genes_df_allcells_CCreg$consequence_simplified <- 
    mut_genes_df_allcells_CCreg$consequence
mut_genes_df_allcells_CCreg$consequence_simplified[
    mut_genes_df_allcells_CCreg$consequence_simplified %in% 
        c("stop_retained", "start_lost", "stop_lost", "stop_gained")] <- "nonsense"
mut_genes_df_allcells_CCreg$consequence_simplified[
    mut_genes_df_allcells_CCreg$consequence_simplified %in% 
        c("splice_donor", "splice_acceptor", "splice_region")] <- "splicing"
dplyr::arrange(mut_genes_df_allcells_CCreg, FDR) %>% 
  dplyr::select(-location) %>% 
   DT::datatable(., options = list(pageLength = 20))

data_frame(PValue_default = mut_genes_df_allcells[["PValue"]],
           PValue_CCreg = mut_genes_df_allcells_CCreg[["PValue"]],
           sig_default = mut_genes_df_allcells[["FDR"]] < 0.1,
           sig_PCreg = mut_genes_df_allcells_CCreg[["FDR"]] < 0.1) %>%
    ggplot(aes(x = -log10(PValue_default), y = -log10(PValue_CCreg),
               shape = sig_default, colour = sig_PCreg)) +
    geom_point() +
    geom_smooth(aes(group = 1), colour = "firebrick") +
    geom_abline(slope = 1, intercept = 0) +
    scale_color_manual(values = c("gray60", "black"))

sites_by_line %>% 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_manual(values = colorRampPalette(brewer.pal(8, "Accent"))(18)) +
    guides(colour = FALSE) +
    ggtitle("Clone tagging variants by consequence class: all lines") +
    xlab("number of variants") + ylab("consequence") +
    theme_bw(16)

sites_by_line %>% group_by(Consequence, donor_short_id) %>%
    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(~donor_short_id) +
#    scale_color_brewer(palette = "Set2") +
    scale_color_manual(values = colorRampPalette(brewer.pal(8, "Accent"))(18)) +
    guides(colour = FALSE) +
    ggtitle("Clone tagging variants by consequence class: all lines") +
    xlab("number of variants") + ylab("consequence") +
    theme_bw(16)

devtools::session_info()