Single-Cell Metabolic Analysis of Sézary Syndrome Heterogeneity

1 Introduction

Sézary syndrome (SS) is an aggressive leukemic variant of cutaneous T-cell lymphoma (CTCL) characterized by metabolic reprogramming. Recent studies show that malignant Sézary cells upregulate oxidative phosphorylation (OXPHOS), unlike most cancers that rely on aerobic glycolysis [web:104][web:105]. This analysis uses scMetabolism to characterize metabolic heterogeneity at single-cell resolution.

Key metabolic features in Sézary syndrome: - Elevated OXPHOS in malignant clones [web:104] - Altered glycolysis and glutaminolysis [web:104] - Fatty acid metabolism dysregulation [web:110][web:113] - Oxidative stress and ROS production [web:110][web:113]

Tool: scMetabolism GitHub

2 load libraries

# Data Processing
library(dplyr)
library(Seurat)
library(tibble)
library(tidyr)
library(stringr)

# Visualization
library(ggplot2)
library(ComplexHeatmap)
library(patchwork)
library(circlize)
library(viridis)

# Metabolic Analysis
library(scMetabolism)
library(rsvd)

# Set theme
theme_set(theme_classic())

3 Load Seurat Object

# Load your integrated Sézary syndrome object
SS <- readRDS("../../1-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_Renamed_Clusters_Cell_state-03-12-2025.rds.rds")

# Set identity to clusters
Idents(SS) <- "seurat_clusters"

# Summary
cat("Total cells:", ncol(SS), "\n")

Total cells: 49305 

cat("Total features:", nrow(SS), "\n")

Total features: 36601 

cat("Clusters:", length(unique(SS$seurat_clusters)), "\n")

Clusters: 14 

cat("Samples:", length(unique(SS$orig.ident)), "\n")

Samples: 9 

4 Prepare Seurat Object for scMetabolism

# 1. Create v3 assay from RNA counts
rna_counts <- GetAssayData(SS, assay = "RNA", layer = "counts")
RNA_v3 <- CreateAssayObject(counts = rna_counts)

# 2. Add the v3 assay with a temporary name
SS[["RNA_v3"]] <- RNA_v3

# 3. Switch default assay away from RNA
DefaultAssay(SS) <- "RNA_v3"

# 4. Delete the old RNA Assay5
SS[["RNA"]] <- NULL

# 5. Create RNA from RNA_v3
SS[["RNA"]] <- SS[["RNA_v3"]]

# 6. Switch default to RNA FIRST (before deleting RNA_v3!)
DefaultAssay(SS) <- "RNA"

# 7. Now delete RNA_v3
SS[["RNA_v3"]] <- NULL

# Verify the setup
Assays(SS)                    # Should show RNA (v3 version)

[1] "ADT"                          "prediction.score.celltype.l1" "prediction.score.celltype.l2" "prediction.score.celltype.l3" "SCT"                         
[6] "RNA"                         

class(SS[["RNA"]])            # Should be "Assay" not "Assay5"

[1] "Assay"
attr(,"package")
[1] "SeuratObject"

DefaultAssay(SS)              # Should be "RNA"

[1] "RNA"

# Run scMetabolism
SS <- sc.metabolism.Seurat(obj = SS,
                           method = "AUCell",
                           imputation = FALSE,
                           ncores = 4,
                           metabolism.type = "KEGG")

Your choice is: KEGG
Start quantify the metabolism activity...

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

# Check results
dim(SS@assays$METABOLISM$score)

[1]    85 49305

head(rownames(SS@assays$METABOLISM$score), 10)

 [1] "Glycolysis / Gluconeogenesis"                "Citrate cycle (TCA cycle)"                   "Pentose phosphate pathway"                  
 [4] "Pentose and glucuronate interconversions"    "Fructose and mannose metabolism"             "Galactose metabolism"                       
 [7] "Ascorbate and aldarate metabolism"           "Starch and sucrose metabolism"               "Amino sugar and nucleotide sugar metabolism"
[10] "Pyruvate metabolism"                        

5 Run scMetabolism Analysis check

# Check results
cat("\nMetabolism assay dimensions:", dim(SS@assays$METABOLISM$score), "\n")

Metabolism assay dimensions: 85 49305 

cat("Pathways quantified:", nrow(SS@assays$METABOLISM$score), "\n")

Pathways quantified: 85 

cat("\nFirst 10 pathways:\n")

First 10 pathways:

print(head(rownames(SS@assays$METABOLISM$score), 10))

 [1] "Glycolysis / Gluconeogenesis"                "Citrate cycle (TCA cycle)"                   "Pentose phosphate pathway"                  
 [4] "Pentose and glucuronate interconversions"    "Fructose and mannose metabolism"             "Galactose metabolism"                       
 [7] "Ascorbate and aldarate metabolism"           "Starch and sucrose metabolism"               "Amino sugar and nucleotide sugar metabolism"
[10] "Pyruvate metabolism"                        

6 Extract Metabolism Scores

# Extract metabolism matrix (official method from GitHub)
metabolism.matrix <- SS@assays$METABOLISM$score

cat("Pathways:", nrow(metabolism.matrix), "| Cells:", ncol(metabolism.matrix), "\n")

Pathways: 85 | Cells: 49305 

print(head(rownames(metabolism.matrix), 10))

 [1] "Glycolysis / Gluconeogenesis"                "Citrate cycle (TCA cycle)"                   "Pentose phosphate pathway"                  
 [4] "Pentose and glucuronate interconversions"    "Fructose and mannose metabolism"             "Galactose metabolism"                       
 [7] "Ascorbate and aldarate metabolism"           "Starch and sucrose metabolism"               "Amino sugar and nucleotide sugar metabolism"
[10] "Pyruvate metabolism"                        

7 Top Scoring Pathways Analysis

# Get metabolism scores
DefaultAssay(SS) <- "METABOLISM"
met_mat <- SS@assays$METABOLISM$score
pathway_means <- rowMeans(met_mat)

# Top 10 pathways
N <- 20
top_paths_global <- names(sort(pathway_means, decreasing = TRUE))[1:N]
cat("Top 20 pathways:\n")

Top 20 pathways:

print(top_paths_global)

 [1] "Oxidative phosphorylation"               "Pyruvate metabolism"                     "Glycolysis / Gluconeogenesis"           
 [4] "One carbon pool by folate"               "Cysteine and methionine metabolism"      "Citrate cycle (TCA cycle)"              
 [7] "Pentose phosphate pathway"               "Propanoate metabolism"                   "Pyrimidine metabolism"                  
[10] "Glyoxylate and dicarboxylate metabolism" "Glutathione metabolism"                  "Purine metabolism"                      
[13] "Fatty acid elongation"                   "Drug metabolism - other enzymes"         "Fructose and mannose metabolism"        
[16] "Vitamin B6 metabolism"                   "Riboflavin metabolism"                   "Biosynthesis of unsaturated fatty acids"
[19] "Phenylalanine metabolism"                "N-Glycan biosynthesis"                  

8 Visualization

# ============================================================================
# ADD ALL PATHWAYS TO METADATA (FIXED)
# ============================================================================
for (pathway in rownames(metabolism.matrix)) {
  SS@meta.data[[pathway]] <- as.numeric(metabolism.matrix[pathway, ])
}

cat("✓ Added", nrow(metabolism.matrix), "pathways to metadata\n")

✓ Added 85 pathways to metadata

# ============================================================================
# UMAP PLOTS - Top 4 Pathways
# ============================================================================
library(patchwork)
library(viridis)

p1 <- FeaturePlot(SS, features = top_paths_global[1], reduction = "umap", 
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[1]) +
  theme(plot.title = element_text(face = "bold"))

p2 <- FeaturePlot(SS, features = top_paths_global[2], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[2]) +
  theme(plot.title = element_text(face = "bold"))

p3 <- FeaturePlot(SS, features = top_paths_global[3], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[3]) +
  theme(plot.title = element_text(face = "bold"))

p4 <- FeaturePlot(SS, features = top_paths_global[4], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[4]) +
  theme(plot.title = element_text(face = "bold"))

(p1 | p2) / (p3 | p4)

# All top 10
all_plots <- list()
for (i in 1:10) {
  all_plots[[i]] <- FeaturePlot(
    SS, 
    features = top_paths_global[i], 
    reduction = "umap",
    pt.size = 0.3,
    cols = c("lightgrey", "red")
  ) + 
    ggtitle(top_paths_global[i]) +
    NoLegend() +
    theme(plot.title = element_text(size = 9, face = "bold"))
}

wrap_plots(all_plots, ncol = 3)

# Key pathways
p_oxphos <- FeaturePlot(SS, "Oxidative phosphorylation", 
                        reduction = "umap", pt.size = 0.5,
                        cols = viridis(100)) +
  ggtitle("Oxidative Phosphorylation (#1)")

p_glyc <- FeaturePlot(SS, "Glycolysis / Gluconeogenesis",
                      reduction = "umap", pt.size = 0.5,
                      cols = viridis(100)) +
  ggtitle("Glycolysis / Gluconeogenesis")

p_oxphos | p_glyc

# DotPlot
DotPlot.metabolism(
  obj = SS,
  pathway = top_paths_global,
  phenotype = "seurat_clusters",
  norm = "y"
) +
  ggtitle("Top 10 Pathways") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

# BoxPlot
BoxPlot.metabolism(
  obj = SS,
  pathway = top_paths_global[1:3],
  phenotype = "seurat_clusters",
  ncol = 1
)

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

# Heatmap
library(ComplexHeatmap)
library(circlize)

clusters <- SS$seurat_clusters
avg_by_cluster <- sapply(levels(factor(clusters)), function(cl) {
  cells <- which(clusters == cl)
  rowMeans(metabolism.matrix[, cells, drop = FALSE])
})
colnames(avg_by_cluster) <- levels(factor(clusters))

avg_top10 <- avg_by_cluster[top_paths_global, ]
avg_top10_scaled <- t(scale(t(avg_top10)))

Heatmap(
  avg_top10_scaled,
  name = "Z-score",
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  column_title = "Top 10 Pathways by Cluster"
)

NA
NA

9 UMAP PLOTS - Top 4 Pathways

# ============================================================================
# UMAP PLOTS - Top 4 Pathways
# ============================================================================
library(patchwork)
library(viridis)

p1 <- FeaturePlot(SS, features = top_paths_global[1], reduction = "umap", 
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[1]) +
  theme(plot.title = element_text(face = "bold"))

p2 <- FeaturePlot(SS, features = top_paths_global[2], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[2]) +
  theme(plot.title = element_text(face = "bold"))

p3 <- FeaturePlot(SS, features = top_paths_global[3], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[3]) +
  theme(plot.title = element_text(face = "bold"))

p4 <- FeaturePlot(SS, features = top_paths_global[4], reduction = "umap",
                  pt.size = 0.5, cols = viridis(100)) +
  ggtitle(top_paths_global[4]) +
  theme(plot.title = element_text(face = "bold"))

(p1 | p2) / (p3 | p4)

# All top 10
all_plots <- list()
for (i in 1:10) {
  all_plots[[i]] <- FeaturePlot(
    SS, 
    features = top_paths_global[i], 
    reduction = "umap",
    pt.size = 0.3,
    cols = c("lightgrey", "red")
  ) + 
    ggtitle(top_paths_global[i]) +
    NoLegend() +
    theme(plot.title = element_text(size = 9, face = "bold"))
}

wrap_plots(all_plots, ncol = 3)

10 UMAP PLOTS - Top 20 Pathways

# All top 10
all_plots <- list()
for (i in 1:20) {
  all_plots[[i]] <- FeaturePlot(
    SS, 
    features = top_paths_global[i], 
    reduction = "umap",
    pt.size = 0.3,
    cols = c("lightgrey", "red")
  ) + 
    ggtitle(top_paths_global[i]) +
    NoLegend() +
    theme(plot.title = element_text(size = 9, face = "bold"))
}

wrap_plots(all_plots, ncol = 3)

11 FeaturePlots - Top 20 Pathways

# Key pathways
p_oxphos <- FeaturePlot(SS, "Oxidative phosphorylation", 
                        reduction = "umap", pt.size = 0.5,
                        cols = viridis(100)) +
  ggtitle("Oxidative Phosphorylation (#1)")

p_glyc <- FeaturePlot(SS, "Glycolysis / Gluconeogenesis",
                      reduction = "umap", pt.size = 0.5,
                      cols = viridis(100)) +
  ggtitle("Glycolysis / Gluconeogenesis")

p_oxphos | p_glyc

12 Dot PLOTS

# DotPlot
DotPlot.metabolism(
  obj = SS,
  pathway = top_paths_global,
  phenotype = "orig.ident",
  norm = "y"
) +
  ggtitle("Top 10 Pathways") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

# DotPlot
DotPlot.metabolism(
  obj = SS,
  pathway = top_paths_global,
  phenotype = "seurat_clusters",
  norm = "y"
) +
  ggtitle("Top 10 Pathways") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

13 Box PLOTS

# BoxPlot
BoxPlot.metabolism(
  obj = SS,
  pathway = top_paths_global[1:3],
  phenotype = "orig.ident",
  ncol = 1
)

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

# BoxPlot
BoxPlot.metabolism(
  obj = SS,
  pathway = top_paths_global[1:3],
  phenotype = "seurat_clusters",
  ncol = 1
)

Please Cite: 
Yingcheng Wu, Qiang Gao, et al. Cancer Discovery. 2021. 
https://pubmed.ncbi.nlm.nih.gov/34417225/   

14 Heatmap

# Heatmap
library(ComplexHeatmap)
library(circlize)

clusters <- SS$orig.ident
avg_by_cluster <- sapply(levels(factor(clusters)), function(cl) {
  cells <- which(clusters == cl)
  rowMeans(metabolism.matrix[, cells, drop = FALSE])
})
colnames(avg_by_cluster) <- levels(factor(clusters))

avg_top10 <- avg_by_cluster[top_paths_global, ]
avg_top10_scaled <- t(scale(t(avg_top10)))

Heatmap(
  avg_top10_scaled,
  name = "Z-score",
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  column_title = "Top 10 Pathways by Cluster"
)

# Heatmap
library(ComplexHeatmap)
library(circlize)

clusters <- SS$seurat_clusters
avg_by_cluster <- sapply(levels(factor(clusters)), function(cl) {
  cells <- which(clusters == cl)
  rowMeans(metabolism.matrix[, cells, drop = FALSE])
})
colnames(avg_by_cluster) <- levels(factor(clusters))

avg_top10 <- avg_by_cluster[top_paths_global, ]
avg_top10_scaled <- t(scale(t(avg_top10)))

Heatmap(
  avg_top10_scaled,
  name = "Z-score",
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  column_title = "Top 10 Pathways by Cluster"
)

NA
NA
NA

15 Heatmap

# Heatmap
library(ComplexHeatmap)
library(circlize)

clusters <- SS$orig.ident
avg_by_cluster <- sapply(levels(factor(clusters)), function(cl) {
  cells <- which(clusters == cl)
  rowMeans(metabolism.matrix[, cells, drop = FALSE])
})
colnames(avg_by_cluster) <- levels(factor(clusters))

avg_top10 <- avg_by_cluster[top_paths_global, ]
avg_top10_scaled <- t(scale(t(avg_top10)))

Heatmap(
  avg_top10_scaled,
  name = "Z-score",
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  column_title = "Top 10 Pathways by Cluster"
)

# Heatmap
library(ComplexHeatmap)
library(circlize)

clusters <- SS$seurat_clusters
avg_by_cluster <- sapply(levels(factor(clusters)), function(cl) {
  cells <- which(clusters == cl)
  rowMeans(metabolism.matrix[, cells, drop = FALSE])
})
colnames(avg_by_cluster) <- levels(factor(clusters))

avg_top10 <- avg_by_cluster[top_paths_global, ]
avg_top10_scaled <- t(scale(t(avg_top10)))

Heatmap(
  avg_top10_scaled,
  name = "Z-score",
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
  column_title = "Top 10 Pathways by Cluster"
)

NA
NA
NA