diff --git a/scRNAseq_mixed/seurat_analysis.R b/scRNAseq_mixed/seurat_analysis.R
index 84cb77eef8c12eb4f0cf46966865e241e54f19cf..9ac16c565f3593a022d437761d0ccbe5ef74ec9b 100644
--- a/scRNAseq_mixed/seurat_analysis.R
+++ b/scRNAseq_mixed/seurat_analysis.R
@@ -11,78 +11,21 @@ suppressPackageStartupMessages(library(ggrepel))
 suppressPackageStartupMessages(library(clusterProfiler))
 suppressPackageStartupMessages(library(RColorBrewer))
 ## load seurat data
-h48_final <- readRDS("data/mixed_h48_final.rds")
-d4_final <- readRDS("data/mixed_D4_final.rds")
-Idents(h48_final) <- "RNA_snn_h.orig.ident_res.1.2"
-Idents(d4_final) <- "RNA_snn_h.orig.ident_res.1.2"
-h48_final@meta.data[['orig.ident']] <- recode_factor(
-  h48_final@meta.data[["orig.ident"]],
-  "Sample3" = "AAV9",
-  "Sample5" = "Spk",
-  "Sample7" = "UT")
-h48_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]] <- recode_factor(
-  h48_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]], 
-  "0" = "Undifferentiated", 
-  "1" = "Astrocytes", 
-  "2" = "Astrocytes", 
-  "3" = "Astrocytes", 
-  "4" = "OPC", 
-  "5" = "Oligodendrocytes", 
-  "6" = "Astrocytes", 
-  "7" = "Astrocytes", 
-  "8" = "Oligodendrocytes",
-  "9" = "Immature.Neurons",
-  "10" = "Oligodendrocytes",
-  "11" = "Immature.Neurons",
-  "12" = "Immature.Astrocytes",
-  "13" = "Immature.Astrocytes",
-  "14" = "Astrocytes",
-  "15" = "Neurons",
-  "16" = "Oligodendrocytes",
-  "17" = "OPC")
-h48_final$RNA_snn_h.orig.ident_res.1.2 <- factor(
-  h48_final$RNA_snn_h.orig.ident_res.1.2,
-  levels = c("Astrocytes",
-             "Neurons",
-             "Oligodendrocytes",
-             "Immature.Astrocytes",
-             "Immature.Neurons",
-             "OPC",
-             "Undifferentiated"))
-d4_final@meta.data[['orig.ident']] <- recode_factor(
-  d4_final@meta.data[["orig.ident"]],
-  "Sample2" = "AAV2",
-  "Sample4" = "AAV9",
-  "Sample6" = "Spk",
-  "Sample8" = "UT")
-d4_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]] <- recode_factor(
-  d4_final@meta.data[["RNA_snn_h.orig.ident_res.1.2"]], 
-  "0" = "Undifferentiated", 
-  "1" = "Oligodendrocytes", 
-  "2" = "Neurons", 
-  "3" = "Astrocytes", 
-  "4" = "Immature.Neurons", 
-  "5" = "OPC", 
-  "6" = "Astrocytes", 
-  "7" = "Oligodendrocytes", 
-  "8" = "Oligodendrocytes",
-  "9" = "Immature.Neurons",
-  "10" = "Oligodendrocytes",
-  "11" = "Oligodendrocytes",
-  "12" = "Astrocytes",
-  "13" = "Immature.Astrocytes",
-  "14" = "OPC",
-  "15" = "Astrocytes",
-  "16" = "Astrocytes")
-d4_final$RNA_snn_h.orig.ident_res.1.2 <- factor(
-  d4_final$RNA_snn_h.orig.ident_res.1.2,
-  levels = c("Astrocytes",
-             "Neurons",
-             "Oligodendrocytes",
-             "Immature.Astrocytes",
-             "Immature.Neurons",
-             "OPC",
-             "Undifferentiated"))
+tp <- c("h48","d4")
+for (t in tp) {
+  sc_obj <- readRDS(paste0("data/mixed_",t,"_final.rds"))
+  Idents(sc_obj) <- "Cell_Type"
+  sc_obj$Cell_type <- factor(
+    sc_obj$Cell_Type,
+    levels = c("Astrocytes",
+               "Neurons",
+               "Oligodendrocytes",
+               "Immature.Astrocytes",
+               "Immature.Neurons",
+               "OPC",
+               "Undifferentiated"))
+  assign(paste0(t,"_final"), sc_obj)
+}
 ggplotColours <- function(n = 6, h = c(0, 360) + 15){
   if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
   hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
@@ -90,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
 ## umap
 DimPlot(object = h48_final,
         reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.1.2",
+        group.by = "Cell_Type",
         label = T,
         repel = T,
         pt.size = 0.2,
@@ -105,7 +48,7 @@ DimPlot(object = h48_final,
                      values = ggplotColours(n = 7))
 
 ## donut plot
-ggplot_object_clusters_labels <- melt(table(h48_final$RNA_snn_h.orig.ident_res.1.2))
+ggplot_object_clusters_labels <- melt(table(h48_final$Cell_Type))
 colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
 totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
 totalB$percentage <- totalB$percentage*100
@@ -159,7 +102,7 @@ calc_helper <- function(object,genes){
 aav9 <- subset(h48_final, orig.ident=="AAV9")
 gfp_percentage_cell_type <- PrctCellExpringGene(aav9,
                                                 genes ="GFP",
-                                                group.by = "RNA_snn_h.orig.ident_res.1.2")
+                                                group.by = "Cell_Type")
 gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100
 gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)]
 colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types")
@@ -201,7 +144,7 @@ h48_aav9_ut <- subset(h48_final, orig.ident == "AAV9" | orig.ident == "UT")
 d4_aav9_ut <- subset(d4_final, orig.ident == "AAV9" | orig.ident == "UT")
 DimPlot(object = h48_aav9_ut,
         reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.1.2",
+        group.by = "Cell_Type",
         split.by = "orig.ident",
         label = F,
         repel = T,
@@ -218,7 +161,7 @@ DimPlot(object = h48_aav9_ut,
                      values = ggplotColours(n = 7))
 DimPlot(object = d4_aav9_ut,
         reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.1.2",
+        group.by = "Cell_Type",
         split.by = "orig.ident",
         label = F,
         repel = T,
@@ -234,12 +177,11 @@ DimPlot(object = d4_aav9_ut,
                               "Undifferentiated"),
                      values = ggplotColours(n = 7))
 samples <- c("AAV9", "UT")
-tp <- c("h48","d4")
 for (t in tp) {
   obj_t <- get(paste0(t,"_final"))
   for (i in samples) {
     obj_i <- subset(obj_t, orig.ident == i)
-    ggplot_object_clusters_labels <- melt(table(obj_i$RNA_snn_h.orig.ident_res.1.2))
+    ggplot_object_clusters_labels <- melt(table(obj_i$Cell_Type))
     colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
     totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
     totalB$percentage <- totalB$percentage*100
@@ -276,15 +218,9 @@ donut_d4_final <- (donut_d4_AAV9 | donut_d4_UT)
 ## markers
 for (t in tp) {
   obj_t <- get(paste0(t,"_final"))
-  neurons <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Neurons" | RNA_snn_h.orig.ident_res.1.2=="Immature.Neurons")
-  astro <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Astrocytes" | RNA_snn_h.orig.ident_res.1.2=="Immature.Astrocytes")
-  oligo <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Oligodendrocytes" | RNA_snn_h.orig.ident_res.1.2=="OPC")
+  neurons <- subset(obj_t, Cell_Type=="Neurons" | Cell_Type=="Immature.Neurons")
+  astro <- subset(obj_t, Cell_Type=="Astrocytes" | Cell_Type=="Immature.Astrocytes")
+  oligo <- subset(obj_t, Cell_Type=="Oligodendrocytes" | Cell_Type=="OPC")
   assign(paste0("neurons_",t,"_obj_x_markers"), neurons)
   assign(paste0("astro_",t,"_obj_x_markers"), astro)
   assign(paste0("oligo_",t,"_obj_x_markers"), oligo)
diff --git a/scRNAseq_organoids/seurat_analysis.R b/scRNAseq_organoids/seurat_analysis.R
index 92b1b45ae66ec6cd6a5b1700e2708633a4b00fbc..f5ecbb2fc11b8dbcf87c134ecf7631be86b234c5 100644
--- a/scRNAseq_organoids/seurat_analysis.R
+++ b/scRNAseq_organoids/seurat_analysis.R
@@ -11,70 +11,21 @@ suppressPackageStartupMessages(library(ggrepel))
 suppressPackageStartupMessages(library(clusterProfiler))
 suppressPackageStartupMessages(library(RColorBrewer))
 ## load seurat data
-h48_final <- readRDS("data/organoids_h48_final.rds")
-d5_final <- readRDS("data/organoids_D5_final.rds")
-Idents(h48_final) <- "RNA_snn_h.orig.ident_res.0.6"
-Idents(d5_final) <- "RNA_snn_h.orig.ident_res.0.6"
-h48_final@meta.data[['orig.ident']] <- recode_factor(
-  h48_final@meta.data[["orig.ident"]],
-  "Sample1" = "AAV2",
-  "Sample3" = "AAV9",
-  "Sample5" = "Spk",
-  "Sample7" = "UT")
-h48_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]] <- recode_factor(
-  h48_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]], 
-  "0" = "Astrocytes", 
-  "1" = "Immature.Astrocytes", 
-  "2" = "Undifferentiated", 
-  "3" = "Neurons", 
-  "4" = "Immature.Neurons", 
-  "5" = "OPC", 
-  "6" = "Astrocytes", 
-  "7" = "Oligodendrocytes", 
-  "8" = "Undifferentiated",
-  "9" = "Astrocytes",
-  "10" = "Undifferentiated",
-  "11" = "Undifferentiated")
-h48_final$RNA_snn_h.orig.ident_res.0.6 <- factor(
-  h48_final$RNA_snn_h.orig.ident_res.0.6,
-  levels = c("Astrocytes",
-             "Neurons",
-             "Oligodendrocytes",
-             "Immature.Astrocytes",
-             "Immature.Neurons",
-             "OPC",
-             "Undifferentiated"))
-d5_final@meta.data[['orig.ident']] <- recode_factor(
-  d5_final@meta.data[["orig.ident"]],
-  "Sample2" = "AAV2",
-  "Sample4" = "AAV9",
-  "Sample6" = "Spk",
-  "Sample8" = "UT")
-d5_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]] <- recode_factor(
-  d5_final@meta.data[["RNA_snn_h.orig.ident_res.0.6"]], 
-  "0" = "Astrocytes", 
-  "1" = "Immature.Neurons", 
-  "2" = "Undifferentiated", 
-  "3" = "Neurons", 
-  "4" = "Oligodendrocytes", 
-  "5" = "Undifferentiated", 
-  "6" = "Undifferentiated", 
-  "7" = "Immature.Astrocytes", 
-  "8" = "OPC",
-  "9" = "Immature.Neurons",
-  "10" = "Astrocytes",
-  "11" = "Undifferentiated",
-  "12" = "Astrocytes",
-  "13" = "Neurons")
-d5_final$RNA_snn_h.orig.ident_res.0.6 <- factor(
-  d5_final$RNA_snn_h.orig.ident_res.0.6,
-  levels = c("Astrocytes",
-             "Neurons",
-             "Oligodendrocytes",
-             "Immature.Astrocytes",
-             "Immature.Neurons",
-             "OPC",
-             "Undifferentiated"))
+tp <- c("h48","d5")
+for (t in tp) {
+  sc_obj <- readRDS(paste0("data/organoids_",t,"_final.rds"))
+  Idents(sc_obj) <- "Cell_Type"
+  sc_obj$Cell_type <- factor(
+    sc_obj$Cell_Type,
+    levels = c("Astrocytes",
+               "Neurons",
+               "Oligodendrocytes",
+               "Immature.Astrocytes",
+               "Immature.Neurons",
+               "OPC",
+               "Undifferentiated"))
+  assign(paste0(t,"_final"), sc_obj)
+}
 ggplotColours <- function(n = 6, h = c(0, 360) + 15){
   if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
   hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
@@ -82,7 +33,7 @@ ggplotColours <- function(n = 6, h = c(0, 360) + 15){
 ## umap
 DimPlot(object = h48_final,
         reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.0.6",
+        group.by = "Cell_Type",
         label = T,
         repel = T,
         pt.size = 0.2,
@@ -97,7 +48,7 @@ DimPlot(object = h48_final,
                      values = ggplotColours(n = 7))
 
 ## donut plot
-ggplot_object_clusters_labels <- melt(table(h48_final$RNA_snn_h.orig.ident_res.0.6))
+ggplot_object_clusters_labels <- melt(table(h48_final$Cell_Type))
 colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
 totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
 totalB$percentage <- totalB$percentage*100
@@ -149,10 +100,10 @@ calc_helper <- function(object,genes){
   }else{return(NA)}
 }
 
-aav9 <- subset(h48_final, orig.ident=="AAV9")
+aav9 <- subset(d5_final, orig.ident=="AAV9")
 gfp_percentage_cell_type <- PrctCellExpringGene(aav9,
                                                 genes ="GFP",
-                                                group.by = "RNA_snn_h.orig.ident_res.0.6")
+                                                group.by = "Cell_Type")
 gfp_percentage_cell_type$Percentage <- gfp_percentage_cell_type$Cell_proportion*100
 gfp_percentage_cell_type <- gfp_percentage_cell_type[,c(4,3)]
 colnames(gfp_percentage_cell_type) <- c("Percentage", "Cell_Types")
@@ -190,49 +141,14 @@ ggplot(data=gfp_percentage_cell_type,
     axis.text.y = element_text(size=18*96/72),
     legend.position = "none",
     aspect.ratio = 1.1)
-h48_aav9_ut <- subset(h48_final, orig.ident == "AAV9" | orig.ident == "UT") 
-d4_aav9_ut <- subset(d4_final, orig.ident == "AAV9" | orig.ident == "UT")
-DimPlot(object = h48_aav9_ut,
-        reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.1.2",
-        split.by = "orig.ident",
-        label = F,
-        repel = T,
-        pt.size = 0.2,
-        label.size = 6,
-        ncol = 2) + 
-  scale_color_manual(labels=c("Astrocytes",
-                              "Neurons",
-                              "Oligodendrocytes",
-                              "Immature Astrocytes",
-                              "Immature Neurons",
-                              "OPC",
-                              "Undifferentiated"),
-                     values = ggplotColours(n = 7))
-DimPlot(object = d4_aav9_ut,
-        reduction = "umap.harmony.orig.ident",
-        group.by = "RNA_snn_h.orig.ident_res.1.2",
-        split.by = "orig.ident",
-        label = F,
-        repel = T,
-        pt.size = 0.2,
-        label.size = 6,
-        ncol = 2) + 
-  scale_color_manual(labels=c("Astrocytes",
-                              "Neurons",
-                              "Oligodendrocytes",
-                              "Immature Astrocytes",
-                              "Immature Neurons",
-                              "OPC",
-                              "Undifferentiated"),
-                     values = ggplotColours(n = 7))
+
 samples <- c("AAV9", "UT")
-tp <- c("h48","d4")
+tp <- c("h48","d5")
 for (t in tp) {
   obj_t <- get(paste0(t,"_final"))
   for (i in samples) {
     obj_i <- subset(obj_t, orig.ident == i)
-    ggplot_object_clusters_labels <- melt(table(obj_i$RNA_snn_h.orig.ident_res.1.2))
+    ggplot_object_clusters_labels <- melt(table(obj_i$Cell_Type))
     colnames(ggplot_object_clusters_labels) <- c("Cell_type", "Count")
     totalB <- ggplot_object_clusters_labels %>% mutate(percentage = Count/sum(Count))
     totalB$percentage <- totalB$percentage*100
@@ -264,41 +180,88 @@ for (t in tp) {
     assign(paste0("donut_",t,"_", i), donut_i)
   }
 }
-donut_h48_final <- (donut_h48_AAV9 | donut_h48_UT)
-donut_d4_final <- (donut_d4_AAV9 | donut_d4_UT)
+donut_final <- (donut_h48_UT | donut_h48_AAV9 | donut_d5_UT | donut_d5_AAV9)
 ## markers
 for (t in tp) {
   obj_t <- get(paste0(t,"_final"))
-  neurons <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Neurons" | RNA_snn_h.orig.ident_res.1.2=="Immature.Neurons")
-  astro <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Astrocytes" | RNA_snn_h.orig.ident_res.1.2=="Immature.Astrocytes")
-  oligo <- subset(
-    obj_t, 
-    RNA_snn_h.orig.ident_res.1.2=="Oligodendrocytes" | RNA_snn_h.orig.ident_res.1.2=="OPC")
+  neurons <- subset(obj_t, Cell_Type=="Neurons" | Cell_Type=="Immature.Neurons")
+  astro <- subset(obj_t, Cell_Type=="Astrocytes" | Cell_Type=="Immature.Astrocytes")
+  oligo <- subset(obj_t, Cell_Type=="Oligodendrocytes" | Cell_Type=="OPC")
   assign(paste0("neurons_",t,"_obj_x_markers"), neurons)
   assign(paste0("astro_",t,"_obj_x_markers"), astro)
   assign(paste0("oligo_",t,"_obj_x_markers"), oligo)
 }
-for (k in ls()[grep("_obj_x_markers", ls())]) {
-  cells <- get(k)
-  degs <- FindMarkers(object = cells, 
-                      ident.1 = "AAV9", 
-                      ident.2 = "UT", 
-                      group.by = "orig.ident",
-                      only.pos = FALSE, 
-                      min.pct = 0, 
-                      test.use = "wilcox", 
-                      logfc.threshold = 0,
-                      min.cells.group = 0)
-  assign(paste0(k,"_aav9_vs_ut_degs"), degs)
+
+cell_types <- c("neurons_h48", "astro_h48", "oligo_h48",
+                "neurons_d5", "astro_d5", "oligo_d5")
+samples_degs <- c("AAV9", "AAV2", "Spk")
+
+for (k in cell_types) {
+  cells <- get(paste0(k,"_obj_x_markers"))
+  for (j in samples_degs) {
+    print(paste("Extracting differentially expressed genes of",j,"vs UT"))
+    degs <- FindMarkers(cells, 
+                        ident.1 = j, 
+                        ident.2 = "UT", 
+                        group.by = "orig.ident", 
+                        only.pos = FALSE, 
+                        min.pct = 0.1, 
+                        test.use = "wilcox", 
+                        logfc.threshold = 0.1)
+    assign(paste(k,j,"vs_UT_degs",sep = "_"), degs, envir = .GlobalEnv)
+    degs_all <- FindMarkers(object = cells, 
+                            ident.1 = j, 
+                            ident.2 = "UT", 
+                            group.by = "orig.ident",
+                            only.pos = FALSE, 
+                            min.pct = 0, 
+                            test.use = "wilcox", 
+                            logfc.threshold = 0,
+                            min.cells.group = 0)
+    assign(paste(k,j,"vs_UT_degs_all",sep = "_"), degs_all, envir = .GlobalEnv)
+    ## volcano plot
+    print(paste("Making volcano plot for DEGs of...", j, " vs UT" ))
+    logfc.thr <- 0
+    adj.pval.thr <- 0.05
+    data <- data.frame(gene = row.names(degs), 
+                       pval = -log10(degs$p_val_adj), 
+                       lfc = degs$avg_log2FC)
+    data <- na.omit(data)
+    data <- mutate(data,color=case_when(data$lfc>logfc.thr & data$pval> -log10(adj.pval.thr)~"Overexpressed",
+                                        data$lfc< -logfc.thr & data$pval> -log10(adj.pval.thr)~"Underexpressed",
+                                        abs(data$lfc)<logfc.thr|data$pval< -log10(adj.pval.thr)~"NonSignificant"))
+    data <- data[order(data$pval, decreasing = TRUE),]
+    highl_up <- head(subset(data, color != "NonSignificant" & lfc > 0), 10)
+    highl_down <- head(subset(data, color != "NonSignificant" & lfc < 0), 10)
+    highl <- rbind(highl_up, highl_down)
+    
+    sig <- subset(degs, degs$p_val_adj<0.05)
+    n_sig <- nrow(sig)
+    
+    vol <- ggplot(data, aes(x = lfc, y = pval, colour = color)) +
+      geom_point(size = 3, alpha = 0.8, na.rm = T) +
+      geom_hline(yintercept = -log10(adj.pval.thr), colour = "darkgray") +
+      geom_label_repel(data = highl, 
+                       aes(label = gene), 
+                       show.legend = FALSE, 
+                       size = 6,
+                       family = "Arial") +
+      xlab(expression(log[2]("Fold Change"))) + 
+      ylab(expression(-log[10]("adjusted p-value"))) +
+      scale_color_manual(name = "Expression",
+                         values = c(Overexpressed = "red",
+                                    Underexpressed = "royalblue3",
+                                    NonSignificant = "black")) +
+      ggtitle(paste("Significant genes = ", n_sig)) +
+      theme_bw()
+    assign(paste(k,j,"vs_UT_volcano",sep = "_"), vol, envir = .GlobalEnv)
+  }
 }
+rm(degs_all)
 ## gsea
 h_gmt <- read.gmt("data/h.all.v7.2.symbols.gmt")
-for (i in ls()[grep("aav9_vs_ut_degs", ls())]) {
-  gene_res <- get(i)
+for (g in ls()[grep("degs_all", ls())]) {
+  gene_res <- get(g)
   logfc_symbol <- as.vector(gene_res$avg_log2FC)
   names(logfc_symbol) <- row.names(gene_res)
   logfc_symbol <- sort(logfc_symbol, decreasing = TRUE)
@@ -307,60 +270,63 @@ for (i in ls()[grep("aav9_vs_ut_degs", ls())]) {
                      nPerm = 10000, 
                      pvalueCutoff = 1)
   contr.gsea.symb <- as.data.frame(contr.gsea)
-  assign(paste0("GSEA_",i), contr.gsea.symb)
+  assign(paste0("GSEA_",g), contr.gsea.symb)
 }
 ## GSEA heatmaps
 type <- c("astro", "neurons", "oligo")
+samples_degs <- c("AAV9", "AAV2", "Spk")
 for (p in type) {
-  hallmark.full <- data.frame()
-  for (ds in ls()[grep(paste0("GSEA_", p), ls())]) {
-    ds.t <- get(ds)
-    ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalues"), drop = FALSE]
-    ds.t.filt$Dataset <- ds
-    if (nrow(hallmark.full) == 0) {
-      hallmark.full <- ds.t.filt
-    } else {
-      hallmark.full <- rbind(hallmark.full, ds.t.filt)
+  for (j in samples_degs) {
+    hallmark.full <- data.frame()
+    for (ds in ls()[grep(paste0("GSEA_.*",p,".*",j), ls())]) {
+      ds.t <- get(ds)
+      ds.t.filt <- ds.t[,c("ID", "setSize", "enrichmentScore", "NES", "pvalue", "p.adjust", "qvalue"), drop = FALSE]
+      ds.t.filt$Dataset <- ds
+      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", "qvalue")]
+    hallmark.full.filt$ID <- gsub(x = hallmark.full.filt$ID, "HALLMARK_", "")
+    hallmark.full.filt$sig <- ifelse(
+      hallmark.full.filt$qvalue <= 0.05 & hallmark.full.filt$qvalue > 0.01, '*', 
+      ifelse(hallmark.full.filt$qvalue <= 0.01 & hallmark.full.filt$qvalue > 0.001, '**', 
+             ifelse(hallmark.full.filt$qvalue <= 0.001 & hallmark.full.filt$qvalue > 0.0001, '***', 
+                    ifelse(hallmark.full.filt$qvalue <= 0.0001, '****', ''))))
+    hallmark.order <- hallmark.full.filt %>% 
+      group_by(ID) %>% 
+      summarise(Pos = sum(NES))
+    hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
+    hallmark.full.filt.tt <- reshape2::melt(reshape2::dcast(hallmark.full.filt, 
+                                                            ID ~ Dataset, 
+                                                            value.var="NES"), id.vars = c("ID"))
+    colnames(hallmark.full.filt.tt) <- c("ID", "Dataset", "NES")
+    hallmark.full.filt.tt$ID <- factor(hallmark.full.filt.tt$ID, levels = hallmark.order.terms$ID)
+    hallmark.full.filt.tt <- merge(hallmark.full.filt.tt, 
+                                   hallmark.full.filt[, c("ID", "Dataset", "sig")], 
+                                   by = c("ID", "Dataset"), 
+                                   all.x = TRUE)
+    labels <- c("Day4", "Day2")
+    p.hallmark <- ggplot(hallmark.full.filt.tt, 
+                         aes(x = Dataset, y = ID)) +
+      geom_tile(aes(fill = NES), colour = "white") +
+      geom_text(aes(label = paste(sig)), size=4*96/72, fontface = "bold") +
+      scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
+                           limits = c(-3.5, 3.5),
+                           na.value = "white") +
+      ylab("") + 
+      xlab("") +
+      coord_fixed(ratio = 0.6) +
+      scale_x_discrete(labels = labels) +
+      theme_bw(base_size = 12) +
+      theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 18*96/72),
+            axis.text.y = element_text(face = "bold", size = 12*96/72),
+            legend.text = element_text(face = "bold", size = 14*96/72),
+            legend.title = element_text(face = "bold", size = 12*96/72))
+    assign(paste("p.hallmark", p,j,sep = "_"), p.hallmark)
   }
-  hallmark.full.filt <- hallmark.full[, c("ID", "NES", "Dataset", "qvalues")]
-  hallmark.full.filt$ID <- gsub(x = hallmark.full.filt$ID, "HALLMARK_", "")
-  hallmark.full.filt$sig <- ifelse(
-    hallmark.full.filt$qvalues <= 0.05 & hallmark.full.filt$qvalues > 0.01, '*', 
-    ifelse(hallmark.full.filt$qvalues <= 0.01 & hallmark.full.filt$qvalues > 0.001, '**', 
-           ifelse(hallmark.full.filt$qvalues <= 0.001 & hallmark.full.filt$qvalues > 0.0001, '***', 
-                  ifelse(hallmark.full.filt$qvalues <= 0.0001, '****', ''))))
-  hallmark.order <- hallmark.full.filt %>% 
-    group_by(ID) %>% 
-    summarise(Pos = sum(NES))
-  hallmark.order.terms <- hallmark.order[order(hallmark.order$Pos, decreasing = FALSE), "ID", drop = FALSE]
-  hallmark.full.filt.tt <- reshape2::melt(reshape2::dcast(hallmark.full.filt, 
-                                                          ID ~ Dataset, 
-                                                          value.var="NES"), id.vars = c("ID"))
-  colnames(hallmark.full.filt.tt) <- c("ID", "Dataset", "NES")
-  hallmark.full.filt.tt$ID <- factor(hallmark.full.filt.tt$ID, levels = hallmark.order.terms$ID)
-  hallmark.full.filt.tt <- merge(hallmark.full.filt.tt, 
-                                 hallmark.full.filt[, c("ID", "Dataset", "sig")], 
-                                 by = c("ID", "Dataset"), 
-                                 all.x = TRUE)
-  labels <- c("Day4", "Day2")
-  p.hallmark <- ggplot(hallmark.full.filt.tt, 
-                       aes(x = Dataset, y = ID)) +
-    geom_tile(aes(fill = NES), colour = "white") +
-    geom_text(aes(label = paste(sig)), size=4*96/72, fontface = "bold") +
-    scale_fill_gradientn(colours = colorRampPalette(rev(brewer.pal(11,"RdBu")))(100),
-                         limits = c(-3.5, 3.5),
-                         na.value = "white") +
-    ylab("") + 
-    xlab("") +
-    coord_fixed(ratio = 0.6) +
-    scale_x_discrete(labels = labels) +
-    theme_bw(base_size = 12) +
-    theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 18*96/72),
-          axis.text.y = element_text(face = "bold", size = 12*96/72),
-          legend.text = element_text(face = "bold", size = 14*96/72),
-          legend.title = element_text(face = "bold", size = 12*96/72))
-  assign(paste0("p.hallmark_", p), p.hallmark)
 }