# HALLMARK GSEA
library('dplyr')
library('ggplot2')
library('reshape2')
library('RColorBrewer')
plot.dir="/DATA/31/molteni/vexas_mouse_3/results/heatmap/"
wdir="/DATA/31/molteni/vexas_mouse_3/results"

#heatmap custom cluster & order (05/06/24)
wdir="/DATA/31/molteni/vexas_mouse_3/results"
plot.dir="/DATA/31/molteni/vexas_mouse_3/results/heatmap/"
Full_final <- readRDS("/DATA/31/molteni/vexas_mouse_3/results/Full/01-seurat/Full_final_annot_filt.rds")
hallmark_ds <- c("A_C")
clusters <- unique(Full_final$celltypes_annot)
hallmark.full <- data.frame()
for (ds in hallmark_ds) {
  for (i in clusters) {
    s <- paste(ds,"markers_cluster", i,sep="_") #qui cambiare 1 con i e inserire ciclo for
    ff <- paste(wdir, "GSEA", "AvsC",s, "tables", paste(s, "NA-gsea-h.all.v7.2.symbols.txt", sep = "-"), sep = "/")
    ds.t <- read.table(ff, header = T, sep = "\t")
    ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalues"), drop = FALSE]
    ds.t.filt$Dataset <- paste(ds,i,sep="_")
    if (nrow(hallmark.full) == 0) {
      hallmark.full <- ds.t.filt
    } else {
      hallmark.full <- rbind(hallmark.full, ds.t.filt)
    }
  }
}

hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset","p.adjust")]

##
all_hallmark <- hallmark.full.filt
all_hallmark$Dataset <- gsub("A_C", "VEXAS-like_WT", all_hallmark$Dataset)

all_hallmark$star <- cut(all_hallmark$p.adjust, breaks=c(-Inf, 0.001, 0.01, 0.05, Inf), label=c("***", "**", "*", ""))  # Create column of significance labels
all_hallmark$ID <- gsub(x = all_hallmark$ID, "HALLMARK_", "")
hallmark.order <- all_hallmark %>% group_by(ID) %>% summarise(Pos = sum(NES))
hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
all_hallmark.nn <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="NES"), id.vars = c("ID"))
colnames(all_hallmark.nn) <- c("ID", "Dataset", "NES")
all_hallmark.pp <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="star"), id.vars = c("ID"))
colnames(all_hallmark.pp) <- c("ID", "Dataset", "STARS")
all_hallmark.tt <- merge(all_hallmark.nn,all_hallmark.pp)
all_hallmark.tt$ID <- factor(all_hallmark.tt$ID, levels = hallmark.order.terms$ID)

level_order <- c("VEXAS-like_WT_Erythroid_prog", "VEXAS-like_WT_Prolif_Myelo_NK_prec", "VEXAS-like_WT_Promyelocyte",
                 "VEXAS-like_WT_Monocyte_biased_prog","VEXAS-like_WT_pDC_biased_prog","VEXAS-like_WT_GMP_2",
                 "VEXAS-like_WT_GMP_1","VEXAS-like_WT_CMP","VEXAS-like_WT_HSC_MPP",
                 "VEXAS-like_WT_LMPP","VEXAS-like_WT_CLP", "VEXAS-like_WT_Prolif_lymphoid_prec", 
                 "VEXAS-like_WT_B_cell_progenitors_1", "VEXAS-like_WT_B_cell_progenitors_2")

all_hallmark.tt$Dataset <- factor(all_hallmark.tt$Dataset, levels =level_order )
p.hallmark <- ggplot(all_hallmark.tt, aes(x = Dataset, y = ID)) +
  theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(axis.text.y=element_text(size=6))+
  geom_tile(aes(fill = NES), colour = "white") +
  geom_text(aes(label = STARS), color="black", size=3) +
  scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
                       limits = c(-5, 5),
                       na.value = "grey"
                       #guide = guide_colourbar(reverse = TRUE)
  ) +
  ylab("") + 
  xlab("") +
  coord_fixed(ratio = 0.4)

ggsave(filename = paste(plot.dir, paste("hallmark_AvsC_ordered_last", "heatmap.pdf", sep = "-"), sep = "/"), 
       plot = p.hallmark,
       width = 10, height = 8, dpi = 600)



#heatmap specific comparison GMP2 vs GMP1 (05/06/24)
wdir="/DATA/31/molteni/vexas_mouse_3/results"
plot.dir="/DATA/31/molteni/vexas_mouse_3/results/heatmap/"
Full_final <- readRDS("/DATA/31/molteni/vexas_mouse_3/results/Full/01-seurat/Full_final_annot_filt.rds")
hallmark_ds <- c("A", "C", "Full")
clusters <- unique(Full_final$celltypes_annot)
hallmark.full <- data.frame()
for (ds in hallmark_ds) {
  #for (i in clusters) {
    #s <- paste("_",ds,sep="") #qui cambiare 1 con i e inserire ciclo for
    ff <- paste(wdir, "GSEA_GMP2_GMP1", paste("GMP2vsGMP1" , ds , sep="_"),"tables", paste(paste("GMP2vsGMP1" , ds , sep="_"),"NA-gsea-h.all.v7.2.symbols.txt", sep = "-"), sep = "/")
    ds.t <- read.table(ff, header = T, sep = "\t")
    ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalues"), drop = FALSE]
    ds.t.filt$Dataset <- paste(ds,"GMP2vsGMP1",sep="_")
    if (nrow(hallmark.full) == 0) {
      hallmark.full <- ds.t.filt
    } else {
      hallmark.full <- rbind(hallmark.full, ds.t.filt)
    #}
  }
}

hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset","p.adjust")]

##
all_hallmark <- hallmark.full.filt

all_hallmark$star <- cut(all_hallmark$p.adjust, breaks=c(-Inf, 0.001, 0.01, 0.05, Inf), label=c("***", "**", "*", ""))  # Create column of significance labels
all_hallmark$ID <- gsub(x = all_hallmark$ID, "HALLMARK_", "")
hallmark.order <- all_hallmark %>% group_by(ID) %>% summarise(Pos = sum(NES))
hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
all_hallmark.nn <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="NES"), id.vars = c("ID"))
colnames(all_hallmark.nn) <- c("ID", "Dataset", "NES")
all_hallmark.pp <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="star"), id.vars = c("ID"))
colnames(all_hallmark.pp) <- c("ID", "Dataset", "STARS")
all_hallmark.tt <- merge(all_hallmark.nn,all_hallmark.pp)
all_hallmark.tt$ID <- factor(all_hallmark.tt$ID, levels = hallmark.order.terms$ID)

level_order <- c("A_GMP2vsGMP1", "C_GMP2vsGMP1", "Full_GMP2vsGMP1")

all_hallmark.tt$Dataset <- factor(all_hallmark.tt$Dataset, levels =level_order )
p.hallmark <- ggplot(all_hallmark.tt, aes(x = Dataset, y = ID)) +
  theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(axis.text.y=element_text(size=6))+
  geom_tile(aes(fill = NES), colour = "white") +
  geom_text(aes(label = STARS), color="black", size=3) +
  scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
                       limits = c(-5, 5),
                       na.value = "grey"
                       #guide = guide_colourbar(reverse = TRUE)
  ) +
  ylab("") + 
  xlab("") +
  coord_fixed(ratio = 0.4)

ggsave(filename = paste(plot.dir, paste("hallmark_GMP2_GMP1", "heatmap.pdf", sep = "-"), sep = "/"), 
       plot = p.hallmark,
       width = 10, height = 8, dpi = 600)

#heatmap specific comparison Bprog1 vs BProg2 (05/06/24)
wdir="/DATA/31/molteni/vexas_mouse_3/results"
plot.dir="/DATA/31/molteni/vexas_mouse_3/results/heatmap/"
Full_final <- readRDS("/DATA/31/molteni/vexas_mouse_3/results/Full/01-seurat/Full_final_annot_filt.rds")
hallmark_ds <- c("A", "C", "Full")
clusters <- unique(Full_final$celltypes_annot)
hallmark.full <- data.frame()
for (ds in hallmark_ds) {
  #for (i in clusters) {
  #s <- paste("_",ds,sep="") #qui cambiare 1 con i e inserire ciclo for
  ff <- paste(wdir, "GSEA_BP1_BP2", paste("BP1vsBP2" , ds , sep="_"),"tables", paste(paste("BP1vsBP2" , ds , sep="_"),"NA-gsea-h.all.v7.2.symbols.txt", sep = "-"), sep = "/")
  ds.t <- read.table(ff, header = T, sep = "\t")
  ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalues"), drop = FALSE]
  ds.t.filt$Dataset <- paste(ds,"BP1vsBP2",sep="_")
  if (nrow(hallmark.full) == 0) {
    hallmark.full <- ds.t.filt
  } else {
    hallmark.full <- rbind(hallmark.full, ds.t.filt)
    #}
  }
}

hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset","p.adjust")]

##
all_hallmark <- hallmark.full.filt

all_hallmark$star <- cut(all_hallmark$p.adjust, breaks=c(-Inf, 0.001, 0.01, 0.05, Inf), label=c("***", "**", "*", ""))  # Create column of significance labels
all_hallmark$ID <- gsub(x = all_hallmark$ID, "HALLMARK_", "")
hallmark.order <- all_hallmark %>% group_by(ID) %>% summarise(Pos = sum(NES))
hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
all_hallmark.nn <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="NES"), id.vars = c("ID"))
colnames(all_hallmark.nn) <- c("ID", "Dataset", "NES")
all_hallmark.pp <- melt(reshape2::dcast(all_hallmark, ID ~ Dataset, value.var="star"), id.vars = c("ID"))
colnames(all_hallmark.pp) <- c("ID", "Dataset", "STARS")
all_hallmark.tt <- merge(all_hallmark.nn,all_hallmark.pp)
all_hallmark.tt$ID <- factor(all_hallmark.tt$ID, levels = hallmark.order.terms$ID)

level_order <- c("A_BP1vsBP2", "C_BP1vsBP2", "Full_BP1vsBP2")

all_hallmark.tt$Dataset <- factor(all_hallmark.tt$Dataset, levels =level_order )
p.hallmark <- ggplot(all_hallmark.tt, aes(x = Dataset, y = ID)) +
  theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(axis.text.y=element_text(size=6))+
  geom_tile(aes(fill = NES), colour = "white") +
  geom_text(aes(label = STARS), color="black", size=3) +
  scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
                       limits = c(-5, 5),
                       na.value = "grey"
                       #guide = guide_colourbar(reverse = TRUE)
  ) +
  ylab("") + 
  xlab("") +
  coord_fixed(ratio = 0.4)

ggsave(filename = paste(plot.dir, paste("hallmark_BP1_BP2", "heatmap.pdf", sep = "-"), sep = "/"), 
       plot = p.hallmark,
       width = 10, height = 8, dpi = 600)