Rewrite visualisation of marker genes

analysis/visualise-marker-genes.Rmd

@@ -14,6 +14,7 @@ library(scater)
 library(scran)
 library(dplyr)
 library(schex)
+library(patchwork)
 
 source(here::here("code", "find-marker-genes.R"))
 ```
@@ -24,16 +25,6 @@ sim <- readRDS(here::here(config$sim_data_folder, "standard_sim_1.rds"))
 ```
 
 
-```{r process-data}
-# FIXME: Remove when the simulations are rerun!
-clusters <- quickCluster(sim)
-sim <- computeSumFactors(sim, clusters = clusters)
-sim <- logNormCounts(sim)
-sim <- runPCA(sim)
-sim <- runTSNE(sim, dimred = "PCA")
-```
-
-
 ```{r find-marker-genes}
 umgs <- find_unique_marker_genes(sim)
 # FIXME: `find_marker_genes` should do this.
@@ -52,11 +43,11 @@ umgs[[1]] %>%
   slice_head(n = 10)
 
 plotTSNE(sim, colour_by = "Group")
-plotTSNE(sim, colour_by = "Gene7047")
+plotTSNE(sim, colour_by = "Gene5425")
 ```
 
 ```{r}
-plotExpression(sim, x = "label", features = "Gene7047")
+plotExpression(sim, x = "label", features = "Gene5425")
 plotExpression(sim, x = "label", features = "Gene3112")
 plotExpression(sim, x = "label", features = "Gene1056")
 plotExpression(sim, x = "label", features = "Gene9300")
@@ -114,3 +105,33 @@ plot_hexbin_feature(
 sim
 
 ```
+
+```{r visuliase-detected-marker-genes}
+plot_top_marker_genes <- function(res, n, clus) {
+  stopifnot(length(clus) == 1)
+  mgs <- get_top_genes(res, n)[[clus]]
+  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, 10, 1)
+plot_top_marker_genes(scran_t_all_res, 10, 1)
+plot_top_marker_genes(deseq2_res, 10, 1)
+
+find_unique_marker_genes(data)[[1]] %>% 
+  arrange(desc(fc)) %>% 
+  slice_head(n = 10) %>% 
+  pull(gene)
+
+get_top_genes(poisson_res, 10)[[1]]
+
+```