Add more concordance analysis

analysis/method-concordance.Rmd

@@ -11,6 +11,8 @@ knitr::opts_chunk$set(echo = TRUE)
 
 ```{r libraries, message = FALSE, warning = FALSE}
 library(dplyr)
+library(ggupset)
+library(ggplot2)
 ```
 
 **Update: now using actual simulations**
@@ -26,20 +28,24 @@ NB: scran::findMarkers has top in any but not in all
 result_test_1 <- readRDS(here::here("results", "sim_1-scran_t_any.rds"))
 result_test_2 <- readRDS(here::here("results", "sim_1-scran_wilcox_any.rds"))
 
+
+
 group_1_1 <- result_test_1$result$Group1
 group_1_2 <- result_test_2$result$Group1
 
 test_genes_1 <- group_1_1 %>% 
   as.data.frame() %>% 
-  filter(Top < 4) %>% 
+  filter(Top < 6) %>% 
   rownames()
 
 test_genes_2 <- group_1_2 %>% 
   as.data.frame() %>% 
-  filter(Top < 4) %>% 
+  filter(Top < 6) %>% 
   rownames()
 
-intersect(test_genes_1, test_genes_2)
+length(test_genes_1)
+length(test_genes_2)
+length(intersect(test_genes_1, test_genes_2))
 ```
 
 No intersection!
@@ -64,12 +70,40 @@ length(test_genes_1)
 length(intersect(de_inds_group_1, test_genes_alt))
 length(intersect(de_inds_group_1, test_genes_1))
 
+test_genes_1
+
+
 ```
 
 Why is there an order of magnitude difference in found vs actual?
 
+Let's think more about what it means to be DE in this context...
+
+* Not all exactly 1000...
+
+```{r de-genes}
 
+group_names <- rep(names(de_genes), lengths(de_genes))
 
+test <- tibble(gene = unlist(de_genes), group_names) %>% 
+  group_by(gene) %>% 
+  summarise(group = list(group_names), .groups = "drop")
+
+test %>%
+  ggplot(aes(x = group)) +
+  geom_bar() +
+  scale_x_upset()
+
+de_genes <- sim_1@metadata$de_inds
+
+```
+
+An upset plot of all the Sim 1 methods (4 so far)
+
+```{r sim-1-upset}
+here::here("results", "sim-1-*")
+
+```