selected_mutations_diffexp_and_ORA.Rmd

---
title: "Mutational analysis: differential expression + ORA"
author: "Camilo Posso"
date: "02/16/2022"
output: 
  html_document:
    code_folding: hide
    toc: true
editor_options: 
  chunk_output_type: inline
---

## Goal


The goal of this markdown is to make the ORA and heatmap plots coming from the
differential expression for selected mutations.


```{r include=FALSE}
library(dplyr)
library(ggplot2)
library(kableExtra)
library(gridExtra)
library(grid)

source("../util/synapseUtil.R")
source("../util/loading_data.R")
source("../util/mutational_analysis_helper.R")
source("../util/make_plots_util.R")


load("../Misc/load.combined.data 3-09-2022.RData")
# load.combined.data()
metadata <- load.metadata()
samples <- global.data %>%
  pull(Barcode.ID) %>% unique()

total_mutations <- WES.data %>%
  pull(Gene) %>% unique()

all_mutation_data <- load_mutational_sample_data()

IDH1_IDH2 <- "IDH1+IDH2"
KRAS_NRAS <- "KRAS+NRAS"

```



```{r}

## Mutations with at least 4 patients having said mutation
WES_reduced <- WES.data %>%
  filter(Barcode.ID %in% samples) %>% 
  group_by(Barcode.ID, Gene) %>%
  dplyr::slice(1) %>%
  ungroup(Barcode.ID, Gene)

mutation_count <- WES_reduced %>%
  group_by(Gene) %>%
  summarise(n = n())

mutations <-  mutation_count %>%
  filter(n > 3) %>%
  pull(Gene)

WES_samples <- WES_reduced %>%
  group_by(Gene) %>%
  mutate(n = n()) %>%
  arrange(-n) %>%
  pull(Barcode.ID) %>%
  unique()

WES_reduced <- WES_reduced %>%
  mutate(value = 1) %>%
  dplyr::select(Barcode.ID, Gene, value)


m.WES <- make.msnset(WES_reduced, "Gene", "Barcode.ID", "value", meta)  
mat <- exprs(m.WES)[, WES_samples]
mat[is.na(mat)] <- 0

```


Of the `r length(total_mutations)` mutations recorded in at least one sample, `r length(mutations)`
of them are present in at least 4 samples. Below we summarize the mutational data for these
`r length(mutations)` mutations using a heatmap. Note the columns represent samples, while the
light beige color denotes mutation.


```{r}
pheatmap(mat[mutations, ], color = c("#16182A", "#F5E9D9"), 
         legend = FALSE, show_colnames = FALSE)
```


## Differential Expression

Focusing on any one of the `r length(mutations)`, say \(\mu\), we can ask. Which phospho-sites 
are differentially expressed in \(\mu\) mutated samples vs non mutated samples? What if we use
global or RNA data? Below we collect this information into a table.


```{r include=FALSE}

m.RNA <- make.msnset(RNA.data, feature.col = "Gene",
                     value.col = "RNA counts", metadata = meta)
m.global <- make.msnset(global.data, feature.col = "Gene",
                        value.col = "LogRatio", metadata = meta)
m.phospho <- make.msnset(phospho.data, feature.col = "SiteID",
                         value.col = "LogRatio", metadata = meta)


datasets <- list(m.global, m.phospho, m.RNA)
names(datasets) <- c("global", "phospho", "RNA")

mutations <- c(mutations, c(IDH1_IDH2, KRAS_NRAS, "NPM1_clinical", "FLT3.ITD"))
## NPM1 mutations are made from clinical summary data used instead of the WES data
## so we remove NPM1 here
NPM1_WES_index <- which(mutations == "NPM1")
mutations <- mutations[-NPM1_WES_index]

differential_expression_results <- lapply(mutations, check_diff_exp_individual, 
                                          all_mutation_data, datasets)

summary_table <- lapply(differential_expression_results, function(xx){xx[[1]]}) %>% 
  do.call("rbind", .)

names(differential_expression_results) = sapply(differential_expression_results, function(xx){
  xx[[1]]$mutation[[1]]
})

saveRDS(differential_expression_results, "./raw_tables/all_mutations_diffexp_results.RDS")
write.table(summary_table, "raw_tables/all_mutations_diffexp_summary.txt", sep = "\t")

```



```{r}
differential_expression_results <- readRDS("./raw_tables/all_mutations_diffexp_results.RDS")
summary_table <- read.table("raw_tables/all_mutations_diffexp_summary.txt", sep = "\t")
colnames(summary_table)[3:5] <- gsub("\\.", " ", colnames(summary_table)[3:5])
selected_mutations <- c("IDH1", "IDH2", "TP53", "NRAS", "KRAS", "FLT3.ITD", "NPM1_clinical", 
                        IDH1_IDH2, KRAS_NRAS)

combined_diffexp <- lapply(selected_mutations, function(mutation){
    dff <- differential_expression_results[[mutation]][[2]] %>%
      mutate(mutation = mutation)
  }) %>% do.call("rbind", .)

# write.table(combined_diffexp, "./raw_tables/selected_mutations_combined_differential_expression_results.txt",
#             sep = "\t", quote = F)
# 
# upload.plot("selected_mutations_combined_differential_expression_results.txt", 
#             parentId = "syn26529351")
#
# combined_diffexp <- read.table("raw_tables/selected_mutations_combined_differential_expression_results.txt", sep = "\t")

```



```{r}
plot_colors <- data.frame("Global" = data_type_colors[[1]],
                 "Phospho" = data_type_colors[[2]],
                 "RNA" = data_type_colors[[3]])
plot_df <- summary_table %>%
  select(-`Number.mutated`) %>%
  pivot_longer(-mutation) %>%
  mutate(name = sub("^.*in ", "", name)) %>%
  mutate(name = case_when(name == "global" ~ "Global",
                          name == "phospho" ~ "Phospho",
                          name == "RNA" ~ "RNA")) %>%
  mutate(individual_bar = paste(mutation, name)) %>%
  dplyr::rename(`Data type` = name) %>%
  arrange(-value) %>%
  mutate(mutation = factor(mutation, levels = unique(mutation)))

p1 <- ggplot(plot_df, aes(y = mutation, x = value, fill = `Data type`)) + 
  geom_bar(stat = 'identity', pos = position_dodge(), width = 0.7) + xlab("Features differentially expressed") +
  scale_fill_manual(values = plot_colors) + scale_x_continuous(trans = "log1p", breaks = c(0, 10, 100, 1000, 4000)) + 
  theme(legend.position='left') + ylab("")


count_df <- summary_table %>%
  select(mutation, Number.mutated) %>%
  mutate(mutation = factor(mutation, levels = unique(plot_df$mutation)))
p2 <- ggplot(count_df, aes(y = mutation, x = Number.mutated/210)) +
  geom_bar(stat = 'identity', width = 0.4) + xlab("Fraction of mutated samples") + ylab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())


lay <- rbind(c(1,1,1,1,1,1,1,2,2,2,2,2), 
             c(1,1,1,1,1,1,1,2,2,2,2,2))
top_grob <- textGrob("Differential expression by mutation", 
                     gp = gpar(fontsize = 20))
p_all <- arrangeGrob(p1, p2, top = top_grob, layout_matrix = lay)


ggsave(filename = "./plots/mutational_differential_expression_barplot_count.pdf", 
       plot = p_all, width = 15, height = 10)
upload.plot("./plots/mutational_differential_expression_barplot_count.pdf", 
            parentId = "syn26529351")
```



```{r}
library(DT)
DT::datatable(summary_table)
```


```{r eval=FALSE, include=FALSE}
### Run this to make all the heatmaps. A few will fail, due to the exactly one feature
### being selected to plot.** or failure to compute distances due to many NA
dataset_names <- names(datasets)

for (i in 1:length(differential_expression_results)){
  mutation = names(differential_expression_results)[[i]]
  diff_exp <- differential_expression_results[[i]][[2]]
  
  
  features <- lapply(1:3, function(j){
    diff_exp %>%
      filter(datatype == dataset_names[[j]]) %>%
      filter(adj.P.Val < 0.05) %>%
      head(20) %>%
      pull(feature)
  })
  
  for (j in 1:3){
    if (length(features[[j]]) > 0){
      try(plot_diffexp_heatmap(mutation, all_mutation_data, datasets[j], features[[j]], 
                               folder = "./plots/differential_expression_heatmaps/", 
                               suffix = dataset_names[[j]]))
    }
  }
}

```



```{r}
## Setup the databases to use in the tests
library(clusterProfiler)
library(msigdbr)

t2g_hallmark <- msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol)
## 50 unique terms
# unique(t2g_hallmark$gs_name)

t2g_oncogenic <- msigdbr(species = "Homo sapiens", category = "C6") %>% 
  dplyr::select(gs_name, gene_symbol)
## 189 unique terms
# unique(t2g_oncogenic$gs_name)

t2g_gobp <- msigdbr(species = "Homo sapiens", category = "C5", subcategory = "BP") %>% 
  dplyr::select(gs_name, gene_symbol) %>%
  mutate(gs_name = gsub("_", " ", gs_name))

```




```{r}
## Here we compute the enrichment results

selected_mutations <- c("IDH1", "IDH2", "TP53", "NRAS", "KRAS", "FLT3.ITD", "NPM1_clinical", 
                        IDH1_IDH2, KRAS_NRAS)


## need same order as datasets here
universes <- list(unique(global.data$Gene), unique(phospho.data$Gene), 
                  unique(RNA.data$Gene))
names(universes) <- names(datasets)

gsea_ora_results <- lapply(selected_mutations, function(mutation){
  print(mutation)
  gsea_all = gsea_global = gsea_rna <- data.frame(Description = c(), 
                                                  GeneRation = c(), 
                                                  BgRatio = c(), 
                                                  pvalue = c(), 
                                                  p.adjust = c(), 
                                                  Count = c(),
                                                  DB = c())
  gsea_global_go <- gsea_global_hallmark <- gsea_global_oncogenic <- gsea_all
  gsea_rna_go <- gsea_rna_hallmark <- gsea_rna_oncogenic <- gsea_all
  
  genes <- differential_expression_results[[mutation]][[2]] %>%
    filter(adj.P.Val < 0.05) %>%
    filter(datatype == "global") %>%
    pull(feature)
  
  if (length(genes) > 0){
    gsea_global_go <- try(enricher(genes, TERM2GENE = t2g_gobp, universe = universes[[1]])@result %>%
                            dplyr::rename(core_enrichment = geneID) %>%
                            compress_enrichment(colname = "p.adjust", descending = FALSE) %>%
                            dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                            mutate(mutation = mutation, data_type = "global", DB = "GO_BP"))
    
    gsea_global_hallmark <- try(enricher(genes, TERM2GENE = t2g_hallmark, universe = universes[[1]])@result %>%
                                  dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                                  mutate(mutation = mutation, data_type = "global", DB = "hallmark"))
    
    gsea_global_oncogenic <- try(enricher(genes, TERM2GENE = t2g_oncogenic, universe = universes[[1]])@result %>%
                                   dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                                   mutate(mutation = mutation, data_type = "global", DB = "oncogenic"))
    
    gsea_global <- rbind(gsea_global_go, gsea_global_hallmark, gsea_global_oncogenic)
  }
  
  genes <- differential_expression_results[[mutation]][[2]] %>%
    filter(adj.P.Val < 0.05) %>%
    filter(datatype == "RNA") %>%
    pull(feature)

  if (length(genes) > 0){
    gsea_rna_go <- try(enricher(genes, TERM2GENE = t2g_gobp, universe = universes[[3]])@result %>%
                         dplyr::rename(core_enrichment = geneID) %>%
                         compress_enrichment(colname = "p.adjust", descending = FALSE) %>%
                         dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                         mutate(mutation = mutation, data_type = "RNA", DB = "GO_BP"))
    
    gsea_rna_hallmark <- try(enricher(genes, TERM2GENE = t2g_hallmark, universe = universes[[3]])@result %>%
                               dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                               mutate(mutation = mutation, data_type = "RNA", DB = "hallmark"))
    
    gsea_rna_oncogenic <- try(enricher(genes, TERM2GENE = t2g_oncogenic, universe = universes[[3]])@result %>%
                                dplyr::select(Description, GeneRatio, BgRatio, pvalue, p.adjust, Count) %>%
                                mutate(mutation = mutation, data_type = "RNA", DB = "oncogenic"))
    
    gsea_rna <- rbind(gsea_rna_go, gsea_rna_hallmark, gsea_rna_oncogenic)
  }
  
  gsea_all <- rbind(gsea_all, gsea_global) %>%
    rbind(gsea_rna) %>%
    mutate(mutation = mutation)
  
  return(gsea_all)
}) 

gsea_ora_results_df <- do.call("rbind", gsea_ora_results) %>%
  filter(!grepl("Error", p.adjust)) %>%
  mutate(p.adjust = as.numeric(p.adjust),
         Count = as.numeric(Count))

write.table(gsea_ora_results_df, "raw_tables/mutation_analysis_ORA_results_combined.txt", sep = "\t", quote = F)
# upload.plot("raw_tables/mutation_analysis_ORA_results_combined.txt", parentId = "syn27220694")

```



```{r, ora plots for the selected mutations}
## Here we make the plots using the enrichment results
to_plot <- gsea_ora_results_df %>%
  filter(p.adjust < 0.05) 

categories_plot <- to_plot %>%
  dplyr::select(mutation, data_type, DB) %>%
  unique()
rownames(categories_plot) <- 1:nrow(categories_plot)

lapply(1:nrow(categories_plot), function(i){
  plot_mutation <- categories_plot$mutation[[i]]
  plot_data_type <- categories_plot$data_type[[i]]
  plot_DB <- categories_plot$DB[[i]]
  plot_title <- paste(plot_DB, plot_mutation, plot_data_type)
  
  plot_path <- "./plots/over_representation_analysis/ORA"
  plot_path <- paste(plot_path, plot_DB, plot_mutation, plot_data_type, sep = "_") %>%
    paste0(".png")
  plot_df <- to_plot %>%
    filter(mutation == plot_mutation) %>%
    filter(data_type == plot_data_type) %>%
    filter(DB == plot_DB) %>%
    dplyr::select(Description, Count, p.adjust) %>%
    dplyr::rename(enrichment = Count,
                  pathway = Description,
                  adj_p_val = p.adjust) %>%
    head(20)
  
  png(file = plot_path, width = 1200, height = 800)
  plot_enrichment_result(plot_df, "Count", enrichment_title = plot_title)
  dev.off()
  
  ## uploading to synapse
  upload.plot(plot_path, parentId = "syn27220694")
  
  "done"
})

```



```{r}


```