DSP_DE_Analysis_Report.qmd

---
title: "DSP DE Analysis Template"
format:
  html:
    code-fold: true
editor: visual

params:
  data.folder: "$path_to_folder"
  QC.folder: "$path_to_folder"
  results.folder: "$path_to_folder"
  run.folder: "$path_to_folder"
---

# DE Analysis - Vessel

```{r Setup}
#| warning: false
#| message: false

library(GeomxTools)
library(dplyr)
library(limma)
library(edgeR)
library(ggplot2)
library(ggrepel)
library(stringr)
library(PCAtools)
library(readxl)
library(gridExtra)
library(grid)
library(DSPWorkflow)

# Source the helper functions
source("DSP_functions.R")

results.folder <- "$path_to_folder"
data.folder <- "$path_to_folder"
QC.folder <- "$path_to_folder"
run.folder <- "$path_to_folder"

DE.results.folder <- paste0(results.folder, run.folder)
DE.results.files <- list.files(DE.results.folder)
```

### Load the Normalized Counts and Annotation

```{r Load Normalized Counts}

normalized.object.list <- get(load(file = paste0(params$QC.folder,
                                       "$normalized.object.RDA")))

```

### Set up the Annotation Colors

```{r Annotation Colors}

# Set up annotation colors for the heatmaps
anno.colors = list(
      Treatment_group = c(C = "slateblue", 
                T1 = "coral", 
                T2 = "seagreen1", 
                T3 = "cadetblue1", 
                T4 = "pink")
      )
```

### Run DE

#### T1 vs. C

```{r DE contrast: T1 vs. C}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
# Order is Condition, Reference
contrast.groups <- c("T1", "C")

contrast.field <- "Treatment_group"

# Define the name of the contrast
contrast.name <- paste0(contrast.groups[[2]], 
                        "_", 
                        contrast.groups[[1]])

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list()

contrast.groups.list[[contrast.field]] = contrast.groups
 


# Filter data for the chosen annotations
object <- normalized.object.list$object

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_object_for_lmm(object = object, 
                            subset.list = contrast.groups.list)

# Create summary table of group counts
summary.table.df <- pData(lmm.input$subset.object) %>%  
  select(c(names(contrast.groups.list)))

summary.table <- table(summary.table.df)

print("Sample Numbers per Annotation Group")
print("-----------------------------------")
print(summary.table)

# Check if the DE results have already been generated
DE.result.file <- grepl(contrast.name, DE.results.files)
if(sum(DE.result.file) == 1) {
  
  # Load the previously generated DE results
  results.df <- as.data.frame(read.csv(paste0(results.folder, run.folder, contrast.name, "_de.results.csv")))
                              
} else {
  # Use the function from the DSPWorkflow package

  # Within slide analysis
  # Listed contrasts are condition, reference 
  results.list <- diffExpr(object = lmm.input$subset.object, 
                        analysis.type = "Within Groups", 
                        region.col = "Treatment_group", 
                        regions = c("T1", "C"), 
                        group.col = "region", 
                        groups = c("region_B"), 
                        n.cores = parallel::detectCores())
  
  # Create the results df
  results.df <- results.list$results
  
  write.results <- FALSE
  if(write.results == TRUE){
    write.csv(results.df, paste0(results.folder, run.folder, contrast.name, "_de.results.csv"))
  }
  
  
}

# Adjust column names
logfc.column <- colnames(results.df[grepl("logFC",colnames(results.df))])
results.df$logfc <- results.df[[logfc.column]]

pval.column <- colnames(results.df[grepl("_pval",colnames(results.df))])
results.df$pval <- results.df[[pval.column]]

adj.pval.column <- colnames(results.df[grepl("adjpval",colnames(results.df))])
results.df$padj <- results.df[[adj.pval.column]]

results.df$gene <- results.df$Gene

# Keep only the necessary columns
results.df <- results.df %>% select(c("gene", 
                                      "logfc", 
                                      "pval", 
                                      "padj"))


# Export the results
export.deg.list <- FALSE
if(export.deg.list == TRUE){ 
  write.csv(results.df, 
            file = paste0(results.folder, run.folder, "DE/DEG_lists/", contrast.name, "_deg_list.csv"), 
            row.names = FALSE)  
}

# Make the volcano plot
volcano.output <- make_volcano(lmm.results = results.df, 
                               title = contrast.name, 
                               legend.title = "Expression in T1", 
                               x.axis.title = "logFC: C (cond) versus T1 (ref)", 
                               fc.limit = 1, 
                               pos.label.limit = 1, 
                               neg.label.limit = -1)

# Export the volcano
export.volcano <- FALSE
if(export.volcano == TRUE){ 
  ggsave(filename = paste0(results.folder, run.folder, "DE/volcano/", contrast.name, "_volcano_plot.png"), 
         width = 14, 
         height = 10)
}

print(volcano.output$volcano.plot)

# Subset the annotation just for the heatmap annotations of interest
annotation.heatmap <- lmm.input$annotation %>% 
  select(names(anno.colors)) %>% 
  arrange(Treatment_group)

heatmap.plot <- make_heatmap(
  normalized.log.counts.df = lmm.input$log.counts, 
  de.results = results.df, 
  top.degs = TRUE, 
  annotation.column = annotation.heatmap,
  annotation.row = NULL, 
  anno.colors = anno.colors, 
  cluster.rows = TRUE, 
  cluster.columns = TRUE, 
  main.title = contrast.name, 
  row.gaps = NULL, 
  column.gaps = NULL, 
  show.colnames = FALSE, 
  show.rownames = TRUE)

print(heatmap.plot)

export.heatmap <- FALSE
if(export.heatmap == TRUE){
  
  ggsave(heatmap.plot, 
         filename = paste0(results.folder, run.folder, "DE/heatmap/", contrast.name, "_heatmap_plot.png"), 
         width = 14, 
         height = 10)
  
}

# GSEA Preranked List

gsea.preranked.df <- gsea_preranked_list(contrast.field = contrast.field, 
                                         contrast.levels = contrast.groups, 
                                         annotation = lmm.input$annotation,
                                         log.counts = lmm.input$log.counts)

export.gsea.input <- FALSE
if(export.gsea.input == TRUE){
  
  write.csv(GSEA.preanked.df, file = paste0(results.folder, run.folder, "GSEA/", contrast.name, "_gsea_preranked_input.csv"), row.names = FALSE)
  
}

```

#### T4 vs. C

```{r DE contrast: T4 vs. C}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
# Order is Condition, Reference
contrast.groups <- c("T4", "C")

contrast.field <- "Treatment_group"

# Define the name of the contrast
contrast.name <- paste0(contrast.groups[[2]], 
                        "_", 
                        contrast.groups[[1]])

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list()

contrast.groups.list[[contrast.field]] = contrast.groups
 


# Filter data for the chosen annotations
object <- normalized.object.list$object

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_object_for_lmm(object = object, 
                            subset.list = contrast.groups.list)

# Create summary table of group counts
summary.table.df <- pData(lmm.input$subset.object) %>%  
  select(c(names(contrast.groups.list)))

summary.table <- table(summary.table.df)

print("Sample Numbers per Annotation Group")
print("-----------------------------------")
print(summary.table)

# Check if the DE results have already been generated
DE.result.file <- grepl(contrast.name, DE.results.files)
if(sum(DE.result.file) == 1) {
  
  # Load the previously generated DE results
  results.df <- as.data.frame(read.csv(paste0(results.folder, run.folder, contrast.name, "_de.results.csv")))
                              
} else {
  # Use the function from the DSPWorkflow package

  # Within slide analysis
  # Listed contrasts are condition, reference 
  results.list <- diffExpr(object = lmm.input$subset.object, 
                        analysis.type = "Within Groups", 
                        region.col = "Treatment_group", 
                        regions = c("T4", "C"), 
                        group.col = "region", 
                        groups = c("region_B"), 
                        n.cores = parallel::detectCores())
  
  # Create the results df
  results.df <- results.list$results
  
  write.results <- TRUE
  if(write.results == TRUE){
    write.csv(results.df, paste0(results.folder, run.folder, contrast.name, "_de.results.csv"))
  }
  
  
}

# Adjust column names
logfc.column <- colnames(results.df[grepl("logFC",colnames(results.df))])
results.df$logfc <- results.df[[logfc.column]]

pval.column <- colnames(results.df[grepl("_pval",colnames(results.df))])
results.df$pval <- results.df[[pval.column]]

adj.pval.column <- colnames(results.df[grepl("adjpval",colnames(results.df))])
results.df$padj <- results.df[[adj.pval.column]]

results.df$gene <- results.df$Gene

# Keep only the necessary columns
results.df <- results.df %>% select(c("gene", 
                                      "logfc", 
                                      "pval", 
                                      "padj"))


# Export the results
export.deg.list <- FALSE
if(export.deg.list == TRUE){ 
  write.csv(results.df, 
            file = paste0(results.folder, run.folder, "DE/DEG_lists/", contrast.name, "_deg_list.csv"), 
            row.names = FALSE)  
}

# Make the volcano plot
volcano.output <- make_volcano(lmm.results = results.df, 
                               title = contrast.name, 
                               legend.title = "Expression in T4", 
                               x.axis.title = "logFC: C (Ref) versus T4 (Cond)", 
                               fc.limit = 1, 
                               pos.label.limit = 1, 
                               neg.label.limit = -1)

# Export the volcano
export.volcano <- FALSE
if(export.volcano == TRUE){ 
  ggsave(filename = paste0(results.folder, run.folder, "DE/volcano/", contrast.name, "_volcano_plot.png"), 
         width = 14, 
         height = 10)
}

print(volcano.output$volcano.plot)

# Subset the annotation just for the heatmap annotations of interest
annotation.heatmap <- lmm.input$annotation %>% 
  select(names(anno.colors)) %>% 
  arrange(Treatment_group)

heatmap.plot <- make_heatmap(
  normalized.log.counts.df = lmm.input$log.counts, 
  de.results = results.df, 
  top.degs = TRUE, 
  annotation.column = annotation.heatmap,
  annotation.row = NULL, 
  anno.colors = anno.colors, 
  cluster.rows = TRUE, 
  cluster.columns = TRUE, 
  main.title = contrast.name, 
  row.gaps = NULL, 
  column.gaps = NULL, 
  show.colnames = FALSE, 
  show.rownames = TRUE)

print(heatmap.plot)

export.heatmap <- FALSE
if(export.heatmap == TRUE){
  
  ggsave(heatmap.plot, 
         filename = paste0(results.folder, run.folder, "DE/heatmap/", contrast.name, "_heatmap_plot.png"), 
         width = 14, 
         height = 10)
  
}



```

### Venn Diagrams for log FC cutoff of 1

```{r Setup Gene Lists log FC 1}
#| warning: false
#| message: false

# Load all DE results
results.files <- list.files(paste0(results.folder, run.folder))

up.results.list <- list()
down.results.list <- list()

logfc.cutoff <- 1

results.files <- results.files[grepl(".csv", results.files)]


for(file in results.files){
  
  # Load the results for a DE contrast
  de.results <- as.data.frame(read.csv(paste0(results.folder, run.folder, file)))
  
  # Filter for DEGs for up and down regulated
  logfc.column <- colnames(de.results)[grep("logFC", colnames(de.results))]
  adj.pval.column <- colnames(de.results)[grep("adjpval", colnames(de.results))]
  
  up.degs.data <- de.results %>% 
    filter(.data[[logfc.column]] > logfc.cutoff) %>% 
    filter(.data[[adj.pval.column]] < 0.05)
  
  down.degs.data <- de.results %>% 
  filter(.data[[logfc.column]] < -(logfc.cutoff)) %>% 
  filter(.data[[adj.pval.column]] < 0.05)
  
  # Gather the DE gene names and store
  up.degs <- up.degs.data$Gene
  down.degs <- down.degs.data$Gene
  
  contrast.name <- gsub("_de.results.csv", "", file)
  contrast.name <- gsub("_", " vs. ", contrast.name)
  
  if(length(up.degs) > 0){
    
    up.results.list[[contrast.name]] <- up.degs
    
  } else {
    
    # For group with no DEGs
    up.results.list[[contrast.name]] <- "no DEGs"
    
  }
  
  if(length(down.degs) > 0){
    
    down.results.list[[contrast.name]] <- down.degs
    
  } else {
    
    # For group with no DEGs
    down.results.list[[contrast.name]] <- "no DEGs"
    
  }
  
}
```

#### Up-regulated Venn Diagram

```{r Create venn diagram for up degs log FC 1}
#| warning: false
#| message: false


venn.up.degs <- venn.diagram(up.results.list, 
                            main = "Up-regulated DEGs each DE Contrast",
                            category.names = c(names(up.results.list)[1],
                                               names(up.results.list)[2], 
                                               names(up.results.list)[3], 
                                               names(up.results.list)[4]), 
                            fill = c("maroon", 
                                     "#E69F00", 
                                     "#56B4E9", 
                                     "#009E73"), 
                            filename = NULL)


grid.draw(venn.up.degs)
```

# **Cell Type Deconvolution**

#### Set up

```{r Cell Type Deconvolution}
#| warning: false
#| message: false

# Get the reference sets
# Example Human Reference
human.skin.ref <- download_profile_matrix(species = "Human",
                                       age_group = "Adult", 
                                       matrixname = "Skin_HCA")

# Example Mouse Reference
human.immune.ref <- download_profile_matrix(species = "Mouse",
                                       age_group = "Adult", 
                                       matrixname = "PeripheralBlood_MCA.RData")

# Derive the probe set background

# Extract the normalized read counts
object <- normalized.object.list$object
q3.normalized.reads <- object@assayData$q_norm

# Extract the gene (feature) names
feature.data <- fData(object)

# Derive the probe background
probe.background <- derive_GeoMx_background(norm = q3.normalized.reads,
                                            probepool = feature.data$Module,
                                            negnames = c("NegProbe-WTX"))


# Combine terms in the reference, this is optional
# Example list of combining terms, all terms must be included even if not combining some
# Combine common cell types
skin.combine.types <- list(
  Macrophages = c("Macrophages.DC.1", 
                  "Macrophages.DC.2", 
                  "Macrophages.DC.3"), 
  Melanocytes = c("Melanocytes.1", 
                  "Melanocytes.2"), 
  Pericytes = c("Pericytes.1", 
                  "Pericytes.2"), 
  Keratinocytes = c("Keratinocytes.1", 
                  "Keratinocytes.2"), 
  T.cells = c("T.cells"), 
  Secretory.reticular.Fibroblast = c("Secretory.reticular.Fibroblast"), 
  Pro.inflammatory.Fibroblast = c("Pro.inflammatory.Fibroblast"), 
  Secretory.papillary.Fibroblast = c("Secretory.papillary.Fibroblast"), 
  Mesenchymal.Fibroblast = c("Mesenchymal.Fibroblast"), 
  Vascular.EC = c("Vascular.EC"), 
  Lymphatic.EC = c("Lymphatic.EC"), 
  Erythrocytes = c("Erythrocytes")
)

# Run spatial decon for each reference set
# Example using human skin reference
spatial.decon.skin = runspatialdecon(object = object,
                      norm_elt = "q_norm",
                      raw_elt = "exprs",
                      X = human.skin.ref,
                      align_genes = TRUE, 
                      cellmerges = skin.combine.types)
```

#### Heatmap

```{r Deconvolution Heatmap}
#| warning: false
#| message: false

anno.colors = list(
      region = c(tumor = "pink",
                 vessel = "seagreen1", 
                 immune = "steelblue1", 
                 'immune/stroma' = "maroon"), 
      segment = c(full_ROI = "darkgoldenrod2", 
                  LANA_pos = "slateblue")
      )

# Subset the annotation just for the heatmap annotations of interest
spatial.annotation <- spatial.input$annotation %>% 
  select(c(names(anno.colors), AOI_ID))

# Generate a mapping df for AOI IDs
spatial.annotation$sample_ID <- rownames(spatial.annotation)
ID.mapping <- spatial.annotation %>% 
  select(sample_ID, AOI_ID)

# Change the annotation row names
rownames(spatial.annotation) <- spatial.annotation$AOI_ID
spatial.annotation <- spatial.annotation %>% 
  select(names(anno.colors))

# Gather the scores for each cell type with each ref set
skin.cell.scores <- as.data.frame(spatial.decon.skin$beta)
immune.cell.scores <- as.data.frame(spatial.decon.immune$beta)
bahman.cell.scores <- as.data.frame(spatial.decon.bahman$beta)

# Add a column for mapping to AOI ID
skin.cell.scores$sample_ID <- rownames(skin.cell.scores)
immune.cell.scores$sample_ID <- rownames(immune.cell.scores)
bahman.cell.scores$sample_ID <- rownames(bahman.cell.scores)

# Map the AOI ID
skin.cell.scores <- merge(skin.cell.scores, ID.mapping, by = "sample_ID")
immune.cell.scores <- merge(immune.cell.scores, ID.mapping, by = "sample_ID")
bahman.cell.scores <- merge(bahman.cell.scores, ID.mapping, by = "sample_ID")

# Replace the rownames and remove extra ID columns
rownames(skin.cell.scores) <- skin.cell.scores$AOI_ID
skin.cell.scores <- skin.cell.scores %>% 
  select(-AOI_ID, -sample_ID)

rownames(immune.cell.scores) <- immune.cell.scores$AOI_ID
immune.cell.scores <- immune.cell.scores %>% 
  select(-AOI_ID, -sample_ID)

rownames(bahman.cell.scores) <- bahman.cell.scores$AOI_ID
bahman.cell.scores <- bahman.cell.scores %>% 
  select(-AOI_ID, -sample_ID)


# Transform for combining with annotation
skin.cell.scores <- as.data.frame(t(skin.cell.scores))
immune.cell.scores <- as.data.frame(t(immune.cell.scores))
bahman.cell.scores <- as.data.frame(t(bahman.cell.scores))

# Check for cell types with all zeros and remove
row.sums <- rowSums(skin.cell.scores)
zero.sum.skin <- rownames(skin.cell.scores[row.sums == 0, ])
skin.cell.scores.cleaned <- skin.cell.scores[!rownames(skin.cell.scores) %in% zero.sum.skin, ]

row.sums <- rowSums(immune.cell.scores)
zero.sum.immune <- rownames(immune.cell.scores[row.sums == 0, ])
immune.cell.scores.cleaned <- immune.cell.scores[!rownames(immune.cell.scores) %in% zero.sum.immune, ]

row.sums <- rowSums(bahman.cell.scores)
zero.sum.immune <- rownames(bahman.cell.scores[row.sums == 0, ])
bahman.cell.scores.cleaned <- bahman.cell.scores[!rownames(bahman.cell.scores) %in% zero.sum.immune, ]


# Make heatmaps
skin.heatmap <- make_heatmap(normalized.log.counts.df = skin.cell.scores.cleaned,
                         annotation.column = spatial.annotation, 
                         anno.colors = anno.colors, 
                         cluster.rows = TRUE, 
                         cluster.columns = TRUE, 
                         main.title = "Spatial Decon Skin Reference", 
                         show.rownames = TRUE)


immune.heatmap <- make_heatmap(normalized.log.counts.df = immune.cell.scores.cleaned,
                         annotation.column = spatial.annotation, 
                         anno.colors = anno.colors, 
                         cluster.rows = TRUE, 
                         cluster.columns = TRUE, 
                         main.title = "Spatial Decon Immune Reference", 
                         show.rownames = TRUE)

bahman.heatmap <- make_heatmap(normalized.log.counts.df = bahman.cell.scores.cleaned,
                         annotation.column = spatial.annotation, 
                         anno.colors = anno.colors, 
                         cluster.rows = TRUE, 
                         cluster.columns = TRUE, 
                         main.title = "Spatial Decon Bahman Reference", 
                         show.rownames = TRUE)

# Export Heatmaps
global.export.decon <- TRUE
export.decon.heatmaps <- global.export.decon
if(export.decon.heatmaps == TRUE){
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_skin_ref_heatmap.png"), 
         skin.heatmap, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_immune_ref_heatmap.png"), 
         immune.heatmap, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_bahman_ref_heatmap.png"), 
         bahman.heatmap, 
         height = 10, 
         width = 14)
  
}

```

#### Proportion Bar Plots

```{r Deconvolution Proportion Bar Plots}
#| warning: false
#| message: false

# Calculate the percentage for each cell type in each AOI
skin.cell.perc <- sweep(skin.cell.scores.cleaned, 
                               2, 
                               colSums(skin.cell.scores.cleaned), 
                               FUN = "/") * 100

immune.cell.perc <- sweep(immune.cell.scores.cleaned, 
                               2, 
                               colSums(immune.cell.scores.cleaned), 
                               FUN = "/") * 100

bahman.cell.perc <- sweep(bahman.cell.scores.cleaned, 
                               2, 
                               colSums(bahman.cell.scores.cleaned), 
                               FUN = "/") * 100

# Transform for combining with annotation
skin.cell.perc.transform <- as.data.frame(t(skin.cell.perc))
skin.cell.perc.transform$AOI_ID <- rownames(skin.cell.perc.transform)

immune.cell.perc.transform <- as.data.frame(t(immune.cell.perc))
immune.cell.perc.transform$AOI_ID <- rownames(immune.cell.perc.transform)

bahman.cell.perc.transform <- as.data.frame(t(bahman.cell.perc))
bahman.cell.perc.transform$AOI_ID <- rownames(bahman.cell.perc.transform)

# Add column for mapping to annotation
spatial.annotation$AOI_ID <- rownames(spatial.annotation)

# Combine the cell type percentage and annotation and melt for use in ggplot
skin.cell.combine <- merge(skin.cell.perc.transform, spatial.annotation,
                           by = "AOI_ID")
skin.cell.melt <- melt(skin.cell.combine, 
                       id.vars = c("AOI_ID", 
                                   "region", 
                                   "segment"), 
                       variable.name = "cell.type"
                       )

immune.cell.combine <- merge(immune.cell.perc.transform, spatial.annotation,
                           by = "AOI_ID")
immune.cell.melt <- melt(immune.cell.combine, 
                       id.vars = c("AOI_ID", 
                                   "region", 
                                   "segment"), 
                       variable.name = "cell.type"
                       )

bahman.cell.combine <- merge(bahman.cell.perc.transform, spatial.annotation,
                           by = "AOI_ID")
bahman.cell.melt <- melt(bahman.cell.combine, 
                       id.vars = c("AOI_ID", 
                                   "region", 
                                   "segment"), 
                       variable.name = "cell.type"
                       )

# Make the bar plots
skin.bar.plot <- ggplot(data = skin.cell.melt, aes(x = AOI_ID, 
                                                   y = value, 
                                                   fill = cell.type)) + 
  geom_bar(stat = "identity", color = "black") + 
  facet_wrap(~region, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  )
  
immune.bar.plot <- ggplot(data = immune.cell.melt, aes(x = AOI_ID, 
                                                   y = value, 
                                                   fill = cell.type)) + 
  geom_bar(stat = "identity", color = "black") + 
  facet_wrap(~region, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  )

bahman.bar.plot <- ggplot(data = bahman.cell.melt, aes(x = AOI_ID, 
                                                   y = value, 
                                                   fill = cell.type)) + 
  geom_bar(stat = "identity", color = "black") + 
  facet_wrap(~region, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  )


# Export Barplots
export.decon.barplots <- global.export.decon
if(export.decon.barplots == TRUE){
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_skin_ref_barplot.png"), 
         skin.bar.plot, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_immune_ref_barplot.png"), 
         immune.bar.plot, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_bahman_ref_barplot.png"), 
         bahman.bar.plot, 
         height = 10, 
         width = 14)
  
}

```

#### Boxplots

```{r Deconvolution Boxplots}
#| warning: false
#| message: false

tumor.vessel.comparison <- list(
  c("tumor", "vessel"))

tumor.immune.comparison <- list(
  c("tumor", "immune"))

vessel.immune.comparison <- list(
  c("vessel", "immune"))

skin.box.plot <- ggplot(data = skin.cell.melt, aes(x = region, 
                                                   y = value, 
                                                   fill = segment)) + 
  geom_boxplot(color = "black") + 
  facet_wrap(~cell.type, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  ) + 
  labs(y = "Percentage of Total") + 
  stat_compare_means(comparisons = tumor.vessel.comparison, 
                     label = "p.signif", 
                     label.y = 20) + 
  stat_compare_means(comparisons = tumor.immune.comparison, 
                     label = "p.signif", 
                     label.y = 40) + 
  stat_compare_means(comparisons = vessel.immune.comparison, 
                     label = "p.signif", 
                     label.y = 30)

immune.box.plot <- ggplot(data = immune.cell.melt, aes(x = region, 
                                                   y = value, 
                                                   fill = segment)) + 
  geom_boxplot(color = "black") + 
  facet_wrap(~cell.type, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  ) + 
  labs(y = "Percentage of Total") + 
  stat_compare_means(comparisons = tumor.vessel.comparison, 
                     label = "p.signif", 
                     label.y = 20) + 
  stat_compare_means(comparisons = tumor.immune.comparison, 
                     label = "p.signif", 
                     label.y = 40) + 
  stat_compare_means(comparisons = vessel.immune.comparison, 
                     label = "p.signif", 
                     label.y = 30)

bahman.box.plot <- ggplot(data = bahman.cell.melt, aes(x = region, 
                                                   y = value, 
                                                   fill = segment)) + 
  geom_boxplot(color = "black") + 
  facet_wrap(~cell.type, scales = "free_x") + 
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank()
  ) + 
  labs(y = "Percentage of Total") + 
  stat_compare_means(comparisons = tumor.vessel.comparison, 
                     label = "p.signif", 
                     label.y = 20) + 
  stat_compare_means(comparisons = tumor.immune.comparison, 
                     label = "p.signif", 
                     label.y = 40) + 
  stat_compare_means(comparisons = vessel.immune.comparison, 
                     label = "p.signif", 
                     label.y = 30)

# Export Boxplots
export.decon.boxplots <- global.export.decon
if(export.decon.boxplots == TRUE){
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_skin_ref_region_boxplot.png"), 
         skin.box.plot, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_immune_ref_region_boxplot.png"), 
         immune.box.plot, 
         height = 10, 
         width = 14)
  
  ggsave(paste0(results.folder, 
                run.folder, 
                "spatial_decon/spatial_decon_bahman_ref_region_boxplot.png"), 
         bahman.box.plot, 
         height = 10, 
         width = 14)
  
}

```