Hi all, I am working on a multi-sample scRNA-seq dataset using Seurat. I merged several Seurat objects (each sample processed independently through QC and DoubletFinder) using the merge() function. After merging, I noticed that the merged object contains multiple assay layers (e.g., multiple SCT models or multiple RNA layers depending on upstream processing).

I found that Seurat provides a function called JoinLayers(), which can combine layers across samples into a single unified layer. However, I am not fully sure when it is appropriate or necessary to use JoinLayers() in a multi-sample workflow.

Before diving into my specific questions, here is a brief outline of my current workflow for processing multiple samples:

obj.list <- lapply(files_list, function(s){
  mat <- Read10X(s)
  sample_name <- sub('_.*', '', s)
  obj <- CreateSeuratObject(counts = mat, project = sample_name)
  obj$sample <- sample_name
  return(obj)
})

QC Process

for (i in seq_along(obj.list)) {
  obj.list[[i]][["percent.mt"]] = PercentageFeatureSet(obj.list[[i]], 
                                                       pattern = "^MT-")
  print(VlnPlot(obj.list[[i]], 
          features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
          ncol = 3) )
}
#### QC for nFeature, nCount, mt #### 
obj.qc.list <- list()
for (i in seq_along(obj.list)) {
  obj.qc.list[[i]] <- subset(obj.list[[i]], 
                             subset =
                               nFeature_RNA > 300 &
                               nFeature_RNA < 6000 &
                               nCount_RNA > 500 &
                               percent.mt < 10)
  print(VlnPlot(obj.qc.list[[i]], 
                features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
                ncol = 3))
}

Merge & Normalization

merged_obj <- merge(x = obj.qc.list[[1]],
                     y = obj.qc.list[-1])
object <- NormalizeData(merged_obj,
                        normalization.method = "LogNormalize",
                        scale.factor = 10000)
object <- FindVariableFeatures(object, 
                               selection.method = "vst", 
                               nfeatures = 2000)
all.genes <- rownames(object)
object <- ScaleData(object, features = all.genes)

Run PCA & UMAP

object <- RunPCA(object = object)
object <- RunUMAP(object, 
                  dims = 1:20,
                  verbose = FALSE)

Harmony

object <- RunHarmony(object, group.by.vars = "orig.ident")

Clustering cells

object <- FindNeighbors(object, dims = 1:20, k.param = 20)
object <- FindClusters(object, algorithm = 1, resolution = 0.8)

Get expression matrix

expr <- GetAssayData(object)  ## This code went wrong, which shows below

Specifically, I would like to understand:

Under what circumstances should JoinLayers() be used after merging Seurat objects?
For example:

• Only when SCT was run separately per sample?

• Only when assays contain multiple layers (e.g., multiple SCT models)?

• Only before integration workflows (CCA/RPCA/Harmony/SCTransform integration)?

What happens if I do not run JoinLayers()?
In my case, skipping JoinLayers() leads to an error when I try to extract the expression matrix:

expr <- GetAssayData(object)

Error in GetAssayData() at SeuratObject/R/seurat.R:1943:3: ! GetAssayData doesn't work for multiple layers in v5 assay. Run rlang::last_trace() to see where the error occurred. Called from: signal_abort(cnd, .file)

So it seems that downstream steps relying on a unified assay layer (e.g., extracting counts/data, etc.) may fail if multiple layers remain.

Is there a recommended best practice for multi-sample workflows regarding JoinLayers()?
For example:

• Should I always run SCT per sample > merge > JoinLayers() > integration?

• Or merge > NormalizeData > VariableFeaturePlot > ScaleData > Harmony > JoinLayers()

• Or is it better to merge raw RNA counts first and then run SCT only once on the merged object?

Any clarification or examples would be greatly appreciated. Thanks!