Draft additional concordance analysis

analysis/concordance.Rmd

@@ -160,7 +160,7 @@ dendrogram_plot <- ggplot() +
             aes(x = x, y = y, label = label, hjust = 0), 
             size = 4
   ) +
-  coord_flip(xlim = c(3, 55)) +
+  coord_flip(xlim = c(3, 57)) +
   scale_y_reverse(limits = c(4, -1.5)) + 
   theme_dendro()
 dendrogram_plot
@@ -253,7 +253,6 @@ saveRDS(is_non_rank_plot, here::here("figures", "raw", "strategy-info.rds"))
 saveRDS(package_plot, here::here("figures", "raw", "package-info.rds"))
 ```
 
-
 ## Dimensionality reduction: logistic PCA
 
 Perform Logistic PCA
@@ -566,3 +565,107 @@ saveRDS(
   here::here("figures", "raw", "logfc-abs-prop-intersect.rds")
 )
 ```
+
+```{r intersection-dendrogram}
+plot_dendrogram_by_dataset <- function(data, data_id) {
+  
+  n <- 10
+  all_concordance_data <- data %>% 
+    filter(data_id == !!data_id) %>% 
+    rowwise() %>% 
+    mutate(mgs = list(get_top_sel_mgs(mgs, n = n)$gene)) %>% 
+    ungroup() %>% 
+    select(pars, method, cluster, data_id, mgs)
+  
+  intersection_data <- all_concordance_data %>% 
+    expand_grid(pars_2 = unique(all_concordance_data$pars)) %>% 
+    left_join(
+      dplyr::rename(all_concordance_data, mgs_2 = mgs), 
+      by = c("pars_2" = "pars", "cluster")
+    ) %>% 
+    rowwise() %>% 
+    mutate(prop_intersect = length(intersect(mgs, mgs_2)) / n) %>% 
+    ungroup() %>% 
+    select(pars, pars_2, cluster, prop_intersect)  %>%
+    group_by(pars, pars_2) %>% 
+    summarise(prop_intersect = mean(prop_intersect), .groups = "drop") %>% 
+    pivot_wider(names_from = pars, values_from = prop_intersect)
+  
+  intersection_mat <- as.matrix(intersection_data[, -1])
+  rownames(intersection_mat) <- intersection_data[[1]]
+  
+  hier_clust <- hclust(dist(intersection_mat))
+  hier_clust$labels <- pars_lookup[hier_clust$labels]
+  
+  hier_clust_data <- dendro_data(hier_clust, type = "rectangle")
+  
+  dendrogram_plot <- ggplot() +
+    geom_segment(data = segment(hier_clust_data), 
+                 aes(x = x, y = y, xend = xend, yend = yend)
+    ) +
+    geom_text(data = label(hier_clust_data), 
+              aes(x = x, y = y, label = label, hjust = 0), 
+              size = 4
+    ) +
+    coord_flip(xlim = c(3, 57)) +
+    scale_y_reverse(limits = c(4, -1.5)) + 
+    theme_dendro()
+  
+  dendrogram_plot
+}
+
+plot_dendrogram_by_dataset(concordance_data, "pbmc3k")
+plot_dendrogram_by_dataset(concordance_data, "lawlor")
+plot_dendrogram_by_dataset(concordance_data, "endothelial")
+```
+
+```{r}
+n_cells_data <- lapply(
+  list.files(here::here("data", "real_data"), full.names = TRUE), 
+  function(x) {
+   data <- readRDS(x)
+   out <- tibble(
+     data_id = tools::file_path_sans_ext(basename(x)),
+     cell_type = unique(data$label),
+     n = as.vector(table(data$label))
+   )
+   rm(data)
+   out
+  }
+) %>% 
+  bind_rows(!!!.)
+```
+
+```{r}
+n <- 10
+all_concordance_data <- concordance_data %>% 
+  rowwise() %>% 
+  mutate(mgs = list(get_top_sel_mgs(mgs, n = n)$gene)) %>% 
+  ungroup() %>% 
+  select(pars, method, cluster, data_id, mgs)
+
+intersect_n_cells_data <- all_concordance_data %>% 
+  expand_grid(pars_2 = unique(all_concordance_data$pars)) %>% 
+  left_join(
+    dplyr::rename(all_concordance_data, mgs_2 = mgs), 
+    by = c("pars_2" = "pars", "cluster")
+  ) %>% 
+  rowwise() %>% 
+  mutate(prop_intersect = length(intersect(mgs, mgs_2)) / n) %>% 
+  ungroup() %>% 
+  select(pars, pars_2, cluster, prop_intersect, data_id.x)  %>%
+  group_by(cluster, data_id.x) %>% 
+  summarise(prop_intersect = mean(prop_intersect), .groups = "drop") %>% 
+  left_join(n_cells_data, by = c("cluster" = "cell_type", "data_id.x" = "data_id"))
+
+intersect_n_cells_data %>% 
+  filter(!is.na(n)) %>% 
+  ggplot(aes(n, prop_intersect)) + 
+  geom_point() + 
+  labs(
+    x = "Number of cells", 
+    y = "Intersect proportion"
+  ) + 
+  theme_bw()
+```
+