Seurat v5 Skill

Guide for single-cell and multimodal data analysis with Seurat v5 in R. This SKILL.md contains the essential quick reference. Detailed workflows are in references/ files, read the relevant one when you need step-by-step code for a specific analysis.

Update check

Throttle to one check per 24 hours per installation; never mutate the skill directory without explicit user consent.

1. If <this-skill-dir>/.last_update exists and is less than 24 hours old, skip this step entirely.

1. Otherwise, fetch the latest tag from upstream:

   git -C <this-skill-dir> ls-remote --tags origin 'v*' 2>/dev/null \
     | awk '{print $2}' | sed 's|refs/tags/||' \
     | sort -V | tail -1

1. Compare with this skill's metadata.version from the frontmatter. If the upstream tag is strictly newer (semver), tell the user one line and ask:

"A newer version of this skill is available: vX.Y.Z → vA.B.C. Want me to git pull?"

If they say yes, run git -C <this-skill-dir> pull --ff-only. Refresh .last_update either way so the prompt doesn't repeat for 24 hours.

1. If upstream is the same or older, refresh .last_update silently and continue.

1. On any failure (offline, not a git checkout — e.g. ClawHub-installed copy, read-only path, no permission), swallow the error silently and continue with the user's task. Do not mention the failure.

Reference Files

Read the relevant reference file when the user's task matches a topic:

Quick Reference

Standard Workflow

obj = CreateSeuratObject(counts = counts, project = "my_project", min.cells = 3, min.features = 200)
obj[["percent.mt"]] = PercentageFeatureSet(obj, pattern = "^MT-")
obj = subset(obj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

# Option A: Log-normalize
obj = NormalizeData(obj)
obj = FindVariableFeatures(obj)
obj = ScaleData(obj)

# Option B: SCTransform (replaces the 3 steps above)
obj = SCTransform(obj)

obj = RunPCA(obj)
obj = FindNeighbors(obj, dims = 1:30)
obj = FindClusters(obj, resolution = 0.5)
obj = RunUMAP(obj, dims = 1:30)
DimPlot(obj, reduction = "umap", label = TRUE)

Differential Expression

markers = FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
markers = FindMarkers(obj, ident.1 = "cluster1", ident.2 = "cluster2")
markers = FindMarkers(obj, ident.1 = "cluster1", test.use = "DESeq2", slot = "counts")

Integration

# v5 layer-based integration
obj[["RNA"]] = split(obj[["RNA"]], f = obj$batch)
obj = NormalizeData(obj)
obj = FindVariableFeatures(obj)
obj = ScaleData(obj)
obj = RunPCA(obj)
obj = IntegrateLayers(obj, method = CCAIntegration, orig.reduction = "pca",
  new.reduction = "integrated.cca")
# Also: RPCAIntegration, HarmonyIntegration, FastMNNIntegration, scVIIntegration
obj = FindNeighbors(obj, reduction = "integrated.cca", dims = 1:30)
obj = FindClusters(obj, resolution = 0.5)
obj = RunUMAP(obj, reduction = "integrated.cca", dims = 1:30)
obj[["RNA"]] = JoinLayers(obj[["RNA"]])

Subsetting

subset(obj, idents = "B")                              # by cluster identity
subset(obj, idents = c("B", "NK"), invert = TRUE)      # exclude clusters
subset(obj, subset = MS4A1 > 2.5)                      # by expression
subset(obj, subset = condition == "treated")            # by metadata
subset(obj, downsample = 100)                           # downsample per cluster

Key Visualization

DimPlot(obj, reduction = "umap", group.by = "celltype", label = TRUE)
FeaturePlot(obj, features = c("CD3D", "MS4A1", "CD8A"))
VlnPlot(obj, features = c("CD3D", "MS4A1"), group.by = "celltype")
DotPlot(obj, features = c("CD3D", "MS4A1", "CD14"), group.by = "celltype")
DoHeatmap(obj, features = top_markers) + NoLegend()
FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

Object Access

Cells(obj)                    # cell barcodes
Features(obj)                 # gene names
Idents(obj)                   # active identities
obj[[]]                       # metadata data.frame
obj$nCount_RNA                # single metadata column
Embeddings(obj, "pca")        # PCA embeddings
obj[["RNA"]]$counts           # raw counts (v5 layer)
DefaultAssay(obj)             # current default assay
Layers(obj)                   # list layers
VariableFeatures(obj)         # HVGs
FetchData(obj, vars = c("UMAP_1", "UMAP_2", "CD3D"))  # mixed data access

Pseudobulk

bulk = AggregateExpression(obj, group.by = c("celltype", "sample"), return.seurat = TRUE)

Multi-Assay (CITE-seq)

obj[["ADT"]] = CreateAssayObject(counts = adt.counts)
obj = NormalizeData(obj, assay = "ADT", normalization.method = "CLR", margin = 2)
DefaultAssay(obj) = "ADT"
FeaturePlot(obj, features = "adt_CD3")

Spatial

obj = Load10X_Spatial(data.dir = "path/to/spaceranger/outs")
SpatialDimPlot(obj)
SpatialFeaturePlot(obj, features = "MS4A1")

Format Conversion

# Seurat to AnnData
library(SeuratDisk)
SaveH5Seurat(obj, filename = "obj.h5Seurat")
Convert("obj.h5Seurat", dest = "h5ad")

# Seurat to SingleCellExperiment
sce = as.SingleCellExperiment(obj)

# SingleCellExperiment to Seurat
obj = as.Seurat(sce)