I am analyzing Visium HD spatial transcriptomics data with cell segmentation. The experiment was designed using two Visium HD slides, each slide containing two capture areas. Each capture area contains a tissue microarray (TMA) with 9 patient cores, so in total there are 4 TMAs (TMA1–TMA4) representing 36 patient samples.

Each TMA was processed separately by the 10x pipeline, and the filtered matrices were loaded into Seurat. After merging the objects and performing clustering, I attempted to run FindAllMarkers() to obtain cluster markers, but the output is empty.

Below is the relevant part of the workflow starting from merging the Seurat objects.

# Merge Seurat objects from four TMAs
# obj1, obj2, obj3, obj4 correspond to TMA1–TMA4
obj <- merge(obj1, y = list(obj2, obj3, obj4),
             add.cell.ids = c("TMA1","TMA2","TMA3","TMA4"))

After merging:

An object of class Seurat
18093 features across ~173000 cells
Active assay: RNA
4 layers present: counts.1 counts.2 counts.3 counts.4

Each layer corresponds to one TMA.

### QC filtering
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")

obj <- subset(
  obj,
  subset = nFeature_RNA > 200 &
           nFeature_RNA < 6000 &
           percent.mt < 20
)

Normalization and dimensionality reduction

objlog <- NormalizeData(obj)
objlog <- FindVariableFeatures(objlog)
objlog <- ScaleData(objlog)
objlog <- RunPCA(objlog)

Check batch structure and run UMAP

objlog <- RunUMAP(objlog, dims = 1:30)
DimPlot(objlog, group.by = c("batch","sample_id"))

Batch correction with Harmony

objlog <- RunHarmony(objlog, group.by.vars = "batch")
objlog <- RunUMAP(objlog, reduction = "harmony", dims = 1:30)
objlog <- FindNeighbors(objlog, reduction = "harmony", dims = 1:30)
objlog <- FindClusters(objlog)

After clustering:

Seurat object
~18000 genes
~172000 cells
layers:
counts.1 counts.2 counts.3 counts.4
data.1 data.2 data.3 data.4
scale.data
reductions:
pca
umap
harmony

Attempt to identify cluster markers

Idents(objlog) <- "seurat_clusters"
markers <- FindAllMarkers(
  objlog,
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25)

The function runs through all clusters:

Calculating cluster 0
Calculating cluster 1
...
Calculating cluster 28

However the resulting markers object is empty.

I also tried using Presto:

library(presto)

markers <- FindAllMarkers(objlog, test.use = "wilcox", only.pos = TRUE, min.pct = 0.25)

but the result is still empty.

Attempt to join layers

Because the assay contains multiple layers, I attempted:

objlog[["RNA"]] <- JoinLayers(objlog)

or

objlog[["RNA"]] <- JoinLayers(objlog, layers = "data")

However this step does not finish even after >20 hours on this dataset (~172k cells).

Questions

1. Why does FindAllMarkers() return no markers in this layered Seurat v5 object?
2. Is it necessary to run JoinLayers() before marker detection in Seurat v5?
3. If so, what is the recommended approach for performing marker detection on large layered objects (~170k cells) without extremely long runtimes?

Any suggestions or best practices for handling marker detection in large Seurat v5 layered datasets would be greatly appreciated.