CPTR9_DE_Analysis.qmd

---
title: "CPTR9_DE_Analysis"
format:
  html:
    code-fold: true
editor: visual
---

```{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 where the results should be exported
results.folder <- "results/"
run.folder <- "7_23_2024/"

```

```{r Load R object}
#| warning: false
#| message: false

# Load the R object containing normalized reads

normalized.object <- readRDS("CPTR9_normalized_dsp_object.RDS")

normalized.object
dim(normalized.object)
```

## Differential Expression

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


# Set up annotation colors for the heatmaps
anno.colors = list(
      cell_type = c(epithelial = "slateblue", 
                myeloid = "coral", 
                other = "green"), 
      region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue")
      
      )

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

```

### DE contrast: Myeloid vs. PanCK

```{r DE contrast: Myeloid vs. PanCK}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
segment.types <- c("Myeloid", "PanCK")

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

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


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

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_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, "DE/", 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 = "segment", 
                        regions = c("Myeloid", "PanCK"), 
                        group.col = "class", 
                        groups = c("cancer"), 
                        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, "DE/", 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 Myeloid", 
                               x.axis.title = "Myeloid (cond) versus PanCK (ref)", 
                               fc.limit = 1, 
                               pos.label.limit = 3, 
                               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)

# Make the heatmap

# Define new annotations if needed
anno.colors = list(region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue"))

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

# Order the annotation
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$region), ]
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$segment), ]

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)

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)
  
}



```

### DE contrast: PanCK tumor vs PanCK healthy epithelium

```{r DE contrast PanCK tumor vs PanCK healthy epithelium}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
segment.types <- c("PanCK")
region.types <- c("tumor", "healthyepithelium")


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

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list(
  "segment" = segment.types, 
  "region" = region.types)


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

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_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, "DE/", 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 = "region", 
                        regions = c("tumor", "healthyepithelium"), 
                        group.col = "segment", 
                        groups = c("PanCK"), 
                        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, "DE/", 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 Tumor", 
                               x.axis.title = "Tumor (cond) versus Healthy Epi (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)

# Make the heatmap

# Define new annotations if needed
anno.colors = list(region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue"))

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

# Order the annotation
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$region), ]
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$segment), ]

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)

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)
  
}



```

### DE contrast: Myeloid Tumor vs Myeloid Interstitial

```{r DE contrast Myeloid Tumor vs Myeloid Interstitial}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
segment.types <- c("Myeloid")
region.types <- c("tumor", "interstitial")


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

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list(
  "segment" = segment.types, 
  "region" = region.types)


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

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_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, "DE/", 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 = "region", 
                        regions = c("tumor", "interstitial"), 
                        group.col = "segment", 
                        groups = c("Myeloid"), 
                        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, "DE/", 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 Tumor", 
                               x.axis.title = "Tumor (cond) versus Intersitial (ref) logFC", 
                               fc.limit = 1, 
                               pos.label.limit = 1, 
                               neg.label.limit = -2.5)

# 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)

# Make the heatmap

# Define new annotations if needed
anno.colors = list(region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue"))

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

# Order the annotation
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$region), ]
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$segment), ]

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)

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)
  
}



```

### DE contrast: Myeloid Cluster vs Myeloid Interstitial

```{r DE contrast Myeloid Cluster vs Myeloid Interstitial}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
segment.types <- c("Myeloid")
region.types <- c("cluster", "interstitial")


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

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list(
  "segment" = segment.types, 
  "region" = region.types)


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

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_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, "DE/", 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 = "region", 
                        regions = c("cluster", "interstitial"), 
                        group.col = "segment", 
                        groups = c("Myeloid"), 
                        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, "DE/", 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 Cluster", 
                               x.axis.title = "Cluster (cond) versus Intersitial (ref) logFC", 
                               fc.limit = 1, 
                               pos.label.limit = 1, 
                               neg.label.limit = -1.5)

# 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)

# Make the heatmap

# Define new annotations if needed
anno.colors = list(region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue"))

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

# Order the annotation
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$region), ]
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$segment), ]

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)

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)
  
}



```

### DE contrast: Myeloid Cluster vs Myeloid Tumor

```{r DE contrast Myeloid Cluster vs Myeloid Tumor}
#| warning: false
#| message: false

# Setup annotation groups for the contrast
segment.types <- c("Myeloid")
region.types <- c("cluster", "tumor")


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

# Gather the column names to be used in defining contrasts
contrast.groups.list <- list(
  "segment" = segment.types, 
  "region" = region.types)


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

# Generate the counts, annotation, and subset object for lmm
lmm.input <- subset_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, "DE/", 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 = "region", 
                        regions = c("cluster", "tumor"), 
                        group.col = "segment", 
                        groups = c("Myeloid"), 
                        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, "DE/", 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 Cluster", 
                               x.axis.title = "Cluster (cond) versus Tumor (ref) logFC", 
                               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)

# Make the heatmap

# Define new annotations if needed
anno.colors = list(region = c(tumor = "pink",
                 interstitial = "seagreen1", 
                 healthyepithelium = "cadetblue1",
                 cluster = "darkorange"), 
      segment = c(PanCK = "lightsalmon", 
                  Myeloid = "darkolivegreen1", 
                  other_cells = "lightblue"))

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

# Order the annotation
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$region), ]
annotation.heatmap <- annotation.heatmap[order(annotation.heatmap$segment), ]

heatmap.plot <- make_heatmap(
  normalized.log.counts.df = lmm.input$log.counts, 
  de.results = results.df, 
  top.degs = FALSE, 
  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)

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)
  
}



```