fgsea Analysis for ManuScript_Feb2026

1 load libraries

2 LOAD DATA & SETUP

# Make sure the file name matches exactly what you have
fg_all <- read.csv("fgsea_all_results.csv", stringsAsFactors = FALSE)

cat("Loaded", nrow(fg_all), "pathways from fgsea_all_results.csv\n")

Loaded 9200 pathways from fgsea_all_results.csv

# Expanded Proliferation Keywords
# STRICT EXCLUSION LIST (Updated)
prolif_terms <- c(
  "CELL_CYCLE", "MITOTIC", "G2M", "E2F", "SPINDLE", 
  "CHROMOSOME", "DNA_REPLICATION", "NUCLEAR_DIVISION",
  "ORGANELLE_FISSION", "KINETOCHORE", "CENTROSOME",
  "REPLICATION", "SEGREGATION", "DIVISION", "M_PHASE", 
  "KINESINS", "MEIOSIS", "OOCYTE", 
  "MICROTUBULE", "CYTOSKELETON", "TRAFFIC", "GOLGI", "CYCLIN",
  "RECOMBINATION", "REPAIR", "REPLICATIVE", "POLO_LIKE"
)

3 DATA PREPARATION FUNCTION

prepare_data <- function(fg_tbl, topN = 3, exclude_prolif = FALSE) {
  if (exclude_prolif) {
    fg_tbl <- fg_tbl %>% filter(!grepl(paste(prolif_terms, collapse = "|"), pathway, ignore.case = TRUE))
  }

  # Get top N up and down per database
  fg_plot <- bind_rows(
    fg_tbl %>% filter(dataset == "hallmark") %>% { bind_rows(slice_max(., NES, n = topN), slice_min(., NES, n = topN)) },
    fg_tbl %>% filter(dataset == "kegg") %>% { bind_rows(slice_max(., NES, n = topN), slice_min(., NES, n = topN)) },
    fg_tbl %>% filter(dataset == "reactome") %>% { bind_rows(slice_max(., NES, n = topN), slice_min(., NES, n = topN)) },
    fg_tbl %>% filter(dataset == "go_bp") %>% { bind_rows(slice_max(., NES, n = topN), slice_min(., NES, n = topN)) }
  )

  # Format Labels
  fg_plot %>%
    mutate(
      db_prefix = case_when(
        dataset == "hallmark" ~ "HALLMARK",
        dataset == "kegg" ~ "KEGG",
        dataset == "reactome" ~ "REACTOME",
        dataset == "go_bp" ~ "GOBP"
      ),
      clean_pathway = gsub("^HALLMARK_|^KEGG_|^REACTOME_|^GOBP_", "", pathway),
      plot_label = paste0(db_prefix, "_", clean_pathway)
    ) %>%
    arrange(NES) %>%
    mutate(plot_label = factor(plot_label, levels = unique(plot_label)))
}

4 PLOTTING FUNCTION

create_plot <- function(data, color_var, color_label, title_text) {
  ggplot(data, aes(x = NES, y = plot_label)) +
    geom_point(aes(shape = dataset, size = leadingEdgeCount, color = !!sym(color_var)), alpha = 0.9) +
    geom_vline(xintercept = 0, linetype = "solid", color = "gray80", linewidth = 0.5) +

    scale_color_gradientn(
      colors = c("red", "orange", "blue"),
      trans = "log10",
      name = color_label,
      guide = guide_colorbar(reverse = TRUE)
    ) +

    scale_shape_manual(
      values = c("hallmark" = 17, "kegg" = 15, "reactome" = 3, "go_bp" = 16),
      guide = "none"
    ) +

    scale_size_continuous(range = c(3, 8), name = "Leading edge genes") +

    theme_minimal() +
    labs(x = "Normalized Enrichment Score (NES)", y = NULL, title = title_text) +
    theme(
      axis.text.y = element_text(size = 14, face = "bold", color = "black"),
      # X-axis labels (Numbers) - BOLD & BIGGER
      axis.text.x = element_text(size = 10, face = "bold", color = "black"),
      # X-axis TITLE (The Text "Normalized Enrichment Score...") - BOLD
      axis.title.x = element_text(size = 14, face = "bold", color = "black", margin = margin(t = 10)),
      
      # Plot Title
      plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
      
      # Legend
      legend.position = "right",
      legend.box = "vertical",
      legend.title = element_text(face = "bold", size = 10),

      panel.grid.major.y = element_line(color = "gray95")
    )
}

5 GENERATE 4 PLOTS

# A) All Pathways (padj)
df_all <- prepare_data(fg_all, exclude_prolif = FALSE)
p1 <- create_plot(df_all, "padj", "FDR (padj)", "Global Pathway Alterations (Malignant vs. Normal CD4+)")
ggsave("Fig1_All_padj.png", p1, width = 16, height = 8, dpi = 300)
ggsave("Fig1_All_padj.pdf", p1, width = 16, height = 8)

# B) All Pathways (p-value)
p2 <- create_plot(df_all, "pval", "P-value", "Nominally Significant Pathway Alterations")
ggsave("Fig2_All_pval.png", p2, width = 16, height = 8, dpi = 300)
ggsave("Fig2_All_pval.pdf", p2, width = 16, height = 8)

# C) Non-Proliferation (padj)
df_no_prolif <- prepare_data(fg_all, exclude_prolif = TRUE)
p3 <- create_plot(df_no_prolif, "padj", "FDR (padj)", "Functional & Immune Signatures (Non-Proliferative)")
ggsave("Fig3_NonProlif_padj.png", p3, width = 16, height = 8, dpi = 300)
ggsave("Fig3_NonProlif_padj.pdf", p3, width = 16, height = 8)

# D) Non-Proliferation (p-value)
p4 <- create_plot(df_no_prolif, "pval", "P-value", "Exploratory Functional Signatures (Non-Proliferative)")
ggsave("Fig4_NonProlif_pval.png", p4, width = 16, height = 8, dpi = 300)
ggsave("Fig4_NonProlif_pval.pdf", p4, width = 16, height = 8)

print("Created 4 figures: Fig1..Fig4 (.png and .pdf)")

[1] "Created 4 figures: Fig1..Fig4 (.png and .pdf)"

p1

p2

p3

p4