Seurat (R) — single-cell RNA-seq

Overview

Seurat (Satija lab) is the dominant R toolkit for single-cell/spatial analysis. Use it when the analysis must live in R (Bioconductor downstream, R-based collaborators, lab convention) — otherwise [[scanpy]]/[[anndata]] in Python is the equivalent and is already well covered in this vault. Seurat v5 introduced the layered Assay5 object and a BPCells-backed on-disk path for very large datasets.

When to use vs scanpy

• Seurat — R/tidyverse/Bioconductor downstream, integration with DESeq2/edgeR pseudobulk, SingleR annotation, established lab pipelines.
• scanpy — Python ecosystem, scVI/deep models ([[scvi-tools]]), very large atlases, AnnData-native tooling.

Convert between them rather than re-running: use sceasy, zellkonverter (SingleCellExperiment ↔ AnnData), or SeuratDisk (.h5Seurat ↔ .h5ad).

Standard workflow

library(Seurat)

# 1. Load (10x output) and create object with basic gene/cell filters
counts <- Read10X(data.dir = "filtered_feature_bc_matrix/")
obj <- CreateSeuratObject(counts, min.cells = 3, min.features = 200)

# 2. QC: mitochondrial fraction, then filter
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")  # "^mt-" for mouse
VlnPlot(obj, c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
obj <- subset(obj, nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 15)

# 3a. Standard normalization path
obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj, nfeatures = 2000)
obj <- ScaleData(obj)

# 3b. OR SCTransform (preferred; replaces NormalizeData+FindVariableFeatures+ScaleData)
# obj <- SCTransform(obj, vars.to.regress = "percent.mt")

# 4. Linear + non-linear dim reduction, graph clustering
obj <- RunPCA(obj)
ElbowPlot(obj)                       # choose #PCs
obj <- FindNeighbors(obj, dims = 1:30)
obj <- FindClusters(obj, resolution = 0.5)   # higher res -> more clusters
obj <- RunUMAP(obj, dims = 1:30)
DimPlot(obj, reduction = "umap", label = TRUE)

# 5. Marker genes
markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
FeaturePlot(obj, c("CD3D", "MS4A1", "LYZ"))

saveRDS(obj, "obj.rds")

Integration (multiple samples / batch correction)

Seurat v5 unifies integration under IntegrateLayers. Split layers by batch first:

obj[["RNA"]] <- split(obj[["RNA"]], f = obj$sample)
obj <- NormalizeData(obj) |> FindVariableFeatures() |> ScaleData() |> RunPCA()
obj <- IntegrateLayers(obj, method = CCAIntegration,        # or HarmonyIntegration,
                       orig.reduction = "pca",              #    RPCAIntegration,
                       new.reduction = "integrated.cca")    #    scVIIntegration
obj <- FindNeighbors(obj, reduction = "integrated.cca", dims = 1:30)
obj <- FindClusters(obj); obj <- RunUMAP(obj, reduction = "integrated.cca", dims = 1:30)
obj <- JoinLayers(obj)   # re-join before DE / pseudobulk

HarmonyIntegration mirrors [[harmonypy]]; RPCA is faster/more conservative than CCA for large or weakly-overlapping datasets.

Gotchas

• v5 layers: counts/data/scale.data are now layers. After integration, call JoinLayers() before differential expression or pseudobulk, or you'll operate per-batch.
• SCTransform + markers: run PrepSCTFindMarkers() before FindMarkers on SCT data.
• Species mito pattern: ^MT- (human) vs ^mt- (mouse); a wrong pattern silently yields percent.mt == 0 and bad QC.
• Pseudobulk DE beats single-cell Wilcoxon for condition contrasts: aggregate with AggregateExpression()/Seurat::PseudobulkExpression() then hand off to [[pydeseq2]] or edgeR. See [[de-summary]] for interpreting results.
• Memory: for >500k cells use the BPCells on-disk backend or switch to [[scanpy]].

Related

R counterpart to [[scanpy]]/[[anndata]]; integrate with [[harmonypy]], [[scvi-tools]] (via scVIIntegration), [[bioconductor-bridge]], and pseudobulk DE via [[pydeseq2]].