Identify Top5 Markers

1. load libraries

2. Load Seurat Object

All_samples_Merged <- readRDS("../../../0-Seurat_RDS_OBJECT_FINAL/Seurat_object_Final_changes/All_samples_Merged_with_STCAT_Annotation_final-5-09-2025.rds")

DefaultAssay(All_samples_Merged) <- "RNA"

All_samples_Merged <- NormalizeData(
    All_samples_Merged,
    normalization.method = "LogNormalize",
    scale.factor = 1e4,
    verbose = TRUE
  )

3.find Top markers

Idents(All_samples_Merged) <- "seurat_clusters"

# ---------------------------------------------------------
# 2️⃣ Find marker genes per cluster
SS_markers <- FindAllMarkers(
  All_samples_Merged,
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25,
  min.pct.diff = 0.2
  
)

library(dplyr)

# Precise blacklist for uninformative genes
blacklist_patterns <- c(
  "^TRAV", "^TRBV", "^TRGV", "^TRDV", "^TRBC", "^TRAC", "^TRDC", "^TRGC", # TCR
  "^IGH", "^IGK", "^IGL", "^IGJ",                                         # Ig genes
  "^RPL", "^RPS",                                                         # ribosomal
  "^MT-",                                                                 # mitochondria
  "^HBA", "^HBB", "^HB[ABZ]",                                             # hemoglobins
  "^NEAT1$", "^MALAT1$",                                                  # optional lncRNAs
  "^XIST$"                              )

blacklist_regex <- paste(blacklist_patterns, collapse = "|")

# Preview which markers will be removed
to_remove <- SS_markers %>%
  filter(grepl(blacklist_regex, gene, ignore.case = TRUE))
message("Rows to remove: ", nrow(to_remove))
head(to_remove$gene)

# Filter markers (keep important metabolic/proliferation genes)
SS_markers_filtered <- SS_markers %>%
  filter(!grepl(blacklist_regex, gene, ignore.case = TRUE))

4. Top5 Markers

library(dplyr)

# ---------------------------------------------------------
# Save filtered markers
write.csv(SS_markers_filtered, file = "SS_markers_filtered.csv", row.names = FALSE)

# ---------------------------------------------------------
# Extract top 25 markers per cluster
top25_markers <- SS_markers_filtered %>%
  filter(p_val_adj < 0.05) %>%  # ensure statistical significance
  group_by(cluster) %>%
  slice_max(order_by = avg_log2FC, n = 25) %>%
  ungroup()

write.csv(top25_markers, file = "SS_markers_top25.csv", row.names = FALSE)

# ---------------------------------------------------------
# Extract top 5 markers per cluster
top5_markers <- SS_markers_filtered %>%
  filter(p_val_adj < 0.05) %>%  # ensure statistical significance
  group_by(cluster) %>%
  slice_max(order_by = avg_log2FC, n = 5) %>%
  ungroup()

write.csv(top5_markers, file = "SS_markers_top5.csv", row.names = FALSE)

message("Filtered markers, top25, and top5 markers saved successfully.")