fiumara_20241115_UploadScript.R

### upload on gitlab

# 1) load R packages
library(Seurat)
library(GEOquery)
library(limma)
library(dplyr)
library(data.table)
library(clusterProfiler)
library(ggplot2)
library(RColorBrewer)
library(qusage)
library(UCell)

# 2) define input and output paths
# set path to the working directory
path <- ""

# 3) load scRNA-seq BM-object

## 3.1) restrict analysis to CD34+ clusters
obj <- readRDS(file = paste0(path, "BM_20240222.rds"))
CD34p_clusters <- as.character(c(2, 8, 14, 15, 16, 19, 20, 26, 27, 30, 33))

obj$seurat <- obj$RNA_snn_h.orig.ident_res.1.8
obj <- SetIdent(object = obj, value = "seurat")

obj <- subset(obj, subset = (seurat %in% CD34p_clusters))
saveRDS(object = obj, file = paste0(path, "BM_CD34p.rds"))

## 3.2) free memory space
l <- ls()
l <- l[!l %in% "path"]
rm(list = l); gc()



# 4) load published scRNA-seq dataset

## 4.1) Ainciburu et al.

### 4.1.1) download tar file and metadata
destfolder <- paste0(path, "GSE180298/"); dir.create(path = destfolder, showWarnings = F, recursive = T)
destfile <- paste0(destfolder, "GSE180298_RAW.tar")
download.file(url = "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE180298&format=file", destfile = destfile)

metadata_files <- c("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE180298&format=file&file=GSE180298%5Felderly%5Fmetadata%2Etxt%2Egz",
                    "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE180298&format=file&file=GSE180298%5Fmds%5Fmetadata%2Etxt%2Egz",
                    "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE180298&format=file&file=GSE180298%5Fyoung%5Fmetadata%2Etxt%2Egz")
download.file(url = metadata_files, destfile = paste0(destfolder, "GSE180298_", c("elderly", "mds", "young"), "_metadata.txt.gz"))

### 4.1.2) untar file
untar(tarfile = destfile, exdir = paste0(destfolder, "GSE180298_RAW/"))
files <- list.files(path = paste0(destfolder, "GSE180298_RAW/"), full.names = T)
md <- list.files(path = destfolder, full.names = T)
md_files <- md[grepl(md, pattern = "metadata")]

### 4.1.3) define seurat object
metadata <- c()
for(i in seq_len(length(md_files))){
  gunzip(md_files[i], remove = FALSE, overwrite = TRUE)
  unzfile <- gsub(md_files[i], pattern = ".gz", replacement = "")
  temp <- read.delim(file = unzfile, header = T, sep = "\t")
  metadata <- rbind(metadata, temp)
}
metadata <- as.data.frame(metadata)
metadata$orig.ident <- strsplit2(rownames(metadata), split = "\\_")[, 2]

counts <- cell_counts <- samples <- c()
for(i in seq_len(length(files))){
  
  split_filename <- strsplit2(files[i], split = "/")
  sampleid <- strsplit2(split_filename[, ncol(split_filename)], split = "\\_")[, 2]
  print(paste0("Processing file ", sampleid, " (", i, " over ", length(files), ")"))
  
  # load h5 object
  temp <- Read10X_h5(filename = files[i], use.names = TRUE, unique.features = TRUE)
  coln <- paste0(strsplit2(colnames(temp), split = "-")[, 1], "_", sampleid)
  colnames(temp) <- coln
  
  # store
  temp_object <- CreateSeuratObject(counts = temp)
  counts_temp <- matrix(data = temp_object@assays$RNA@counts, ncol = ncol(temp_object))
  rownames(counts_temp) <- rownames(temp_object)
  colnames(counts_temp) <- colnames(temp_object)
  
  # store info
  cell_counts <- c(cell_counts, ncol(counts_temp))
  samples <- c(samples, sampleid)
  
  if(is.null(nrow(counts))){
    counts <- counts_temp
    gs <- rownames(counts)
  }else{
    gs <- intersect(gs, rownames(counts_temp))
    counts <- cbind(counts[match(gs, rownames(counts)), ], 
                    counts_temp[match(gs, rownames(counts_temp)), ])
  }
}
names(cell_counts) <- samples

md <- metadata[rownames(metadata) %in% colnames(counts),]; dim(md)
m <- match(rownames(md), colnames(counts)); table(is.na(m))
counts <- counts[, m]
GSE180298 <- CreateSeuratObject(counts = counts, meta.data = metadata)

### 4.1.4) select only elderly donors
GSE180298 <- SetIdent(object = GSE180298, value = "orig.ident")
elderly_sampleid <- unique(GSE180298$orig.ident)[grepl(unique(GSE180298$orig.ident), pattern = "elderly")]
GSE180298 <- subset(GSE180298, subset = (orig.ident %in% elderly_sampleid))
saveRDS(object = GSE180298, file = paste0(path, "GSE180298_elderly.rds")); gc()

## 4.1.5) free memory space
l <- ls()
l <- l[!l %in% "path"]
rm(list = l); gc()


## 4.2) Wu et al.

### 4.2.1) download tar file
destfolder <- paste0(path, "GSE196052/"); dir.create(path = destfolder, showWarnings = F, recursive = T)
destfile <- paste0(destfolder, "GSE196052_RAW.tar")
download.file(url = "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE196052&format=file", destfile = destfile)

metadata_files <- c("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE196052&format=file&file=GSE196052%5FdataCount%5FCD34%2Ecsv%2Egz")
download.file(url = metadata_files, destfile = paste0(destfolder, "GSE196052_dataCount_CD34.csv.gz"))

### 4.2.2) untar file
untar(tarfile = destfile, exdir = paste0(destfolder, "GSE196052_RAW/"))
files_annotation <- list.files(path = paste0(destfolder, "GSE196052_RAW/"), full.names = T)
files <- list.files(path = destfolder, full.names = T)

### 4.2.3) define seurat object
md <- read.delim(file = files[grepl(files, pattern = "SraRun")], header = T, sep = ",")
table(md$subject_id.status, md$tissue.cell_type)
mdfiles <- files_annotation[grepl(files_annotation, pattern = "GSM")]
metadata <- c()
for(i in seq_len(length(mdfiles))){
  x <- mdfiles[i]
  temp <- read.delim(gzfile(x), header = T, sep = ",") 
  metadata <- rbind(metadata, temp)
}

f <- files[grepl(files, pattern = "CD34.csv.gz")]
counts_cd34 <- fread(f) %>% as.data.frame()
rn <- as.character(counts_cd34[, 1])
counts_cd34 <- counts_cd34[, -1]
rownames(counts_cd34) <- rn

cd34_cells <- colnames(counts_cd34)
cd34_cells <- gsub(x = cd34_cells, pattern = "CD34_", replacement = "")
wu_annot <- read.csv(file = paste0(destfolder, "CD34_metaDatatSNECellType_ALiceManual.csv"), header = T)
m <- match(cd34_cells, wu_annot$orig.ident); table(is.na(m))
status <- wu_annot$group[m]
sampleid <- wu_annot$subject[m]
mdt <- data.frame(sampleid, status)
rownames(mdt) <- colnames(counts_cd34)

GSE196052 <- CreateSeuratObject(counts = counts_cd34, meta.data = mdt)
GSE196052 <- SetIdent(object = GSE196052, value = "orig.ident")
saveRDS(object = GSE196052, file = paste0(path, "GSE196052.rds"))

## 4.2.4) free memory space
l <- ls()
l <- l[!l %in% "path"]
rm(list = l); gc()




# 5) Define single seurat object

## 5.1) load datasets
fiumara <- readRDS(paste0(path, "BM_CD34p.rds"))
GSE180298 <- readRDS(paste0(path, "GSE180298_elderly.rds"))
GSE196052 <- readRDS(paste0(path, "GSE196052.rds"))

## 5.2) define common features
genes <- table(c(rownames(fiumara), rownames(GSE180298), rownames(GSE196052)))
common_genes <- names(genes)[genes == 3]; length(common_genes)

## 5.3) add sampleinfo
fiumara$source <- "fiumara"
fiumara$sample_info <- paste0("fiumara_", fiumara$orig.ident)
GSE180298$source <- "GSE180298"
GSE180298$sample_info <- paste0("GSE180298_", GSE180298$orig.ident)
GSE196052$source <- "GSE196052"
GSE196052$sample_info <- paste0("GSE196052_", GSE196052$sampleid)

## 5.4) extract counts relative to these genes and define Seurat object
counts <- cbind(fiumara@assays$RNA@counts[common_genes,],
                GSE180298@assays$RNA@counts[common_genes,],
                GSE196052@assays$RNA@counts[common_genes,])

vars <- c("source", "sample_info", "nCount_RNA", "nFeature_RNA")
md <- rbind(fiumara@meta.data[, vars],
            GSE180298@meta.data[, vars],
            GSE196052@meta.data[, vars])

obj <- CreateSeuratObject(counts = counts, meta.data = md)
obj$percent.mt <- (colSums(obj@assays$RNA@counts[grepl(rownames(obj), pattern = "^MT-"),])/colSums(obj@assays$RNA@counts))*100
obj$percent.rb <- (colSums(obj@assays$RNA@counts[grepl(rownames(obj), pattern = "^RPL"),])/colSums(obj@assays$RNA@counts))*100
saveRDS(object = obj, file = paste0(path, "combined.rds")); gc()

## 5.5) free memory space
l <- ls()
l <- l[!l %in% c("path", "obj")]
rm(list = l); gc()



# 6) removing low quality cells, normalization, scaling, integration and clustering

## 6.1) cell filtering
obj$keep <- (obj$nFeature_RNA > 200) & (obj$percent.mt < 25)
table(obj$keep, obj$source)
obj <- subset(obj, subset = (keep %in% TRUE))
saveRDS(obj, file = paste0(path, "combined_filtered.rds"))

## 6.2) normalization
varfeatures <- 1000
obj <- NormalizeData(object = obj)
obj <- FindVariableFeatures(obj, selection.method = "vst",
                            nfeatures = varfeatures,
                            verbose = T)

## 6.3) scaling
reg_vars = c("percent.mt", "nCount_RNA")
obj <- ScaleData(object = obj, vars.to.regress = reg_vars,  display.progress = T, features = rownames(obj))
saveRDS(obj, file = paste0(path, "combined_normscaled.rds"))

## 6.4) dimensionality reduction
max_pca <- 100
obj <- RunPCA(object = obj, features = VariableFeatures(object = obj), npcs = max_pca, reduction.name="pca", reduction.key="PC_")

explvar <- ((obj@reductions$pca@stdev^2)/sum((obj@reductions$pca@stdev^2)))*100
delta <- explvar - c(explvar[-1], 0)
opt_delta <- length(delta[delta > 1e-2])
opt_explvar <- min(which(cumsum(explvar) > 80))
opt <- min(opt_delta, opt_explvar)

obj <- RunUMAP(object = obj, seed.use = 123, reduction = "pca", dims = 1:opt)

## 6.5) harmony dataset integration
integration.var <- c("source", "sample_info")
obj <- RunHarmony(object = obj, 
                  group.by.vars = integration.var, 
                  max.iter.harmony = 30, 
                  plot_convergence = FALSE, 
                  reduction.save = "harmony")
obj <- RunUMAP(object = obj, 
               seed.use = 123, 
               dims = 1:opt, 
               reduction = "harmony", 
               reduction.name = "harmony_umap", 
               reduction.key = "UMAPh_", 
               return.model = TRUE)

## 6.6) find neighboring cells
obj <- FindNeighbors(object = obj, 
                     dims = 1:opt, 
                     force.recalc = T, 
                     reduction = "harmony", 
                     graph.name = c("RNA_nn_h.iv", "RNA_snn_h.iv"))

## 6.7) cell clustering
clu_res <- seq(0.1, 1, by = 0.1)
for(res in clu_res){
  obj <- FindClusters(object = obj, 
                      algorithm = 1, 
                      resolution = as.numeric(res), 
                      graph.name = "RNA_snn_h.iv")
}

## 6.8) assign celltype to each cluster
obj$seurat <- obj$RNA_snn_h.iv_res.0.5
obj$seurat_annotation <- NA
obj$seurat_annotation[obj$seurat %in% "0"] <- "HSC/MPP"
obj$seurat_annotation[obj$seurat %in% c("1", "6", "11")] <- "Mature Ery"
obj$seurat_annotation[obj$seurat %in% "7"] <- "VEXAS Ery/CMP"
obj$seurat_annotation[obj$seurat %in% c("8", "22")] <- "Immature Ery"
obj$seurat_annotation[obj$seurat %in% "2"] <- "MPP/CMP"
obj$seurat_annotation[obj$seurat %in% "12"] <- "PreBNK"
obj$seurat_annotation[obj$seurat %in% "3"] <- "CMP/GMP"
obj$seurat_annotation[obj$seurat %in% "4"] <- "GP"
obj$seurat_annotation[obj$seurat %in% "13"] <- "MDP"
obj$seurat_annotation[obj$seurat %in% "5"] <- "MyeloLympho/CMP"
obj$seurat_annotation[obj$seurat %in% c("9", "21")] <- "VEXAS Immature Ery"
obj$seurat_annotation[obj$seurat %in% "10"] <- "Undefined"
obj$seurat_annotation[obj$seurat %in% "14"] <- "BEM"
obj$seurat_annotation[obj$seurat %in% "16"] <- "Monocyte Progenitors"
obj$seurat_annotation[obj$seurat %in% c("15", "19")] <- "MLP"
obj$seurat_annotation[obj$seurat %in% "17"] <- "MEP"
obj$seurat_annotation[obj$seurat %in% "18"] <- "VEXAS Mature Ery"
obj$seurat_annotation[obj$seurat %in% "20"] <- "VEXAS MPP-Ery"
obj$seurat_annotation[obj$seurat %in% "23"] <- "VEXAS Myelo/CMP"

annotated_cols <- c("HSC/MPP" = '#FB0207',
                    "MyeloLympho/CMP" = '#c6dbef',
                    "MPP/CMP" = "#7fcdbb",
                    "CMP/GMP" = '#9ecae1',
                    "Monocyte Progenitors" = '#66CCFF',
                    "MDP" = '#0F80FF',
                    "GP" = "#08519c",
                    "PreBNK" = '#118040',
                    "MLP" = "#FECC66",
                    "BEM" = "#bdbdbd",
                    "MEP" = '#f768a1',
                    "Immature Ery" = '#B17DFC',
                    "Mature Ery" = "#800080",
                    "VEXAS MPP-Ery" = '#fde0dd',
                    "VEXAS Ery/CMP" = '#fcc5c0',
                    "VEXAS Immature Ery" = '#fa9fb5',
                    "VEXAS Mature Ery" = "#dd3497",
                    "VEXAS Myelo/CMP" = "#ccebc5",
                    "Undefined" = '#d9d9d9')
levs <- names(annotated_cols)
obj$seurat_annotation <- factor(obj$seurat_annotation, levels = levs)

## 6.9) define status and vexas-mutation variables
obj$cell_barcode <- strsplit2(rownames(obj@meta.data), split = "\\_")[, 2] 
GSE196052_annot <- read.csv(file = paste0(path, "GSE196052/CD34_metaDatatSNECellType_ALiceManual.csv"), header = T)
GSE196052_cases <- GSE196052_annot[GSE196052_annot$group %in% "PT",]

### 6.9.1) status
obj$status <- "HD"
obj$status[(obj$source %in% "GSE196052") & (obj$cell_barcode %in% GSE196052_cases$orig.ident)] <- "PT"
obj$status[(obj$source %in% "fiumara")] <- "PT"
table(obj$source, obj$status)

### 6.9.2) VEXAS mutation
GSE196052_pt2upn <- paste0("GSE196052_PT", 1:9)
names(GSE196052_pt2upn) <- paste0("GSE196052_UPN", c(6, 11, 1, 10, 13, 14, 15, 16, 17))
obj$sample_info_upn <- obj$sample_info
for(i in seq_len(length(GSE196052_pt2upn))){
  obj$sample_info_upn[obj$sample_info_upn %in% GSE196052_pt2upn[i]] <- names(GSE196052_pt2upn)[i]
}
table(obj$source, obj$sample_info_upn)
table(obj$sample_info, obj$sample_info_upn)

obj$vexas_mutation <- "HD"
obj$vexas_mutation[obj$sample_info_upn %in% c(paste0("fiumara_BM-0", 1),
                                              paste0("GSE196052_UPN", c(1, 10, 11, 13, 16, 17)))] <- "THR"
obj$vexas_mutation[obj$sample_info_upn %in% c(paste0("fiumara_BM-0", c(2, 3, 8)),
                                              paste0("GSE196052_UPN", c(6)))] <- "VAL"
obj$vexas_mutation[obj$sample_info_upn %in% c(paste0("fiumara_BM-0", c(4, 9)),
                                              paste0("GSE196052_UPN", c(14, 15)))] <- "LEU"
saveRDS(obj, file = paste0(path, "combined_annotated.rds"))

table(obj$status, obj$vexas_mutation)
table(obj$source, obj$vexas_mutation)
table(obj$vexas_mutation, obj$sample_info_upn)

## 6.10) free memory space
l <- ls()
l <- l[!l %in% c("path", "obj")]
rm(list = l); gc()



# 7) Celltype-wise VEXAS vs HD (DE and GSEA)
outpath <- paste0(path, "DE_GSEA/"); dir.create(path = outpath, showWarnings = F, recursive = T)

## 7.1) define function to run clusterProfiler GSEA
gsea_run <- function(marks, gmt){
  
  # load gmt file
  gmt.obj <- clusterProfiler::read.gmt(gmt)
  
  # order DE results by logFC
  genes <- marks$avg_log2FC
  names(genes) <- marks$gene_name
  genes <- genes[order(genes, decreasing = T)]
  genes <- genes[!duplicated(names(genes))]
  
  # run GSEA
  gsea <- GSEA(geneList = genes, TERM2GENE = gmt.obj, nPerm = 10000, pvalueCutoff = 1)
  return(gsea)
}

## 7.2) download hallmarks gene sets
hallmarks_gsea <- c("https://data.broadinstitute.org/gsea-msigdb/msigdb/release/7.4/h.all.v7.4.symbols.gmt")
download.file(url = hallmarks_gsea, destfile = paste0(outpath, "h.all.v7.4.symbols.gmt"))

## 7.3) VEXAS vs HD

### 7.3.1) DE analysis
annclusters <- names(annotated_cols)
mincells <- 10
for(i in seq_len(length(annclusters))){
  cl <- annclusters[i]
  cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
  temp <- subset(obj, subset = (seurat_annotation %in% cl))
  temp <- SetIdent(object = temp, value = "status")
  
  if(all(table(temp$status) >= mincells)){
    
    de <- FindMarkers(temp, 
                      ident.1 = "PT", 
                      ident.2 = "HD", 
                      test.use = "wilcox",
                      min.pct = 0.1, 
                      logfc.threshold = 0)
    marks <- de[order(de$p_val_adj, decreasing = F),]
    marks$gene_name <- rownames(marks)
    
    write.table(x = marks, file = paste0(outpath, "de_", cl_id, ".txt"), sep = '\t', row.names = F)
    
  }
}

### 7.3.2) GSEA Hallmarks
for(i in seq_len(length(annclusters))){
  cl <- annclusters[i]
  cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
  defile <- paste0(outpath, 'de_', cl_id, ".txt")
  
  if(file.exists(defile)){
    marks <- read.table(file = defile, sep = "\t", header = T)
    gsea <- gsea_run(marks, gmt = hallmarks_gsea)
    write.table(x = gsea, file = paste0(outpath, 'gsea_', cl_id, ".txt"), sep = '\t', row.names = F)
  }
  
}


## 7.4) VEXAS_MUT vs HD

### 7.4.1) DE analysis
annclusters <- names(annotated_cols)
mincells <- 10

mut <- c("LEU", "THR", "VAL")

for(m in mut){
  
  levs <- c(m, "HD")
  sub <- subset(obj, subset = (vexas_mutation %in% levs))
  
  for(i in seq_len(length(annclusters))){
    
    cl <- annclusters[i]
    cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
    temp <- subset(sub, subset = (seurat_annotation %in% cl))
    temp <- SetIdent(object = temp, value = "vexas_mutation")
    
    if(all(table(temp$vexas_mutation) >= mincells)){
      
      de <- FindMarkers(temp, 
                        ident.1 = m, 
                        ident.2 = "HD", 
                        test.use = "wilcox",
                        min.pct = 0.1, 
                        logfc.threshold = 0)
      marks <- de[order(de$p_val_adj, decreasing = F),]
      marks$gene_name <- rownames(marks)
      
      write.table(x = marks, file = paste0(outpath, "de_", cl_id, "_", m, "vHD.txt"), sep = '\t', row.names = F)
      
    }
  }
}

### 7.4.2) GSEA Hallmarks
for(m in mut){  
  for(i in seq_len(length(annclusters))){
    cl <- annclusters[i]
    cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
    defile <- paste0(outpath, "de_", cl_id, "_", m, "vHD.txt")
    
    if(file.exists(defile)){
      marks <- read.table(file = defile, sep = "\t", header = T)
      gsea <- gsea_run(marks, gmt = hallmarks_gsea)
      write.table(x = gsea, file = paste0(outpath, 'gsea_', cl_id, "_", m, "vHD.txt"), sep = '\t', row.names = F)
    }
  }
}


## 7.5) free memory space
l <- ls()
l <- l[!l %in% c("path", "obj")]
rm(list = l); gc()



# 8) UCell module scores and wilcoxon test

## 8.1) load marker gene set
vexas_50 <- qusage::read.gmt(file = paste0(path, "xenograft_signatures/custom_vexas_50.gmt"))
vexas_signature <- vexas_50[[1]]

## 8.2) compute UCell module scores
names(vexas_signature) <- "VEXAS_Xenograft_sig50"
ncol <- ncol(obj@meta.data)
obj <- AddModuleScore_UCell(obj, features = vexas_signature)
colnames(obj@meta.data) <- c(colnames(obj@meta.data)[seq_len(ncol)], names(vexas_signature))

## 8.3) Celltype-wise wilcoxon test: VEXAS vs HD
x <- melt(data = obj@meta.data, id.vars = c("status", "seurat_annotation"), measure.vars = c("VEXAS_Xenograft_sig50"))
w <- x %>% 
  dplyr::group_by(variable, seurat_annotation) %>% 
  dplyr::summarise(pvalue = wilcox.test(x = value[status == "PT"], y = value[status == "HD"])$p.value)
w$p.adjust <- p.adjust(p = w$pvalue)
w <- w[order(w$p.adjust, decreasing = F),]
write.table(x = w, file = paste0(path, "xenograft_signatures/wilcoxon.txt"), sep = "\t", row.names = F, col.names = T, quote = F)




# 9) figures
figpath <- paste0(path, "figures/"); dir.create(path = figpath, showWarnings = F, recursive = T)

## 9.1) figure 5d: Annotated UMAP
obj <- SetIdent(object = obj, value = "seurat_annotation")
g <- DimPlot(obj, reduction = "harmony_umap",
             label.box = T, label = T, label.color = T, label.size = 2) +
  scale_color_manual(values = annotated_cols, limits = levs) +
  scale_fill_manual(values = annotated_cols, limits = levs) +
  ggtitle("Cluster Annotation") +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "right",
        legend.text = element_text(size=7))
ggsave(g, filename = paste0(figpath, "Fig5D_UMAP_AnnotatedClusters.png"),
       width = 10, height = 7, limitsize = FALSE)

## 9.2) figure 5e/5f: GSEA Hallmarks

### 9.2.1) load results
res_VEXASvHD <- c()
for(i in seq_len(length(annclusters))){
  cl <- annclusters[i]
  cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
  file <- paste0(outpath, 'gsea_', cl_id, ".txt")
  
  if(file.exists(file)){
    x <- read.delim(file = file, header = T)
    x <- x[order(x$p.adjust, -x$NES),]
    res_VEXASvHD <- rbind(res_VEXASvHD,
                          data.frame(x, celltype = cl, test = "VEXASvOLD"))
  }
}

res_MUTvHD <- c()
for(m in mut){  
  for(i in seq_len(length(annclusters))){
    cl <- annclusters[i]
    cl_id <- gsub(x = cl, pattern = "/", replacement = "-")
    file <- paste0(outpath, 'gsea_', cl_id, "_", m, "vHD.txt")
    
    if(file.exists(file)){
      id <- paste0(m, "vHD")
      x <- read.delim(file = file, header = T)
      x <- x[order(x$p.adjust, -x$NES),]
      res_MUTvHD <- rbind(res_MUTvHD,
                          data.frame(x, celltype = cl, test = id))
    }
  }
}

res <- rbind(res_VEXASvHD, res_MUTvHD)
res$significance_asterisk <- ""
res$significance_asterisk[res$p.adjust < 0.05] <- "*"
res$significance_asterisk[res$p.adjust < 0.01] <- "**"
res$significance_asterisk[res$p.adjust < 0.001] <- "***"
res <- res %>% tidyr::complete(ID, celltype, test) %>% as.data.frame()
res$ID <- gsub(x = res$ID, pattern = "HALLMARK_", replacement = "")

### 9.2.2) plot
g <- res %>%
  ggplot() + 
  theme_bw() +
  facet_grid(. ~ test) +
  geom_tile(aes(x = celltype,  y = ID, fill = NES)) + 
  geom_text(aes(x = celltype,  y = ID, label = significance_asterisk), size = 2) + 
  theme(plot.title = element_text(hjust = 0.5, size = 10),
        axis.text.x = element_text(angle = 45, , vjust = 1, hjust = 1),
        legend.position = "top",
        strip.background =element_rect(fill="white")) + 
  ylab("") + xlab("") +
  scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
                       limits = c(-4, 4),
                       na.value = "grey")
ggsave(g, filename = paste0(figpath, "Fig5EF_GSEA_Celltype_Hallmarks.png"), 
       width = length(unique(res$celltype))*0.3*length(unique(res_complete$test)), 
       height = length(unique(res$ID))*0.2, limitsize = FALSE)


## 9.3) figure 5g: Monocyte xenograft signature CD34+

### 9.3.1) load wilcoxon test results
w <- read.table(file = paste0(path, "xenograft_signatures/wilcoxon.txt"), sep = "\t", header = T)
w$significance_asterisk <- ""
w$significance_asterisk[w$p.adjust < 0.05] <- "*"
w$significance_asterisk[w$p.adjust < 0.01] <- "**"
w$significance_asterisk[w$p.adjust < 0.001] <- "***"

### 9.3.2) plot
g <- obj@meta.data %>%
  ggplot() + 
  theme_classic() +
  geom_violin(aes(x = seurat_annotation, y = VEXAS_Xenograft_sig50, fill = status), scale = "width") + 
  geom_text(data = w, aes(x = seurat_annotation, y = 0.7, label = significance_asterisk)) +
  theme(plot.title = element_text(hjust = 0.5, size = 10),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = "right") + 
  ylab("UCell Module score") + xlab("") +
  ggtitle("VEXAS xenograft signature (n = 50)") + 
  theme(axis.text.x = element_text(angle = 45, , vjust = 1, hjust = 1)) + 
  scale_fill_manual(values = adjustcolor(col = c("#F8766D", "#02818a"), alpha.f = 0.8), name = "")
ggsave(g, 
       filename = paste0(figpath, "Fig5G_UCell_MonocyteXenograft_WilcoxonTest.png"), 
       width = length(unique(obj$seurat_annotation))*5*0.1, 
       height = 5, 
       limitsize = FALSE)