GSEA tile plot script

Script/Figure_5b.R

+# Loads the necessary libraries for data analysis and visualization
+library(reshape2)
+library(ggplot2)
+library(dplyr)
+library(RColorBrewer)
+library(gplots)
+library(pheatmap)
+library(Seurat)
+
+# Loads the GSEA results from an RDS file
+gsea_res <- readRDS("/beegfs/scratch/ric.gregori/ric.gregori/Infertility/object/gsea_intracluster_OAT_NOA_vs_FER.rds")
+
+# Removes comparisons between NOA and OAT
+index<- c(3, 6 ,9 ,12, 15, 18, 21, 24, 27)
+gsea_res<-gsea_res[-index] 
+
+# Divides the results into two groups: OAT vs FER and NOA vs FER
+gsea_OAT_vs_FER<- gsea_res[-c(2,4,6,8,10,12,14,16,18)]
+gsea_NOA_vs_FER <- gsea_res[-c(1,3,5,7,9,11,13,15,17)]
+# Combines the two groups
+gsea_res <- append(gsea_OAT_vs_FER, gsea_NOA_vs_FER)
+
+# Initializes empty list and dataframe
+gsea_df_melt <- list()
+counter <- 1
+
+gsea_df <- NULL
+
+# Extracts GSEA results for each comparison
+for(i in names(gsea_res)){
+  print(i)
+  if(is.null(gsea_res[[i]])){
+    next()
+  }
+  else if (nrow(gsea_res[[i]]@result) == 0){
+    next
+  }
+  else{
+    gsea_df[[i]] <- gsea_res[[i]]@result
+  }
+}
+
+
+# Defines the correct order for the classifications
+right_classification <- c("CD4+ Tfh_OAT_vs_FER", "CD4+ Th2_OAT_vs_FER", "CD4+ Th1/Th17_OAT_vs_FER", "CD4+ CD8+_OAT_vs_FER", "CD4+ CD8+ TEM_OAT_vs_FER", "CD8+ TCM(Gzk+)_OAT_vs_FER", "CD8+ TE(TOX+)_OAT_vs_FER", "CD8+ TE (GzK+)_OAT_vs_FER", "CD8+ TEM_OAT_vs_FER","CD4+ Tfh_NOA_vs_FER", "CD4+ Th2_NOA_vs_FER", "CD4+ Th1/Th17_NOA_vs_FER", "CD4+ CD8+_NOA_vs_FER", "CD4+ CD8+ TEM_NOA_vs_FER", "CD8+ TCM(Gzk+)_NOA_vs_FER", "CD8+ TE(TOX+)_NOA_vs_FER", "CD8+ TE (GzK+)_NOA_vs_FER", "CD8+ TEM_NOA_vs_FER")
+
+
+# Gets the indices to reorder the dataframe
+order_indices <- match(right_classification, names(gsea_df))
+
+# Reorder the dataframe
+gsea_df<- gsea_df[order_indices]
+
+# Extracts information from each dataframe
+for (i in names(gsea_df)) {
+  cluster <- unlist(strsplit(x = i, split = "_"))[1]
+  comparison <- unlist(sub(x = i, "^.+?_", ""))
+  gsea_df[[i]]$cluster <- cluster
+  gsea_df[[i]]$comparison <- comparison
+  gsea_df[[i]]$ID <- gsub(gsea_df[[i]]$ID, pattern = "HALLMARK_", replacement = "")
+  gsea_df_melt[[i]] <- melt(data = gsea_df[[i]], id.vars = c("ID","p.adjust","comparison", "cluster"),  measure.vars = c("NES"))
+  gsea_df_melt[[i]]$value <- round(gsea_df_melt[[i]]$value, digits = 2)
+  counter <- counter + 1
+}
+
+# Combines all dataframes into a single dataframe
+total <- bind_rows(gsea_df, .id = "Classification")
+total_melt <- bind_rows(gsea_df_melt, .id = "Classification")
+
+#min(total_melt$value, na.rm = T)#-2.4
+#max(total_melt$value, na.rm = T)#2.37
+
+# Defines the levels of the factors for ordering variables
+total_melt$ID <- factor(total_melt$ID, levels = term_right)
+total_melt$Classification <- factor(total_melt$Classification, levels = unique(right_classification))
+
+# Transforms the data into a format suitable for visualization
+total_cast_melt <- melt(dcast(data = total_melt, formula = ID ~ Classification, value.var = 'value'), id.vars = 'ID' ) 
+
+# Extracts the cluster and comparison from the column labels
+total_cast_melt$cluster <- stringr::str_split_fixed(total_cast_melt$variable, "_", 2)[,1]
+total_cast_melt$comparison <- stringr::str_split_fixed(total_cast_melt$variable, "_", 2)[,2]
+
+# Defines the GSEA terms and their groups
+
+term_right<-c("NOTCH_SIGNALING",
+              "MTORC1_SIGNALING",
+              "KRAS_SIGNALING_DN",
+              "WNT_BETA_CATENIN_SIGNALING",
+              "TGF_BETA_SIGNALING",
+              "PI3K_AKT_MTOR_SIGNALING",
+              "EPITHELIAL_MESENCHYMAL_TRANSITION",
+              "MYOGENESIS",
+              "APICAL_JUNCTION",
+              
+              "COAGULATION",
+              "COMPLEMENT",
+              "ALLOGRAFT_REJECTION",
+              "IL2_STAT5_SIGNALING",
+              "IL6_JAK_STAT3_SIGNALING",
+              "INFLAMMATORY_RESPONSE",
+              "INTERFERON_GAMMA_RESPONSE",
+              "INTERFERON_ALPHA_RESPONSE",
+              "TNFA_SIGNALING_VIA_NFKB",
+              
+              "G2M_CHECKPOINT",
+              "MITOTIC_SPINDLE",
+              "MYC_TARGETS_V2",
+              "MYC_TARGETS_V1",
+              
+              "APOPTOSIS",
+              "UV_RESPONSE_DN",
+              "UV_RESPONSE_UP",
+              "P53_PATHWAY",
+              "UNFOLDED_PROTEIN_RESPONSE",
+              "PROTEIN_SECRETION",
+              
+              "GLYCOLYSIS",
+              "HYPOXIA",
+              "CHOLESTEROL_HOMEOSTASIS",
+              "ESTROGEN_RESPONSE_EARLY",
+              "ESTROGEN_RESPONSE_LATE",
+              "FATTY_ACID_METABOLISM",
+              "BILE_ACID_METABOLISM")
+
+
+group_letter <- c(rep("A_CellSignaling", 9), rep("B_ImmuneResponse", 9), rep("C_CellCycle", 4), rep("D_CellStress", 6), rep("E_Metabolism", 7))
+terms_group <- data.frame(term_right = term_right, group_letter = group_letter)
+
+# Merges the data with the information about the GSEA term groups
+total_cast_melt <- merge(total_cast_melt, terms_group, by.x = "ID", by.y = "term_right", all.x = T)
+# Define the order of comparison
+total_cast_melt$comparison <- factor(total_cast_melt$comparison, levels = c("OAT_vs_FER", "NOA_vs_FER"))
+# Define the order of clusters
+total_cast_melt$cluster <- factor(total_cast_melt$cluster, levels = c("CD4+ Tfh","CD4+ Th2","CD4+ Th1/Th17","CD4+ CD8+","CD4+ CD8+ TEM","CD8+ TCM(Gzk+)","CD8+ TE(TOX+)","CD8+ TE (GzK+)","CD8+ TEM"))
+# Define the order of groups (pathways)
+total_cast_melt$group_letter <- factor(total_cast_melt$group_letter, levels = c("E_Metabolism","D_CellStress","C_CellCycle","B_ImmuneResponse","A_CellSignaling"))
+
+#sort(unique(total_cast_melt$group_letter))
+
+# Create the GSEA tile plot using ggplot2
+gsea_tile_plot <-ggplot(data = total_cast_melt, mapping = aes(x = cluster, y = ID))+  
+  facet_grid(group_letter~comparison, scales = "free", space = "free") + 
+  geom_tile(aes(fill=value), width=0.80, height=0.80, color="white",lwd = 0.1,linetype = 5, size=0.20) + labs(fill = "NES") + 
+  ggtitle(label = "GSEA hallmarks", subtitle = " Intracluster -  padj < 0.1")+
+  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=10),
+        strip.text.x = element_text(size = 16),
+        strip.text.y = element_text(size =0),
+        strip.background = element_blank(),
+        plot.title = element_text(hjust = 0, size = 18),
+        plot.subtitle = element_text(hjust = 0, size = 12),
+        axis.ticks = element_blank(),
+        axis.text.y = element_text(size = 9),
+        axis.title = element_blank(),
+        panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        panel.spacing.x = unit(.3, "cm"),
+        panel.spacing.y = unit(.3,"cm"),
+        panel.background = element_rect(fill = "white")) +
+  scale_fill_gradient2(na.value = '#ECECED', 
+                       limits=c(-2.5, 2.5), breaks=seq(-3,3,by=1),
+                       low='#005295',high='#770000', mid = 'white') # +
+
+# Save the plot to a file
+ggsave(filename = "GSEA_tile_plot.png", width = 7, height = 6, plot =gsea_tile_plot)
\ No newline at end of file