fgsea Analysis for ManuScript_Feb2026-top10 significant pathways only

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

# 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", "CHECKPOINTS",
  "TRANSCRIPTION", "S_PHASE", "ANAPHASE", "TELOPHASE", "PROPHASE", 
  "CYTOKINESIS", "SPINDLE_ASSEMBLY", "SPINDLE_CHECKPOINT", 
  "MITOTIC_SPINDLE", "MITOTIC_CHECKPOINT", "MITOTIC_G1", 
  "MITOTIC_S_PHASE", "MITOTIC_G2M"
)

3 DATA PREPARATION FUNCTION

prepare_data <- function(fg_tbl, topN = 10, exclude_prolif = FALSE, 
                         sig_col = "padj", sig_cutoff = 0.05) {
  
  # 1. Filter by Significance FIRST
  # We select pathways where the chosen column (padj or pval) is less than the cutoff
  fg_tbl <- fg_tbl %>% filter(!!sym(sig_col) < sig_cutoff)
  
  # 2. Filter out proliferation (optional)
  if (exclude_prolif) {
    fg_tbl <- fg_tbl %>% filter(!grepl(paste(prolif_terms, collapse = "|"), pathway, ignore.case = TRUE))
  }
  
  # Check if empty after filtering
  if(nrow(fg_tbl) == 0) {
    warning("No pathways met the significance criteria!")
    return(NULL)
  }

  # 3. 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)) }
  )
  
  if(nrow(fg_plot) == 0) return(NULL)

  # 4. 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) {
  if(is.null(data)) return(NULL)
  
  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"),
      axis.text.x = element_text(size = 10, face = "bold", color = "black"),
      axis.title.x = element_text(size = 14, face = "bold", color = "black", margin = margin(t = 10)),
      plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
      legend.position = "right",
      legend.box = "vertical",
      legend.title = element_text(face = "bold", size = 10),
      panel.grid.major.y = element_line(color = "gray95")
    )
}

5 Generate Plots (Significance Filter Applied)

# A) All Pathways (ONLY significant by padj < 0.05)
df_all_padj <- prepare_data(fg_all, exclude_prolif = FALSE, sig_col = "padj", sig_cutoff = 0.05)
if(!is.null(df_all_padj)) {
  p1 <- create_plot(df_all_padj, "padj", "FDR (padj)", "Significant Pathways (FDR < 0.05)")
  print(p1)
  ggsave("Fig1_Sig_All_padj.png", p1, width = 20, height = 18, dpi = 300)
}

# B) All Pathways (ONLY significant by pval < 0.05)
df_all_pval <- prepare_data(fg_all, exclude_prolif = FALSE, sig_col = "pval", sig_cutoff = 0.05)
if(!is.null(df_all_pval)) {
  p2 <- create_plot(df_all_pval, "pval", "P-value", "Significant Pathways (Nominal P < 0.05)")
  print(p2)
  ggsave("Fig2_Sig_All_pval.png", p2, width = 20, height = 18, dpi = 300)
}

# C) Non-Proliferation (ONLY significant by padj < 0.05)
df_no_prolif_padj <- prepare_data(fg_all, exclude_prolif = TRUE, sig_col = "padj", sig_cutoff = 0.05)
if(!is.null(df_no_prolif_padj)) {
  p3 <- create_plot(df_no_prolif_padj, "padj", "FDR (padj)", "Significant Non-Prolif Pathways (FDR < 0.05)")
  print(p3)
  ggsave("Fig3_Sig_NonProlif_padj.png", p3, width = 20, height = 18, dpi = 300)
} else {
  message("No Non-Proliferation pathways met FDR < 0.05 criteria.")
}

# D) Non-Proliferation (ONLY significant by pval < 0.05)
df_no_prolif_pval <- prepare_data(fg_all, exclude_prolif = TRUE, sig_col = "pval", sig_cutoff = 0.05)
if(!is.null(df_no_prolif_pval)) {
  p4 <- create_plot(df_no_prolif_pval, "pval", "P-value", "Significant Non-Prolif Pathways (Nominal P < 0.05)")
  print(p4)
  ggsave("Fig4_Sig_NonProlif_pval.png", p4, width = 20, height = 18, dpi = 300)
}