Cluster Enrichment Analysis - Validation of Cluster Annotations

1 Purpose

This script performs functional enrichment analysis (GO, KEGG, Reactome, Hallmark) on the top 100 marker genes per cluster (Supplementary Table S6) to independently validate the cluster names you assigned manually. If the enriched pathways align with your proposed name (e.g., “NK-like” cluster enriched for “natural killer cell mediated cytotoxicity”), the annotation is supported.

2 load libraries

3 Load Top100 Marker Gene Table

markers <- read_excel("Supplementary_Table_S6.xlsx")

# Expected columns: pval, avglog2FC, pct.1, pct.2, pvaladj, cluster, gene
markers <- markers %>%
  rename_with(tolower) %>%
  mutate(cluster = as.character(cluster))

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

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

4 Map Gene Symbols to Entrez IDs

gene_map <- bitr(unique(markers$gene),
                  fromType = "SYMBOL",
                  toType = "ENTREZID",
                  OrgDb = org.Hs.eg.db)

markers_mapped <- markers %>%
  left_join(gene_map, by = c("gene" = "SYMBOL")) %>%
  filter(!is.na(ENTREZID))

unmapped_genes <- setdiff(unique(markers$gene), gene_map$SYMBOL)
print(unmapped_genes)

5 Prepare Background Gene Universe

universe_entrez <- unique(markers_mapped$ENTREZID)

6 Define Enrichment Function

run_cluster_enrichment <- function(cluster_id, marker_df) {

  sub <- marker_df %>% filter(cluster == cluster_id)
  entrez_ids <- unique(sub$ENTREZID)

  message("Cluster ", cluster_id, ": ", length(entrez_ids), " mapped genes")

  results <- list()

  results$GO_BP <- tryCatch(
    enrichGO(gene = entrez_ids, universe = universe_entrez,
             OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP",
             pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2,
             readable = TRUE),
    error = function(e) NULL)

  results$KEGG <- tryCatch(
    enrichKEGG(gene = entrez_ids, universe = universe_entrez,
               organism = "hsa", pAdjustMethod = "BH",
               pvalueCutoff = 0.05, qvalueCutoff = 0.2),
    error = function(e) NULL)

  results$Reactome <- tryCatch(
    enrichPathway(gene = entrez_ids, universe = universe_entrez,
                   organism = "human", pAdjustMethod = "BH",
                   pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE),
    error = function(e) NULL)

  hallmark_sets <- msigdbr(species = "Homo sapiens", category = "H") %>%
    dplyr::select(gs_name, entrez_gene)

  results$Hallmark <- tryCatch(
    enricher(gene = entrez_ids, universe = universe_entrez,
             TERM2GENE = hallmark_sets, pAdjustMethod = "BH",
             pvalueCutoff = 0.05, qvalueCutoff = 0.2),
    error = function(e) NULL)

  return(results)
}

7 Run Enrichment for Every Cluster

all_results <- list()

for (cl in clusters_list) {
  all_results[[cl]] <- tryCatch(
    run_cluster_enrichment(cl, markers_mapped),
    error = function(e) {
      message("Cluster ", cl, " failed: ", e$message)
      NULL
    }
  )
}

8 Export Enrichment Tables to Excel

wb <- createWorkbook()

for (cl in names(all_results)) {
  res <- all_results[[cl]]
  if (is.null(res)) next

  for (ont in names(res)) {
    obj <- res[[ont]]
    if (is.null(obj) || nrow(as.data.frame(obj)) == 0) next

    sheet_name <- substr(paste0("C", cl, "_", ont), 1, 31)
    addWorksheet(wb, sheet_name)
    writeData(wb, sheet_name, as.data.frame(obj))
  }
}

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

9 Visualize Top Pathways per Cluster

for (cl in names(all_results)) {
  res <- all_results[[cl]]
  if (is.null(res)) next

  cat("\n## Cluster", cl, "\n")

  for (ont in c("GO_BP", "KEGG", "Reactome", "Hallmark")) {
    obj <- res[[ont]]
    if (is.null(obj) || nrow(as.data.frame(obj)) == 0) next

    cat("\n###", ont, "- Cluster", cl, "\n")
    p <- dotplot(obj, showCategory = 10) +
      ggtitle(paste("Cluster", cl, "-", ont)) +
      theme(axis.text.y = element_text(size = 9))
    print(p)
    ggsave(filename = paste0("cluster", cl, "_", ont, "_dotplot.png"),
           plot = p, width = 8, height = 6, dpi = 300)
  }
}

## Cluster 0 

### KEGG - Cluster 0 

## Cluster 1 

## Cluster 2 

## Cluster 3 

## Cluster 4 

## Cluster 5 

### Reactome - Cluster 5 

## Cluster 6 

## Cluster 7 

### GO_BP - Cluster 7 

### KEGG - Cluster 7 

### Reactome - Cluster 7 

### Hallmark - Cluster 7 

## Cluster 8 

## Cluster 9 

## Cluster 10 

## Cluster 11 

### GO_BP - Cluster 11 

### Reactome - Cluster 11 

### Hallmark - Cluster 11 

## Cluster 12 

### GO_BP - Cluster 12 

### Hallmark - Cluster 12 

## Cluster 13 

### GO_BP - Cluster 13 

### KEGG - Cluster 13 

### Reactome - Cluster 13 

### Hallmark - Cluster 13 

10 Cross-Check: Annotation vs Enrichment

Manual cross-check: proposed name vs. expected enriched terms

Cluster	Proposed_Name	Expected_Enriched_Terms
0	MHC-II high aberrant state	antigen processing and presentation, MHC class II, interferon response (STAT1)
1	NK-like cytotoxic	natural killer cell mediated cytotoxicity, NK receptor signalling, KIR family activity
2	Th2-like	Th2 cytokine production, IL-13 signalling, chemokine (CCL17) response
3	Naive/CD4 reference	T cell differentiation, naive T cell, TCR signalling
4	Migratory/Adhesion state	cell adhesion, extracellular matrix interaction, cell migration (ENPP2/autotaxin signalling)
5	Stem-like	stem cell population maintenance, chromatin remodeling (HMGA2), T-cell lineage commitment (RUNX1)
6	Th2-like (Activated)	T cell activation, costimulation (4-1BB/TNFRSF9), Th2 cytokine production
7	Cycling (G2/M)	mitotic nuclear division, kinetochore organization, G2/M cell cycle checkpoint
8	Metabolically reprogrammed	lipid biosynthesis, metabolic reprogramming, lymphocyte activation (CD38)
9	GZMA-cytotoxic	T cell mediated cytotoxicity, granzyme-mediated apoptotic signalling
10	Central memory/CD4 reference	T cell quiescence, central memory differentiation
11	Pro-inflammatory	inflammatory response, S100 alarmin signalling, TNF-family receptor signalling
12	GZMB-high inflammatory	granzyme-mediated cytotoxicity, inflammatory chemokine signalling, GM-CSF response
13	IFN stimulated	type I interferon signalling, response to virus, ISGylation

11 Notes

• Enrichment used each cluster’s top 100 marker genes (Supplementary Table S6).
• Background universe restricted to the union of all top100 genes across clusters.
• Clusters 3 and 10 are healthy CD4 T-cell reference populations, not malignant states.
• tryCatch() wraps every enrichment call so one failed database (e.g., no significant KEGG terms for a cluster) won’t stop the whole loop. ```