Add real data heatmaps

analysis/visualise-marker-genes.Rmd

@@ -13,17 +13,83 @@ library(SingleCellExperiment)
 library(scater)
 library(scran)
 library(dplyr)
-library(schex)
+library(tidyr)
+library(pheatmap)
 library(patchwork)
 
 source(here::here("code", "find-marker-genes.R"))
+source(here::here("code", "top-genes.R"))
+source(here::here("code", "analysis-utils.R"))
 ```
 
 ```{r load-data}
 config <- yaml::read_yaml(here::here("config.yaml"))
-sim <- readRDS(here::here(config$sim_data_folder, "standard_sim_1.rds"))
+res_paths <- here::here(list.files(config$results_folder, full.names = TRUE))
+
+sel_genes <- list()
+for (i in seq_along(res_paths)) {
+  print(i)
+  res <- readRDS(res_paths[[i]])
+  # FIXME: Some of the runs error.
+  if (!(length(res$result) == 0)) {
+    sel_genes[[i]] <- reformat_found_mgs(res)
+  }
+  rm(res)
+}
+
+res_names <- substr(basename(res_paths), 1, nchar(basename(res_paths)) - 4)
+selected_genes <- tibble(res_name = res_names, sel_genes = sel_genes)
+
+mgs_pars <- retrive_simulation_parameters() %>% 
+  mutate(res_name = substr(file_name, 1, nchar(file_name) - 4)) %>% 
+  left_join(selected_genes, by = "res_name")
+
+mgs_data <- mgs_pars %>% 
+  dplyr::rename(sel_mgs = sel_genes) %>% 
+  # HACK
+  rowwise() %>% 
+  mutate(group_id = list(names(sel_mgs))) %>% 
+  unchop(c(sel_mgs, group_id))
+
+sim <- readRDS(here::here("data", "sim_data", "standard_sim_1.rds"))
 ```
 
+```{r}
+top_mgs <- mgs_data %>% 
+  filter(sim_label == "standard_sim" & rep == 1 & group_id == "Group1") %>% 
+  filter(pars == "seurat_DESeq2") %>% 
+  rowwise() %>% 
+  mutate(top_genes = list(
+    sel_mgs %>% 
+      arrange(p_value) %>% 
+      slice(1:10) %>% 
+      pull(gene)
+    ) 
+  ) %>% 
+  ungroup() %>% 
+  pull(top_genes) %>% 
+  purrr::pluck(1)
+
+mgs_sce <- sim[top_mgs, ]
+plot_data <- tibble(group = colLabels(mgs_sce), 
+                    as.data.frame(t(logcounts(mgs_sce))))
+
+heatmap_data <- plot_data %>% 
+  group_by(group) %>% 
+  summarise(across(everything(), mean)) %>% 
+  tibble::column_to_rownames("group") %>% 
+  as.matrix() %>% 
+  t()
+
+pheatmap(heatmap_data)
+pheatmap(heatmap_data,  breaks=seq(0, 2000, length.out=101))
+```
+
+```{r}
+res <- readRDS(here::here("results", "standard_sim_1-scran_t_all.rds"))
+test <- getMarkerEffects(res$result[[1]])[1:10, ]
+pheatmap(getMarkerEffects(res$result[[1]])[1:10, ])
+```
 
 ```{r find-marker-genes}
 umgs <- find_unique_marker_genes(sim)
@@ -48,10 +114,9 @@ plotTSNE(sim, colour_by = "Gene5425")
 
 ```{r}
 plotExpression(sim, x = "label", features = "Gene5425")
-plotExpression(sim, x = "label", features = "Gene3112")
-plotExpression(sim, x = "label", features = "Gene1056")
-plotExpression(sim, x = "label", features = "Gene9300")
-```
+plotExpression(sim, x = "label", features = "Gene4904")
+plotExpression(sim, x = "label", features = "Gene9834")
+```   
 
 ```{r}
 scran_mgs <- findMarkers(sim, pval.type = "all")
@@ -73,7 +138,7 @@ top_scran_mgs <- scran_mgs[[1]] %>%
 Much better looking...
 
 ```{r}
-plotExpression(sim, x = "label", features = "Gene4350")
+plotExpression(sim, x = "label", features = "Gene3991")
 plotExpression(sim, x = "label", features = "Gene3348")
 plotExpression(sim, x = "label", features = "Gene360")
 plotExpression(sim, x = "label", features = "Gene3452")
@@ -107,9 +172,10 @@ sim
 ```
 
 ```{r visuliase-detected-marker-genes}
-plot_top_marker_genes <- function(res, n, clus) {
+plot_top_marker_genes <- function(res, data, n, clus) {
   stopifnot(length(clus) == 1)
-  mgs <- get_top_genes(res, n)[[clus]]
+  mgs <- reformat_found_mgs(res, n)[[clus]]
+  mgs <- mgs$gene
   plots <- lapply(mgs, function(x) {
   plotTSNE(data, colour_by = x) +
     ggtitle(x) + 
@@ -123,9 +189,9 @@ poisson_res <- readRDS(here::here("results", "standard_sim_1-seurat_poisson.rds"
 deseq2_res <- readRDS(here::here("results", "standard_sim_1-seurat_DESeq2.rds"))
 scran_t_all_res <- readRDS(here::here("results", "standard_sim_1-scran_t_all.rds"))
 
-plot_top_marker_genes(poisson_res, 10, 1)
-plot_top_marker_genes(scran_t_all_res, 10, 1)
-plot_top_marker_genes(deseq2_res, 10, 1)
+plot_top_marker_genes(poisson_res, data, 10, 1)
+plot_top_marker_genes(scran_t_all_res, data, 10, 1)
+plot_top_marker_genes(deseq2_res, data, 10, 1)
 
 find_unique_marker_genes(data)[[1]] %>% 
   arrange(desc(fc)) %>%