2V_02_Fig2_gsea_deconvolution.Rmd

---
title: "2V_02_Fig2_gsea_deconvolution"
author: "DanielZucha"
date: "2024-05-11"
output: html_document
---


Hello, 

In this markdown we prepare the visual representations for the figure 2, summarizing the biological processes and cell types present in the lesioned areas. To develop a general understanding of altered processed, we submitted the regional markers genes for gene enrichment analysis using the online [MetaScape tool](https://metascape.org/gp/index.html#/main/step1). The enrichment was for performed for each lesional area (area + time) separately. The MetaScape allows a network for visualisation purposes, drawing the general scheme of events. To assess cell type distribution, we used RCTD algorithm ([spacexr package](https://github.com/dmcable/spacexr)) developed by Cable et al. Being a reference-based deconvolution algorithm, we reanalyzed the publicly available [Zeng et al 2023](http://www.ncbi.nlm.nih.gov/pubmed/37063847) single-cell stroke [dataset](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE227651) as it matched our experimental design and cell types of all plausible classes were included, i.e. neurons, glial cells, vascular or peripheral immune cells.


```{r libraries}
# daily libraries
library(tidyverse)
library(dplyr)
library(Seurat)
library(ggplot2)
library(openxlsx)
library(magrittr)
library(stringr)
library(patchwork)

# graphics
library(circlize)
library(grid)
library(magick)
library(ggrepel)

# sourcing custom functions and color palettes
source("code/supporting_functions_MCAO.R")
```

```{r working space and lists}
# ws
if(!dir.exists("results/PNAS/Fig02")){dir.create("results/PNAS/Fig02")
  ws <- "results/PNAS/Fig02"}else{ws <- "results/PNAS/Fig02"}

# lists
plot.list <- list()
result.list <- list()
```

load data
```{r load data}
# spatial seurat processed with a standardized seurat pipeline
spatial.seurat <- readRDS(file = file.path("data", "seurat_spatial_integrated_ready_rctdZeng2023.rds"))

# misc 
sections <- c("Ctrl", "1DPI", "3DPI", "7DPI")
ctx_regions <- c("CTX1-4", "CTX5", "CTX6", "ISD1c", "ISD1p", "ISD3c", "ISD3p", "ISD7c", "ISD7p")
```

# Functional Annotation
load the summarizing GO parent terms
```{r GO parents}
result.list[["Enrichments"]] <- readRDS(file = file.path("ws", "Metascape_GO_parent_terms.Rds"))
result.list[["Enrichments"]] %>% glimpse
```

```{r GO summaries as module scores}
# get the gene lists of GO summary parent terms
result.list[["GOsummaries"]] <- list()
for (i in result.list[["Enrichments"]][["Process_Description"]]) {
  # select the gene symbols of every process
  temp <-  result.list[["Enrichments"]] %>%
    filter(Process_Description %in% i) %>%
    pull(Symbols) %>%
    strsplit(",")
  
  # append the results
  result.list[["GOsummaries"]] %<>% append(temp)
}
names(result.list[["GOsummaries"]]) <-
  result.list[["Enrichments"]] %>% pull(Process_Description)

# turn the GO summaries into module scores.
spatial.seurat %<>% AddModuleScore(features = result.list[["GOsummaries"]],
                                   name = names(result.list[["GOsummaries"]]),
                                   assay = "SCT")

# remove environmental variables
rm(i, temp)
```

Visualize selected processes
```{r SpatialPlot of the GO summaries}
features <-
  spatial.seurat@meta.data %>% 
  colnames() %>% 
  stringr::str_subset(pattern = "^(1DPI_|3DPI_|7DPI_)")
  
plot.list[["SpatialPlots_GOsummaries"]] <-
  features %>% lapply(\(x) {
    # range value for the feature
    limits <- c(min(spatial.seurat@meta.data[[x]]),
                max(spatial.seurat@meta.data[[x]]))
    
    # # unify the color scheme such that middle value is 0
    # color_values <- c(
    #   limits[1],
    #   limits[1] + (limits[2] - limits[1]) * 1 / 6,
    #   limits[1] + (limits[2] - limits[1]) * 2 / 6,
    #   0,
    #   # Midpoint corresponding to 0
    #   limits[1] + (limits[2] - limits[1]) * 4 / 6,
    #   limits[1] + (limits[2] - limits[1]) * 5 / 6,
    #   limits[2]
    # )
    
    #Plot the spatial data
    plot <- spatial.seurat %>% SpatialPlot(
      features = x,
      images = sections,
      crop = TRUE,
      ncol = 4,
      pt.size.factor = 2.5,
      stroke = 0,
      alpha = 1,
      image.alpha = 0,
      combine = TRUE
    ) &
      theme_mk &
      remove_grid &
      NoAxes() &
      theme(legend.position = "none",
            panel.border = element_blank()) &
      scale_fill_gradientn(
        colours = col.list$gradient_blue_yellow_red,
        limits = limits
      ) &
      scale_alpha_continuous(range = c(0.5, 1)) &
      labs(title = NULL, subtitle = x)
    
    return(plot)
  })

names(plot.list[["SpatialPlots_GOsummaries"]]) <- features
print(plot.list[["SpatialPlots_GOsummaries"]][["7DPI_Gliogenesis53"]])

png_save_show(
  plot = plot.list[["SpatialPlots_GOsummaries"]][["7DPI_Gliogenesis53"]], 
  file = file.path(ws, "Spatial_GO_Gliogenesis.png"), 
  dpi = 1000,
  height = 70, 
  width = 100
)

```

The GO term networks were prepared natively with MetaScape web interface.

# Deconvolution

First, we start by looking at the neuronal content in the cortex. The lesions have significantly less neuronal content.
```{r neuronal content in ctx}
# neuronal loss
{
  temp <- spatial.seurat@meta.data %>% filter(
    BrainAreas %in% c("Cortex", "Lesion")
  ) %>% select(Condition, BrainAreas, RCTD_Neurons)

  t.test(RCTD_Neurons ~ BrainAreas, data = temp)
  
  plot.list[["DensityPlot_NeuronalContent"]] <- ggplot(temp, aes(.data[["RCTD_Neurons"]], fill = .data[["BrainAreas"]])) + 
    geom_density(alpha = 0.8, color = "white", linewidth = 0.5) +
    scale_x_reverse() +
    annotate(geom = "text", x = 0.2, y = 8.4, label = expression("P < 2 x 10"^-16), hjust = 0, vjust = 0.5, size = 2.45, lineheight = 0.7, color = "black") +
    annotate(geom = "text", x = 0.375, y = 9, label = "Lesioned cortex", hjust = 0.5, vjust = 0.5, size = 2.45, lineheight = 0.7, color = "black") +
    annotate(geom = "text", x = 0.875, y = 9, label = "Intact cortex", hjust = 0.5, vjust = 0.5, size = 2.45, lineheight = 0.7, color = "black") +
    geom_vline(xintercept = 0.75, color = "black", linetype = "dotted", linewidth = 0.5) +
    scale_fill_manual(values = col.list$BrainAreas) + 
    theme_mk +
    remove_grid +
    theme(legend.position = c(0.75, 0.75), 
          legend.justification = c(0, 0.5)) +
    labs(x = "Per-spot neuronal proportion", y = "Density", fill = NULL, title = "Cortical Neuron Content")
  print(plot.list[["DensityPlot_NeuronalContent"]])
  png_save_show(plot = plot.list[["DensityPlot_NeuronalContent"]], 
              file = file.path(ws, "DensityPlot_NeuronalContent.png"), 
              dpi = 1000, 
              height = 50, 
              width = 120)
}
```

Glial content
```{r glial content}
# glial rise
{
  temp <- spatial.seurat@meta.data %>% 
    filter(BrainAreas %in% c("Cortex", "Lesion")) %>% 
    select(Condition, BrainAreas, RCTD_Microglia, RCTD_Astrocytes, RCTD_OLs, RCTD_OPCs) %>% 
    mutate(Glia = RCTD_Microglia + RCTD_Astrocytes + RCTD_OLs + RCTD_OPCs)
  
  t.test(Glia ~ BrainAreas, data = temp)
  
  plot.list[["DensityPlot_GlialContent"]] <- ggplot(temp, aes(.data[["Glia"]], fill = .data[["BrainAreas"]])) + 
    geom_density(alpha = 0.8, color = "white", linewidth = 0.5) +
    annotate(geom = "text", x = 0.75, y = 6, label = expression("P < 2 x 10"^-16), hjust = 0, vjust = 0.5, size = 2.45, lineheight = 0.7, color = "black") +
    scale_fill_manual(values = col.list$BrainAreas) + 
    theme_mk +
    remove_grid +
    theme(legend.position = c(0.75, 0.75), 
          legend.justification = c(0, 0.5)) +
    labs(x = "Per-spot glial proportion", y = "Density", fill = NULL, title = "Cortical Glial Content")
  print(plot.list[["DensityPlot_GlialContent"]])
  png_save_show(plot = plot.list[["DensityPlot_GlialContent"]], 
              file = file.path(ws, "DensityPlot_GlialContent.png"), 
              dpi = 1000, 
              height = 70, 
              width = 90)
}
```

Summarizing stacked barplot to see the full overview of cell types
```{r stacked barplot of rois}
# celltypes of interest
coi <- spatial.seurat@meta.data %>% colnames %>% stringr::str_subset(pattern = "^RCTD_")

# df of celltypes of interest (coi) in the regions of interest (roi)
spatial.seurat@meta.data[["DetailedRegionAnnoShort"]] %>% levels

# are regions of interest defined?
if(!exists("roi")){
  roi <- c(
  "CTX1-4", "CTX5", "CTX6", # intact cortex
  "ISD1c", "ISD1p", "ISD3c", "ISD3p", "ISD7c", "ISD7p") # lesion
  }

df <- spatial.seurat@meta.data %>% 
  filter(DetailedRegionAnnoShort %in% roi) %>% 
  select(DetailedRegionAnnoShort, coi)

# table of means for coi per roi
mean_table <- df %>% 
  group_by(DetailedRegionAnnoShort) %>% 
  summarize(across(starts_with("RCTD_"), \(x) mean(x, na.rm = TRUE)), .groups = "drop") %>% 
  column_to_rownames(var = "DetailedRegionAnnoShort") %>% 
  "*"(100) %>% 
  round(1) %>% 
  rownames_to_column(var = "DetailedRegionAnnoShort")


# Stacked bar plot
{
  long_data <- mean_table %>% 
    pivot_longer(
      cols = starts_with("RCTD_"), 
      names_to = "CellType", 
      values_to = "MeanValue"
    ) %>% mutate(CellType = stringr::str_remove(CellType, "RCTD_"))
  
  # Create the plot
  main_plot <- ggplot(long_data, aes(x = DetailedRegionAnnoShort, y = MeanValue, fill = CellType)) +
    geom_bar(stat = "identity", alpha = 0.9, width = 0.75) + # Stacked bar plot
    geom_text(aes(label = ifelse(MeanValue > 10, paste0(round(MeanValue, 1), "%"), "")), # display when >10%
              position = position_stack(vjust = 0.5), 
              size = 2.45, 
              color = "black") +
    scale_fill_manual(values = col.list$zeng.cols, 
                      guide = guide_legend(nrow = 3)) +
    labs(
      x = NULL,
      y = "Average Spot Composition (%)",
      fill = "Cell Type"
    ) + 
    theme_mk +
    theme(legend.position = "bottom") +
    remove_grid
  
  
  ann_plot <- ggplot(long_data, aes(x = DetailedRegionAnnoShort, y = MeanValue)) +
    geom_bar(stat = "identity", fill = NA) +
    ylim(0, 1) + 
    annotate("text", x = 2, y = 0.5, label = "Intact Cortex", vjust = 0.5, hjust = 0.5, size = 2.45) + 
    annotate("segment", x = 0.5, xend = 3.45, y = 0.1, yend = 0.1, color = "grey", linewidth = 2) + 
    annotate("text", x = 6.5, y = 0.5, label = "Lesion", vjust = 0.5, hjust = 0.5, size = 2.45) + 
    annotate("segment", x = 3.55, xend = 9.5, y = 0.1, yend = 0.1, color = "#AF3039", linewidth = 2) + 
    theme_void() & 
      NoAxes()
  
  plot.list[["StackedBarPlot_RCTD_ROIcomposition"]] <- (ann_plot / main_plot) + plot_layout(heights = c(1,10))

  # Display the plot
  print(plot.list[["StackedBarPlot_RCTD_ROIcomposition"]])
  png_save_show(plot = plot.list[["StackedBarPlot_RCTD_ROIcomposition"]],
              file = file.path(ws, "StackedBarPlot_RCTD_ROIcomposition.png"),
              dpi = 1000,
              height = 120,
              width = 150)
}
```

Close up on the lesion areas - pie chart
```{r close up on the lesion areas}
library(STdeconvolve)

temp <- lapply(sections, \(section){
  
  # df of spot positions
  pos <- spatial.seurat@images[[section]]@coordinates[, c("imagerow", "imagecol")] %>% 
    rename(
      x = imagecol, 
      y = imagerow
    ) %>% 
    mutate(y = y %>% "*"(-1)) # turn the y coordinates upside-down
  
  # define spots in the area of interest
  if (section == "7DPI") {
    division_factor_x <- 2
    division_factor_y <- 3
  } else if(section == "Ctrl"){
    division_factor_x <- 3
    division_factor_y <- 3
  } else {
    division_factor_x <- 4
    division_factor_y <- 4
  }
  
  x.threshold <-
    (max(pos$x) - min(pos$x)) %>% "/"(division_factor_x) %>% "+"(min(pos$x)) %>% round(0)
  y.threshold <-
    (max(pos$y) - min(pos$y)) %>% "/"(division_factor_y) %>% "+"(min(pos$y)) %>% round(0)
  
  pos %<>% filter(x > x.threshold,
                  y > y.threshold)
  
  # for 7DPI section, we also trip top spots
  if(section == "7DPI"){
    top_pos_to_trim <- 
      (max(pos$y) - min(pos$y)) %>% "/"(5)
    top_y_threshold <- 
      max(pos$y) %>% "-"(top_pos_to_trim) %>% round(0)
    
    pos %<>% filter(y < top_y_threshold)
  }
  
  # names of spots to keep
  spots.to.plot <- rownames(pos)
  
  # keep metadata for the area of interest
  md <- spatial.seurat@meta.data %>%
    rownames_to_column("barcodes") %>%
    filter(barcodes %in% spots.to.plot) %>%
    column_to_rownames("barcodes")
  
  # prep data to plot
  theta <- md %>% select(starts_with("RCTD_")) %>% dplyr::rename_with(~ sub("^RCTD_", "", .))
  groups <- md %>% pull(DetailedRegionAnnoShort)
  group_cols <- col.list$cols_mono_short[match(x = groups, table = names(col.list$cols_mono_short))]
  topic_cols <- col.list$zeng.cols[match(x = colnames(theta), table = names(col.list$zeng.cols))]
  
  # plot the pie charts
  pie_chart <- vizAllTopics(
    theta = theta, 
    pos = pos, 
    groups = groups, 
    group_cols = group_cols, 
    topicCols = topic_cols,
    plotTitle = section, 
    lwd = 0.15, 
    r = 60, 
    showLegend = F
  )
  
  return(pie_chart)
}); names(temp) <- paste0("PieChart_RCTD_", sections)

# patchwork the individual piecharts together
plot.list[["PieChart_Lesion_CloseUp"]] <- patchwork::wrap_plots(temp, ncol = 4)
print(plot.list[["PieChart_Lesion_CloseUp"]])

rm(temp)

png_save_show(plot = plot.list[["PieChart_Lesion_CloseUp"]], 
              file = file.path(ws, "PieChart_Lesion_CloseUp.png"), 
              dpi = 1000, 
              height = 100, 
              width = 360)
```

Expression of selected chemokines
```{r chemokine expression}
chemokines <- c("Ccl12", "Ccl3", "Ccl4")

spatial.seurat %<>% AddModuleScore(features = list(chemokines), name = "chemokines")

feature <- "chemokines1"
plot.list[["Spatial_Chemokines"]] <- 
  spatial.seurat %>%
  SpatialPlot(
    features = feature,
    images = sections,
    image.alpha = 0.05,
    stroke = 0,
    pt.size.factor = 2.2,
    crop = TRUE,
    ncol = 4
  ) &
  theme_mk &
  remove_grid &
  NoAxes() &
  theme(legend.position = "none",
        panel.border = element_blank()) &
  scale_alpha_continuous(range = c(0.2, 1)) &
  labs(title = NULL, subtitle = NULL) &
  viridis::scale_fill_viridis(option = 'rocket',
                              direction = -1,
                              limits = c(
                                min(spatial.seurat@meta.data[[feature]]),
                                max(spatial.seurat@meta.data[[feature]])
                              ))

print(plot.list[["Spatial_Chemokines"]])

png_save_show(
  plot = plot.list[["Spatial_Chemokines"]], 
  file = file.path(ws, "Spatial_Chemokines.png"), 
  dpi = 1000,
  height = 70,
  width = 100
)

```

Correlation genes vs Cell type proportions
```{r load gene vs cell correlations}
gene_cell_correlations <- read_all_sheets(file = file.path(ws, "Correlations_gexp_CellProportion_in_lesions.xlsx"))
```

```{r Microglia correlated 3DPI}
# define timepoint and celltype to visualize
temp_df <- gene_cell_correlations[["3DPI"]] %>% 
  select(RCTD_Microglia, GeneSymbols) %>% 
  filter(!is.na(RCTD_Microglia)) %>% 
  arrange(desc(RCTD_Microglia)) %>% 
  mutate(Percentile = row_number() / n() * 100)

# top positive and negative correlations
top_pos <- temp_df %>% 
  filter(GeneSymbols %in% c(
    "Ctsd", "Ctsz", "Ctsb",
    "C1qa", "C1qb", "C1qc", 
    "Cst3", 
    "Spp1", "Lpl", "Abhd12",
    "Ly86", "Hexb", "Trem2", "Cd9"))
# top_neg <- temp %>% top_n(-5, RCTD_Microglia)
# top_genes <- bind_rows(top_pos, top_neg)

# Create the rank plot
plot.list[["RankPlot_Genes_vs_Microglia_3DPI"]] <- 
  ggplot(temp_df, aes(x = Percentile, y = RCTD_Microglia, color = RCTD_Microglia)) +
  geom_point() + 
  geom_label_repel(data = top_pos, aes(label = GeneSymbols), 
                   size = 2.45,
                   color = "white", 
                   fill = col.list$zeng.cols["Microglia"], 
                   box.padding = 0.3, 
                   point.padding = 0.2, 
                   segment.color = 'grey50',
                   max.overlaps = Inf,
                   min.segment.length = 0,
                   force = 3, 
                   direction = "both", # Move labels horizontally
                   nudge_x = 25 , # Ensure labels start at least at x = 10
                   segment.size = 0.2) + 
  labs(title = NULL, 
       subtitle = "Microglia-correlated genes\n3DPI", 
       x = "Gene Percentile",
       y = expression(rho[Pearson])) + 
  scale_color_gradient2(low = "grey40", mid = "#FFFFFF", high = "#D01B1B", midpoint = 0) +
  theme_mk + 
  remove_grid + 
  theme(legend.position = "none")

print(plot.list[["RankPlot_Genes_vs_Microglia_3DPI"]])

png_save_show(
  plot = plot.list[["RankPlot_Genes_vs_Microglia_3DPI"]], 
  file = file.path(ws, "RankPlot_correlation_genes_microglia_3DPI.png"), 
  show_plot = FALSE, 
  width = 110, 
  height = 60
)
```

```{r Astrocytes correlated genes 7DPI}
# define timepoint and celltype to visualize
temp_df <- gene_cell_correlations[["7DPI"]] %>% 
  select(RCTD_Astrocytes, GeneSymbols) %>% 
  filter(!is.na(RCTD_Astrocytes)) %>% 
  arrange(desc(RCTD_Astrocytes)) %>% 
  mutate(Percentile = row_number() / n() * 100)

# top positive and negative correlations
top_pos <- temp_df %>% top_n(n = 20, wt = RCTD_Astrocytes)

top_pos <- temp_df %>% 
  filter(GeneSymbols %in% c(
    "Clu", "Aldoc", "Gfap",
    "Mt1", "Mt2", "Mt3", 
    "Aqp4", "Fxyd1", "Prdx6", 
    "Atp1a2", "Vim"
))
# top_neg <- temp %>% top_n(-5, RCTD_Microglia)
# top_genes <- bind_rows(top_pos, top_neg)

# Create the rank plot
plot.list[["RankPlot_Genes_vs_Astrocytes_7DPI"]] <- 
  ggplot(temp_df, aes(x = Percentile, y = RCTD_Astrocytes, color = RCTD_Astrocytes)) +
  geom_point() + 
  geom_label_repel(data = top_pos, aes(label = GeneSymbols), 
                   size = 2.45,
                   color = "white", 
                   fill = col.list$zeng.cols["Astrocytes"], 
                   box.padding = 0.3, 
                   point.padding = 0.2, 
                   segment.color = 'grey50',
                   max.overlaps = Inf,
                   min.segment.length = 0,
                   force = 3, 
                   direction = "both", # Move labels horizontally
                   nudge_x = 25 , # Ensure labels start at least at x = 10
                   segment.size = 0.2) + 
  labs(title = NULL, 
       subtitle = "Astrocytes-correlated genes\n7DPI", 
       x = "Gene Percentile",
       y = expression(rho[Pearson])) + 
  scale_color_gradient2(low = "grey40", mid = "#FFFFFF", high = "#02818A", midpoint = 0) +
  theme_mk + 
  remove_grid + 
  theme(legend.position = "none")

print(plot.list[["RankPlot_Genes_vs_Astrocytes_7DPI"]])

png_save_show(
  plot = plot.list[["RankPlot_Genes_vs_Astrocytes_7DPI"]], 
  file = file.path(ws, "RankPlot_correlation_genes_astrocytes_7DPI.png"), 
  show_plot = FALSE, 
  width = 110, 
  height = 60
)

```

```{r Oligodendrocytes correlated genes 7DPI}
# define timepoint and celltype to visualize
temp_df <- gene_cell_correlations[["7DPI"]] %>% 
  select(RCTD_OLs, GeneSymbols) %>% 
  filter(!is.na(RCTD_OLs)) %>% 
  arrange(desc(RCTD_OLs)) %>% 
  mutate(Percentile = row_number() / n() * 100)

# top positive and negative correlations
top_pos <- temp_df %>% top_n(n = 20, wt = RCTD_OLs)

top_pos <- temp_df %>% 
  filter(GeneSymbols %in% c(
    "Plp1", "Mbp", "Mobp", 
    "Cldn11", "Il33", "Klk6", 
    "Serpina3n", "C4b", "Trf", 
    "Car2", "Cnp", "Tubb4a"
))
# top_neg <- temp %>% top_n(-5, RCTD_Microglia)
# top_genes <- bind_rows(top_pos, top_neg)

# Create the rank plot
plot.list[["RankPlot_Genes_vs_OLs_7DPI"]] <- 
  ggplot(temp_df, aes(x = Percentile, y = RCTD_OLs, color = RCTD_OLs)) +
  geom_point() + 
  geom_label_repel(data = top_pos, aes(label = GeneSymbols), 
                   size = 2.45,
                   color = "white", 
                   fill = col.list$zeng.cols["OLs"], 
                   box.padding = 0.3, 
                   point.padding = 0.2, 
                   segment.color = 'grey50',
                   max.overlaps = Inf,
                   min.segment.length = 0,
                   force = 3, 
                   direction = "both", # Move labels horizontally
                   nudge_x = 25 , # Ensure labels start at least at x = 10
                   segment.size = 0.2) + 
  labs(title = NULL, 
       subtitle = "Oligodendrocytes-correlated genes\n7DPI", 
       x = "Gene Percentile",
       y = expression(rho[Pearson])) + 
  scale_color_gradient2(low = "grey40", mid = "#FFFFFF", high = "#8C96C6", midpoint = 0) +
  theme_mk + 
  remove_grid + 
  theme(legend.position = "none")

print(plot.list[["RankPlot_Genes_vs_OLs_7DPI"]])

png_save_show(
  plot = plot.list[["RankPlot_Genes_vs_OLs_7DPI"]], 
  file = file.path(ws, "RankPlot_correlation_genes_oligodendrocytes_7DPI.png"), 
  show_plot = FALSE, 
  width = 110, 
  height = 60
)

```

Individual cell type populations in cortex vs lesion
```{r vlnplot for individual coi}
# choose celltypes of interest to plot
coi_to_plot <- c(
    "RCTD_Neurons",
    "RCTD_Astrocytes",
    "RCTD_Microglia",
    "RCTD_OLs",
    "RCTD_OPCs",
    "RCTD_PeripheralMyeloidCells"
  )

# regions of interest
roi <- c(
  "CTX1-4", "CTX5", "CTX6", # intact cortex
  "ISD1c", "ISD1p", "ISD3c", "ISD3p", "ISD7c", "ISD7p" # lesion
  )

plot.list[["VlnPlot_COI_in_ROI"]] <- spatial.seurat %>% enh_vlnplot(
  feature = coi_to_plot, 
  grouping = "BrainAreas", 
  colors = col.list$BrainAreas, 
  compare_means = T,
  stat_test = "t.test", 
  ref.group = "Cortex",
  idents = roi, 
  ncol = 6, 
  combine = T
) & 
  theme_mk & 
  remove_grid & 
  xlab(NULL) & 
  NoLegend()

print(plot.list[["VlnPlot_COI_in_ROI"]])


# change the plot names and keep just one legend
plot.list[["VlnPlot_COI_in_ROI"]] <- lapply(seq_along(plot.list[["VlnPlot_COI_in_ROI"]]), \(plot){
  
  # define the plot title
  title <- coi_to_plot[plot] %>% stringr::str_remove(pattern = "^RCTD_")
  if(title == "PeripheralMyeloidCells"){title <- "Peripheral Myeloid Cells"}
  
  # define the plot output
  ## give the first plot a y axis label and new title
  if(plot == 1){
    plot.list[["VlnPlot_COI_in_ROI"]][[plot]] <- plot.list[["VlnPlot_COI_in_ROI"]][[plot]] +
      ggtitle(title) + 
      ylab("Per-spot proportion")
  } else {
    # for the rest just fix the title
   plot.list[["VlnPlot_COI_in_ROI"]][[plot]] <- plot.list[["VlnPlot_COI_in_ROI"]][[plot]] +
     ggtitle(title)
  }
})

plot.list[["VlnPlot_COI_in_ROI"]] %<>% wrap_plots(ncol = length(.))
  
print(plot.list[["VlnPlot_COI_in_ROI"]])
```

Each cell population spatially
```{r SpatialPlot of RCTD}
temp <- lapply("RCTD_Microglia", ## RCTD results
               \(celltype_level) {
                 celltype <- celltype_level %>% str_remove(pattern = "^RCTD_")
                 plot <- spatial.seurat %>% SpatialFeaturePlot(
                   features = celltype_level,
                   images = sections,
                   crop = T,
                   ncol = 4,
                   pt.size.factor = 2.5,
                   stroke = 0,
                   alpha = 1,
                   image.alpha = 0.05,
                   combine = T
                 ) &
                   theme_mk & 
                   remove_grid &
                   NoAxes() & 
                   NoLegend() &
                   theme(legend.position = "none", 
                         panel.border = element_blank()) & 
                   scale_alpha_continuous(range = c(0.2, 1)) &
                   viridis::scale_fill_viridis(option = 'rocket',
                                               direction = -1,
                                               limits = c(0, 1), 
                                               breaks = c(0, 0.5, 1)) &
                   labs(title = NULL)
                 
                 plot <-
                   plot + plot_annotation(title = celltype)
               })
names(temp) <- paste0("SpatialPlot_", coi %>% str_remove(pattern = "^RCTD_"))
plot.list %<>% append(values = temp)

print(plot.list$SpatialPlot_Microglia)

# Display the plot
for(plot in names(plot.list)){
  png_save_show(
    plot = plot.list[[plot]],
    file = file.path(ws, paste0(plot, ".png")),
    dpi = 1000,
    height = 50,
    width = 150
  )
  dev.off()
}
```

# wrap plot

```{r }
plot.list[["wrap_gliogenesis_chemokines_microglia"]] <- 
  patchwork::wrap_plots(
    list(
      plot.list[["SpatialPlots_GOsummaries"]][["7DPI_Gliogenesis53"]], 
      temp[["SpatialPlot_Microglia"]], 
      plot.list[["Spatial_Chemokines"]]),
    ncol = 1
  )
print(plot.list[["wrap_gliogenesis_chemokines_microglia"]])

png_save_show(
  plot = plot.list[["wrap_gliogenesis_chemokines_microglia"]],
  file = file.path(ws, "wrap_gliogenesis_chemokines_microglia.png"),
  dpi = 1000,
  height = 150,
  width = 200
)



```