FGSEA Analysis of Malignant CD4 T-cell States (Positive Markers Only) in Sezary Syndrome

Overview

This R Markdown performs FGSEA-based pathway characterization of malignant CD4+ T-cell clusters using positive markers only (Version A).

Rationale: for cluster annotation, the relevant question is “what biological programs define this cluster?”, not “which pathways are increased or decreased relative to other clusters?”. Using only genes enriched within each cluster (FindAllMarkers(only.pos = TRUE)), ranked by an unsigned statistic combining effect size and significance, keeps FGSEA focused on defining/up-regulated programs and substantially reduces the “unbalanced positive/negative gene-level statistic” warnings that occur when a signed ranking is dominated by one-sided pathways (e.g. interferon-stimulated genes, which are almost universally positive in an IFN-high cluster).

The workflow:

1. Identify cluster-specific positive markers using Seurat FindAllMarkers(only.pos = TRUE), retaining all tested positive genes (no top-N filtering, since FGSEA requires a full ranked list).
2. Rank genes per cluster using: \[ Score = avg\_log2FC \times \sqrt{-log10(p_{adj})} \] with an explicit floor on p_val_adj to avoid infinite ranks.
3. Run FGSEA against Hallmark, Reactome, KEGG, and GO:BP gene sets.
4. Export enrichment tables including NES, p-value, adjusted p-value, and leading-edge genes – the genes actually driving each enrichment signal.
5. Sanity-check known clusters (13 = IFN, 7 = G2M) and visualize top pathways.

The resulting evidence chain for cluster naming is: FindAllMarkers(only.pos=TRUE) -> defining marker genes -> FGSEA on those genes -> defining pathways -> leading-edge genes link the pathway back to the specific markers driving it. This is the primary evidence to cite in the manuscript’s cluster annotation table.

1. Load packages

library(Seurat)
library(dplyr)
library(fgsea)
library(msigdbr)
library(ggplot2)
library(tibble)
library(stringr)

set.seed(123)

2. Load Seurat object

All_samples_Merged <- readRDS("../../../1-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_Renamed_Clusters_Cell_state-03-12-2025.rds.rds")

All_samples_Merged
DefaultAssay(All_samples_Merged)

Idents(All_samples_Merged) <- "seurat_clusters"

3. Differential expression analysis (positive markers only)

All positive genes are retained for FGSEA – no top-N filtering, since FGSEA requires a full ranked list rather than a short marker set.

markers_pos <- FindAllMarkers(
  object = All_samples_Merged,
  only.pos = TRUE,
  min.pct = 0.1,
  logfc.threshold = 0
)

head(markers_pos)

4. Generate ranking statistic

The ranking combines expression effect size and statistical confidence, without a sign term since only positive markers are retained:

\[ Score = avg\_log2FC \times \sqrt{-log10(p_{adj})} \]

Fix applied: in large single-cell datasets, p_val_adj is frequently exactly 0 due to very large cell numbers, which makes -log10(p_val_adj) = Inf for many genes simultaneously. A numerical floor is applied to p_val_adj before the log transform, preserving relative differences between highly significant genes rather than tying them all to one collapsed maximum value.

p_floor <- min(markers_pos$p_val_adj[markers_pos$p_val_adj > 0], na.rm = TRUE) / 10

markers_pos <- markers_pos %>%
  mutate(
    p_val_adj_floored = ifelse(p_val_adj == 0, p_floor, p_val_adj),
    fgsea_rank = avg_log2FC * sqrt(-log10(p_val_adj_floored))
  )

sum(is.infinite(markers_pos$fgsea_rank))

[1] 0

sum(markers_pos$p_val_adj == 0)

[1] 5175

5. Prepare pathway databases

library(msigdbr)

msigdbr_species <- "Homo sapiens"

# Hallmark
hallmark_sets <- msigdbr(species = msigdbr_species, collection = "H")
hallmark_sets <- split(hallmark_sets$gene_symbol, hallmark_sets$gs_name)

# Reactome
reactome_sets <- msigdbr(species = msigdbr_species, collection = "C2",
                          subcollection = "CP:REACTOME")
reactome_sets <- split(reactome_sets$gene_symbol, reactome_sets$gs_name)

# KEGG (legacy, equivalent to the old "CP:KEGG")
kegg_sets <- msigdbr(species = msigdbr_species, collection = "C2",
                      subcollection = "CP:KEGG_LEGACY")
kegg_sets <- split(kegg_sets$gene_symbol, kegg_sets$gs_name)

# GO Biological Process
go_sets <- msigdbr(species = msigdbr_species, collection = "C5",
                    subcollection = "GO:BP")
go_sets <- split(go_sets$gene_symbol, go_sets$gs_name)

length(hallmark_sets)

[1] 50

length(reactome_sets)

[1] 1839

length(kegg_sets)

[1] 186

length(go_sets)

[1] 7538

Note: if your installed msigdbr version is < 10.0, replace collection = with category = and subcollection = with subcategory = throughout this chunk, and use "CP:KEGG" instead of "CP:KEGG_LEGACY".

6. FGSEA function

run_fgsea_cluster <- function(cluster_id, markers, pathways) {

  message("Running cluster: ", cluster_id)

  cluster_markers <- markers %>% dplyr::filter(cluster == cluster_id)

  ranks <- cluster_markers$fgsea_rank
  names(ranks) <- cluster_markers$gene

  ranks <- ranks[!duplicated(names(ranks))]
  ranks <- ranks[!is.na(ranks) & !is.infinite(ranks)]
  ranks <- sort(ranks, decreasing = TRUE)

  fgsea_result <- fgsea(
    pathways = pathways,
    stats = ranks,
    minSize = 15,
    maxSize = 500,
    eps = 0,
    nPermSimple = 10000
  )

  fgsea_result <- as.data.frame(fgsea_result)

  fgsea_result$leadingEdge_genes <- sapply(
    fgsea_result$leadingEdge,
    function(x) paste(x, collapse = ", ")
  )

  fgsea_result$cluster <- cluster_id
  fgsea_result$n_genes_ranked <- length(ranks)

  return(fgsea_result)
}

7. Define clusters

clusters <- sort(unique(markers_pos$cluster))
clusters

 [1] 0  1  2  3  4  5  6  7  8  9  10 11 12 13
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13

8. Run FGSEA across all four databases

hallmark_results <- bind_rows(lapply(clusters, run_fgsea_cluster,
                                      markers = markers_pos, pathways = hallmark_sets))
hallmark_results$database <- "Hallmark"

reactome_results <- bind_rows(lapply(clusters, run_fgsea_cluster,
                                      markers = markers_pos, pathways = reactome_sets))
reactome_results$database <- "Reactome"

kegg_results <- bind_rows(lapply(clusters, run_fgsea_cluster,
                                  markers = markers_pos, pathways = kegg_sets))
kegg_results$database <- "KEGG"

go_results <- bind_rows(lapply(clusters, run_fgsea_cluster,
                                markers = markers_pos, pathways = go_sets))
go_results$database <- "GO_BP"

9. Combine all databases

all_fgsea_results <- bind_rows(hallmark_results, reactome_results, kegg_results, go_results)

all_fgsea_results <- all_fgsea_results %>%
  dplyr::select(cluster, database, pathway, NES, pval, padj, size,
                n_genes_ranked, leadingEdge_genes) %>%
  dplyr::arrange(cluster, padj)

head(all_fgsea_results)

10. Export complete FGSEA table

write.csv(all_fgsea_results, "FGSEA_PositiveMarkers_All_Pathways_with_LeadingEdge_Genes.csv",
          row.names = FALSE)

11. Significant pathways

significant_fgsea <- all_fgsea_results %>% dplyr::filter(padj < 0.05)
write.csv(significant_fgsea, "FGSEA_PositiveMarkers_Significant_Pathways.csv", row.names = FALSE)

12. Top pathways per cluster

top_pathways <- all_fgsea_results %>%
  dplyr::group_by(cluster) %>%
  dplyr::arrange(padj) %>%
  dplyr::slice_head(n = 10)

write.csv(top_pathways, "FGSEA_PositiveMarkers_Top10_Pathways_Per_Cluster.csv", row.names = FALSE)
top_pathways

13. Sanity check on known clusters

Expected: Cluster 13 should show Hallmark INTERFERON_ALPHA_RESPONSE / INTERFERON_GAMMA_RESPONSE with leading-edge genes IFIT2, OASL, ISG15, IFIT3, MX1; Cluster 7 should show G2M_CHECKPOINT, E2F_TARGETS, MITOTIC_SPINDLE.

hallmark_results %>%
  dplyr::filter(cluster == 13) %>%
  dplyr::arrange(padj) %>%
  dplyr::select(pathway, NES, padj, leadingEdge_genes) %>%
  head(10)

hallmark_results %>%
  dplyr::filter(cluster == 7) %>%
  dplyr::arrange(padj) %>%
  dplyr::select(pathway, NES, padj, leadingEdge_genes) %>%
  head(10)

14. Plot Hallmark enrichment per cluster

plot_data <- hallmark_results %>%
  dplyr::filter(padj < 0.05) %>%
  dplyr::group_by(cluster) %>%
  dplyr::slice_max(order_by = NES, n = 5)

if (nrow(plot_data) > 0) {
  ggplot(plot_data, aes(x = reorder(pathway, NES), y = NES, fill = NES > 0)) +
    geom_col() +
    facet_wrap(~cluster, scales = "free") +
    coord_flip() +
    scale_fill_manual(values = c("TRUE" = "#B2182B", "FALSE" = "#4393C3"), guide = "none") +
    labs(x = NULL, y = "Normalized Enrichment Score (NES)",
         title = "FGSEA on positive markers only (Version A)") +
    theme_classic(base_size = 10)
} else {
  cat("No Hallmark pathways passed padj < 0.05 in any cluster.\n")
}

15. Save R objects

save(markers_pos, all_fgsea_results,
     hallmark_results, reactome_results, kegg_results, go_results,
     file = "FGSEA_SS_Malignant_CD4_PositiveMarkers_results.RData")

16. Notes on using this output for cluster naming justification

• FGSEA_PositiveMarkers_All_Pathways_with_LeadingEdge_Genes.csv is the primary file to cite from, since it links each pathway’s NES/FDR directly to the specific leading-edge genes driving the signal (e.g. “driven by IFIT2, IFIT3, OASL”).
• Report NES magnitude alongside FDR, not significance alone: a larger NES means the pathway’s genes are more strongly concentrated toward the top of that cluster’s positive-marker ranking.
• This positive-marker FGSEA approach and your existing ORA (enrichGO/enrichKEGG/enrichPathway) approach answer closely related but distinct questions and should be reported as complementary: ORA asks “are my top-N marker genes overrepresented in a pathway?”; FGSEA asks “is a pathway’s gene set skewed toward the top of my full positive-marker ranking for this cluster?” Agreement between the two provides stronger support for a cluster name than either method alone.
• The evidence chain to report per cluster: marker genes (from FindAllMarkers(only.pos=TRUE)) -> defining pathway (from FGSEA) -> leading-edge genes (confirming which markers actually drive the pathway signal). Example: Cluster 13 marker genes IFIT2/OASL/ISG15 -> Hallmark INTERFERON_ALPHA_RESPONSE/Reactome interferon signaling -> annotation “IFN-stimulated malignant state”.

17. Session Info

sessionInfo()