diff --git a/README.md b/README.md
index 9be7b78..1376d27 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,9 @@
-# kirchmair_2023
+## Kirchmair et al., Frontiers in Immunology 2023
 
-Data analyses for Kirchmair et al., 2023
+Data analyses for '<sup>13</sup>C tracer analysis reveals the landscape of metabolic checkpoints in human CD8<sup>+</sup> T cell differentiation and exhaustion' (Kirchmair et al., Frontiers in Immunology 2023)
 
 
+ 
 ## Environment setup
 ```bash
 # source lib/make_env.sh # initial code to make conda envs
@@ -17,7 +18,7 @@ Rscript -e 'devtools::install_github("AlexanderKirchmair/DeLuciatoR")' # version
 
 mkdir logs
 ```
-
+ 
 Set up [NGSCheckMate-1.0.0](https://github.com/parklab/NGSCheckMate)
 ```bash
 cd lib
@@ -36,8 +37,8 @@ cd ..
 ```bash
 conda activate cd8
 ```
-
-Memory differentiation samples: [GSE234099](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE234099)
+ 
+Memory differentiation samples ([GSE234099](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE234099)):
 ```bash
 mkdir -p data/rnaseq/MEM/00_RAW
 accs=$(awk 'NR>1 {print $2 "-" $1}' "tables/GSE234099.txt")
@@ -47,8 +48,8 @@ do
   while [ $(qstat -s pr | grep -w -c "DOWNLOAD") -gt 3 ]; do sleep 3; done
 done
 ```
-
-Exhaustion samples: [GSE234100](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE234100)
+ 
+Exhaustion samples ([GSE234100](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE234100)):
 ```bash
 mkdir -p data/rnaseq/EXH/00_RAW
 accs=$(awk 'NR>1 {print $2 "-" $1}' "tables/GSE234100.txt")
@@ -58,7 +59,7 @@ do
   while [ $(qstat -s pr | grep -w -c "DOWNLOAD") -gt 3 ]; do sleep 3; done
 done
 ```
-
+ 
 ### Preprocessing
 
 Trimming:
@@ -76,7 +77,7 @@ do
 done
 ```
 
-Read alignment and quantification using the [nf-core/rnaseq-3.4](https://nf-co.re/rnaseq/3.4) pipeline (set genome paths in 'lib/run_rnaseq.sh'):
+Read alignment and quantification using the [nf-core/rnaseq-3.4](https://nf-co.re/rnaseq/3.4) pipeline (set genome paths in `lib/run_rnaseq.sh`):
 ```bash
 bash -i lib/run_rnaseq.sh 'tables/samplesheet_mem.csv' 'data/rnaseq/MEM/02_NF_results'
 mv .nextflow.log logs/mem.nextflow.log
@@ -103,17 +104,21 @@ mv r_script.r.Rout data/rnaseq/EXH/samplecheck/
 Rscript lib/plot_NGSCheckMate.R
 ```
 
+Gene sets were prepared by running `Rscript lib/prepare_genesets.R`.
 
-## Metabolomics data
-13C metabolomics data: 'data/metabolomics'
 
+## Metabolomics data
+<sup>13</sup>C metabolomics data: `data/metabolomics`
+ 
 
 ## Seahorse data
-Seahorse data: 'data/seahorse'
-
+Seahorse data: `data/seahorse`
+ 
 
 ## Analysis
-Gene sets were prepared by running 'Rscript lib/prepare_genesets.R'.
+
+The main analyses can be reproduced by rendering the the .Rmd files:
+
 ```bash
 conda activate cd8
 Rscript -e "rmarkdown::render('analyses/01-RNA-Differentiation.Rmd')"
@@ -122,12 +127,16 @@ Rscript -e "rmarkdown::render('analyses/03-RNA-Exhaustion.Rmd')"
 Rscript -e "rmarkdown::render('analyses/04-13C-Exhaustion.Rmd')"
 Rscript -e "rmarkdown::render('analyses/05-RNA-Exhaustion-Public.Rmd')"
 Rscript -e "rmarkdown::render('analyses/06-RNA-Mitochondria.Rmd')"
+Rscript -e "rmarkdown::render('analyses/07-Public-Dataset-Comparison.Rmd')"
 ```
 
 
-## Results (figures and tables)
+## Results
+
+To reproduce the final figures and tables, run:
+
 ```bash
 conda activate cd8
-Rscript -e "rmarkdown::render('analyses/07-Results.Rmd')"
+Rscript -e "rmarkdown::render('analyses/08-Results.Rmd')"
 ```
 
diff --git a/analyses/07-Public-Dataset-Comparison.Rmd b/analyses/07-Public-Dataset-Comparison.Rmd
new file mode 100644
index 0000000..e6ba332
--- /dev/null
+++ b/analyses/07-Public-Dataset-Comparison.Rmd
@@ -0,0 +1,64 @@
+---
+title: "Public-Dataset-Comparison"
+author: "Alexander Kirchmair"
+params:
+  data:    ../data/public
+---
+
+```{r setup, include=FALSE}
+library(Seurat)
+library(DESeq2)
+library(limma)
+library(datamisc)
+library(dplyr)
+```
+
+
+Subset markers
+```{r}
+RNA <- list()
+RNAmem <- readRDS(fp("../data", "RNAmem.rds"))
+RNAexh <- readRDS(fp("../data", "RNAexh.rds"))
+
+RNA$counts <- cjoin(RNAmem$counts, RNAexh$counts)
+RNA$design <- full_join(mutate(RNAmem$design, Celltype = as.character(Celltype), Donor = as.character(Donor)),
+                        mutate(RNAexh$design, Exp = NULL, Celltype = as.character(Celltype), Donor = as.character(as.numeric(Donor)+3)))
+rownames(RNA$design) <- c(rownames(RNAmem$design), rownames(RNAexh$design))
+
+RNA$markers <- getmarkers(as.matrix(RNA$counts)[,rownames(RNA$design)], RNA$design, group = "Celltype", formula = ~ group, log2FC > 1 & padj < 0.05)
+RNA$markers <- lapply(RNA$markers, function(x){ x[1:min(length(x),200)] })
+```
+
+
+Comparison to public bulk RNA sequencing data
+```{r}
+public <- list()
+public$bulk <- read.csv("../data/public/cd8_subset_profiles.csv", row.names = 1)
+public$gsva_bulk <- runGSVA(public$bulk, genesets = RNA$markers)
+colnames(public$gsva_bulk) <- paste0(colnames(public$gsva_bulk), " (public)")
+```
+
+
+Comparison to public single-cell RNA sequencing data
+```{r}
+if (!file.exists(fp(params$data, "CD8_Tcellmap.rds"))){
+  download.file("https://singlecell.mdanderson.org/TCM/download/CD8",
+                fp(params$data, "CD8_Tcellmap.rds"), method = "wget")
+}
+tcellmap <- readRDS(fp(params$data, "CD8_Tcellmap.rds"))
+
+Idents(tcellmap) <- tcellmap$cell.type
+tcellmap <- NormalizeData(tcellmap)
+public$sc <- AverageExpression(tcellmap, group.by = "cell.type", assays = "RNA")[["RNA"]]
+colnames(public$sc) <- make.names(colnames(public$sc))
+public$sc <- public$sc[order(rowMeans(public$sc), decreasing = TRUE),]
+public$gsva_sc <- runGSVA(public$sc[1:15000,], genesets = RNA$markers)
+
+```
+
+
+```{r}
+saveRDS(public, "../data/public.rds")
+```
+
+
diff --git a/analyses/07-Results.Rmd b/analyses/08-Results.Rmd
similarity index 86%
rename from analyses/07-Results.Rmd
rename to analyses/08-Results.Rmd
index 131719b..782726e 100644
--- a/analyses/07-Results.Rmd
+++ b/analyses/08-Results.Rmd
@@ -6,6 +6,7 @@ params:
   figures:  ../results
   tables:  ../results
   format:   pdf
+  label:    FALSE
 ---
 
 ```{r setup, include=FALSE}
@@ -173,7 +174,7 @@ hm <- RNAmem$heatmap_data |>
                           name = "z-expr.",
                           right_annotation = ha,
                           row_split = rsplit,
-                          row_gap = grid::unit(7, "pt"),
+                          row_gap = grid::unit(5, "pt"),
                           row_title_rot = 0,
                           column_title_gp = grid::gpar(fontsize = titlesize),
                           row_title_gp = grid::gpar(fontsize = titlesize),
@@ -185,18 +186,6 @@ hm <- RNAmem$heatmap_data |>
 
 hm_anno_colors <- pals::polychrome(length(hm_annotation)) |> colorspace::darken(amount = 0.7) |> setNames(names(hm_annotation))
 
-hmdraw <- function(){
-  draw(hm, heatmap_legend_side = "left")
-  for(i in seq(hm_annotation)) {
-      decorate_annotation(annotation = "pathways",
-                          slice = i,
-                          code = {
-                            grid.rect(x = 0, width = unit(1, "mm"), gp = gpar(fill = hm_anno_colors[i], col = NA), just = "left")
-                            grid.text(paste(hm_annotation[[i]], collapse = "\n"), gp = gpar(col = hm_anno_colors[i], fontsize = fontsize), x = unit(5, "mm"), just = "left")
-                          })
-  }
-}
-
 E <- grid.grabExpr(hmdraw())
 
 ```
@@ -218,6 +207,9 @@ DEE
 
 Fig1 <-  wrap_plots(p, guides = "collect", tag_level = "new", design = layout)
 Fig1 <- Fig1 + plot_annotation(tag_levels = 'A',  theme = theme(plot.tag = element_text(size = 8))) & theme_fig
+if (params$label == TRUE){
+  Fig1 <- plot_grid(Fig1, textGrob("Figure 1", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 Fig1 %>% saveplot(file.path(params$figures, "Figure1"), dev = params$format, width =  full_width, height =  full_height*0.8)
 
 ```
@@ -236,22 +228,22 @@ C13cells <- subset(C13mem, Sampletype == "cells")
 Plots
 ```{r}
 
-dir.create(fp(params$figures, "isoplots"), showWarnings = FALSE)
-
-C13cells |>
-  isoplot(summarise_by = Celltype, assay = "mid_clean",
-          col_cex = 0.5, linewidth = 2, nacolor = "white", colorscale = c("white", "#1702a1"),
-          dir = fp(params$figures, "isoplots"), title = "ID", height = 10, title_size = 120, fontsize = 100,
-          dev = "pdf", label = FALSE, cumulative = TRUE, legend = FALSE)
-
-# plot legend
-leg_plot <- C13cells |>
-  isoplot(summarise_by = Celltype, assay = "mid_clean", mets = "lac", fontsize = 120, colorscale = c("white", "#1702a1"))
-leg_plot <- leg_plot + theme_basic() + theme(plot.background = element_blank(), panel.background = element_blank(), legend.background = element_blank()) +
-  guides(fill = guide_colourbar(barheight = grid::unit(200, "mm"), barwidth = grid::unit(20, "mm"), raster = FALSE, nbin = 20,
-                                ticks.colour = "black", frame.colour = "black", ticks.linewidth = 0.5, frame.linewidth = 0.5))
-cowplot::get_legend(leg_plot) |>
-  cowplot::ggsave2(filename = fp(params$figures, "isoplots", "legend.pdf"), dev = "pdf", width = 10, height = 20)
+# dir.create(fp(params$figures, "isoplots"), showWarnings = FALSE)
+# 
+# C13cells |>
+#   isoplot(summarise_by = Celltype, assay = "mid_clean",
+#           col_cex = 0.5, linewidth = 2, nacolor = "white", colorscale = c("white", "#1702a1"),
+#           dir = fp(params$figures, "isoplots"), title = "ID", height = 10, title_size = 120, fontsize = 100,
+#           dev = "pdf", label = FALSE, cumulative = TRUE, legend = FALSE)
+# 
+# # plot legend
+# leg_plot <- C13cells |>
+#   isoplot(summarise_by = Celltype, assay = "mid_clean", mets = "lac", fontsize = 120, colorscale = c("white", "#1702a1"))
+# leg_plot <- leg_plot + theme_basic() + theme(plot.background = element_blank(), panel.background = element_blank(), legend.background = element_blank()) +
+#   guides(fill = guide_colourbar(barheight = grid::unit(200, "mm"), barwidth = grid::unit(20, "mm"), raster = FALSE, nbin = 20,
+#                                 ticks.colour = "black", frame.colour = "black", ticks.linewidth = 0.5, frame.linewidth = 0.5))
+# cowplot::get_legend(leg_plot) |>
+#   cowplot::ggsave2(filename = fp(params$figures, "isoplots", "legend.pdf"), dev = "pdf", width = 10, height = 20)
 
 ```
 
@@ -376,6 +368,9 @@ p <- list(A, B, C, D, E, F, G, H, I)
 
 Fig3 <-  wrap_plots(p, guides = "collect", tag_level = "new", ncol = 3)
 Fig3 <- Fig3 + plot_annotation(tag_levels = 'A',  theme = theme(plot.tag = element_text(size = 8))) & theme_fig
+if (params$label == TRUE){
+  Fig3 <- plot_grid(Fig3, textGrob("Figure 3", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 Fig3 %>% saveplot(file.path(params$figures, "Figure3"), dev = params$format, width =  full_width, height =  full_height*0.7)
 ```
 
@@ -471,10 +466,13 @@ CE
 
 Fig4 <-  wrap_plots(p, guides = "keep", tag_level = "new", design = layout, widths = c(1,1))
 Fig4 <- Fig4 + plot_annotation(tag_levels = 'A',  theme = theme(plot.tag = element_text(size = 8))) & theme_fig
+if (params$format == "jpg"){
+  Fig4 <- plot_grid(Fig4, textGrob("Figure 4", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 Fig4 %>% saveplot(file.path(params$figures, "Figure4"), dev = params$format, width =  full_width, height =  full_height*0.6)
 
 # plot legend with color gradient separately as png:
-# cowplot::get_legend(C) |> cowplot::plot_grid() |> saveplot(file.path(params$figures, "Figure4legend"), dev = "png", width =  1000, height =  1500)
+cowplot::get_legend(C) |> cowplot::plot_grid() |> saveplot(file.path(params$figures, "Figure4legend"), dev = "png", width =  1000, height =  1500)
 ```
 
 
@@ -489,6 +487,8 @@ RNAexh <- readRDS(fp(params$data, "RNAexh.rds"))
 C13mem_cells <- readRDS(fp(params$data, "C13mem.rds")) |> subset(Sampletype == "cells" & Tracer == "13C")
 C13exh_cells <- readRDS(fp(params$data, "C13exh.rds")) |> subset(Sampletype == "cells" & Tracer == "13C")
 
+public <- readRDS(fp(params$data, "public.rds"))
+
 genesets <- read.delim("../tables/genesets.tsv.gz")
 markers <- readLines("../tables/markers.txt")
 
@@ -508,10 +508,39 @@ expr_mean <- log2(expr + 1) |> summarise_cols(coldata = df, by = Celltype, FUN =
 
 ### Supplementary Figure 2: Differentiation markers.
 ```{r}
-FigS2 <- expr_mean[intersect(markers, rownames(expr)), cols[1:4]] |>
+A <- public$gsva_bulk |> cxheatmap(scale = "row") |> draw() |> grid.grabExpr()
+B <- public$gsva_sc |> cxheatmap(scale = "row") |> draw() |> grid.grabExpr()
+C <- expr_mean[intersect(markers, rownames(expr)), cols[1:4]] |>
   cxheatmap(scale = "", title = "z-expr.", cluster_cols = FALSE,
-            fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
-FigS2 %>% saveplot(file.path(params$figures, "FigureS2"), dev = params$format, width = full_width, height =  full_height*0.8)
+            fontsize = 18, rowcex = 0.7, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE) |> draw() |> grid.grabExpr()
+
+df <- (RNAmem$ora_cluster[c("UDDD","DUDD","DDUD","DDDU")] %L>% mutate(., cluster=.name)) |> Reduce(f=rbind)
+top <- RNAmem$ora_cluster[c("UDDD","DUDD","DDUD","DDDU")] %L>% subset(., padj <= 0.1) %L>% rownames()
+pw <- (top %L>% function(x) x[1:min(length(x),3)]) |> unlist() |> unique()
+df <- subset(df, ID %in% pw)
+df$cluster[df$cluster == "UDDD"] <- "TN"
+df$cluster[df$cluster == "DUDD"] <- "TSCM"
+df$cluster[df$cluster == "DDUD"] <- "TCM"
+df$cluster[df$cluster == "DDDU"] <- "TEM"
+df$cluster <- factor(df$cluster, ordered = TRUE, levels = c("TN","TSCM","TCM","TEM"))
+D <- ggplot(df, aes(cluster, ID, size=Count, color=-log10(padj))) +
+  theme_basic(base_size = 15) +
+  theme(panel.grid.major = element_line(color = "#dddddd")) +
+  geom_point() + xlab("") + ylab("")
+
+
+layout <- "
+AC
+BC
+DD
+"
+
+FigS2 <-  wrap_plots(list(A,B,C,D), guides = "keep", tag_level = "new", design = layout, heights = c(1,1,0.7))
+FigS2 <- FigS2 + plot_annotation(tag_levels = 'A',  theme = theme(plot.tag = element_text(size = 8))) & theme_fig
+if (params$label == TRUE){
+  FigS2 <- plot_grid(FigS2, textGrob("Figure S2", x = unit(0.5, "npc"), y = unit(0.55, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
+FigS2 %>% saveplot(file.path(params$figures, "FigureS2"), dev = params$format, width =  full_width, height =  full_height*0.9)
 
 ```
 
@@ -521,6 +550,11 @@ FigS2 %>% saveplot(file.path(params$figures, "FigureS2"), dev = params$format, w
 genes <- subset(genesets, term == "MetabolicAtlas glycolysis_gluconeogenesis")$gene
 FigS3 <- expr_mean[genes, cols[1:4]] |> cxheatmap(column_title = "MetabolicAtlas Glycolysis/Gluconeogenesis", title = "z-expr.", cluster_cols = FALSE,
                                                   fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
+
+if (params$label == TRUE){
+  FigS3 <- FigS3 |> draw() |> grid.grabExpr()
+  FigS3 <- plot_grid(FigS3, textGrob("Figure S3", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS3 %>% saveplot(file.path(params$figures, "FigureS3"), dev = params$format, width = full_width, height =  full_height*0.8)
 
 ```
@@ -586,6 +620,9 @@ EF
 
 FigS4 <-  wrap_plots(p, guides = "keep", tag_level = "new", design = layout, widths = c(1,1))
 FigS4 <- FigS4 + plot_annotation(tag_levels = 'A') & th
+if (params$format == "jpg"){
+  FigS4 <- plot_grid(FigS4, textGrob("Figure S4", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS4 %>% saveplot(file.path(params$figures, "FigureS4"), dev = params$format, width =  full_width, height =  full_height)
 
 ```
@@ -595,7 +632,12 @@ FigS4 %>% saveplot(file.path(params$figures, "FigureS4"), dev = params$format, w
 ```{r}
 genes <- genesets |> subset(term %in% c("MitoCarta tca cycle")) |> pull(gene)
 FigS5 <- expr_mean[genes, cols[1:4]] |> cxheatmap(column_title = "MitoCarta TCA cycle", title = "z-expr.", cluster_cols = FALSE,
-                                                  fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
+                                                  fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE) 
+
+if (params$label == TRUE){
+  FigS5 <- FigS5 |> draw() |> grid.grabExpr()
+  FigS5 <- plot_grid(FigS5, textGrob("Figure S5", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS5 %>% saveplot(file.path(params$figures, "FigureS5"), dev = params$format, width = full_width, height =  full_height*0.8)
 
 ```
@@ -614,10 +656,15 @@ genes <- list(
 df <- stack(genes) |> dplyr::rename(genes = values, shuttle = ind)
 FigS6 <- expr_mean[unlist(genes), cols[1:4]] |> cxheatmap(column_title = "NADH shuttles", title = "z-expr.", row_split = df$shuttle,
                                                           cluster_rows = FALSE, cluster_cols = FALSE,
-                                                          row_title_gp = grid::gpar(fontsize = 18),
+                                                          row_title_gp = grid::gpar(fontsize = 15),
                                                           row_gap = grid::unit(7, "pt"),
                                                           fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
-FigS6 %>% saveplot(fp(params$figures, "FigureS6"), dev = params$format, width = full_width, height =  full_height*0.8)
+
+if (params$label == TRUE){
+  FigS6 <- FigS6 |> draw() |> grid.grabExpr()
+  FigS6 <- plot_grid(FigS6, textGrob("Figure S6", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
+FigS6 %>% saveplot(fp(params$figures, "FigureS6"), dev = params$format, width = full_width, height =  full_height*0.9)
 
 ```
 
@@ -629,6 +676,11 @@ exp <- c(rep("differentiation", 4), rep("exhaustion", 2))
 FigS7 <- expr_mean[intersect(genes, rownames(expr_mean)), cols] |>
   cxheatmap(column_title = "GO L-amino acid transmembrane transporter activity", column_split = exp, title = "z-expr.", cluster_cols = FALSE,
             fontsize = 18, rowcex = 0.8, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
+
+if (params$label == TRUE){
+  FigS7 <- FigS7 |> draw() |> grid.grabExpr()
+  FigS7 <- plot_grid(FigS7, textGrob("Figure S7", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS7 %>% saveplot(fp(params$figures, "FigureS7"), dev = params$format, width = full_width, height =  full_height*0.8)
 
 ```
@@ -640,6 +692,11 @@ genes <- genesets |> subset(term %in% c("MitoCarta oxphos subunits")) |> pull(ge
 FigS8 <- expr_mean[intersect(genes, rownames(expr_mean)), cols[1:4]] |>
   cxheatmap(column_title = "OXPHOS subunits", title = "z-expr.", cluster_cols = FALSE,
             fontsize = 18, rowcex = 0.6, colcex = 1.2, column_names_rot = 0, column_names_centered = TRUE)
+
+if (params$label == TRUE){
+  FigS8 <- FigS8 |> draw() |> grid.grabExpr()
+  FigS8 <- plot_grid(FigS8, textGrob("Figure S8", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS8 %>% saveplot(fp(params$figures, "FigureS8"), dev = params$format, width = full_width, height =  full_height*0.8)
 
 ```
@@ -683,11 +740,15 @@ th <- theme(axis.text = txt,
 
 FigS9 <-  wrap_plots(A, wrap_plots(B, plot_spacer(), nrow = 1, widths = c(1, 0.3)), guides = "keep", tag_level = "new", ncol = 1)
 FigS9 <- FigS9 + plot_annotation(tag_levels = 'A') & th
+if (params$label == TRUE){
+  FigS9 <- plot_grid(FigS9, textGrob("Figure S9", x = unit(0.5, "npc"), y = unit(0.5, "npc")), ncol = 1, rel_heights = c(0.95,0.05))
+}
 FigS9 %>% saveplot(fp(params$figures, "FigureS9"), dev = params$format, width =  full_width, height =  full_height*0.6)
 ```
 
 
 
+
 # TABLES ----------------------------------------------------------------------------------------------------------------------------------------------------
 
 Data
diff --git a/env/cd8.yml b/env/cd8.yml
index 22d6392..6c9e3ef 100644
--- a/env/cd8.yml
+++ b/env/cd8.yml
@@ -24,15 +24,18 @@ dependencies:
   - bioconductor-annotate=1.72.0=r41hdfd78af_0
   - bioconductor-annotationdbi=1.56.2=r41hdfd78af_0
   - bioconductor-apeglm=1.16.0=r41hc247a5b_2
+  - bioconductor-beachmat=2.10.0=r41hc247a5b_2
   - bioconductor-biobase=2.54.0=r41hc0cfd56_2
   - bioconductor-biocfilecache=2.2.0=r41hdfd78af_0
   - bioconductor-biocgenerics=0.40.0=r41hdfd78af_0
   - bioconductor-biocparallel=1.28.3=r41hc247a5b_1
+  - bioconductor-biocsingular=1.10.0=r41hc247a5b_2
   - bioconductor-biomart=2.50.0=r41hdfd78af_0
   - bioconductor-biostrings=2.62.0=r41hc0cfd56_2
   - bioconductor-clusterprofiler=4.2.0=r41hdfd78af_0
   - bioconductor-complexheatmap=2.10.0=r41hdfd78af_0
   - bioconductor-delayedarray=0.20.0=r41hc0cfd56_2
+  - bioconductor-delayedmatrixstats=1.16.0=r41hdfd78af_0
   - bioconductor-deseq2=1.34.0=r41hc247a5b_3
   - bioconductor-do.db=2.9=r41hdfd78af_13
   - bioconductor-dose=3.20.0=r41hdfd78af_0
@@ -47,6 +50,10 @@ dependencies:
   - bioconductor-ggtree=3.2.0=r41hdfd78af_0
   - bioconductor-go.db=3.14.0=r41hdfd78af_1
   - bioconductor-gosemsim=2.20.0=r41hc247a5b_2
+  - bioconductor-graph=1.72.0=r41hc0cfd56_2
+  - bioconductor-gseabase=1.56.0=r41hdfd78af_0
+  - bioconductor-gsva=1.42.0=r41hc0cfd56_2
+  - bioconductor-hdf5array=1.22.1=r41hc0cfd56_1
   - bioconductor-ihw=1.22.0=r41hdfd78af_0
   - bioconductor-iranges=2.28.0=r41hc0cfd56_2
   - bioconductor-isocorrector=1.12.0=r41hdfd78af_0
@@ -56,7 +63,13 @@ dependencies:
   - bioconductor-matrixgenerics=1.6.0=r41hdfd78af_0
   - bioconductor-org.hs.eg.db=3.14.0=r41hdfd78af_1
   - bioconductor-qvalue=2.26.0=r41hdfd78af_0
+  - bioconductor-rhdf5=2.38.1=r41hbe1951d_0
+  - bioconductor-rhdf5filters=1.6.0=r41hc247a5b_2
+  - bioconductor-rhdf5lib=1.16.0=r41hc0cfd56_2
   - bioconductor-s4vectors=0.32.4=r41hc0cfd56_0
+  - bioconductor-scaledmatrix=1.2.0=r41hdfd78af_0
+  - bioconductor-singlecellexperiment=1.16.0=r41hdfd78af_0
+  - bioconductor-sparsematrixstats=1.6.0=r41hc247a5b_2
   - bioconductor-summarizedexperiment=1.24.0=r41hdfd78af_0
   - bioconductor-tcgabiolinks=2.22.1=r41hdfd78af_0
   - bioconductor-tcgabiolinksgui.data=1.14.1=r41hdfd78af_0
@@ -78,7 +91,7 @@ dependencies:
   - charset-normalizer=3.2.0=pyhd8ed1ab_0
   - comm=0.1.3=pyhd8ed1ab_0
   - coreutils=9.3=h0b41bf4_0
-  - curl=7.86.0=h2283fc2_1
+  - curl=7.86.0=h7bff187_1
   - cutadapt=4.0=py39hbf8eff0_0
   - debugpy=1.6.7=py39h227be39_0
   - decorator=5.1.1=pyhd8ed1ab_0
@@ -136,12 +149,12 @@ dependencies:
   - jupyterlab_server=2.24.0=pyhd8ed1ab_0
   - kernel-headers_linux-64=2.6.32=he073ed8_16
   - keyutils=1.6.1=h166bdaf_0
-  - krb5=1.19.3=h08a2579_0
+  - krb5=1.19.3=h3790be6_0
   - ld_impl_linux-64=2.40=h41732ed_0
   - lerc=4.0.0=h27087fc_0
   - libblas=3.9.0=17_linux64_openblas
   - libcblas=3.9.0=17_linux64_openblas
-  - libcurl=7.86.0=h2283fc2_1
+  - libcurl=7.86.0=h7bff187_1
   - libdeflate=1.17=h0b41bf4_0
   - libedit=3.1.20191231=he28a2e2_2
   - libev=4.33=h516909a_1
@@ -151,20 +164,20 @@ dependencies:
   - libgcc-ng=13.1.0=he5830b7_0
   - libgfortran-ng=13.1.0=h69a702a_0
   - libgfortran5=13.1.0=h15d22d2_0
-  - libgit2=1.6.4=h747ad27_0
+  - libgit2=1.5.1=ha98c156_0
   - libglib=2.76.4=hebfc3b9_0
   - libglu=9.0.0=he1b5a44_1001
   - libgomp=13.1.0=he5830b7_0
   - libiconv=1.17=h166bdaf_0
   - liblapack=3.9.0=17_linux64_openblas
-  - libnghttp2=1.52.0=h61bc06f_0
+  - libnghttp2=1.51.0=hdcd2b5c_0
   - libnsl=2.0.0=h7f98852_0
   - libopenblas=0.3.23=pthreads_h80387f5_0
   - libpng=1.6.39=h753d276_0
   - libsanitizer=13.1.0=hfd8a6a1_0
   - libsodium=1.0.18=h36c2ea0_1
   - libsqlite=3.42.0=h2797004_0
-  - libssh2=1.11.0=h0841786_0
+  - libssh2=1.10.0=haa6b8db_3
   - libstdcxx-devel_linux-64=13.1.0=he3cc6c4_0
   - libstdcxx-ng=13.1.0=hfd8a6a1_0
   - libtiff=4.5.0=h6adf6a1_2
@@ -188,7 +201,7 @@ dependencies:
   - notebook-shim=0.2.3=pyhd8ed1ab_0
   - numpy=1.25.1=py39h6183b62_0
   - openjdk=11.0.1=h516909a_1016
-  - openssl=3.1.2=hd590300_0
+  - openssl=1.1.1v=hd590300_0
   - overrides=7.3.1=pyhd8ed1ab_0
   - packaging=23.1=pyhd8ed1ab_0
   - pandoc=2.19.2=h32600fe_2
@@ -214,7 +227,7 @@ dependencies:
   - pycparser=2.21=pyhd8ed1ab_0
   - pygments=2.15.1=pyhd8ed1ab_0
   - pysocks=1.7.1=pyha2e5f31_6
-  - python=3.9.16=h2782a2a_0_cpython
+  - python=3.9.15=h47a2c10_0_cpython
   - python-dateutil=2.8.2=pyhd8ed1ab_0
   - python-fastjsonschema=2.18.0=pyhd8ed1ab_0
   - python-isal=1.2.0=py39hd1e30aa_0
@@ -318,7 +331,7 @@ dependencies:
   - r-future.apply=1.11.0=r41hc72bb7e_0
   - r-gargle=1.4.0=r41h785f33e_0
   - r-generics=0.1.3=r41hc72bb7e_1
-  - r-gert=1.9.2=r41h492061e_1
+  - r-gert=1.9.2=r41hf3f2ec2_0
   - r-getoptlong=1.0.5=r41hc72bb7e_1
   - r-ggforce=0.4.1=r41h7525677_1
   - r-ggfun=0.0.9=r41hc72bb7e_0
@@ -414,7 +427,7 @@ dependencies:
   - r-nloptr=2.0.3=r41hb13c81a_1
   - r-nnet=7.3_19=r41h57805ef_0
   - r-numderiv=2016.8_1.1=r41hc72bb7e_4
-  - r-openssl=2.0.6=r41habfbb5e_0
+  - r-openssl=2.0.5=r41hb1dc35e_0
   - r-openxlsx=4.2.5.2=r41h38f115c_0
   - r-pals=1.7=r41hc72bb7e_1
   - r-parallelly=1.36.0=r41hc72bb7e_0
@@ -486,6 +499,7 @@ dependencies:
   - r-rsqlite=2.3.1=r41h38f115c_0
   - r-rstatix=0.7.2=r41hc72bb7e_0
   - r-rstudioapi=0.14=r41hc72bb7e_1
+  - r-rsvd=1.0.5=r41hc72bb7e_0
   - r-rvcheck=0.1.8=r41hc72bb7e_1
   - r-rversions=2.1.2=r41hc72bb7e_1
   - r-rvest=1.0.3=r41hc72bb7e_1
diff --git a/lib/R-functions.R b/lib/R-functions.R
index eac3852..5703428 100644
--- a/lib/R-functions.R
+++ b/lib/R-functions.R
@@ -258,3 +258,27 @@ getLimits <- function (x, clip = TRUE, expand = 1, negative = TRUE){
 }
 
 
+getmarkers <- function(data, design, group = "Celltype", formula = ~ group, ...){
+  ct <- unique(design[[group]])
+  markers <- lapply(setNames(ct, ct), function(ctx){
+    design$group <- ifelse(design[[group]] == ctx, ctx, "other")
+    res <- runDESeq2(data, design, formula = formula, contrasts = list(de = c("group", ctx, "other")))
+    subset(res$results[[1]], ...)$gene
+  })
+  markers
+}
+
+
+hmdraw <- function(){
+  draw(hm, heatmap_legend_side = "left")
+  for(i in seq(hm_annotation)) {
+    decorate_annotation(annotation = "pathways",
+                        slice = i,
+                        code = {
+                          grid.rect(x = 0, width = unit(1, "mm"), gp = gpar(fill = hm_anno_colors[i], col = NA), just = "left")
+                          grid.text(paste(hm_annotation[[i]], collapse = "\n"), gp = gpar(col = hm_anno_colors[i], fontsize = fontsize), x = unit(5, "mm"), just = "left")
+                        })
+  }
+}
+
+
diff --git a/lib/make_env.sh b/lib/make_env.sh
index c3be58b..e0032e9 100644
--- a/lib/make_env.sh
+++ b/lib/make_env.sh
@@ -2,7 +2,7 @@
 # Main conda env
 envname="cd8"
 installconda="nextflow=21.10.0 cutadapt=4.0 r-base=4.1.3 r-essentials \
-              bioconductor-isocorrector bioconductor-org.hs.eg.db bioconductor-clusterprofiler bioconductor-deseq2 bioconductor-ihw \
+              bioconductor-isocorrector bioconductor-org.hs.eg.db bioconductor-clusterprofiler bioconductor-deseq2 bioconductor-ihw bioconductor-gsva \
               bioconductor-geoquery bioconductor-complexheatmap bioconductor-biomart bioconductor-limma bioconductor-tximport bioconductor-tcgabiolinks \
               r-nlme r-glmmtmb r-contrast r-multcomp r-betareg r-msigdbr r-r.utils r-dendsort r-pals \
               r-ggpubr r-dplyr r-tibble r-openxlsx r-ggplot2 r-ggrepel r-devtools r-crayon"