Figure 1

Heatmap of significant genes per cluster

Load data and libraries

##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(cowplot)
library(ggrepel)
library(gridtext)
library(png)
library(grid)
#library(scatterpie)
library(patchwork)
library(openxlsx)
library(readxl)

# BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)

source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")

#########
# PATHS #
#########
input_dir <- "../../03_clustering_st_data/"
epi_clus <- "^5$|^6$|^7|^8" # res 0.7

#############
# LOAD DATA #
#############
DEGs_table <- read_csv(paste0("../../results/05_DGE_clusters_st_data/","DGEs_clusters_wilcox.0.7.csv"))
DATA <- readRDS(paste0("../../results/03_clustering_st_data/","seuratObj_clustered.RDS"))


sample_id <- c("P020", "P045", "P050", "P057",
               "P008", "P031", "P080", "P044", "P026", "P105", 
               "P001", "P004", "P014", "P018", "P087", "P118",
               "P021", "P024", "P067", "P081", "P117" ) %>% set_names()

Identify marker genes to seperate epithelial clusters

#######################################
# FILTER BY P-VAL logFC AND pct.diff #
######################################
filter_top.fun <- function(nested_df, n, var, min_pct = 0){
  if(is.na(n[4])){n[4] = min_pct}
  res <- nested_df %>% # res <- top20_epi$data[[3]] %>%
    #mutate(p_val_adj = 1.000-.$p_val_adj) %>%
    filter(p_val_adj < 0.05) %>%
    arrange(desc(avg_log2FC)) %>%
    filter(pct.1 > n[4]) %>%
    top_n(., n=n[1], abs(!!sym(var[1]))) %>%
    arrange(desc(avg_log2FC)) %>%
    #top_n(., n=n[2], Avg_exp) %>%
    top_n(., n=n[2], abs(!!sym(var[2]))) %>%
    top_n(., n=n[3], abs(!!sym(var[3]))) 
  print(res)
  return(res[1:n[3],])
}

sel_epi <- list(
  #Sup_1 = list(n=c(-40, 30, 5),var=c("p_val_adj", "pct.diff", "Avg_exp", "log.pct.diff")),
  "Superficial" = list(n=c(-40, 50, 5),var=c("p_val_adj", "pct.diff", "log.pct.diff")),
  #Sup_1 = list(n=c(-100, 70, 5, 0.5),var=c("pct.diff", "p_val_adj", "log.pct.diff")),
  "Upper IM" = list(n=c(-40, 50, 5),var=c("p_val_adj", "pct.diff", "log.pct.diff")),
  #Sup_2 = list(n=c(20, 5, 5),var=c("avg_log2FC", "log.pct.diff", "log.pct.diff")),
  "Lower IM" = list(n=c(25, 7, 5, .5),var=c("avg_log2FC", "log.pct.diff", "log.pct.diff")),
  "Basal" = list(n=c(10, 5, 5, .8),var=c("avg_log2FC", "log.pct.diff", "log.pct.diff")) )

# Identify the top genes that have a high difference in expression between the clusters
top20_epi <- DEGs_table %>%
  filter(grepl(epi_clus, .$clus)) %>%
  arrange(match(Clusters, names(sel_epi))) %>%
  mutate(Direction = ifelse(avg_log2FC > 0, "UP", "DOWN")) %>%
  filter(Direction == "UP") %>%
  #filter(!(pct.2 > 0.2)) %>%
  nest(., .by = "Clusters") %>%
  mutate(data = pmap(., ~filter_top.fun(..2, n=sel_epi[[..1]]$n, var=sel_epi[[..1]]$var)) ) %>%
  unnest(., cols=c("data"))

## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 7.79e-167       3.11 0.733 0.245 2.00e-162 5_all   FLG     -0.978         1.58
## 2 1.89e-179       2.44 0.721 0.213 4.85e-175 5_all   RNAS…   -0.934         1.76
## 3 1.16e-172       1.82 0.649 0.146 2.98e-168 5_all   DUOX…   -0.795         2.15
## 4 5.19e-164       1.82 0.733 0.246 1.33e-159 5_all   HMOX1   -0.979         1.58
## 5 6.36e-168       1.47 0.617 0.129 1.63e-163 5_all   DHRS9   -0.746         2.26
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 4.36e-142      1.01  0.863 0.353 1.12e-137 6_all   ENDOU   -1.22          1.29
## 2 3.17e-123      0.939 0.83  0.339 8.14e-119 6_all   FAM3D   -1.17          1.29
## 3 1.88e-125      0.906 0.841 0.353 4.84e-121 6_all   CYP4…   -1.19          1.25
## 4 8.67e-113      0.826 0.813 0.34  2.23e-108 6_all   CDA     -1.15          1.26
## 5 9.46e- 94      0.756 0.705 0.294 2.43e- 89 6_all   ACPP    -0.999         1.26
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##      p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene   pct.diff log.pct.diff
##      <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>     <dbl>        <dbl>
## 1 2.06e-83      0.714 0.895 0.537  5.29e-79 7_all   NDRG4     -1.43        0.737
## 2 1.00e-69      0.684 0.898 0.568  2.57e-65 7_all   C10or…    -1.47        0.661
## 3 2.04e-78      0.649 0.935 0.617  5.25e-74 7_all   NOTCH3    -1.55        0.600
## 4 5.16e-67      0.625 0.911 0.615  1.33e-62 7_all   DDR1      -1.53        0.567
## 5 1.80e-59      0.611 0.897 0.59   4.63e-55 7_all   RHOV      -1.49        0.604
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 4.10e-225       2.53 0.92  0.667 1.05e-220 8_all   IGFB…    -1.59        0.464
## 2 4.79e-172       2.31 0.814 0.44  1.23e-167 8_all   APOD     -1.25        0.888
## 3 3.16e-234       2.13 0.946 0.703 8.12e-230 8_all   DCN      -1.65        0.428
## 4 4.97e-218       2.12 0.924 0.668 1.28e-213 8_all   VIM      -1.59        0.468
## 5 6.84e-132       2.10 0.864 0.65  1.76e-127 8_all   SPARC    -1.51        0.411
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>

Identify marker genes to seperate Submucosal clusters

#######################################
# FILTER BY P-VAL logFC AND pct.diff #
######################################
filter_top.fun <- function(nested_df, n, var, filt_pct2=0.9){
  filt <- TRUE
  #if(is.na(n[4])){n[4] = pct}
  res <- nested_df %>% # res <- top20_sub$data[[7]] %>%
    #mutate(p_val_adj = 1.000-.$p_val_adj) %>%
    #filter(Avg_diff > 0) %>%
    #filter(pct.2 < n[4]) %>%
    top_n(., n=n[1], abs(!!sym(var[1]))) %>%
    {if(filt) filter(., between(pct.2, 0.01,filt_pct2 )) else .} %>%
    top_n(., n=n[2], abs(!!sym(var[2]))) %>%
    
    
    #{if(filt) filter(., between(pct.1, 0.8,0.3) | between(pct.2, 0.1,0.4)) else .} %>%
    #arrange(desc(Avg_exp_3)) %>%
    
    top_n(., n=n[3], abs(!!sym(var[3]))) #%>%
    #arrange(desc(!!sym(var[3])))
  print(res)
  return(res)
}

sel_s <- list(
  #"3" = list(n=c(-60, 30, 5),var=c("p_val_adj", "pct.1", "pct.1"), filt_pct2=.25),
  "1" = list(n=c(-60, 30, 5),var=c("p_val_adj", "pct.1", "pct.2"), filt_pct2=.25),
  #"3" = list(filt=F, n=c(-150, 100, 5),var=c("p_val_adj", "pct.diff", "log.pct.diff")),
  "4" = list( n=c(-100, 50, 5),var=c("p_val_adj", "pct.diff", "log.pct.diff"), filt_pct2=.9),
  "0" = list( n=c(-100, 50, 5),var=c("p_val_adj", "avg_log2FC", "log.pct.diff"), filt_pct2=.4),
  #"0" = list( n=c(-60, 30, 5),var=c( "pct.2", "avg_log2FC", "log.pct.diff"), filt_pct2=.4),
  #"0" = list( n=c( -50, 30, 5, 0.2),var=c( "p_val_adj", "avg_log2FC", "log.pct.diff"), filt=T),
  "3" = list(n=c(40, -50, 5),var=c("avg_log2FC", "pct.diff", "log.pct.diff"), filt_pct2=.9),
  "2" = list( n=c(-90, -20, -5),var=c("p_val_adj", "pct.2", "pct.2"), filt_pct2=.5),
  "9" = list( n=c(-90, 50, 5),var=c("p_val_adj", "pct.diff", "log.pct.diff"), filt_pct2=.9),
  "10" = list( n=c(-90, 5, 5),var=c("p_val_adj", "pct.1", "pct.1"), filt_pct2=.15)
  )

# Identify the top genes that have a high difference in expression between the clusters
top20_sub <- DEGs_table %>%
  filter(Regulation == "UP") %>%
  filter(!(grepl(epi_clus, .$clus))) %>%
  filter(!(grepl("^11$|^12$", .$clus))) %>%
  mutate(Direction = ifelse(avg_log2FC > 0, "UP", "DOWN")) %>%
  
  nest(., .by = "Clusters") %>%
  arrange(match(Clusters, names(sel_s))) %>%
  mutate(data = pmap(., ~filter_top.fun(..2, n=sel_s[[..1]]$n, var=sel_s[[..1]]$var, filt=sel_s[[..1]]$filt ))) %>%
  unnest(., cols=c("data"))

## # A tibble: 5 × 12
##      p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene   pct.diff log.pct.diff
##      <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>     <dbl>        <dbl>
## 1 1.98e-83      0.633 0.557 0.241  5.09e-79 1_all   FXYD3    -0.798         1.21
## 2 3.85e-80      0.614 0.56  0.246  9.89e-76 1_all   TRIM29   -0.806         1.19
## 3 6.06e-53      0.401 0.44  0.202  1.56e-48 1_all   C19or…   -0.642         1.12
## 4 9.46e-52      0.449 0.394 0.174  2.43e-47 1_all   SPINT2   -0.568         1.18
## 5 3.05e-51      0.443 0.398 0.179  7.83e-47 1_all   DSG1     -0.577         1.15
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 7.93e-237      1.62  0.704 0.222 2.04e-232 4_all   ACKR1   -0.926         1.67
## 2 5.10e-160      1.10  0.398 0.084 1.31e-155 4_all   TSPA…   -0.482         2.24
## 3 1.46e-137      1.10  0.471 0.135 3.76e-133 4_all   RAMP3   -0.606         1.80
## 4 9.27e-105      0.844 0.368 0.103 2.38e-100 4_all   IL33    -0.471         1.84
## 5 4.94e- 79      0.745 0.349 0.114 1.27e- 74 4_all   ITM2A   -0.463         1.61
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##      p_val avg_log2FC pct.1 pct.2  p_val_adj cluster gene  pct.diff log.pct.diff
##      <dbl>      <dbl> <dbl> <dbl>      <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 1.39e-10     0.0260 0.065 0.143 0.00000357 0_all   NDUF…   -0.208        -1.14
## 2 7.72e-10     0.0369 0.081 0.163 0.0000198  0_all   RGL2    -0.244        -1.01
## 3 2.28e- 9     0.0101 0.063 0.135 0.0000587  0_all   RBM42   -0.198        -1.10
## 4 3.38e- 9     0.0681 0.073 0.148 0.0000869  0_all   SNRPC   -0.221        -1.02
## 5 5.03e- 9     0.0149 0.063 0.133 0.000129   0_all   WARS    -0.196        -1.08
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##      p_val avg_log2FC pct.1 pct.2  p_val_adj cluster gene  pct.diff log.pct.diff
##      <dbl>      <dbl> <dbl> <dbl>      <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 1.07e-11      0.396 0.149 0.288    2.76e-7 3_all   TAX1…   -0.437       -0.951
## 2 9.47e-11      0.364 0.123 0.24     2.43e-6 3_all   PCM1    -0.363       -0.964
## 3 5.58e- 9      0.343 0.069 0.146    1.43e-4 3_all   UPF3A   -0.215       -1.08 
## 4 2.55e- 8      0.414 0.089 0.173    6.55e-4 3_all   USP8    -0.262       -0.959
## 5 8.00e- 7      0.341 0.056 0.114    2.05e-2 3_all   ZFC3…   -0.17        -1.03 
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##      p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##      <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 5.94e-27      0.590 0.235 0.113  1.53e-22 2_all   OGN     -0.348        1.06 
## 2 4.14e-11      0.345 0.154 0.088  1.06e- 6 2_all   PGBD5   -0.242        0.807
## 3 5.72e-11      0.364 0.13  0.071  1.47e- 6 2_all   PLK5    -0.201        0.873
## 4 9.83e-11      0.344 0.19  0.119  2.52e- 6 2_all   WNT16   -0.309        0.675
## 5 1.05e-10      0.248 0.182 0.11   2.70e- 6 2_all   CADM3   -0.292        0.726
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 3.67e-262      2.76  0.72  0.129 9.43e-258 9_all   DES     -0.849         2.48
## 2 8.82e-248      2.76  0.65  0.108 2.27e-243 9_all   ACTG2   -0.758         2.59
## 3 4.12e-137      1.49  0.297 0.033 1.06e-132 9_all   KCNM…   -0.33          3.17
## 4 2.00e-104      1.56  0.329 0.057 5.14e-100 9_all   PCP4    -0.386         2.53
## 5 1.63e- 65      0.946 0.237 0.045 4.18e- 61 9_all   SMOC2   -0.282         2.40
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>
## # A tibble: 5 × 12
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  pct.diff log.pct.diff
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>    <dbl>        <dbl>
## 1 1.18e-291       2.17 0.728 0.138 3.04e-292 10_all  MUST…   -0.866         2.40
## 2 1.39e-222       1.51 0.573 0.104 3.57e-218 10_all  MCAM    -0.677         2.46
## 3 3.44e-158       1.16 0.472 0.094 8.83e-154 10_all  FRZB    -0.566         2.33
## 4 1.49e-141       1.15 0.458 0.099 3.82e-137 10_all  ITGA7   -0.557         2.21
## 5 1.43e-131       1.20 0.486 0.121 3.67e-127 10_all  CCDC3   -0.607         2.01
## # ℹ 3 more variables: clus <dbl>, Regulation <chr>, Direction <chr>

clus_2 <- DEGs_table %>%
  filter(Regulation == "DOWN") %>%
  filter(!(grepl(epi_clus, .$clus))) %>%
  
    filter(Regulation == "DOWN") %>%
    filter((grepl("^2$", .$clus))) %>%
    #top_n(., n=-10, abs(pct.2)) %>%
    filter(., between(pct.2, 0.3,0.8 )) %>%
    filter(., between(pct.1, 0.1,0.3 ))
    # #slice_min(order_by = p_val_adj, n = 20) %>%


c( top20_sub$gene[1:15], clus_2$gene, top20_sub$gene[21:35]) 

##  [1] "FXYD3"    "TRIM29"   "C19orf33" "SPINT2"   "DSG1"     "ACKR1"   
##  [7] "TSPAN7"   "RAMP3"    "IL33"     "ITM2A"    "NDUFB6"   "RGL2"    
## [13] "RBM42"    "SNRPC"    "WARS"     "DSG3"     "HOPX"     "CALML3"  
## [19] "KRT16"    "A2ML1"    "OGN"      "PGBD5"    "PLK5"     "WNT16"   
## [25] "CADM3"    "DES"      "ACTG2"    "KCNMA1"   "PCP4"     "SMOC2"   
## [31] "MUSTN1"   "MCAM"     "FRZB"     "ITGA7"    "CCDC3"

clus_2_genes <- c("CD151", "C12orf57", "NPDC1", "TGFBR2", "PECAM1")

# Housekeeping genes
# clus_2_genes <- c("RPL17", "TUBA","EF1A", "B2M", "RPSD", "GAPDH")

# probably not going to use
Idents(DATA) <- "layers"
Idents(DATA) <- "layers"
DATA %>%
filter(., grepl(epi_clus, DATA$Clusters)) %>%
  RidgePlot(., features = top20_epi$gene[1:2], ncol = 2)

# epithelium
genes <- top20_epi$gene
col <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3")
# Submucosa
genes <- c( top20_sub$gene[1:10],clus_2_genes , top20_sub$gene[21:35]) 
col <- c("#00A9FF","#377EB8","#CD9600","#7CAE00","#e0e067","#FF61CC","#FF9DA7")

df <- DATA %>%
  filter(., grepl(epi_clus, DATA$Clusters)) %>% 
  #filter(., !(grepl("^5$|^6$|^7|^8|^11$|^12$", DATA$Clusters))) %>% 
  mutate(., FetchData(., vars = genes)) %>%
  as_tibble() %>%
  select(.cell, sp_annot, orig.ident, groups, layers, any_of(genes)) %>%
  pivot_longer(cols = any_of(genes), names_to = "Genes", values_to = "Expression") %>%
  mutate(Genes = factor(.$Genes, levels=rev(genes))) %>%
  arrange(layers) %>%
  nest(data= -layers) %>%
  mutate(col = col)


p <- df %>%
  mutate(plot = pmap(., 
        ~ggplot(..2, aes(x = Expression, y = Genes, fill = ..3)) +
                geom_density_ridges() +
                scale_fill_identity() +
                xlim(c(0,4)) +
                theme_ridges() + labs(title = ..1 )  +
                theme(legend.position = "none",
                      plot.margin = unit(c(0,-0,0,0),units = "pt"),
                      axis.title.y = element_blank())
  ))

(p_ <- plot_grid(plotlist = p$plot, ncol=length(col)) )
ggsave("./Figures/Marker_genes_ridges_epi.pdf", p_, width = 12, height = 5)
ggsave("./Figures/Marker_genes_ridges_sub.pdf", p_, width = 21, height = 7)

library(ComplexHeatmap)
library(circlize)

#### GENES THAT ARE UNIQELY OVERLAPING ####
getVennOverlap <- function(lsvenn) {
  
  ItemsList <- gplots::venn(lsvenn, show.plot = FALSE)
  print(lengths(attributes(ItemsList)$intersections))
  return(attributes(ItemsList)$intersections)
}

################
# PLOT HEATMAP #
################
heatmap.fun <- function(genes, cluster){
###############
# GET MATRIX #
###############
DAT <- filter(DATA, grepl(paste0("^",cluster,"$",collapse="|"), DATA$layers))
VariableFeatures(DAT) <- genes
DAT <- Seurat::ScaleData(DAT)
matx <- DAT@assays$RNA$scale.data %>% 
  as.matrix()
  
  
##############
# ANNOTATION #
##############
annot_col <<- DAT@meta.data %>%
  #mutate(layers = factor(.$layers, levels = cluster)) %>%
  dplyr::select(1:7) 

# set row order
matx <- matx[genes,]
#print( dput(rownames(matx)))
#col_order <<- rownames(arrange(annot_col, layers) )

clus_cols <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3",
               "#00A9FF","#377EB8","#CD9600","#7CAE00","#e0e067","#FF61CC","#FF9DA7","#BAB0AC","#A65628")
clus <- c("Superficial","Upper IM","Lower IM","Basal","1","4","0","3","2","9","10","11","12") #
clus_cols <- set_names(clus_cols, clus)
clus_cols <- clus_cols[unique(annot_col$layers)]

ID_cols <- c("#4E79A7","#F28E2B","#E15759","#76B7B2","#59A14F","#EDC948","#9C755F","#B07AA1","#FF9DA7",
             "#E78AC3","#7FC97F","#BEAED4","#BAB0AC","#FDC086","#FFFF99","#386CB0","#F0027F","#BF5B17",
             "#66C2A5","#FC8D62","#8DA0CB")
ID_cols <- set_names(ID_cols, sample_id)
ID_cols <- ID_cols[unique(annot_col$orig.ident)]

group_cols <- c("#56B4E9","#009E73","#CC79A7","#FC8D62")
gr <- c('L1', 'L2', 'L3', 'L4') %>% set_names()
group_cols <- set_names(group_cols, gr)
group_cols <- group_cols[unique(annot_col$groups)]

################
# DRAW HEATMAP #
################

# qwartile 
q <- quantile(matx, c(0.05, 0.95))
c(min(matx), max(matx))

# average expression:
set.seed(123)
H <- Heatmap(matx, name = " ",
             col = circlize::colorRamp2(c(q[1]*2, q[1], 0, q[2], q[2]*2), 
                                        c("#440154FF","#6A51A3","#21908CFF","#FDE725FF","#FF7F00")
                                        #rev(c("#D7191C","#FDAE61", "white","#ABD9E9","#4575B4"))
                                        ),
             #row_km = 6, #column_km = 2, # kmeans change every time you run it
             column_split =  factor(as.character(annot_col$layers), levels = cluster),
             #row_split = 10,  # hierarchical static
             row_order = rownames(matx),
             column_order = rownames(arrange(annot_col, layers) ),
             show_column_names = FALSE,
             show_row_dend = FALSE,
             show_column_dend = FALSE,
             #cluster_columns = T,
             use_raster = F, # removes white vertical lines in the heatmap
             
             # text
             row_title = gt_render("", padding = unit(c(0, 0, 0, 0), "pt")),
             
             #column_labels = gt_render(column_labels, padding = unit(c(0, 0, 0, 0), "pt")),
             column_title_gp = grid::gpar(fontsize = 8),
             #column_labels_gp =  grid::gpar(fontsize = 8),
             
             column_names_gp = grid::gpar(fontsize = 8),
             row_names_gp = grid::gpar(fontsize = 8),
            
             
             # annotation
             # right_annotation = right_anno_row, left_annotation = left_anno_row,
             top_annotation =
               columnAnnotation(Cluster=annot_col$layers,
                                #show_legend = FALSE,
                                show_annotation_name = F,
                                annotation_legend_param =
                                  list(grid_height = unit(.2, "mm"),
                                       grid_width = unit(2, "mm"), title = "",
                                       labels_gp = gpar(fontsize = 7),
                                       title_gp = gpar(fontsize = 8)),
                                simple_anno_size = unit(.1, "cm"),
                                #gap = unit(1, "cm"),
                                col=list( Cluster = clus_cols )
                                ),

             # legend
             heatmap_legend_param = list(legend_height = unit(20, "mm"), grid_width = unit(2, "mm"), 
                                         labels_gp = gpar(fontsize = 7))
             #annotation_legend_param = list(size = unit(2, "mm"))

             ) #colorRampPalette(c(col))(10) 

H <- draw(H, merge_legend = TRUE)
return(H)
}

# Heatmap global options:
ht_opt$COLUMN_ANNO_PADDING = unit(.05, "cm")
ht_opt$HEATMAP_LEGEND_PADDING = unit(0, "cm")
ht_opt$TITLE_PADDING = unit(.05, "cm")
ht_opt$DIMNAME_PADDING = unit(.05, "cm")

lvl_layers <- c("Superficial", "Upper IM", "Lower IM", "Basal","1","4","0","3","2","9","10","11","12")

# all clusters together: 
# genes <- c(top20_epi$gene, top20_sub$gene)
# H <- heatmap.fun(genes = genes, lvl_layers)

##############
# EPITHELIUM #
##############
genes <- top20_epi$gene
H_epi <- heatmap.fun(genes = genes, lvl_layers[1:4])

# save file
H_grob <- grid.grabExpr(draw(H_epi, merge_legend = TRUE))  # show_heatmap_legend = FALSE,
ggsave("./Figures/01/heatmap_clusters_epi.png", H_grob, width = 6, height = 3, limitsize = F)

#############
# SUBMUCOSA #
#############
genes <- c( top20_sub$gene[1:10], clus_2_genes, top20_sub$gene[21:35]) 
H_sub <- heatmap.fun(genes = genes, lvl_layers[5:11])

# save file
H_grob <- grid.grabExpr(draw(H_sub, merge_legend = TRUE))  # show_heatmap_legend = FALSE,
ggsave("./Figures/01/heatmap_clusters_sub.png", H_grob, width = 6, height = 4, limitsize = F)
# NB! when working with the heatmap, make sure output is printed to the console

#######################
# DEFINE MARKER GENES #
#######################
Cell_marker <- c("MS4A1","CD79A", "POU2AF1", "MZB1", "FCRL5", "CD38",
                 "CD247","CD8A","CD3D","CD3G","CD4","CXCR4",  
                 "NCAM1", "GZMA", "GNLY",  "FCER1A", "CD1A", "LTB",  
                 "COL3A1","COL1A1", "APOD","LAMA2", "MYH11","ICAM2","PECAM1","MCAM",   # FLT1 #"COL3A1","COL1A1", "MYH11","LAMA2"
                 "ATG9B", "KPRP","ALOX12", "PRSS3", "BICDL2", "MT1X")
names <- c("B cell", "B cell","B cell", "Plasma cell","Plasma cell", "Plasma cell", 
           "T cell", "T cell","T cell", "T cell", "T cell", "T cell", 
           "NK cells", "NK cells", "NK cells","Myeloid", "Myeloid", "Myeloid", 
           "Fibroblast", "Fibroblast", "Fibroblast","Fibroblast", 
           "Endothelial","Endothelial", "Endothelial", "Endothelial", "Keratinocyte supra", "Keratinocyte supra", 
           "Keratinocyte supra","Keratinocyte supra", "Keratinocyte basal", "Keratinocyte basal")
cell_type <- set_names(names, Cell_marker)
clus_lvl <- rev(c("5", "6", "7", "8", "1","4","0","3","2","9","10","11","12")) 

clus_lvl <- set_names(seq_along(clus_lvl), clus_lvl)
gene_lvl <- set_names(seq_along(Cell_marker), Cell_marker)

################
# FETCH GENES #
################
cell_type <- set_names(names, Cell_marker)
df <- DATA %>%
  mutate(., FetchData(., vars = c(Cell_marker)) ) %>%
  as_tibble() %>%
  select(., .cell, Clusters, any_of(Cell_marker)) %>%
  pivot_longer(., cols = -c(".cell", "Clusters"), 
               names_to = "marker", values_to = "values") %>%
  mutate(., cell = cell_type[as.character(.$marker)]) %>%
  mutate( marker_id = paste0(.$cell," (",.$marker,")")) %>%
  #filter(., Clusters != "11"| Clusters != "12") %>%
  mutate(Clusters = factor(.$Clusters, levels=names(clus_lvl))) %>%
  mutate(marker = factor(.$marker, levels=Cell_marker)) %>%
  group_by(Clusters, marker, cell) %>%
  summarise(Avg = mean(values),
            Pct = sum(values > 0) / length(values) * 100, .groups="drop") %>%
  mutate(., ymin = clus_lvl[as.character(.$Clusters)]-0.5,
            ymax = clus_lvl[as.character(.$Clusters)]+0.5) %>%
  mutate(., xmin = gene_lvl[as.character(.$marker)]-0.5,
            xmax = gene_lvl[as.character(.$marker)]+0.5)

########################
# MARKER GENES DOTPLOT #
########################
library(ggnewscale)
cell_col <- c("#FFD92F","#FFFFCC","#8DA0CB","#eb6062","#A6D854","#E78AC3","#66C2A5","#FC8D62","#FED9A6",
         "#377EB8","#4DAF4A","#B3B3B3","#984EA3","#FF7F00","#FFFF33","#A65628","#F781BF")

clus_col <- rev(c("#FBAAB1","#FFDAB8","#F1D2FF","#E6B3E9","#92DBFF","#ABC9E1", "#E5C264","#BBD99B", "#E1E2A4", "#FF61CC","#FA9DA7", "#B3B3B3", "#A65628"))
lab <- c("B cell", "Plasma cell", "T cell", "NK cells", "Myeloid", "Fibroblast", 
"Endothelial", "Keratinocytes")

ym <- max(df$ymax)
xm <- length(cell_type)+3.5

clus_lvl <- rev(c("5","6","7","8","1","4","0","3","2","9","10","11","12"))
# DATA %>%
#   mutate(Clusters = factor(.$Clusters, levels = clus_lvl)) %>%
# ggplot(., aes(x=Clusters, y=nCount_RNA, fill=Clusters)) + geom_violin() + scale_fill_manual(values = clus_col)

(A <- ggplot(df, aes(x=marker, y=Clusters)) +
  geom_point(aes(size = Pct, fill = Avg), color="white", shape=21) +
  scale_fill_gradientn(colours = viridisLite::magma(100),
                       guide = guide_colorbar(ticks.colour = "white",
                                              frame.colour = "white",
                                              barwidth = .5, barheight = 4),
                       name = "Average\nexpression") +
  #facet_grid(~ cell, scales = "free_x") +
  # Cell type colour bar
  new_scale_fill() +
  geom_rect(aes(ymin=max(ymax), ymax=max(ymax)+.4,
                xmin=xmin, xmax=xmax,
                fill = cell),data=df,alpha = 0.1,show.legend=F) +
  annotate("text", x = c(1.8, 5, 10, 14, 17, 20.5, 24.6, 29.5 ), y = 13.7, label = lab) +
  geom_rect(aes(ymin=min(ymin), ymax=min(ymin)+.3,
                xmin=xmin, xmax=xmax,
                fill = cell),data=df,alpha = 0.1,show.legend=F) +
  scale_fill_manual(values = set_names(cell_col[1:length(unique(cell_type))], unique(cell_type))) +
    
  # add extra borders for the colour bars
  # annotate(x = c(.5,xm,-.5,-.5,.5,-.5), xend=c(.5,xm,xm,xm,xm,-.5),
  #            y=c(1,1,.5,ym,ym+.5,.5), yend=c(ym+.5,ym+.5,.5,ym,ym+.5,ym),
  #            geom="segment",colour = "black", linewidth = .6, alpha = 1) +
  # Cluster bar annnotation
  new_scale_fill() +
  geom_rect(aes(ymin=ymin, ymax=ymax,
                xmin=0.4, xmax=-.5,
                fill = Clusters),data=df,alpha = 1,show.legend=F) +  
  scale_fill_manual(values = clus_col) +
  coord_cartesian(clip="off", xlim=c(.5,xm),ylim=c(.5,ym),expand = F) +
  scale_size("% detected", range = c(0,6)) +

  ylab("Cluster") + xlab("") +
  theme_bw() +
  guides(size = guide_legend(override.aes = list(color = "black"), keywidth = .6, keyheight = .6)) +
  theme(axis.text.x = element_text(size=8, angle=45, hjust=1, color="black"),
        axis.text.y = element_text(size=9, color="black",hjust=.9 ),
        axis.title = element_text(size=9),
        axis.ticks.y = element_blank(),
        panel.border = element_blank(),
        legend.title = element_text(size = 9),
        legend.margin=margin(0,0,-0,0),
        plot.margin = unit(c(1,0,-.4,.1),units = "cm") #trbl
        ) 
)

# dev.new(width=8, height=3.5, noRStudioGD = TRUE)
# ggsave("./Figures/01/marker-gene-dotplot_new.pdf", A, width = 8, height = 3.5) # cell_marker

Not used!

library(edgeR)
pseudo.bulk <- function(DATA, group_var, degs=FALSE){
  # create ID per cluster group
  DATA <- mutate(DATA, ID_clus = paste0(DATA$orig.ident, "_", DATA[[group_var]]))
  
  # get the count matrix for all cells
  sparse_mtrx <- DATA@assays$RNA@counts
  
  # Compute pseudobulk
  mm <- Matrix::sparse.model.matrix(~0 + DATA$ID_clus)
  pseudobulk <- sparse_mtrx %*% mm
  colnames(pseudobulk) <- str_replace(colnames(pseudobulk), "DATA\\$ID_clus", "")
  
  # group
  gr <- str_replace(colnames(pseudobulk), "P\\d\\d\\d_", "")
  gr <- set_names(gr, colnames(pseudobulk))
  
  if(length(degs) > 1){
    # filter genes based on DEG analysis
    matx <- as.matrix(pseudobulk)
    matx <- matx[degs,]
  }else{
    # filter genes based on groups
  dge.list <- DGEList(counts = pseudobulk, group = factor(gr))
  keep <- filterByExpr(dge.list)
  matx <- dge.list[keep, , keep.lib.sizes = FALSE]}
  
  return(matx)
}

# remove specified genes:
remove_genes <- function(x, gene_name) x[!(grepl(gene_name, rownames(x[["RNA"]]))), ]

sample_id <- c("P020","P045","P050","P057",
               "P008","P026","P031","P044","P080","P105", 
               "P001","P004","P014","P018","P087","P108","P118",
               "P021","P024","P067","P081","P117")
groups <- c( "L1","L1","L1","L1",
             "L2","L2","L2","L2","L2","L2",
             "L3","L3","L3","L3","L3","L3","L3",
             "L4","L4","L4","L4","L4")
id_gr <- set_names(groups, sample_id)
######################
# GET DEGS  # 
######################
pseudo_DEGs <- readRDS(paste0("../results/06_DGE_condition_st_data/","Pseudobulk_across_DEGs.RDS"))
pseudo_DEGs <- pseudo_DEGs %>%
  #filter(layers != "11" | layers != "12" ) %>% unnest() %>%
  mutate(test = map(test, ~filter(.x, PValue <= 0.001) )) %>%
  unnest(test) %>% .$symbol %>% unique()
######################
# NORMALIZE TRX DATA # 
######################
matx <- pseudo.bulk(DATA, "layers", degs = pseudo_DEGs)
matx <- pseudo.bulk(DATA, "orig.ident")
#### CPM-log2 ####
# counts per million (CPM) log2 transfomed counts:
matrix_trans <- cpm( matx, log=T)

########
# UMAP # 
########
lvl_layers <- c("Superficial","Upper IM","Lower IM","Basal","3","4","0","2","1","9","10","11","12")
col <- c("#E41A1C","#FF7F00","#C77CFF","#984EA3","#00A9FF","#377EB8","#CD9600","#7CAE00","#e0e067","#FF61CC","#FF9DA7","#999999","#A65628")

lvl_layers <- sample_id
col <-  c("#4E79A7", "#F28E2B", "#E15759", "#76B7B2", "#59A14F", "#EDC948", "#9C755F", "#B07AA1", "#FF9DA7", "#E78AC3", "#7FC97F", "#BEAED4", "#BAB0AC", "#FDC086", "#FFFF99", "#386CB0", "#F0027F", "#BF5B17", "#666666", "#66C2A5", "#FC8D62", "#8DA0CB")

#### UMAP analysis ####
set.seed(1);umap<-uwot::umap(t(matrix_trans), 
                 n_neighbors=15, init="spectral", scale=T)
umap_df <- tibble("UMAP 1"=umap[,1],"UMAP 2"=umap[,2], 
                  "u_ID"=rownames(umap), "group"=gr[rownames(umap)], "ID"=str_extract(rownames(umap), "P\\d\\d\\d") ) %>%
            mutate(group = factor(.$group, levels = lvl_layers)) %>%
            mutate(g = id_gr[.$ID]) 


txt_df <- umap_df %>% 
  filter(., grepl("^P", .$txt)) #%>%
  #mutate("UMAP 2" = .$`UMAP 2`+ c(0, .1, rep(0,ncol(.))))

#### UMAP plotting by groups ####
UMAP.fun <- function(umap_df, group){
  # dev.new(width=8.5, height=8)
  p <- ggplot(umap_df, aes(x=`UMAP 1`, y=`UMAP 2`, fill={{group}}))+ 
    #geom_jitter( shape=21, size=3, color="white", width=.5, height=.5) +  # Tassos used jitter
    geom_point( shape=21, size=3, alpha = .7, stroke =1) +  
    scale_fill_manual(values=col) + 
    # geom_text(data=txt_df, aes(x=`UMAP 1`, y=`UMAP 2`, label=txt), size=3, vjust = -.5, nudge_x = 0.03, color="gray51") +
    theme_void() +
    xlab("UMAP 1") + ylab("UMAP 2") +
    theme(panel.border = element_rect(colour = "black", fill=NA),
          legend.title = element_text(size = 8),
          legend.text = element_text(size = 8),
          legend.position = "bottom",
          plot.margin = unit(c(1,.5,1,.5), "lines"),
          axis.title.x.bottom = element_text(margin = margin(b = 1)),
          axis.title.y.left = element_text(margin = margin(l = 0), angle = 90),
          axis.ticks=element_line(color="black"), 
          axis.title=element_blank()) 
  return(p)
}

(p <- UMAP.fun(umap_df, g))
(p_T <- UMAP.fun(umap_df, Tissue_gr_v3))

# dev.new(width=8, height=3.5, noRStudioGD = TRUE)
(p_ <- plot_grid(p_L, p_T, ncol = 2 ) )