Sézary Syndrome Single-Cell RNA-Seq Marker Analysis-15-09-2025-Part2

load libraries————————————

1. Load and Subset Normal CD4 T Cells

Layers(All_samples_Merged@assays$RNA)

[1] "data"   "counts"

2. Run FindAllMarkers on RNA Assay

# Set assay to RNA
DefaultAssay(All_samples_Merged) <- "RNA"

# Check if RNA assay has normalized data
if (!"data" %in% slotNames(All_samples_Merged[["RNA"]]) || 
    ncol(GetAssayData(All_samples_Merged[["RNA"]], slot = "data")) == 0) {
  
  message("RNA assay does not contain normalized data. Running NormalizeData...")
  All_samples_Merged <- NormalizeData(
    All_samples_Merged,
    assay = "RNA",
    normalization.method = "LogNormalize",
    scale.factor = 10000
  )
  
} else {
  message("RNA assay already contains normalized data. Skipping normalization.")
}

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Idents(All_samples_Merged) <- "seurat_clusters"

# Find markers for all clusters
all_markers <- FindAllMarkers(
    All_samples_Merged,
    logfc.threshold = 0.2,
    test.use = "wilcox",
    min.pct = 0.1,
    min.diff.pct = 0.2,
    only.pos = TRUE,
    max.cells.per.ident = 50,
    assay = "RNA"
)

# Save markers
write.csv(all_markers, "AllMarkers_clusters-17-09-2025-part2-nbi-parameters.csv", row.names = FALSE)

cat("Full marker list saved to AllMarkers_clusters-17-09-2025-part2-nbi-parameters.csv\n")

Full marker list saved to AllMarkers_clusters-17-09-2025-part2-nbi-parameters.csv

3. Top 25 and Top 5 markers per cluster

# Top 25 markers per cluster
all_markers %>%
    group_by(cluster) %>%
    top_n(-25, p_val_adj) %>% 
    # In case of tied p-values, further select the top 25 genes by fold-change
    top_n(25, avg_log2FC) -> 
  top25
head(top25)

write.csv(top25, "Top25Markers_clusters-part2-nbi.csv", row.names = FALSE)

cat("Top 25 marker list saved to Top25Markers_clusters-part2-nbi.csv\n")

Top 25 marker list saved to Top25Markers_clusters-part2-nbi.csv

# Top 5 markers per cluster
all_markers %>%
    group_by(cluster) %>%
    top_n(-5, p_val_adj) %>% 
    # In case of tied p-values, further select the top 25 genes by fold-change
    top_n(5, avg_log2FC) -> 
  top5
head(top5)

write.csv(top25, "Top5Markers_clusters-part2-nbi.csv", row.names = FALSE)

cat("Top 5 marker list saved to Top5Markers_clusters-part2-nbi.csv\n")

Top 5 marker list saved to Top5Markers_clusters-part2-nbi.csv

3. Visualize Markers

DefaultAssay(All_samples_Merged) <- "RNA"

# Unique genes from top5
all_genes <- top5$gene

# 1️⃣ Check overlap with RNA assay
valid_genes <- intersect(unique(all_genes), rownames(All_samples_Merged[["RNA"]]))
invalid_genes <- setdiff(unique(all_genes), rownames(All_samples_Merged[["RNA"]]))

cat("Valid genes:", length(valid_genes), "\n")

Valid genes: 67 

cat("Missing genes:", invalid_genes, "\n\n")

Missing genes:  

# 2️⃣ Check duplicates in top5 list
gene_counts <- table(all_genes)
duplicate_genes <- names(gene_counts[gene_counts > 1])

if(length(duplicate_genes) > 0){
  cat("Genes appearing more than once in top5 list:\n")
  print(duplicate_genes)
} else {
  cat("No duplicates found in top5 list.\n")
}

Genes appearing more than once in top5 list:
[1] "EGFL6"       "RPS4Y1"      "TRAV38-2DV8"

# Scale only the top5 marker genes
All_samples_Merged <- ScaleData(All_samples_Merged, features = unique(top5$gene))

  |                                                                                                            
  |                                                                                                      |   0%
  |                                                                                                            
  |======================================================================================================| 100%

 # DotPlot for top 5 markers
DotPlot(All_samples_Merged, features = unique(top5$gene), cols =
c("grey", "firebrick"), dot.scale = 8, assay = "RNA") + RotatedAxis()