fgsea on Top100 Marker Genes per Cluster

Hallmark + GO:BP + KEGG + Reactome, ranked by avg_log2FC, validated via NES dotplots

1 Install Missing Packages (run once, outside knit)

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install(c("fgsea", "org.Hs.eg.db"), update = FALSE, ask = FALSE)

install.packages(c("readxl", "dplyr", "tidyr", "ggplot2", "msigdbr",
                    "tibble", "stringr"))

2 Load Libraries

library(readxl)
library(dplyr)
library(tidyr)
library(tibble)
library(fgsea)
library(msigdbr)
library(ggplot2)
library(stringr)
library(purrr)
library(openxlsx)

3 Load Top100 Marker Table

markers <- read_excel("../Supplementary_Table_S6.xlsx") %>%
  rename_with(tolower) %>%
  mutate(cluster = as.character(cluster))

# Fix known outdated/renamed symbols
symbol_updates <- c("QARS" = "QARS1", "CARS" = "CARS1", "WARS" = "WARS1")
markers <- markers %>%
  mutate(gene = ifelse(gene %in% names(symbol_updates), symbol_updates[gene], gene)) %>%
  filter(gene != "46083.0")

clusters_list <- as.character(sort(as.numeric(unique(markers$cluster))))
cat("Clusters found (numeric order):", paste(clusters_list, collapse = ", "), "\n")

Clusters found (numeric order): 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 

4 Load Gene Set Databases (Hallmark + GO:BP)

hallmark_sets <- msigdbr(species = "Homo sapiens", collection = "H") %>%
  distinct(gs_name, gene_symbol) %>%
  split(x = .$gene_symbol, f = .$gs_name)

go_bp_sets <- msigdbr(species = "Homo sapiens", collection = "C5", subcollection = "GO:BP") %>%
  distinct(gs_name, gene_symbol) %>%
  split(x = .$gene_symbol, f = .$gs_name)

kegg_df <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:KEGG_LEGACY")
if (nrow(kegg_df) == 0) {
  kegg_df <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:KEGG")
}
kegg_sets <- kegg_df %>%
  distinct(gs_name, gene_symbol) %>%
  split(x = .$gene_symbol, f = .$gs_name)

reactome_sets <- msigdbr(species = "Homo sapiens", collection = "C2", subcollection = "CP:REACTOME") %>%
  distinct(gs_name, gene_symbol) %>%
  split(x = .$gene_symbol, f = .$gs_name)

5 Build Ranked Top100 List per Cluster and Run fgsea

run_top100_fgsea <- function(cluster_id, marker_df, pathway_sets, min_size = 3, max_size = 100) {

  sub <- marker_df %>%
    filter(cluster == cluster_id) %>%
    distinct(gene, avg_log2fc) %>%
    group_by(gene) %>%
    summarise(avg_log2fc = mean(avg_log2fc), .groups = "drop") %>%
    arrange(desc(avg_log2fc))

  ranks <- setNames(sub$avg_log2fc, sub$gene)

  set.seed(42)
  ranks <- ranks + rnorm(length(ranks), mean = 0, sd = 1e-6)
  ranks <- sort(ranks, decreasing = TRUE)

  res <- fgsea(pathways = pathway_sets,
               stats = ranks,
               minSize = min_size,
               maxSize = max_size,
               eps = 0,
               scoreType = "pos")

  res <- res %>% arrange(padj)
  list(ranks = ranks, result = res)
}

fgsea_hallmark <- list()
fgsea_gobp <- list()
fgsea_kegg <- list()
fgsea_reactome <- list()

for (cl in clusters_list) {
  message("Running top100 fgsea for cluster ", cl)

  fgsea_hallmark[[cl]] <- tryCatch(
    run_top100_fgsea(cl, markers, hallmark_sets),
    error = function(e) { message("  Hallmark failed: ", e$message); NULL })

  fgsea_gobp[[cl]] <- tryCatch(
    run_top100_fgsea(cl, markers, go_bp_sets),
    error = function(e) { message("  GO:BP failed: ", e$message); NULL })

  fgsea_kegg[[cl]] <- tryCatch(
    run_top100_fgsea(cl, markers, kegg_sets),
    error = function(e) { message("  KEGG failed: ", e$message); NULL })

  fgsea_reactome[[cl]] <- tryCatch(
    run_top100_fgsea(cl, markers, reactome_sets),
    error = function(e) { message("  Reactome failed: ", e$message); NULL })
}

6 Export Results to Excel

library(openxlsx)
wb <- createWorkbook()

export_map <- list(Hallmark = fgsea_hallmark, GOBP = fgsea_gobp,
                    KEGG = fgsea_kegg, Reactome = fgsea_reactome)

for (cl in clusters_list) {
  for (db in names(export_map)) {
    obj <- export_map[[db]][[cl]]
    if (is.null(obj)) next

    df <- obj$result %>%
      mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")) %>%
      as.data.frame()
    if (nrow(df) == 0) next

    sheet_name <- substr(paste0("C", cl, "_", db), 1, 31)
    addWorksheet(wb, sheet_name)
    writeData(wb, sheet_name, df)
  }
}

saveWorkbook(wb, "Cluster_Top100_fgsea_Results.xlsx", overwrite = TRUE)

7 Combine All Clusters into One Long Dataframe (for faceted dotplot)

build_long_df <- function(fgsea_list, db_name) {
  purrr::map_dfr(names(fgsea_list), function(cl) {
    obj <- fgsea_list[[cl]]
    if (is.null(obj) || nrow(obj$result) == 0) return(NULL)
    obj$result %>%
      mutate(cluster = cl, database = db_name) %>%
      select(cluster, database, pathway, NES, padj, size)
  })
}

library(purrr)
long_hallmark <- build_long_df(fgsea_hallmark, "Hallmark")
long_gobp <- build_long_df(fgsea_gobp, "GOBP")
long_kegg <- build_long_df(fgsea_kegg, "KEGG")
long_reactome <- build_long_df(fgsea_reactome, "Reactome")

all_long <- bind_rows(long_hallmark, long_gobp, long_kegg, long_reactome) %>%
  mutate(cluster = factor(cluster, levels = clusters_list))

write.csv(all_long, "Cluster_Top100_fgsea_long.csv", row.names = FALSE)

8 Dotplot per Cluster Based on NES (Top Pathways)

plot_cluster_dotplot <- function(cluster_id, long_df, n_top = 10, sig_cutoff = 0.05) {

  df <- long_df %>%
    filter(cluster == cluster_id, padj < sig_cutoff) %>%
    arrange(desc(abs(NES))) %>%
    slice_head(n = n_top) %>%
    mutate(pathway = str_replace_all(pathway, "_", " "),
           pathway = str_wrap(pathway, width = 40))

  if (nrow(df) == 0) {
    message("No significant pathways for cluster ", cluster_id)
    return(NULL)
  }

  p <- ggplot(df, aes(x = NES, y = reorder(pathway, NES),
                       size = size, color = -log10(padj))) +
    geom_point() +
    scale_color_gradient(low = "#F4A582", high = "#B2182B", name = "-log10(padj)") +
    scale_size_continuous(name = "Gene set size", range = c(3, 9)) +
    geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
    labs(title = paste0("Cluster ", cluster_id, " - Top100 fgsea (NES-ranked)"),
         x = "Normalized Enrichment Score (NES)", y = NULL) +
    theme_minimal(base_size = 11) +
    theme(axis.text.y = element_text(size = 9))

  print(p)
  ggsave(filename = paste0("cluster", cluster_id, "_top100_NES_dotplot.png"),
         plot = p, width = 9, height = 6, dpi = 300)
}

for (cl in clusters_list) {
  cat("\n## Cluster", cl, "- NES Dotplot\n")
  plot_cluster_dotplot(cl, all_long)
}

8.1 Cluster 0 - NES Dotplot

8.2 Cluster 1 - NES Dotplot

8.3 Cluster 2 - NES Dotplot

8.4 Cluster 3 - NES Dotplot

8.5 Cluster 4 - NES Dotplot

8.6 Cluster 5 - NES Dotplot

8.7 Cluster 6 - NES Dotplot

8.8 Cluster 7 - NES Dotplot

8.9 Cluster 8 - NES Dotplot

8.10 Cluster 9 - NES Dotplot

8.11 Cluster 10 - NES Dotplot

8.12 Cluster 11 - NES Dotplot

8.13 Cluster 12 - NES Dotplot

8.14 Cluster 13 - NES Dotplot

9 Cross-Check: Proposed Name vs Top NES Pathway

proposed_names <- c(
  "0" = "MHC-II high aberrant state", "1" = "NK-like cytotoxic", "2" = "Th2-like",
  "3" = "Naive/CD4 reference", "4" = "Migratory/Adhesion state", "5" = "Stem-like",
  "6" = "Th2-like (Activated)", "7" = "Cycling (G2/M)", "8" = "Metabolically reprogrammed",
  "9" = "GZMA-cytotoxic", "10" = "Central memory/CD4 reference", "11" = "Pro-inflammatory",
  "12" = "GZMB-high inflammatory", "13" = "IFN stimulated"
)

validation_table <- all_long %>%
  filter(padj < 0.05) %>%
  group_by(cluster) %>%
  arrange(desc(abs(NES))) %>%
  slice_head(n = 3) %>%
  summarise(Top3_Pathways_NES = paste0(pathway, " (NES=", round(NES, 2), ")", collapse = "; "),
            .groups = "drop") %>%
  mutate(Proposed_Name = proposed_names[as.character(cluster)]) %>%
  select(cluster, Proposed_Name, Top3_Pathways_NES)

knitr::kable(validation_table, caption = "Proposed name vs. top NES-ranked pathways (top100 fgsea)")

Proposed name vs. top NES-ranked pathways (top100 fgsea)

cluster	Proposed_Name	Top3_Pathways_NES
1	NK-like cytotoxic	GOBP_DEFENSE_RESPONSE_TO_OTHER_ORGANISM (NES=2.83); REACTOME_IMMUNOREGULATORY_INTERACTIONS_BETWEEN_A_LYMPHOID_AND_A_NON_LYMPHOID_CELL (NES=2.18); KEGG_CHEMOKINE_SIGNALING_PATHWAY (NES=2.15)
12	GZMB-high inflammatory	GOBP_RESPONSE_TO_TUMOR_NECROSIS_FACTOR (NES=2.77); GOBP_POSITIVE_REGULATION_OF_LOCOMOTION (NES=2.76); GOBP_POSITIVE_REGULATION_OF_CHEMOTAXIS (NES=2.64)
13	IFN stimulated	GOBP_ANTIVIRAL_INNATE_IMMUNE_RESPONSE (NES=2.33); KEGG_RIG_I_LIKE_RECEPTOR_SIGNALING_PATHWAY (NES=2.24); GOBP_DEFENSE_RESPONSE_TO_OTHER_ORGANISM (NES=2.22)

10 Notes

• fgsea was run using only the top100 marker genes per cluster, ranked by avg_log2FC, per your request. This differs from standard fgsea practice (which uses the full ranked transcriptome), so results should be interpreted as a directional validation signal rather than a fully powered enrichment test.
• maxSize = 100 ensures no gene set larger than the ranked list itself is tested, avoiding meaningless overlaps.
• A tiny jitter (sd = 1e-6) breaks exact ties in avg_log2FC, required by fgsea’s ranking algorithm.
• Dotplots show the top 10 significant pathways (padj < 0.05) per cluster, sized by gene set overlap and colored by significance, x-axis positioned by NES (positive = enriched toward top-ranked/most upregulated genes in that cluster).
• If a cluster shows “No significant pathways” in the dotplot loop, this typically means the top100 genes for that cluster are too sparse or heterogeneous relative to Hallmark/GO:BP gene set sizes – worth checking with a relaxed padj cutoff (e.g., 0.25) before concluding the name lacks support.