Updated analysis for Figure 1 and 2.

analysis/deprecated/molkart.mesmer_seg.misty_analysis.Rmd

@@ -0,0 +1,229 @@
+---
+title: "molkart.misty_analysis2"
+author: "FloWuenne"
+date: "2024-01-14"
+output: workflowr::wflow_html
+editor_options:
+  chunk_output_type: console
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+
+library(tidyverse)
+library(mistyR)
+library(rlist)
+library(FNN)
+library(future)
+library(cowplot)
+library(igraph)
+library(ClusterR)
+library(viridis)
+
+source("./code/functions.R")
+```
+
+
+This markdown will perform Misty analysis on the cell typing based on Mesmer segmentation masks.
+
+```{r}
+plan(multisession, workers = 8)
+
+## Custom functions
+interaction_communities_info <- function(misty.results, concat.views, view,
+                                         trim = 0, trim.measure = "gain.R2", 
+                                         cutoff = 1, resolution = 1) {
+
+  inv <- sign((stringr::str_detect(trim.measure, "gain") |
+                 stringr::str_detect(trim.measure, "RMSE", negate = TRUE)) - 0.5)
+
+  targets <- misty.results$improvements.stats %>%
+    dplyr::filter(
+      measure == trim.measure,
+      inv * mean >= inv * trim
+    ) %>%
+    dplyr::pull(target)
+
+  view.wide <- misty.results$importances.aggregated %>%
+    filter(view == !!view) %>%
+    pivot_wider(
+      names_from = "Target", values_from = "Importance",
+      id_cols = -c(view, nsamples)
+    ) %>% mutate(across(-c(Predictor,all_of(targets)), \(x) x = NA))
+
+
+  mistarget <- setdiff(view.wide$Predictor, colnames(view.wide)[-1])
+  mispred <- setdiff(colnames(view.wide)[-1], view.wide$Predictor)
+
+  if(length(mispred) != 0){
+    view.wide.aug <- view.wide %>% add_row(Predictor = mispred)
+  } else {
+    view.wide.aug <- view.wide
+  }
+
+  if(length(mistarget) != 0){
+    view.wide.aug <- view.wide.aug %>% 
+      bind_cols(mistarget %>% 
+                  map_dfc(~tibble(!!.x := NA)))
+  }
+
+  A <- view.wide.aug %>%
+    column_to_rownames("Predictor") %>%
+    as.matrix()
+  A[A < cutoff | is.na(A)] <- 0
+
+  ## !!! Was buggy
+  G <- graph.adjacency(A[,rownames(A)], mode = "plus", weighted = TRUE) %>%
+    set.vertex.attribute("name", value = names(V(.))) %>%
+    delete.vertices(which(degree(.) == 0))
+
+  Gdir <- graph.adjacency(A[,rownames(A)], "directed", weighted = TRUE) %>%
+    set.vertex.attribute("name", value = names(V(.))) %>%
+    delete.vertices(which(degree(.) == 0))
+
+  C <- cluster_leiden(G, resolution_parameter = resolution, n_iterations = -1)
+
+  mem <- membership(C)
+
+  Gdir <- set_vertex_attr(Gdir, "community", names(mem), as.numeric(mem))
+
+  # careful here the first argument is the predictor and the second the target, 
+  # it might need to come from different view
+  corrs <- as_edgelist(Gdir) %>% apply(1, \(x) cor(
+    concat.views[[view]][, x[1]],
+    concat.views[["intraview"]][, x[2]]
+  )) %>% replace_na(0)
+
+  Gdir <- set_edge_attr(Gdir, "cor", value = corrs)
+  return(Gdir)
+}
+
+cellular_neighborhoods <- function(sample.cells, sample.pos, n, k){
+  misty.views  <- create_initial_view(sample.cells) %>% add_paraview(sample.pos, family = "constant", l = n)
+  clust <- KMeans_rcpp(misty.views[[paste0("paraview.",n)]], k)
+  return(clust$clusters) 
+}
+```
+
+# Introduction
+
+In this markdown we will utilize [MistyR](https://saezlab.github.io/mistyR/) to perform global spatial analysis on the cell-type encodings for our molecular Cartography data.
+
+Make sure to have the latest development version (15.01.2024) : https://github.com/jtanevski/mistyR
+
+## Analyze data using MistyR with low levels cell phenotypes
+
+```{r}
+size_param <- 125
+all.data <- read_tsv("./output/molkart/molkart.misty_celltype_table.mesmer_seg.lowres.tsv")
+
+samples <- all.data %>%
+  pull(sample_ID) %>%
+  unique()
+
+cts <- all.data %>%
+  pull(misty_cts) %>%
+  unique()
+
+cts.names <- make.names(cts, allow_ = FALSE)
+
+## Count number of cells per type
+ct_numbers <- all.data %>%
+  group_by(sample_ID, misty_cts) %>%
+  summarise(n = n()) %>%
+  pivot_wider(names_from = misty_cts, values_from = n) %>%
+  column_to_rownames("sample_ID") %>%
+  as.matrix()
+
+samples %>% walk(\(sample){
+
+  sample.cells <- all.data %>%
+    filter(sample_ID == sample) %>%
+    pull(misty_cts) %>%
+    map(~ .x == cts) %>%
+    list.rbind() %>%
+    `colnames<-`(cts.names) %>%
+    as_tibble()
+
+  sample.pos <- all.data %>%
+    filter(sample_ID == sample) %>%
+    select(X_centroid, Y_centroid)
+
+  l <- size_param / 0.138
+
+  misty.views.cts <- create_initial_view(sample.cells) %>%
+    add_paraview(sample.pos, l) %>%
+    rename_view(paste0("paraview.", l), "paraview") %>%
+    select_markers("intraview", where(~ sd(.) != 0))
+
+  db_name <- paste("results_cts.mesmer_seg.lowres.",size_param,".sqm",sep="")
+
+  run_misty(misty.views.cts, sample, db_name, bypass.intra = TRUE)
+})
+```
+
+
+```{r}
+l <- size_param / 0.138
+db_name <- paste("results_cts.mesmer_seg.lowres.",size_param,".sqm",sep="")
+groups <- samples %>% str_extract("(?<=sample_).+(?=_r)") %>% unique()
+
+misty.results.g <- groups %>% map(~ collect_results(db_name, .x))
+#misty.results.g <- groups %>% map(~ collect_results(paste("results_cts_",as.character(size_param),".sqm",sep=""), .x,)) ##
+names(misty.results.g) <- groups
+
+outdir <- paste("mesmer_seg.misty_figures_",size_param,sep="")
+dir.create(outdir) ##
+
+misty.results.g %>% iwalk(\(misty.results, cond){
+  plot.list <- list()
+  plot_improvement_stats(misty.results, "gain.R2")
+  plot.list <- list.append(plot.list, last_plot())
+  plot_interaction_heatmap(misty.results, "paraview", cutoff = 0.3, clean = TRUE, trim = 5)
+  plot.list <- list.append(plot.list, last_plot())
+  plot_grid(plotlist = plot.list, ncol = 2)
+  ggsave(paste0(outdir,"/", cond, "_stats.pdf"), width = 10, height = 10)
+  #ggsave(paste0("misty_figures_",size_param,"/", cond, "_stats.pdf"), width = 10, height = 10)
+})
+
+dir.create("graphs", recursive = TRUE)
+
+# names(misty.results.g) %>% walk(\(g){
+#   concat.views <- samples %>% str_subset(g) %>% map(\(sample){
+#     sample.cells <- all.data %>%
+#       filter(sample_ID == sample) %>%
+#       pull(misty_cts) %>%
+#       map(~ .x == cts) %>%
+#       list.rbind() %>%
+#       `colnames<-`(cts.names) %>%
+#       as_tibble()
+# 
+#     sample.pos <- all.data %>%
+#       filter(sample_ID == sample) %>%
+#       select(X_centroid, Y_centroid)
+# 
+#     misty.views.cts <- create_initial_view(sample.cells) %>%
+#       add_paraview(sample.pos, l) %>%
+#       rename_view(paste0("paraview.", l), "paraview") %>%
+#       select_markers("intraview", where(~ sd(.) != 0))
+# 
+#   }) %>% list.clean() %>% reduce(\(acc, nxt){
+#     list(
+#       intraview = bind_rows(acc[["intraview"]], nxt[["intraview"]]),
+#       paraview = bind_rows(acc[["paraview"]], nxt[["paraview"]])
+#     )
+#   }, .init = list(intraview = NULL, paraview = NULL))
+# 
+#   out <- interaction_communities_info(misty.results.g[[g]], concat.views,
+#                                       view = "paraview", trim = 5, cutoff = 0.5)
+#   graph_name <- paste0("graphs/para.",size_param,".", g, ".graphml")
+#   write_graph(out, graph_name , "graphml")
+#   #write_graph(out, paste0("graphs/para.",size_param,".", g, ".graphml"), "graphml")
+#  })
+```
+
+```{r}
+# ## Save misty results in R object for easier faster loading
+# saveRDS(misty.results.g,
+#         file = paste0("./output/molkart/misty_results.lowres.",size_param,".rds"))
+```