Unable to run DEA analysis on objects that been through SCTransform with percent.mt and cell cycle score as vars.to.regress

[theresiaakh]

Hi. I am using Seurat 5.3.1 for analysis of my scRNAseq data on pbmcs. I am unable to run differential expression on my object when I put it through SCTransform followed by integration by IntegrateLayers.

This is my workflow:
Filtering and SCTransform:
for (i in seq_along(seurat_list)) {
sample_name <- names(seurat_list)[i]
#Filter
sample_filtered <- subset(
seurat_list[[i]],
subset = nFeature_RNA > 250 &
nFeature_RNA < 7000 &
nCount_RNA > 500 &
nCount_RNA < 25000 &
percent.mt < 15
)

sample_filtered <- NormalizeData(sample_filtered, verbose = FALSE)
sample_filtered <- FindVariableFeatures(sample_filtered, verbose = FALSE)
sample_filtered <- ScaleData(sample_filtered, verbose = FALSE)
#Cell cycle scoring
sample_filtered <- CellCycleScoring(
sample_filtered,
s.features = cc.genes$s.genes,
g2m.features = cc.genes$g2m.genes,
set.ident = TRUE
)

#SCTransform with both mitochondrial AND cell cycle regression
sample_sct <- SCTransform(
sample_filtered,
vars.to.regress = c("percent.mt", "nCount_RNA", "S.Score", "G2M.Score"),
verbose = TRUE
)

seurat_list_sct_cellcycle[[sample_name]] <- sample_sct
}

Integration:
pbmc_merged <- merge(
x = seurat_list_sct_cellcycle[[1]],
y = seurat_list_sct_cellcycle[-1],
add.cell.ids = names(seurat_list_sct_cellcycle),
project = "CART_MM_PBMC"
)

pbmc_merged <- FindVariableFeatures(pbmc_merged, nfeatures = 2000, verbose = FALSE)
pbmc_merged <- RunPCA(pbmc_merged)

pbmc_merged <- FindNeighbors(pbmc_merged, dims = 1:30, reduction = "pca")
pbmc_merged <- FindClusters(pbmc_merged, resolution = 2, cluster.name = "unintegrated_clusters")
pbmc_merged <- RunUMAP(pbmc_merged, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

pbmc_integrated <- IntegrateLayers(
object = pbmc_merged,
method = HarmonyIntegration,
normalization.method = "SCT",
verbose = FALSE
)
pbmc_integrated <- FindNeighbors(
pbmc_integrated,
reduction = "harmony", #
dims = 1:30
)

pbmc_integrated <- FindClusters(
pbmc_integrated,
resolution = 0.5,
cluster.name = "integrated_clusters"
)

pbmc_integrated <- RunUMAP(
pbmc_integrated,
reduction = "harmony", # Use harmony reduction
dims = 1:30,
reduction.name = "umap.harmony"
)

When I later try to do PrepSCTFindMarkers, I get this message:
Minimum UMI unchanged. Skipping re-correction.

When I then proceed with FIndMarkers I get this message:
Error in FindMarkers.SCTAssay(object = data.use, latent.vars = latent.vars, :
Object contains multiple models with unequal library sizes. Run PrepSCTFindMarkers() before running FindMarkers().

I have to change the default assay to RNA instead of SCT to be able to run the differential expression.

I also tried doing the integration accordin gto the vignette with RPCA instead of harmony, but I run in to the same issue.

Thank you so much!