FGSEA Analysis of Malignant CD4 T-cell States (all) in Sezary Syndrome

Overview

This R Markdown performs FGSEA-based pathway characterization of malignant CD4+ T-cell clusters, using the entire ranked differential expression output per cluster (not a top-N marker subset). This is a methodological correction relative to an earlier top100-based FGSEA attempt: FGSEA’s permutation-based null distribution requires a full ranked gene list to be statistically valid. Restricting FGSEA to a short marker list is the reason that earlier run was underpowered; over-representation analysis (ORA) remains the correct tool for short, fixed marker-gene lists, while FGSEA is the correct tool here because the entire transcriptome per cluster is being ranked.

The workflow:

1. Identify cluster-specific differential expression using Seurat FindAllMarkers, retaining all tested genes (no logfc.threshold/top-N filtering).
2. Rank genes per cluster using a composite statistic combining effect size and statistical confidence, with an explicit floor on p-values 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.

1. Load packages

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

set.seed(123)

2. Load Seurat object

Load your integrated malignant CD4+ T-cell Seurat object.

# Example:
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"

Idents(All_samples_Merged)

3. Differential expression analysis

All genes are retained for FGSEA.

No top-100 filtering is applied.

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

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


head(markers_all)

NA

4. Generate ranking statistic

The ranking combines expression effect size and statistical confidence:

\[ Score = sign(avg\_log2FC) \times |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. Replacing all of those with the single maximum finite rank (as in the original script) collapses many genes to an identical tied value, destroying the relative ordering FGSEA depends on. Instead, a numerical floor is applied to p_val_adj before the log transform, preserving relative differences between highly significant genes.

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

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

sum(is.infinite(markers_all$fgsea_rank))

[1] 0

sum(markers_all$p_val_adj == 0)

[1] 8266

5. Prepare pathway databases

library(msigdbr)

msigdbr_species <- "Homo sapiens"


# -------------------------
# Hallmark pathways
# -------------------------

hallmark_sets <- msigdbr(
  species = msigdbr_species,
  collection = "H"
)

hallmark_sets <- split(
  hallmark_sets$gene_symbol,
  hallmark_sets$gs_name
)



# -------------------------
# Reactome pathways
# -------------------------

reactome_sets <- msigdbr(
  species = msigdbr_species,
  collection = "C2",
  subcollection = "CP:REACTOME"
)

reactome_sets <- split(
  reactome_sets$gene_symbol,
  reactome_sets$gs_name
)



# -------------------------
# KEGG pathways
# -------------------------

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
)


# Check pathway numbers

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, the category/subcategory arguments have been renamed to collection/subcollection. If the chunk above errors with “unused argument”, replace category = with collection = and subcategory = with subcollection = throughout this section.

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
  )

  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_all$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_all, pathways = hallmark_sets))
hallmark_results$database <- "Hallmark"

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

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

go_results <- bind_rows(lapply(clusters, run_fgsea_cluster,
                                markers = markers_all, 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_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_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_Top10_Pathways_Per_Cluster.csv", row.names = FALSE)
top_pathways

13. Plot Hallmark enrichment per cluster

plot_data <- hallmark_results %>%
  dplyr::filter(padj < 0.05) %>%
  dplyr::group_by(cluster) %>%
  dplyr::slice_max(order_by = abs(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)") +
    theme_classic(base_size = 10)
} else {
  cat("No Hallmark pathways passed padj < 0.05 in any cluster.\n")
}

14. Save R objects

save(markers_all, all_fgsea_results, hallmark_results, reactome_results,
     kegg_results, go_results, file = "FGSEA_SS_Malignant_CD4_results.RData")

15. Notes on using this output for cluster naming justification

• FGSEA_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 sign and magnitude, not just significance: a large positive NES means the pathway’s genes are concentrated toward the top (up-regulated) end of that cluster’s ranking relative to all other clusters.
• This FGSEA approach and your existing ORA (enrichGO/enrichKEGG/enrichPathway) approach answer different questions and should be reported as complementary, not interchangeable: ORA asks “are my top marker genes overrepresented in a pathway?”; FGSEA asks “is a pathway’s gene set skewed toward the top of my full ranked list for this cluster?” Reporting both, when they agree, is stronger evidence than either alone.

16. Session Info

sessionInfo()

R version 4.5.3 (2026-03-11)
Platform: x86_64-redhat-linux-gnu
Running under: Rocky Linux 9.7 (Blue Onyx)

Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS-OPENMP;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8      
 [8] LC_NAME=C              LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: Europe/Paris
tzcode source: system (glibc)

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] msigdbr_26.1.0              fgsea_1.36.2                future_1.70.0               igraph_2.3.2                pheatmap_1.0.13             scales_1.4.0               
 [7] ggrepel_0.9.8               ggridges_0.5.7              viridis_0.6.5               viridisLite_0.4.3           RColorBrewer_1.1-3          SeuratWrappers_0.4.0       
[13] monocle3_1.4.26             SingleCellExperiment_1.32.0 SummarizedExperiment_1.40.0 GenomicRanges_1.62.1        Seqinfo_1.0.0               IRanges_2.44.0             
[19] S4Vectors_0.48.1            MatrixGenerics_1.22.0       matrixStats_1.5.0           Biobase_2.70.0              BiocGenerics_0.56.0         generics_0.1.4             
[25] CytoTRACE2_1.1.0            RSpectra_0.16-2             Rfast_2.1.5.2               RcppParallel_5.1.11-2       zigg_0.0.2                  Rcpp_1.1.1-1.1             
[31] plyr_1.8.9                  Matrix_1.7-5                magrittr_2.0.5              HiClimR_2.2.1               doParallel_1.0.17           iterators_1.0.14           
[37] foreach_1.5.2               data.table_1.18.4           patchwork_1.3.2             SCpubr_3.0.1                cowplot_1.2.0               lubridate_1.9.5            
[43] forcats_1.0.1               stringr_1.6.0               dplyr_1.2.1                 purrr_1.2.2                 readr_2.2.0                 tidyr_1.3.2                
[49] tibble_3.3.1                tidyverse_2.0.0             plotly_4.12.0               ggplot2_4.0.3               Seurat_5.5.0                SeuratObject_5.4.0         
[55] sp_2.2-1                   

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.23       splines_4.5.3          later_1.4.8            R.oo_1.27.1            polyclip_1.10-7        fastDummies_1.7.6      lifecycle_1.0.5       
  [8] Rdpack_2.6.6           globals_0.19.1         lattice_0.22-9         MASS_7.3-65            sass_0.4.10            limma_3.66.0           rmarkdown_2.31        
 [15] jquerylib_0.1.4        yaml_2.3.12            remotes_2.5.0          httpuv_1.6.17          otel_0.2.0             sctransform_0.4.3      spam_2.11-4           
 [22] spatstat.sparse_3.2-0  reticulate_1.46.0      pbapply_1.7-4          minqa_1.2.8            abind_1.4-8            Rtsne_0.17             presto_1.0.0          
 [29] R.utils_2.13.0         irlba_2.3.7            listenv_0.10.1         spatstat.utils_3.2-3   goftest_1.2-3          spatstat.random_3.5-0  fitdistrplus_1.2-6    
 [36] parallelly_1.47.0      ncdf4_1.24             codetools_0.2-20       DelayedArray_0.36.1    tidyselect_1.2.1       farver_2.1.2           lme4_2.0-1            
 [43] spatstat.explore_3.8-1 jsonlite_2.0.0         progressr_0.19.0       survival_3.8-6         systemfonts_1.3.2      tools_4.5.3            ragg_1.5.2            
 [50] ica_1.0-3              glue_1.8.1             gridExtra_2.3          SparseArray_1.10.10    xfun_0.58              withr_3.0.2            BiocManager_1.30.27   
 [57] fastmap_1.2.0          boot_1.3-32            digest_0.6.39          rsvd_1.0.5             timechange_0.4.0       R6_2.6.1               mime_0.13             
 [64] textshaping_1.0.5      scattermore_1.2        tensor_1.5.1           dichromat_2.0-0.1      spatstat.data_3.1-9    R.methodsS3_1.8.2      utf8_1.2.6            
 [71] httr_1.4.8             htmlwidgets_1.6.4      S4Arrays_1.10.1        uwot_0.2.4             pkgconfig_2.0.3        gtable_0.3.6           rsconnect_1.10.0      
 [78] lmtest_0.9-40          S7_0.2.2               XVector_0.50.0         htmltools_0.5.9        dotCall64_1.2          png_0.1-9              spatstat.univar_3.2-0 
 [85] reformulas_0.4.4       knitr_1.51             rstudioapi_0.19.0      tzdb_0.5.0             reshape2_1.4.5         curl_7.1.0             nlme_3.1-169          
 [92] nloptr_2.2.1           cachem_1.1.0           zoo_1.8-15             KernSmooth_2.23-26     miniUI_0.1.2           pillar_1.11.1          grid_4.5.3            
 [99] vctrs_0.7.3            RANN_2.6.2             promises_1.5.0         xtable_1.8-8           cluster_2.1.8.2        evaluate_1.0.5         cli_3.6.6             
[106] compiler_4.5.3         rlang_1.2.0            future.apply_1.20.2    labeling_0.4.3         stringi_1.8.7          BiocParallel_1.44.0    deldir_2.0-4          
[113] babelgene_22.9         assertthat_0.2.1       lazyeval_0.2.3         spatstat.geom_3.8-1    RcppHNSW_0.7.0         hms_1.1.4              statmod_1.5.2         
[120] shiny_1.13.0           rbibutils_2.4.1        ROCR_1.0-12            bslib_0.11.0           fastmatch_1.1-8