Tutorial 01 — scRNA-seq: QC & Preprocessing

Load a two-sample dataset (ifnb), filter cells, normalize, pick HVGs, scale

Author

Single Cell RNA-seq Workshop

Running the code in this tutorial

Every code chunk here is tagged with #| eval: false, so the published page shows the code without running it. To run it yourself, open the downloaded .qmd in RStudio:

Run one chunk: click the green ▶ (Run Current Chunk) at the chunk’s top-right, or press Ctrl/Cmd + Enter. This runs the chunk regardless of its eval setting — the easiest way to work through the tutorial interactively.
Run a chunk on render: change that chunk’s #| eval: false to #| eval: true (or delete the line) so it executes when the document is rendered.
Render the whole tutorial: click the Render button in RStudio (or run quarto render in a terminal) to execute the chunks top-to-bottom and knit a finished HTML. To run everything on render, set eval: true once in the YAML header at the top.

About this tutorial

The first hands-on tutorial walking the standard Seurat scRNA-seq workflow on a laptop-friendly dataset. We use the ifnb dataset (Kang et al. 2017) — peripheral blood mononuclear cells (PBMCs) from eight lupus patients, split into a control sample and an interferon-β-stimulated sample. Two samples, ~24,000 cells, no 10x click-through.

Note

Companion lecture: Lecture 01 — Pre-processing: QC, doublets, normalization, HVGs · Companion book chapter: Chapter 1 — Preprocessing — the long-form prose treatment of this tutorial’s material, with cross-references to the prerequisite appendices. · HPC version (Talapas): Talapas analysis pipeline — 01_qc_preprocessing.R

Why two samples from the start? Because Tutorial 05 will teach integration — and integration only makes sense if you’ve already loaded, QC-filtered, and clustered samples that come from different conditions / batches. We carry both STIM and CTRL cells through every step.

Tip

Want to run this on the cluster, or scale up later? The Talapas HPC track runs this same ifnb workflow as R + SLURM scripts — see Tutorial — SLURM Basics. To go all the way from raw FASTQs on a larger real dataset, see the optional Full Talapas run from raw FASTQs (a self-paced bonus track).

How to use this page

The rendered HTML shows the code but does not execute it (chunks default to eval: false). To actually run the analysis:

Download the .qmd source: Tutorial_01_QC_Preprocessing.qmd. If your browser saves the file as Tutorial_01_QC_Preprocessing.qmd.txt, drop the trailing .txt so the filename ends in .qmd, then open it in RStudio.
Install the ifnb dataset once (one R command — see Setup).
Work through the chunks, flipping eval: true as you go.

Dataset — `ifnb` (Kang et al. 2017)

What it is. Peripheral blood mononuclear cells (PBMCs) from eight lupus patients. Each donor’s PBMCs were split into two aliquots:

CTRL — control (no stimulation), ~12,000 cells
STIM — treated with interferon-β for 6 hours, ~12,000 cells

Together ~24,000 annotated singlets across two samples on the 10x Chromium 3’ v1 assay (the muscData build of ifnb is larger than the older SeuratData ifnb you may have seen quoted at ~14k). Cell-type labels (e.g. CD14+ Monocytes, CD4 T cells, B cells, NK cells, Dendritic cells) are bundled in the metadata so you can sanity-check clustering and integration without doing your own annotation.

Reference: Kang et al. (2017) Nat Biotechnol 36: 89–94. doi:10.1038/nbt.4042
Origin: GEO accession GSE96583
Distributed via: the Bioconductor muscData package as Kang18_8vs8() — cached locally through ExperimentHub on first download
Used in: Tutorials 01 → 04 of this workshop, and the official Seurat integration vignette

Note

Why this dataset for teaching. It’s small enough to run on a laptop, it has a clean two-condition structure (so integration is meaningful), and the cell types are PBMCs — well-known, well-marked, and easy to reason about. The interferon response is a strong, biologically real signal that drives both batch-style separation (use case for integration) and condition-level differential expression (use case for pseudobulk DE in Tutorial 06).

Disk-space budget

Item	Size
`ifnb` raw counts (compressed cache)	~25 MB
`ifnb` Seurat object in memory	~600 MB
`.rds` snapshots saved between tutorials	~500 MB each

Plan for ~1.5 GB total including intermediate .rds files for Tutorials 01–08.

Data dictionary — what’s in the loaded `ifnb` Seurat object

Slot	Class	Contents
`ifnb@assays$RNA@counts`	`dgCMatrix`	Raw UMI counts. ~35,000 genes × ~24,000 cells (sparse).
`ifnb@meta.data$stim`	factor	`"CTRL"` or `"STIM"` — the sample-of-origin label. This becomes the integration covariate.
`ifnb@meta.data$seurat_annotations`	factor	Author-curated cell-type labels from the original paper. Useful as ground truth in Tutorial 03.
`ifnb@meta.data$nCount_RNA`, `nFeature_RNA`	numeric	Auto-populated by `CreateSeuratObject` — UMIs and unique genes per cell.

Learning goals

By the end of this tutorial you will be able to:

Fetch a published scRNA-seq dataset from Bioconductor’s muscData package and assemble it as a Seurat object
Inspect a multi-sample Seurat object and split it by sample
Add per-cell QC metrics (percent.mt, percent.rb) and visualize them per sample
Choose filtering thresholds that are defensible for your dataset
Run NormalizeData → FindVariableFeatures → ScaleData and understand what each does
Save the QC’d object so Tutorial 02 can pick it up

Common errors / things that bite

library(muscData) errors with “there is no package called ‘muscData’” — it should already be installed from the Tutorial 00 setup; if it’s missing, re-run that setup block.

scDblFinder returns a different result than the workshop — scDblFinder is stochastic. The workshop uses set.seed(2026) immediately before calling it; if you set the seed earlier or not at all, the doublet calls will differ. The 5–10% range is what to expect; exact cell IDs will vary.

SoupX errors with “no raw matrix found” (or subscript out of bounds from load10X) — this is correct behaviour. The Kang18_8vs8() download only ships the filtered counts. SoupX needs both. The SoupX block (code) is an illustrative template for your own 10x data and now stops with a clear message if you run it with the placeholder path. Run it for real on data that has a Cell Ranger raw_feature_bc_matrix/ folder.

Warnings (not errors) during doublet detection — when Step 6 runs you’ll see Layer 'data'/'scale.data' is empty and 'normalizeCounts'/'librarySizeFactors' is deprecated. Both are harmless: doublet detection runs on raw counts before normalization (so those layers are legitimately empty), and the deprecation notes come from scran/scuttle inside scDblFinder, not your code. As long as you see Creating … artificial doublets, it’s working.

Plotting warnings before normalization (harmless) — the QC plots (code) and doublet plots (code) draw per-cell QC metadata (nFeature_RNA, nCount_RNA, …) before normalization, so Seurat prints Default search for "data" layer … utilizing "counts" layer instead. The QC values are metadata, so the layer it uses is irrelevant. Likewise the HVG plot (code) uses a log10 axis and prints log-10 transformation introduced infinite values for genes with zero mean expression (they’re not variable features). The tutorial wraps these plots in suppressWarnings() so you shouldn’t see them, but they’re safe either way.

Gene names look like SYMBOL_ENSG… / SYMBOL-ENSG… — muscData’s Kang18_8vs8() stores each feature as SYMBOL_ENSEMBLID (e.g. ISG15_ENSG00000187608). The loader strips the trailing Ensembl-ID suffix so only the gene symbol remains (ISG15), which is what canonical markers, Azimuth, and org.Hs.eg.db match on. It strips only the Ensembl suffix — not everything after the first dash — so real dashed symbols (HLA-A, MT-CO1) are preserved, and it reports how many names it simplified.

Solutions / instructor copy

The Think about it prompts on this page have inline answers (click “Show answer” / “Click for answer” to expand). For the full runnable instructor copy with every chunk pre-evaluated and outputs shown:

For students: the same .qmd you downloaded already contains every solution inline — just expand the collapsed answers
For instructors: render with quarto render --profile solutions from the Exercise_Folder/ directory. The fully-evaluated HTML is written to docs/Exercise_Folder/_solutions/Tutorial_01_QC_Preprocessing.html. See Exercise_Folder/_quarto-solutions.yml for the build profile

Setup

Set up your project folder first — read this once

These tutorials assume a simple project layout. Make a project folder, put a scripts/ subfolder inside it, and save every downloaded tutorial .qmd into scripts/:

my_scrnaseq_project/
├── scripts/      ← put the downloaded Tutorial_*.qmd files here, and run them from here
├── data/         ← auto-created on first run: the .rds objects passed between tutorials
└── output/       ← auto-created on first run: each module's figures (Mod1/, Mod2/, …)

Run each tutorial from scripts/ (open the .qmd in RStudio from there). The setup chunk below then creates data/ and output/ as siblings of scripts/ — i.e. at the project root — using the paths ../data and ../output/Mod1. (The .. is “one level up from scripts/”.) Two practical rules:

data/ is the hand-off. Tutorial 01 writes ../data/ifnb_preprocessed.rds; Tutorial 02 reads it back from the same ../data/. Keep all the .qmd files together in scripts/ so they share one project data/.
Can’t find a file? The setup chunk prints the exact absolute path (This module writes its figures & tables to: …). If it’s not where you expect, check getwd() — it should be your scripts/ folder. See the downloadable README (linked on the Materials page) for the full layout.

# Packages are installed once in Tutorial 00 (Setup) — here we just load them.
library(Seurat)
library(muscData)
library(SingleCellExperiment)
library(tidyverse)
library(patchwork)   # combine + annotate multi-panel Seurat plots

set.seed(2026)

# --- Figure saving --------------------------------------------------------
# save_fig() writes every figure as a .png (for viewing) AND a .svg (vector,
# editable in Illustrator / Inkscape for a manuscript). Pass the .png path; the
# .svg is written alongside with the same basename. (.svg uses the 'svglite'
# package, installed in Tutorial 00.)
save_fig <- function(filename, plot, width, height, dpi = 300, ...) {
  ggplot2::ggsave(filename, plot, width = width, height = height, dpi = dpi, ...)
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch(ggplot2::ggsave(svg_path, plot, width = width, height = height, ...),
           error = function(e) message("  (could not write ", basename(svg_path),
                                       " - install 'svglite'? ", conditionMessage(e), ")"))
}

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables.
# Every figure/table chunk below writes a file named  Mod1_C<chunk>_<name>
# into ../output/Mod1/ so it can be cross-referenced from the rest of the site.
#   Mod1 = Module 1 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod1"
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create("../data",  showWarnings = FALSE)   # shared .rds hand-off folder (sibling of scripts/, created one level up)
message("This module writes its figures & tables to: ", normalizePath(out_dir))

Note

muscData (installed in the Tutorial 00 setup) ships the data through Bioconductor’s ExperimentHub: the first call to Kang18_8vs8() downloads ~25 MB into your local cache, and subsequent calls are instant.

Tip

Can’t install muscData? Download the preassembled ifnb_raw.rds instead and replace the load block below with ifnb <- readRDS("../data/ifnb_raw.rds"). It is the exact object Step 1 produces, so you pick up right at QC. Details on the Datasets page.

Step 1 — Load the `ifnb` dataset

# Fetch the Kang et al. 2017 PBMC dataset (CTRL + IFN-beta STIM) from
# Bioconductor's ExperimentHub. The first call downloads ~25 MB and caches it;
# every subsequent call is a no-op.
sce <- Kang18_8vs8()

# Keep only true singlets that have an author-assigned cell-type label (drop the
# multiplexed multiplets and unlabeled droplets). This leaves ~24,000 cells — the
# muscData build is larger than the historical SeuratData `ifnb` (~14k).
sce <- sce[, sce$multiplets == "singlet" & !is.na(sce$cell)]

# Gene (feature) name cleanup. muscData's Kang18_8vs8 stores each feature as
# SYMBOL_ENSEMBLID (e.g. "ISG15_ENSG00000187608"). We keep only the readable
# SYMBOL so canonical markers ("CD14", "MS4A1", ...) and reference tools (Azimuth,
# org.Hs.eg.db) match by name. Without this, Seurat would also rewrite the '_' to
# '-' and warn ("Feature names cannot have underscores ...").
#
# We strip ONLY a trailing "_ENSEMBLID" suffix, on the original underscore names.
# We deliberately do NOT cut at the first dash — real symbols contain dashes
# (HLA-A, MT-CO1, HLA-DRB1) and must be preserved. If a feature has no symbol
# (blank / pure Ensembl ID) we keep its original name so nothing is dropped, and
# make.unique() de-duplicates any symbols that collapse to the same name.
cm   <- counts(sce)
orig <- rownames(cm)
sym  <- sub("_ENS[A-Z]*[0-9]+(\\.[0-9]+)?$", "", orig)   # SYMBOL_ENSG... -> SYMBOL
blank <- is.na(sym) | sym == "" | grepl("^ENS[A-Z]*[0-9]+", sym)
sym[blank] <- orig[blank]
rownames(cm) <- make.unique(gsub("_", "-", sym))         # Seurat-safe, unique
n_changed <- sum(rownames(cm) != orig)
if (n_changed > 0) {
  i <- which(rownames(cm) != orig)[1]
  message(n_changed, " of ", length(orig), " feature names simplified to gene ",
          "symbols (e.g. '", orig[i], "' -> '", rownames(cm)[i], "')")
}

# Convert to a Seurat object, mirroring the metadata column names the rest of
# this tutorial series expects:
#   - `stim` as a factor with levels c("CTRL", "STIM")
#   - `seurat_annotations` for the author-curated cell-type labels
ifnb <- CreateSeuratObject(
  counts    = cm,
  meta.data = as.data.frame(colData(sce))
)
ifnb$stim <- factor(toupper(ifnb$stim), levels = c("CTRL", "STIM"))
ifnb$seurat_annotations <- factor(ifnb$cell)

ifnb

# What samples are in here?
table(ifnb$stim)

# Author-curated cell-type labels (we'll use these as ground truth in Tut 03)
table(ifnb$seurat_annotations)

# --- Tables out: per-sample and per-cell-type cell counts ------------------
tibble::enframe(table(ifnb$stim),
                name = "sample", value = "n_cells") |>
  dplyr::mutate(n_cells = as.integer(n_cells)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C2_cells_per_sample.csv"))

tibble::enframe(table(ifnb$seurat_annotations),
                name = "cell_type", value = "n_cells") |>
  dplyr::mutate(n_cells = as.integer(n_cells)) |>
  dplyr::arrange(dplyr::desc(n_cells)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C2_cells_per_celltype.csv"))

Think about it

Why do real integration tutorials almost always have at least two samples?
What’s the difference between a “sample” (stim) and a “cell type” (seurat_annotations) in this dataset?

Show answers

Integration aligns the shared biological signal between samples while preserving any real differences between them. With one sample there’s nothing to align — clustering on PCs alone is enough. The ifnb dataset has CTRL and STIM as two samples drawn from the same donors, which is the perfect minimal setup to teach the alignment step.
stim is the technical / experimental covariate — which condition the cell came from. seurat_annotations is the biological identity — what kind of cell it is. Integration aims to remove stim-driven structure from the embedding without erasing seurat_annotations-driven structure.

Step 2 — Inspect dimensions, split by sample

The ifnb object already arrives as a single Seurat object containing both conditions. For a two-sample QC workflow, it’s useful to look at QC metrics per sample rather than globally — that way you’ll catch a sample-specific quality issue (e.g. one sample with much lower depth) before it confounds the analysis.

# Record the starting matrix dimensions
dim_raw <- dim(ifnb)
dim_raw

# --- Table out: starting matrix dimensions (features x barcodes) -----------
tibble::tibble(
  stage    = "Step 1 — loaded (singlets, annotated)",
  features = dim_raw[1],
  barcodes = dim_raw[2]
) |>
  readr::write_csv(file.path(out_dir, "Mod1_C3_starting_dimensions.csv"))

Think about it

The ifnb object holds both CTRL and STIM cells in a single counts matrix. What information distinguishes a CTRL barcode from a STIM barcode?
Why might it matter to split visualizations by stim even though clustering happens on the joint matrix?

Show answers

The stim column in meta.data. Each barcode also carries the -1 / -2 GEM-well suffix from Cell Ranger, but you should never rely on barcode suffixes for sample membership — the metadata is the source of truth.
Per-sample plots reveal batch-level QC issues (lower median depth in one sample, different mitochondrial fractions, different cell counts per cluster) that can be hidden by a global histogram. If one sample’s distribution is shifted, you may need different filter thresholds for it — or you may need to investigate before going further.

Step 3 — Add QC metrics

ifnb[["percent.mt"]] <- PercentageFeatureSet(ifnb, pattern = "^MT-")
ifnb[["percent.rb"]] <- PercentageFeatureSet(ifnb, pattern = "^RP[SL]")

head(ifnb@meta.data)

# --- Table out: per-sample summary of the QC metrics just computed ---------
ifnb@meta.data |>
  dplyr::group_by(stim) |>
  dplyr::summarise(
    n_cells           = dplyr::n(),
    median_nFeature   = median(nFeature_RNA),
    median_nCount     = median(nCount_RNA),
    median_percent_mt = round(median(percent.mt), 2),
    median_percent_rb = round(median(percent.rb), 2),
    .groups = "drop"
  ) |>
  readr::write_csv(file.path(out_dir, "Mod1_C4_qc_metric_summary.csv"))

Note

The ifnb raw matrix uses HUGO gene symbols, so ^MT- matches the 13 protein-coding mitochondrial genes (MT-ND1, MT-CO1, …) and ^RP[SL] matches the ribosomal protein genes (RPS*, RPL*). For mouse data the conventions are lowercase: ^mt- and ^Rp[sl].

Think about it

What biological signal does percent.mt capture?
Why might percent.rb (ribosomal-protein fraction) be useful?

Show answers

Dying or stressed cells leak cytoplasmic RNA while retaining mitochondrial RNA, so their mitochondrial fraction rises. High percent.mt flags apoptotic / compromised cells.
Ribosomal-protein genes are highly expressed and broadly variable across cell states. A very low or very high percent.rb per cluster is sometimes a tell-tale sign of a metabolic / activation difference rather than a cell-type difference — useful to separate “this cluster is biologically distinct” from “this cluster is just translationally quiet.”

Step 4 — Visualize QC, per sample

# Violin panel: four QC metrics, split by sample, fully labelled.
# These features are per-cell QC METADATA plotted before normalization, so Seurat
# warns that the "data" layer is empty and falls back to "counts". That's
# irrelevant here (metadata isn't read from an expression layer), so we silence
# the expected note.
p_qc_vln <- suppressWarnings(
  VlnPlot(ifnb,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
        group.by = "stim",
        ncol = 4,
        pt.size = 0) &
  xlab("Sample (stim condition)") &
  theme(plot.title = element_text(size = 11)))
p_qc_vln <- p_qc_vln +
  patchwork::plot_annotation(
    title    = "Per-cell QC metrics by sample — ifnb (CTRL vs IFN-β STIM)",
    subtitle = "Genes/cell, UMIs/cell, mitochondrial % and ribosomal % distributions",
    caption  = "Module 1 · QC & Preprocessing"
  )
p_qc_vln

# Scatter: UMIs vs unique genes, coloured by sample, with a linear trend
p_qc_scatter <- suppressWarnings(FeatureScatter(ifnb,
               feature1 = "nCount_RNA",
               feature2 = "nFeature_RNA",
               group.by = "stim")) +
  geom_smooth(method = "lm") +
  labs(
    title    = "Sequencing depth vs library complexity",
    subtitle = "Each point is a cell; the line is a per-group linear fit",
    x        = "nCount_RNA (total UMIs per cell)",
    y        = "nFeature_RNA (unique genes per cell)",
    colour   = "Sample"
  )
p_qc_scatter

# --- Figures out -----------------------------------------------------------
save_fig(file.path(out_dir, "Mod1_C5_qc_violins.png"),  p_qc_vln,
       width = 12, height = 4, dpi = 300)
save_fig(file.path(out_dir, "Mod1_C5_qc_scatter.png"),  p_qc_scatter,
       width = 6,  height = 5, dpi = 300)

Reading the output

Violin panel — one violin per sample (x-axis = CTRL vs STIM) for each of the four QC metrics; the width of a violin at a given height is how many cells sit at that value, so you’re reading the shape of each per-cell distribution. Healthy cells form the fat central body; the thin upper tail on nFeature_RNA/nCount_RNA is where doublets hide, and the upper tail on percent.mt is where dying cells sit. Scatter — each point is one cell: x = total UMIs (nCount_RNA), y = unique genes (nFeature_RNA). The two scale together along a tight curve; points far below the trend are low-complexity/dying cells, points far above at the top-right are doublet-like. The per-sample fit lines let you spot a sample with systematically different depth. (scater::isOutlier() is a data-driven alternative to eyeballing these thresholds.)

Think about it

Compare the violins between CTRL and STIM. Are the medians similar? Does either sample have a much fatter upper tail on nFeature_RNA?
What would a doublet look like on the nCount_RNA vs nFeature_RNA scatter? A dying cell?

Show answers

In ifnb the two samples are very similar in median QC metrics — that’s a feature, not a bug, because the donors are matched. A noticeably fatter upper tail in one sample would suggest more doublets there (perhaps a higher loading concentration). A noticeably lower median nFeature_RNA would suggest sequencing depth or cell-quality differences.
Doublets sit high on both axes, often above the main trend line. Dying cells are low on both and tend to have high percent.mt (color the scatter by percent.mt or facet by it).

Step 5 — Optional: ambient-RNA correction with SoupX

In a 10x droplet experiment, lysed-cell debris in the bulk solution leaks ambient mRNA into every droplet. After Cell Ranger, this shows up as a low background of “wrong” transcripts in every cell — most visibly when a cell expresses a marker for a different cell type at low level. SoupX estimates the ambient profile and corrects each cell’s counts. It runs before QC filtering and replaces the counts matrix with a cleaned one.

Warning

SoupX needs both the raw and filtered Cell Ranger matrices (raw_feature_bc_matrix/ and filtered_feature_bc_matrix/) because it learns the soup profile from empty droplets. Kang18_8vs8() only ships the filtered counts, so you cannot run SoupX on ifnb directly. The code below is the canonical pattern — use it on your own 10x data, or in the optional Full Talapas run from raw FASTQs, where the raw matrix is available.

Optional reference pattern — left #| eval: false because it needs a raw Cell Ranger matrix (not shipped with ifnb). You don’t need to change anything here; use it on your own 10x data.

# Pattern for any 10x dataset where you have the raw Cell Ranger output.
# This block does NOT run on ifnb (which ships only *filtered* counts) and will
# NOT run as-is: edit CR_OUTS to point at a real Cell Ranger 'outs/' directory
# that contains BOTH raw_feature_bc_matrix/ and filtered_feature_bc_matrix/.
library(SoupX)
library(DropletUtils)

# 1. Point at the Cell Ranger 'outs/' directory (containing raw_feature_bc_matrix/
#    and filtered_feature_bc_matrix/)
CR_OUTS <- "path/to/cellranger/outs"   # <-- EDIT to your own Cell Ranger outs/

# Fail early with a clear message rather than a cryptic "subscript out of bounds"
# if the path is still the placeholder or is missing the two required matrices.
if (!dir.exists(file.path(CR_OUTS, "raw_feature_bc_matrix")) ||
    !dir.exists(file.path(CR_OUTS, "filtered_feature_bc_matrix"))) {
  stop("SoupX needs a real Cell Ranger 'outs/' folder containing BOTH ",
       "raw_feature_bc_matrix/ and filtered_feature_bc_matrix/.\n",
       "  '", CR_OUTS, "' is not valid — ifnb ships only the filtered counts, so ",
       "this block is a template for your OWN 10x data, not a step you run on ifnb.")
}

sc <- load10X(CR_OUTS)

# 2. (Recommended) provide preliminary clusters so SoupX can refine the soup
#    profile per cluster — compute them with a quick Seurat workflow first:
# sc <- setClusters(sc, setNames(seu$seurat_clusters, colnames(seu)))
# sc <- setDR(sc, Embeddings(seu, "umap"))

# 3. Estimate the ambient contamination fraction
sc  <- autoEstCont(sc)
rho <- sc$fit$rhoEst
message("Estimated ambient fraction: ", round(rho, 3))

# 4. Adjust the counts (round to integers so downstream UMI tools are happy)
adj <- adjustCounts(sc, roundToInt = TRUE)

# 5. Use `adj` as the input to CreateSeuratObject() in Step 1 instead of the raw
#    filtered matrix; everything else in this tutorial then runs on the cleaned counts.

Think about it

Why does SoupX need the raw (pre-filter) matrix? What information does it pull from there?
When is SoupX overkill? When is it essential?

Show answers

The raw matrix contains droplets that captured no real cell — they hold only ambient RNA. SoupX uses those empty droplets to estimate the soup’s expression profile, then subtracts a fraction of that profile from every called cell.
Overkill: small contamination (rho ≲ 0.05), well-prepared samples, no marker bleed-through. Essential: complex tissues with high cell death / large cell-type differences (tumor microenvironment, mucosal tissues), or any time you see a marker like HBB or IGKC showing up faintly in cells that shouldn’t express it.

Step 6 — Doublet detection with `scDblFinder`

A doublet is one barcode’s worth of counts that came from two cells getting trapped in the same droplet. They show up as cells with anomalously high nCount_RNA and nFeature_RNA and a mixed expression profile (e.g. T-cell markers AND B-cell markers in the same cell). The simple nFeature_RNA < 2500 filter applied later (code) catches some, but doublets that pair two different cell types of similar size will pass that filter — you need a model.

scDblFinder simulates artificial doublets, trains a classifier to distinguish them from singletons, and returns a per-cell score and class. Critically, doublets form within a sample, not across samples, so we run it per sample.

library(scDblFinder)

set.seed(2026)

# Run scDblFinder per sample, then merge the results back into the Seurat object.
# NOTE on the messages you'll see here (all harmless, NOT errors):
#   - "Layer 'data'/'scale.data' is empty" — doublet detection runs BEFORE
#     normalization (Step 8), so those layers aren't populated yet; scDblFinder
#     works on raw counts and normalizes internally. We wrap the conversion in
#     suppressWarnings() to quiet those expected notes.
#   - "'normalizeCounts'/'librarySizeFactors' is deprecated" — these come from
#     scran/scuttle *inside* scDblFinder, not from this code; safe to ignore.
ifnb_list <- SplitObject(ifnb, split.by = "stim")
ifnb_list <- lapply(ifnb_list, function(x) {
  sce <- suppressWarnings(as.SingleCellExperiment(x))
  sce <- scDblFinder(sce, returnType = "sce")
  x$scDblFinder.score <- colData(sce)$scDblFinder.score
  x$scDblFinder.class <- colData(sce)$scDblFinder.class    # "singlet" or "doublet"
  x
})
ifnb <- merge(ifnb_list[[1]], y = ifnb_list[-1])

# `merge()` in Seurat v5 leaves the RNA assay with per-sample split layers
# (counts.CTRL / counts.STIM, etc.). Downstream tools like FindMarkers expect
# a single joined layer, so rejoin them now.
ifnb <- JoinLayers(ifnb)

# How many doublets were called per sample?
table(sample = ifnb$stim, class = ifnb$scDblFinder.class)

# --- Table out: doublet calls per sample -----------------------------------
as.data.frame(table(sample = ifnb$stim, class = ifnb$scDblFinder.class)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C7_doublet_counts.csv"))

# Inspect: doublets should sit at the high end of nCount/nFeature.
# (Same pre-normalization metadata plot as Step 4 — silence the expected
# "data layer empty, using counts" note.)
p_dbl <- suppressWarnings(
  VlnPlot(ifnb,
        features = c("nFeature_RNA", "nCount_RNA"),
        group.by = "scDblFinder.class", pt.size = 0) &
  xlab("scDblFinder class") &
  theme(plot.title = element_text(size = 11)))
p_dbl <- p_dbl +
  patchwork::plot_annotation(
    title    = "Library size and complexity for called singlets vs doublets",
    subtitle = "Doublets are expected to sit at the high end of both metrics",
    caption  = "Module 1 · QC & Preprocessing"
  )
p_dbl

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod1_C7_doublet_violins.png"), p_dbl,
       width = 7, height = 4, dpi = 300)

Reading the output

The table() is a 2×2 of sample × class: the key number is the doublet fraction per sample (doublet / (singlet + doublet)), which should land around 5–10% for a normal 10x run — much higher points to chip overloading. In the violins (x-axis = singlet vs doublet), the doublet group should sit visibly higher on both nFeature_RNA and nCount_RNA; that separation is the sanity check that the classifier keyed on the expected signal. scDblFinder.score (0–1) is the per-cell confidence if you want to threshold more or less aggressively than the default class call.

Think about it

Why run scDblFinder per sample instead of on the merged object?
scDblFinder marks ~5–10% of cells as doublets in a typical 10x run. If your run shows 25%, what likely happened?
Some pipelines run doublet detection after clustering and remove doublet-enriched clusters wholesale. Why might that be more conservative than removing per-cell doublet calls?

Click for answer

Doublets are made when two cells from the same droplet load are co-encapsulated. Cells from different samples (CTRL and STIM here) were prepared in separate channels, so a STIM and a CTRL cell can never be a doublet. Running per-sample also lets each sample’s doublet-rate be modeled correctly.
Either you over-loaded the chip (the loading concentration was too high — typical for “nominal 10k cells, recovered 18k”) or there’s a quality problem (debris, large cell aggregates) making the simulated-doublet classifier confused. Inspect nCount_RNA/nFeature_RNA distributions and the ambient fraction.
Per-cell doublet calls have false positives. If you trust them and remove individual cells, you may also be removing legitimate cycling cells, large activated cells, or megakaryocytes (which look doublet-like by sheer transcriptional output). Removing doublet-enriched clusters is more conservative — only act on a cluster where the doublet fraction is clearly elevated and the markers don’t make biological sense.

Step 7 — Filter low-quality cells

ifnb <- subset(
  ifnb,
  subset = nFeature_RNA > 200 &
           nFeature_RNA < 2500 &
           percent.mt   < 5 &
           scDblFinder.class == "singlet"
)
ifnb
table(ifnb$stim)

dim_after_filter <- dim(ifnb)

# --- Table out: cells surviving the QC filter, per sample ------------------
tibble::enframe(table(ifnb$stim),
                name = "sample", value = "n_cells_after_filter") |>
  dplyr::mutate(n_cells_after_filter = as.integer(n_cells_after_filter)) |>
  readr::write_csv(file.path(out_dir, "Mod1_C8_cells_after_filter.csv"))

Think about it

Did the CTRL or STIM sample lose more cells in this filter? If they lost noticeably different fractions, what would that suggest?
How would you set thresholds empirically rather than by default?

Show answers

If one sample loses a much larger fraction at percent.mt < 5 it usually means its cells were under more stress at dissociation or stimulation. For ifnb, IFN-β stimulation can mildly elevate percent.mt in some donors — flag the sample, but don’t necessarily relax the threshold without thinking.
Look at the per-sample violin / histogram shape; choose thresholds that trim the clear tails (not the mode). For an automated rule, use scater::isOutlier() or median ± k·MAD — and consider running it per sample so a globally lopsided distribution doesn’t push the cutoff in the wrong direction for one sample.

Step 8 — Normalize

ifnb <- NormalizeData(ifnb)
# Equivalent explicit form:
# NormalizeData(ifnb,
#               normalization.method = "LogNormalize",
#               scale.factor = 10000)

dim_after_normalize <- dim(ifnb)

Think about it

What does LogNormalize do algebraically?
When would you prefer SCTransform() instead?

Show answers

For each cell: divide each gene’s count by the cell’s total count, multiply by scale.factor (default 10 000), then take log1p. In formula form: log(1 + (count_gc / total_c) * 10000).
SCTransform models the mean–variance relationship per gene with a regularized negative binomial, which is usually better for low-count datasets and more faithful downstream variance. It also replaces the NormalizeData → FindVariableFeatures → ScaleData chain in a single call. For ifnb either approach works fine.

Step 9 — Highly variable features

ifnb <- FindVariableFeatures(
  ifnb,
  selection.method = "vst",
  nfeatures = 2000
)

top10 <- head(VariableFeatures(ifnb), 10)
plot1 <- VariableFeaturePlot(ifnb)
p_hvg <- LabelPoints(plot = plot1, points = top10, repel = TRUE) +
  labs(
    title    = "Highly variable genes (vst) — top 2,000 selected",
    subtitle = "Top 10 most variable genes labelled; many are interferon-stimulated genes",
    x        = "Average expression (log scale)",
    y        = "Standardized variance",
    colour   = "Selected as HVG"
  )
# VariableFeaturePlot draws mean expression on a log10 x-axis. Genes with zero
# mean (undetected across all retained cells) become -Inf and ggplot drops them
# with "log-10 transformation introduced infinite values" — harmless, since those
# genes are not variable features. Silence it at draw and save time.
suppressWarnings(print(p_hvg))

# --- Outputs: HVG plot + the labelled top-10 genes as a table --------------
suppressWarnings(
  save_fig(file.path(out_dir, "Mod1_C10_variable_features.png"), p_hvg,
         width = 7, height = 5, dpi = 300))

tibble::tibble(rank = seq_along(top10), gene = top10) |>
  readr::write_csv(file.path(out_dir, "Mod1_C10_top10_hvgs.csv"))

Reading the output

Each dot is a gene: x-axis = average expression (log scale), y-axis = standardized variance (how much more variable the gene is than expected for its mean). The ~2,000 genes coloured/highlighted as variable are the ones that rise above the fitted mean–variance curve — these are the only genes PCA will use in Tutorial 02. The labelled top-10 are the most variable of all; in ifnb they’re dominated by interferon-stimulated genes (ISG15, IFI6, IFIT1…) because the STIM half of the data switches them on and the CTRL half doesn’t. Seeing ISGs at the top here is your first preview of the batch-style split you’ll visualize in Tutorial 02 and correct in Tutorial 05.

Think about it

Which genes are at the top of top10? Many will be interferon-stimulated genes (ISGs) like ISG15, IFI6, IFIT1 — why?
What happens if you set nfeatures = 500? Or nfeatures = 10000?

Show answers

Because half the dataset (the STIM sample) has IFN-β–induced ISG expression while the other half (CTRL) does not, ISGs have enormous variance across cells and rise to the top of the HVG list. This is exactly the signal we’ll see drive batch-style separation in Tutorial 02 — and it’s exactly what integration in Tutorial 05 will need to handle without erasing.
Too few (500) risks missing subtle cell-type signals and is brittle to noisy gene choices. Too many (10 000) pulls in low-information genes and dilutes the signal. 1 500–3 000 is a common sweet spot; 2 000 is a good default.

Step 10 — Scale

ifnb <- ScaleData(ifnb)

By default ScaleData only scales the variable features identified in Step 9 (~2,000 HVGs) — that’s all PCA needs, and it keeps the scale.data slot (and the saved .rds) small. If you later need scaled values for a specific non-HVG gene (e.g. to draw it on a heatmap), re-run ScaleData(ifnb, features = <those genes>) at that point.

Think about it

What does ScaleData do to each gene numerically?
Why is scaling necessary before PCA?
Could you vars.to.regress = "stim" here and skip integration in Tutorial 05?

Show answers

Centers each gene’s expression to mean 0 and scales to unit variance across cells (z-scoring).
PCA finds directions of maximum variance. Without scaling, highly-expressed genes (which inherently have large raw variance) would dominate PCs, regardless of whether they are biologically informative.
You could, but regressing on a binary covariate is a brittle form of batch correction — it removes a linear effect uniformly across all cells, but real batch effects are usually non-linear and cell-type-specific. Modern integration methods (Harmony, anchors, scVI) are much better tools. We use ScaleData here as the standard pre-PCA step and reserve the integration question for Tutorial 05.

Counts matrix progression — summary

dim_summary <- tibble::tribble(
  ~step,                                                    ~features_after,        ~barcodes_after,
  "Step 1 — Kang18_8vs8 (singlets, raw)",                    dim_raw[1],             dim_raw[2],
  "Step 7 — subset() QC filter (nFeature_RNA, percent.mt)",  dim_after_filter[1],    dim_after_filter[2],
  "Step 8 — NormalizeData (no dimension change)",            dim_after_normalize[1], dim_after_normalize[2]
)
dim_summary

# --- Table out: counts-matrix progression summary --------------------------
readr::write_csv(dim_summary, file.path(out_dir, "Mod1_C12_dimension_progression.csv"))

You should see roughly:

Step	Features after	Barcodes after
Step 1 — `Kang18_8vs8` (singlets, annotated)	~35,635	~24,500
Step 7 — `subset()` QC filter (`nFeature_RNA`, `percent.mt`)	~35,635	~23,000
Step 8 — `NormalizeData` (no dimension change)	~35,635	~23,000

(Numbers will vary by a few hundred depending on package versions.)

Save and continue to Tutorial 02

# Make sure data/ exists
dir.create("../data", showWarnings = FALSE)

saveRDS(ifnb, file = "../data/ifnb_preprocessed.rds")

Tip

The saved .rds is the starting point for Tutorial 02 — Dimensionality reduction & clustering. We deliberately keep both CTRL and STIM cells in the same object — Tutorial 02 will run PCA + clustering on the unintegrated joint data so you can see the batch effect before fixing it in Tutorial 05.

Credits

Dataset: Kang et al. (2017) Nat Biotechnol 36: 89–94, doi:10.1038/nbt.4042 — fetched via Bioconductor’s muscData::Kang18_8vs8()
Reference vignettes: Seurat integration tutorial · Seurat PBMC 3k tutorial
Inspiration for the multi-step QC pacing: scNotebooks (Sengupta lab)

About this tutorial

Dataset — ifnb (Kang et al. 2017)

Disk-space budget

Data dictionary — what’s in the loaded ifnb Seurat object

Learning goals

Setup

Step 1 — Load the ifnb dataset

Step 2 — Inspect dimensions, split by sample

Step 3 — Add QC metrics

Step 4 — Visualize QC, per sample

Step 5 — Optional: ambient-RNA correction with SoupX

Step 6 — Doublet detection with scDblFinder

Step 7 — Filter low-quality cells

Step 8 — Normalize

Step 9 — Highly variable features

Step 10 — Scale

Counts matrix progression — summary

Save and continue to Tutorial 02

Credits

Dataset — `ifnb` (Kang et al. 2017)

Data dictionary — what’s in the loaded `ifnb` Seurat object

Step 1 — Load the `ifnb` dataset

Step 6 — Doublet detection with `scDblFinder`