Single-Cell RNA-Seq Analysis with Seurat v5 in R: Quality Control, Normalization, Clustering, and Marker Identification

Single-cell RNA sequencing (scRNA-seq) has transformed molecular biology from bulk tissue averages into single-cell resolution atlases. By capturing the transcriptomes of tens of thousands of individual cells simultaneously, scRNA-seq reveals previously unobservable cellular heterogeneity, identifies rare stem or circulating progenitor populations, resolves dynamic differentiation trajectories, and maps the microenvironments of complex diseases such as solid tumors and neurodegenerative disorders.

The R package Seurat (version 5), maintained by the Satija Lab at the New York Genome Center, is the most widely adopted analytical framework for single-cell genomics. Seurat v5 introduces a modular Assay5 structure with split layers, dramatically accelerating analysis on large-scale datasets while supporting seamless integration of multi-modal, spatial, and cross-batch experiments.

In this publication-grade tutorial, we walk step-by-step through a complete end-to-end scRNA-seq pipeline in R—from raw 10x Genomics matrices to publication-ready UMAP visualizations and cell-type annotations.

If your research extends into target structure validation and drug discovery, explore our companion guides on R Environment Setup, foundational Single-Cell RNA-Seq Theory, and our specialized Molecular Dynamics Simulation Services.

1. Introduction & Real-World Biological Context

In standard bulk RNA-seq, transcripts from millions of heterogeneous cells are pooled and sequenced together. The resulting expression values represent an ensemble average:

Y_bar_g = (1 / N) * sum(Y_gi)

While effective for broad comparative condition testing (e.g., healthy liver versus treated liver), bulk sequencing conceals vital biological phenomena:

• Masked Subpopulations: A 5-fold upregulation in a rare 2% immune subset (such as FoxP3+ regulatory T cells or tumor-initiating stem cells) appears as negligible background noise in bulk sequencing.
• Transcriptional Gradients: Continuous differentiation trajectories (e.g., hematopoiesis) are flattened into artificial discrete averages.
• Cell-to-Cell Stochasticity: Bimodal gene expression (where 50% of cells express a transcript at high levels and 50% are off) cannot be distinguished from uniform moderate expression across all cells.

scRNA-seq solves this challenge by tagging individual cellular transcripts with cell-specific barcodes and Unique Molecular Identifiers (UMIs), producing an exact count matrix spanning G distinct genomic features (genes) across C individual cell barcodes.

2. Prerequisites & Environment Setup

Seurat v5 relies on modern R (version ≥ 4.3.0) and high-performance Bioconductor linear modeling packages such as glmGamPoi. Run the following installation commands in your R console or RStudio environment.

# 1. Enable BiocManager and CRAN repositories
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager", repos = "https://cloud.r-project.org")

# 2. Install core CRAN dependencies
install.packages(c(
  "Seurat",
  "SeuratObject",
  "tidyverse",
  "patchwork",
  "Matrix",
  "scales",
  "harmony"
), repos = "https://cloud.r-project.org")

# 3. Install Bioconductor statistical packages
BiocManager::install(c(
  "glmGamPoi",      # Fast negative binomial GLM fitting for SCTransform v2
  "SingleR",        # Automated reference-based cell type annotation
  "celldex",        # Curated reference transcriptomic databases
  "scDblFinder",    # Fast computational doublet detection
  "ComplexHeatmap"  # High-density publication heatmaps
), update = FALSE, ask = FALSE)

# 4. Verify version integrity
packageVersion("Seurat") # Must output >= 5.0.0

3. Input Data Format & Preprocessing

Standard droplet-based single-cell platforms (e.g., 10x Genomics Chromium) process raw sequencing reads (FASTQ) through aligners like Cell Ranger to output three core files representing the sparse expression matrix:

1. barcodes.tsv.gz: The list of all identified cell cellular barcodes (e.g., AAACCCAAGCGTATGG-1).
2. features.tsv.gz: The list of Ensembl gene IDs, official gene symbols, and feature types (Gene Expression, Antibody Capture).
3. matrix.mtx.gz: Coordinate-format sparse count matrix specifying (Gene Index, Cell Index, Raw UMI Count).

Inspecting Raw File Anatomy

# Verify output directory structure
ls -lh filtered_feature_bc_matrix/
# Output:
# barcodes.tsv.gz  (150 KB)
# features.tsv.gz  (280 KB)
# matrix.mtx.gz    (18 MB)

Because 95–98% of values in single-cell count matrices are zeros (due to low per-cell mRNA capture efficiency and biological transcriptional bursting), Seurat reads these matrices as compressed sparse column matrices (dgCMatrix), using less than 5% of the memory of a dense floating-point matrix.

4. The Step-by-Step Computational Workflow

[10x Raw Counts] ---> [CreateSeuratObject] ---> [Mitochondrial QC] ---> [scDblFinder (Doublets)]
                                                                                  |
[Cell Annotation] <--- [FindAllMarkers] <--- [UMAP / Clustering] <--- [SCTransform v2]

Step 1: Initialize the Seurat v5 Object

library(Seurat)
library(tidyverse)
library(patchwork)
library(SingleR)
library(celldex)
library(scDblFinder)

# Set global seed for reproducible stochastic operations
set.seed(42)

# Load 10x Genomics filtered matrix directory
raw_counts <- Read10X(data.dir = "filtered_feature_bc_matrix/")

# Initialize Seurat v5 object
# min.cells: Include genes detected in at least 3 cells
# min.features: Include cells with at least 200 detected genes
pbmc <- CreateSeuratObject(
  counts = raw_counts,
  project = "PBMC_Immune_Atlas",
  min.cells = 3,
  min.features = 200
)

# Inspect Seurat v5 Assay5 layer structure
print(pbmc)
# Notice: Assay 'RNA' contains count layer: counts

Step 2: Quality Control (QC) Metrics & Filtering

Damaged, dying, or lysed cells exhibit leaky cell membranes where cytoplasmic mRNA escapes, leaving behind disproportionately high ratios of mitochondrial transcripts. Conversely, empty droplets contain ambient RNA, while homotypic/heterotypic doublets exhibit abnormally elevated total UMI counts.

We calculate:

• Mitochondrial Percentage (percent.mt): Ratio of reads mapping to genes prefixed with MT- (human) or mt- (mouse).
• Ribosomal Percentage (percent.ribo): Ratio of reads mapping to ribosomal proteins (RPS / RPL).

# Calculate percentage of mitochondrial and ribosomal genes
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc[["percent.ribo"]] <- PercentageFeatureSet(pbmc, pattern = "^RP[SL]")

# Visualize QC distributions via Violin Plots
qc_vln <- VlnPlot(
  pbmc,
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.ribo"),
  ncol = 4,
  pt.size = 0.1
) & theme(axis.title.x = element_blank())

ggsave("QC_PreFilter_Violin.png", qc_vln, width = 12, height = 4, dpi = 300)

# Filter low-quality barcodes and apoptotic cells
# Adjust thresholds according to experimental tissue type:
pbmc_filtered <- subset(
  pbmc,
  subset = nFeature_RNA > 500 &
           nFeature_RNA < 4500 &
           nCount_RNA > 1000 &
           nCount_RNA < 25000 &
           percent.mt < 10
)

cat("Remaining cells after quality filtering:", ncol(pbmc_filtered), "\n")

Step 3: In Silico Doublet Removal with scDblFinder

Droplet encapsulation occasionally captures two cells in a single microfluidic droplet. scDblFinder generates artificial doublets in silico by pooling real cell profiles, modeling non-linear doublet manifolds:

library(SingleCellExperiment)

# Convert Seurat object to SingleCellExperiment
sce <- as.SingleCellExperiment(pbmc_filtered)

# Execute doublet simulation and classification
sce <- scDblFinder(sce)

# Transfer doublet classification back to Seurat metadata
pbmc_filtered$scDblFinder_class <- sce$scDblFinder.class
pbmc_filtered$scDblFinder_score <- sce$scDblFinder.score

# Retain only confirmed singlet droplets
pbmc_clean <- subset(pbmc_filtered, subset = scDblFinder_class == "singlet")

cat("Confirmed singlets:", ncol(pbmc_clean), "\n")

Step 4: Normalization & Variance Stabilization via SCTransform v2

Traditional scRNA-seq workflows apply global scaling (LogNormalize), dividing each gene count by total cell counts, multiplying by 10,000, and taking the natural log. However, LogNormalize fails to eliminate technical sequencing depth confounding for high-variance genes.

Seurat v5 defaults to SCTransform v2, which uses regularized negative binomial regression with generalized linear modeling (glmGamPoi), stabilizing variance across sequencing depths:

# Execute SCTransform v2 with regression against mitochondrial variance
pbmc_clean <- SCTransform(
  pbmc_clean,
  vst.flavor = "v2",
  vars.to.regress = c("percent.mt"),
  verbose = FALSE
)

# Active assay automatically switches to 'SCT'
DefaultAssay(pbmc_clean) <- "SCT"

Step 5: Linear Dimensionality Reduction (PCA)

We project thousands of variable genes into orthogonal principal components (PCs) that capture the major biological axes of transcriptional variance:

# Run Principal Component Analysis on top 3,000 SCT variable features
pbmc_clean <- RunPCA(pbmc_clean, npcs = 50, verbose = FALSE)

# Generate Elbow Plot to objectively determine the dimensional cutoff
elbow_p <- ElbowPlot(pbmc_clean, ndims = 50) +
  geom_vline(xintercept = 30, linetype = "dashed", color = "red") +
  labs(title = "Elbow Plot: Selection of Significant PCs")

ggsave("PCA_ElbowPlot.png", elbow_p, width = 6, height = 4, dpi = 300)

Examine the plot where the standard deviation of eigenvalues reaches an inflection point (plateaus). For standard PBMC or tissue datasets, 30 dimensions captures genuine biological signal while excluding high-order noise.

Step 6: Graph-Based Clustering (Louvain / Leiden)

Seurat clusters cells using a shared nearest neighbor (SNN) graph. It computes cell-cell Euclidean distances in 30-dimensional PCA space, forms a K-nearest neighbor graph (KNN), refines edges via Jaccard similarity, and optimizes modularity using the Louvain or Leiden algorithm:

# 1. Build SNN Graph using first 30 PCs
pbmc_clean <- FindNeighbors(pbmc_clean, dims = 1:30, reduction = "pca")

# 2. Identify Clusters across multiple resolution parameters
# Resolution 0.4–1.2 typically matches biological tissue hierarchy
pbmc_clean <- FindClusters(pbmc_clean, resolution = c(0.4, 0.6, 0.8), algorithm = 1)

# Set resolution 0.6 as primary clustering level
Idents(pbmc_clean) <- "SCT_snn_res.0.6"
cat("Identified", length(unique(Idents(pbmc_clean))), "distinct clusters.\n")

Step 7: Non-Linear Dimensionality Reduction (UMAP)

Uniform Manifold Approximation and Projection (UMAP) compresses the 30-dimensional neighborhood graph into two visual coordinates while preserving local and global topological structure:

# Compute 2D UMAP projection
pbmc_clean <- RunUMAP(pbmc_clean, dims = 1:30, reduction = "pca")

# Visualize clusters on UMAP
umap_plot <- DimPlot(
  pbmc_clean,
  reduction = "umap",
  label = TRUE,
  label.size = 5,
  repel = TRUE,
  pt.size = 0.5
) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "right") +
  labs(title = "PBMC scRNA-seq: Unsupervised Louvain Clusters")

ggsave("UMAP_Clusters.png", umap_plot, width = 8, height = 6, dpi = 300)

Step 8: Differential Marker Gene Identification

To biologically identify what cell type each cluster represents, we run differential expression comparing each cluster against all other remaining cells using the non-parametric Wilcoxon rank-sum test:

# Before running differential expression in Seurat v5, prep the SCT assay
pbmc_clean <- PrepSCTFindMarkers(pbmc_clean)

# Identify positive marker genes for every cluster
all_markers <- FindAllMarkers(
  pbmc_clean,
  assay = "SCT",
  only.pos = TRUE,
  min.pct = 0.25,             # Gene must be detected in >= 25% of cluster cells
  logfc.threshold = 0.50,     # Minimum 1.41-fold change (log2(1.41) ~ 0.50)
  test.use = "wilcox"
)

# Filter for statistically significant markers
top_markers <- all_markers %>%
  filter(p_val_adj < 0.01) %>%
  group_by(cluster) %>%
  slice_max(order_by = avg_log2FC, n = 5)

write.csv(top_markers, "Top5_Cluster_Markers.csv", row.names = FALSE)

Step 9: Reference-Based Automated Annotation via SingleR

While manual marker inspection remains the gold standard, SingleR provides automated, unbiased cell typing by comparing single-cell expression profiles against curated reference atlases (e.g., Human Primary Cell Atlas):

# Load curated Human Primary Cell Atlas reference
ref_hpca <- celldex::HumanPrimaryCellAtlasData()

# Extract normalized log-counts from Seurat object
query_counts <- GetAssayData(pbmc_clean, assay = "SCT", layer = "data")

# Run SingleR classification
singler_results <- SingleR(
  test = query_counts,
  ref = ref_hpca,
  labels = ref_hpca$label.main,
  clusters = Idents(pbmc_clean) # Annotate per cluster for robust inference
)

# Inspect cluster annotation mappings
annotation_map <- data.frame(
  Cluster = rownames(singler_results),
  CellType = singler_results$pruned.labels
)
print(annotation_map)

# Assign biological names back to Seurat clusters
new_cluster_ids <- annotation_map$CellType
names(new_cluster_ids) <- levels(pbmc_clean)
pbmc_annotated <- RenameIdents(pbmc_clean, new_cluster_ids)
pbmc_annotated$CellType <- Idents(pbmc_annotated)

5. Data Visualization & Result Interpretation

Publication-ready figures in single-cell biology require clear multi-panel compositions showing marker specificity across clusters.

5.1 Canonical Marker Profiling (DotPlot & FeaturePlot)

# Canonical immune cell lineage markers:
# CD3D/CD3E: T cells | CD4: Helper T | CD8A: Cytotoxic T | MS4A1: B cells
# CD14/FCGR3A: Monocytes | NKG7: NK cells | PPBP: Platelets
canonical_genes <- c("CD3D", "CD3E", "CD4", "CD8A", "MS4A1", "CD14", "FCGR3A", "NKG7", "PPBP")

# 1. Expression Dot Plot
dot_fig <- DotPlot(
  pbmc_annotated,
  features = canonical_genes,
  cols = c("lightgrey", "#dc2626")
) +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Lineage Specificity of Canonical Immune Markers",
       x = "Gene Symbol",
       y = "Annotated Cell Identity")

# 2. Gradient FeaturePlot on UMAP
feature_fig <- FeaturePlot(
  pbmc_annotated,
  features = c("CD3E", "MS4A1", "CD14", "NKG7"),
  cols = c("#f1f5f9", "#3b82f6", "#1d4ed8"),
  ncol = 2
)

# 3. Final Annotated UMAP
annot_umap <- DimPlot(
  pbmc_annotated,
  reduction = "umap",
  group.by = "CellType",
  label = TRUE,
  repel = TRUE,
  pt.size = 0.6
) +
  scale_color_brewer(palette = "Set2") +
  theme_minimal(base_size = 13) +
  labs(title = "Annotated Human PBMC Cellular Atlas")

# Assemble final multi-panel figure
final_panel <- (annot_umap | dot_fig) / feature_fig +
  plot_annotation(tag_levels = 'A')

ggsave("Figure_scRNA_Publication_Atlas.pdf", final_panel, width = 14, height = 10, dpi = 300)

5.2 How to Interpret the Output Figures

• Dot Size vs. Color Intensity: In the DotPlot, dot diameter represents the percentage of cells in that cluster expressing the gene (pct.exp), while color intensity indicates average expression levels. A high-fidelity marker exhibits both high diameter (>70%) and dark red intensity exclusively within its target cluster.
• UMAP Separation: Clear spatial isolation between lineages (e.g., B cells vs. Monocytes) demonstrates strong orthogonal gene expression. In contrast, continuous gradients (e.g., Naive CD4+ to Memory CD4+ T cells) reflect transitional transcriptional programs.
• Platelet or Erythrocyte Contamination: Clusters showing high PPBP (platelet) or HBB (hemoglobin) expression with low gene diversity represent ambient sequencing artifacts that can be safely gated out.

6. Common Errors & Troubleshooting

7. Frequently Asked Questions (FAQ)

How many cells should be sequenced per biological sample?

A target recovery of 5,000 to 10,000 cells per sample provides sufficient statistical power to detect sub-populations comprising 1–2% of the total tissue volume.

When should batch integration algorithms like Harmony be applied?

When samples are multiplexed across different sequencing runs, donors, or collection dates, running IntegrateLayers(object, method = HarmonyIntegration) corrects technical batch offsets while conserving biological phenotypes.

Can Seurat v5 analyze spatial transcriptomics data?

Yes, Seurat v5 natively ingests 10x Genomics Visium, Visium HD, and Xenium spatial coordinates alongside standard single-cell droplet matrices.

Why does SCTransform v2 recommend against regressing out cell cycle scores?

Regressing cell cycle genes indiscriminately can erase genuine biological phenotypes in proliferating tissues, such as activated germinal center B cells or expanding cancer clones.

Continue Your Computational Biology Journey

Deepen your bioinformatics expertise with our curated learning materials:

• Learn the foundations in our Introduction to Statistics in R.
• Compare bulk and single-cell approaches with our guide on Differential Gene Expression Analysis.
• Validate drug target binding sites using GROMACS Protein-Ligand Simulation Protocol or submit targets to our Molecular Dynamics Services.