19 Marker Gene Identification

Marker genes are identified for every cluster at every clustering resolution computed on the weighted nearest neighbor (WNN) graph in the Cross-Modal Integration chapter (wsnn_res.0.1 through wsnn_res.2, in steps of 0.1). Testing across the full resolution sweep — rather than a single chosen resolution — allows the marker sets to be compared side by side and supports the choice of a final resolution based on how cleanly clusters separate biologically, not just on clustree topology alone.

For each resolution, a Wilcoxon rank-sum test is run per cluster (cluster vs. all other cells) using presto::wilcoxauc, which is a fast implementation of the same test used by Seurat::FindAllMarkers(test.use = "wilcox"). The AUC statistic reflects how well a gene’s expression discriminates a cluster from the rest of the dataset: an AUC near 1 (or near 0) indicates a gene that is strongly up- (or down-) regulated in that cluster relative to the rest, while an AUC near 0.5 indicates no discriminative power.

19.1 Run marker genes detection across all resolutions

Every wsnn_res.* column in the metadata is treated as a separate clustering to test. Resolutions that collapse to a single cluster are skipped, since a one-cluster grouping has no marker genes to detect against. For each remaining resolution, Idents is set to that resolution’s cluster labels, normalised RNA expression data is pulled from the data layer, and presto::wilcoxauc is run per gene per cluster. Results from every resolution are combined into a single long-format table, tagged with the resolution they came from in the annotation column.

library(presto)
library(dplyr)
library(tibble)

DefaultAssay(combined) <- "RNA"
combined <- JoinLayers(combined, assay = "RNA")

resolutions <- colnames(combined@meta.data)[grep("wsnn_res", colnames(combined@meta.data))]

genes_list <- data.frame()

for (res in resolutions) {
  cluster_count <- length(levels(factor(combined@meta.data[[res]])))
  if (cluster_count == 1) next
  
  Idents(combined) <- combined@meta.data[[res]]
  counts_combined <- GetAssayData(combined, layer = "data")
  cluster_labels <- combined@meta.data[[res]]
  
  markers <- presto::wilcoxauc(counts_combined, cluster_labels) |>
    as.data.frame() |>
    rownames_to_column("gene") |>
    mutate(annotation = res)
  
  genes_list <- rbind(genes_list, markers)
}

19.2 Run marker peaks detection across all resolutions

library(presto)
library(dplyr)
library(tibble)

DefaultAssay(combined) <- "ATAC"


resolutions <- colnames(combined@meta.data)[grep("wsnn_res", colnames(combined@meta.data))]

peaks_list <- data.frame()

for (res in resolutions) {
  cluster_count <- length(levels(factor(combined@meta.data[[res]])))
  if (cluster_count == 1) next
  
  Idents(combined) <- combined@meta.data[[res]]
  counts_combined <- GetAssayData(combined, layer = "data")
  cluster_labels <- combined@meta.data[[res]]
  
  markers <- presto::wilcoxauc(counts_combined, cluster_labels) |>
    as.data.frame() |>
    rownames_to_column("peak") |>
    mutate(annotation = res)
  
  peaks_list <- rbind(peaks_list, markers)
}

19.3 Save marker genes and peaks in parquet format/table

The combined marker table across all resolutions is written to parquet for downstream use (e.g. resolution selection, cluster annotation, or loading into the ShinyCellModular explorer).

library(arrow)
write_parquet(genes_list, "data/markers_genes_list_wsnn_all_resolutions.parquet")
write_parquet(peaks_list, "data/markers_peaks_list_wsnn_all_resolutions.parquet")