Remove additional unneeded analysis

analysis/visualise-marker-genes.Rmd

----
-title: "Visualise marker genes"
-output: html_document
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE)
-```
-
-```{r libraries}
-library(ggplot2)
-library(SingleCellExperiment)
-library(scater)
-library(scran)
-library(dplyr)
-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"))
-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)
-# FIXME: `find_marker_genes` should do this.
-umgs <- lapply(umgs, function(x) {
-  abs_fc <- ifelse(x$fc > 1, x$fc, 1 / x$fc)
-  direction <- ifelse(x$fc > 1, "up", "down")
-  cbind(x, abs_fc, direction)
-})
-```
-
-Concentrate on group 1 for now.
-
-```{r plot-mgs}
-umgs[[1]] %>% 
-  arrange(desc(abs_fc)) %>% 
-  slice_head(n = 10)
-
-plotTSNE(sim, colour_by = "Group")
-plotTSNE(sim, colour_by = "Gene5425")
-```
-
-```{r}
-plotExpression(sim, x = "label", features = "Gene5425")
-plotExpression(sim, x = "label", features = "Gene4904")
-plotExpression(sim, x = "label", features = "Gene9834")
-```   
-
-```{r}
-scran_mgs <- findMarkers(sim, pval.type = "all")
-genes <- rownames(scran_mgs[[1]])
-
-scran_mgs[[1]] %>% 
-  as_tibble() %>% 
-  bind_cols(gene = genes) %>% 
-  arrange(FDR)
-
-top_scran_mgs <- scran_mgs[[1]] %>% 
-  as_tibble() %>% 
-  bind_cols(gene = genes) %>% 
-  arrange(FDR) %>% 
-  head(n = 10) %>% 
-  pull(gene)
-```  
-
-Much better looking...
-
-```{r}
-plotExpression(sim, x = "label", features = "Gene3991")
-plotExpression(sim, x = "label", features = "Gene3348")
-plotExpression(sim, x = "label", features = "Gene360")
-plotExpression(sim, x = "label", features = "Gene3452")
-```
-
-```{r}
-umgs[[1]] %>% 
-  as_tibble() %>% 
-  filter(gene %in% top_scran_mgs)
-```
-
-```{r}
-umgs[[1]] %>% 
-  ggplot(aes(x = abs_fc)) +
-  geom_histogram(bins = 30) + 
-  theme_bw()
-```
-
-```{r schex}
-hex_sim <- make_hexbin(sim, nbins = 40, dimension_reduction = "TSNE")
-
-plot_hexbin_feature(
-  hex_sim, 
-  type = "logcounts", 
-  feature = "Gene7047", 
-  action = "median", 
-)
-
-sim
-
-```
-
-```{r visuliase-detected-marker-genes}
-plot_top_marker_genes <- function(res, data, n, clus) {
-  stopifnot(length(clus) == 1)
-  mgs <- reformat_found_mgs(res, n)[[clus]]
-  mgs <- mgs$gene
-  plots <- lapply(mgs, function(x) {
-  plotTSNE(data, colour_by = x) +
-    ggtitle(x) + 
-    theme(legend.position = "none")
-  })
-  wrap_plots(plots)
-}
-
-data <- readRDS(here::here("data", "sim_data", "standard_sim_1.rds"))
-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, 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)) %>% 
-  slice_head(n = 10) %>% 
-  pull(gene)
-
-get_top_genes(poisson_res, 10)[[1]]
-
-```