L01 — Pre-processing: QC, doublets, normalization, HVGs

Where this lecture fits

Goals of this lecture

  • See where pre-processing sits in the canonical pipeline
  • Learn the three core QC metrics, what they mean, and what counts as normal
  • Recognize and handle ambient-RNA contamination and droplet doublets
  • Pick a normalization + feature-selection strategy that fits your data

Note

Why this matters.

  • Every downstream geometric step (PCA, clustering, UMAP) treats the post-QC matrix as ground truth.

  • Garbage in → confidently-wrong clusters out.

  • There is no “fix the clustering” — you have to fix the QC.

What is the same about scRNA-seq

  • Same core measurement: count of reads / UMIs per gene
  • Same statistical concerns: normalization, batch effects, multiple testing, replication
  • Many downstream questions (DE, enrichment, co-expression) are the same questions at a new resolution — and reuse bulk tools (DESeq2, limma, WGCNA) after pseudobulk aggregation (Lecture 06)

What is genuinely new about scRNA-seq

  • Sparsity: ~90–95% zeros in typical droplet data — partly real biology, partly technical dropout
  • Cells are not replicates: cells from one subject are not statistically independent — using n = cells is the most common mistake in the field
  • Cell type is a latent variable you have to recover from the data
  • Integration across samples is a first-class problem, not an afterthought

Sparsity

  • Bulk: nearly every gene is detected in every sample (~5% zeros)

  • Smart-seq2: ~50–60% zeros — moderate dropout but rich per-cell signal

  • 10x 3’ v3 droplet: ~90–95% zeros — most (cell, gene) pairs are zero

  • The implication for analysis: every step has to handle counts that are mostly zero.

  • That’s why we have specialized normalisation (SCTransform), HVG selection, and graph-based clustering instead of the bulk RNA-seq toolkit.

Bulk vs. single-cell — the tools

Bulk vs. single-cell toolchain comparison

Raw data processing (which we don’t do)

Raw data processing (which we don’t do)

  • Common tools:
    • Cell Ranger (10x, official)
    • STARsolo (STAR-based, fast, flexible)
    • alevin-fry / salmon (selective alignment, very fast)
    • kb-python (kallisto | bustools) (pseudoalignment)
  • Final product: a barcode(cell) × feature(gene) count matrix
  • We will start with the count matrix

The downstream workflow at a glance

  • That cell × gene matrix is the input to every analytical step the rest of the workshop covers
  • The eight steps fall naturally into four chunks, each with a dedicated lecture:
    • Pre-processing (Lec 01 / Tut 01): QC → doublets → normalisation → HVGs
    • Geometry (Lec 02 + 05 / Tut 02 + 05): scale → PCA → neighbour graph → Leiden → t-SNE/UMAP → integration
    • Annotation (Lec 03, Lec 04 / Tut 03, Tut 04): markers, manual annotation, reference-based annotation (Azimuth)
    • Differential expression & function (Lec 06–07 / Tut 06–07): pseudobulk DESeq2 → functional enrichment
    • Differential abundance (Lec 08)
  • *Bonus modules** extend the same conceptual toolkit to other data modalities

Tip

A 4-step mental compression. From here on out, every chunk of the pipeline (a) cleans something, (b) puts things on a sensible scale, (c) finds structure, (d) names that structure with biology.

Generic scRNA-seq pipeline at a glance

Specific Seurat workflow

Quality control of the feature by barcode matrix

Why Quality Control?



A droplet that “looks like a cell” in your matrix may actually be

  • an empty droplet with ambient soup

  • a dying cell leaking cytoplasm

  • two (or more) cells in one droplet

  • a real cell

QC is how we triage them.

Why Quality Control?

Quality control — raster reference

Three core metrics

  • nFeature_RNA: genes detected per cell — proxy for capture quality
  • nCount_RNA: total UMI counts per cell — proxy for sequencing depth in that cell
  • percent.mt: % reads mapping to mitochondrial genes — proxy for stress / dying cells (cytoplasmic mRNA leaks before mitochondrial mRNA)

Normal ranges — a starting cheat sheet

Note

Typical “healthy” thresholds for 10x 3’ v3 on dissociated tissue. These are starting points, not rules — always look at the distributions first.

Metric Typical lower Typical upper Notes
nFeature_RNA (genes/cell) 200–500 2,500–6,000 Upper bound is dataset-specific; doublets push past it
nCount_RNA (UMIs/cell) 500–1,000 ~50,000 Very high counts often = doublets, very low = empty
percent.mt — 5–10% (PBMC), 15–20% (tumor / fragile tissue) Nuclei (snRNA) should be < ~1–2%

Warning

Watch-outs.

  • Brain / muscle / cardiomyocytes naturally have high mitochondrial content — a 5% threshold will throw away your biology.
  • snRNA-seq should have almost no mitochondrial reads; if it does, your nuclei prep contaminated cytoplasm.
  • A bimodal nCount_RNA distribution often means two cell sizes (e.g., lymphocytes + neutrophils), not doublets — don’t gate it away.

What QC distributions actually look like

  • Each violin shows what each “category” of barcode looks like in a typical 10x v3 dataset
  • Real cells form the central bulk you want to keep;
    • doublets sit at the high end;
    • dying cells sit at the low-feature / high-mt corner;
    • empty droplets pile up at very low feature counts
  • The correct cutoffs aren’t universal — they depend on your tissue and protocol

Droplets can also be empty — or contaminated

Ambient RNA — raster reference

  • Many barcodes correspond to empty droplets with only ambient RNA — filter with EmptyDrops or the Cell Ranger knee-point
  • Even real cells are contaminated by ambient RNA soup; estimate and subtract with SoupX, DecontX, or CellBender

Doublet detection

Doublets

  • Droplets sometimes capture two cells
  • Run doublet detection on each sample before integration
  • Expected doublet rate scales with loading: ~0.8–1% per 1,000 cells loaded on 10x Chromium → ~8% at 10k cells, ~16% at 20k

Doublets

  • Simulate artificial doublets by averaging random pairs of real cells
  • Project real + simulated cells together in reduced dimensions
  • Real cells that sit close to many simulated doublets are flagged as likely doublets

Doublets

Tool Ecosystem Notes
scDblFinder R / Bioc Fast, recommended default
DoubletFinder R Needs pK sweep
Scrublet Python Classic choice in scanpy flows

Warning

Doublet detection runs per sample, not on a merged object.

  • Inter-sample mixtures are not doublets — they’re integration noise
  • Pretending they are will silently drop real cells.

Note

Seurat v5 note.

  • After merge() in Seurat v5, the RNA assay stores split layers (counts.CTRL, counts.STIM, …).
  • Run JoinLayers() after merging the per-sample objects so that downstream steps see a single layer.
  • The tutorial does this automatically after the doublet-detection merge.

Doublet Score Distribution

Doublet Score Distribution



  • Doublet detectors output a continuous score per cell — high = “looks like a simulated doublet”
  • In a healthy run the distribution is bimodal: a tall left peak (real cells) and a small right shoulder (likely doublets)
  • The threshold is set at the expected doublet rate implied by your loading concentration (e.g. 8% at 10k cells loaded on Chromium)
  • If the distribution is unimodal, your sample either has no detectable doublets (rare in a 10k-cell run) or your simulator parameters need tuning

Normalization

Why normalize?



  • Cells are sequenced to wildly different depths (factor-of-10 spread is common).
  • Without correction, a “highly expressed” gene in a deeply sequenced cell looks the same as a moderately expressed one in a shallow cell.
  • Normalization puts cells on a comparable scale before downstream math.

Two dominant approaches

  • Log-normalization (NormalizeData in Seurat) — fast, simple, the default in most tutorials
  • Variance-stabilizing (SCTransform in Seurat) — models the negative-binomial mean–variance relationship; often better for low-count cells but slower
# Seurat — log-normalization
seu <- NormalizeData(seu) |>
  FindVariableFeatures(nfeatures = 2000) |>
  ScaleData()

# Seurat — SCTransform (alternative, replaces the three calls above)
seu <- SCTransform(seu, verbose = FALSE)

Tip

When to reach for SCTransform: noisy or very low-depth data (~3,000 UMIs/cell median), or when you can see depth artefacts in your UMAP. For typical 10x v3 data, log-normalization is fine.

Warning

Don’t normalize twice. SCTransform replaces NormalizeData + ScaleData. Running both on the same object stacks transforms and breaks downstream assumptions.

Log-norm vs SCTransform

.

  • Log-normalization removes depth but leaves a mean–variance trend for medium/high-expression genes
  • SCTransform’s regularized NB-GLM residuals stay flat across expression levels

Feature selection (HVGs)

Why Highly Variable Genes (HVGs)?

  • Use the ~1,000–3,000 most variable genes for downstream geometry
  • Dramatically reduces noise from ubiquitously-expressed housekeeping genes
  • Functions
    • Seurat: FindVariableFeatures()
    • scanpy: sc.pp.highly_variable_genes()

Note

Typical settings: nfeatures = 2000 is the Seurat default and a defensible choice for most datasets. Bump to 3,000–5,000 for atlas-scale or highly heterogeneous data; drop to 1,000 only if your dataset is small and unimodal.

Warning

Pitfalls in HVG selection. - Including mitochondrial / ribosomal genes as HVGs lets technical variation drive your PCs — Many pipelines exclude these patterns by default. - If you regress out cell-cycle scores, do it before HVG selection, otherwise cycling genes dominate the variable list.

What the HVG (variable feature) plot shows

What the HVG (variable feature) plot shows

  • Every gene gets a (mean, variance) pair. The grey curve is the expected mean–variance relationship for technical / housekeeping genes (roughly Poisson with overdispersion).
  • Genes that sit well above the trend are over-dispersed — their variance can’t be explained by sampling noise alone, so they’re carrying biological signal. Those are your HVGs.
  • FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2000) returns the top 2,000 such genes, ranked by standardized residual variance. VariableFeaturePlot() is exactly this picture.
  • Marker labels in the corner are the exact reason you do this: in a PBMC run the top HVGs include the genes you’d already expect (LYZ, S100A9, CD74, MS4A1, …) — you should recognise them.

Let’s Get Started

The ifnb dataset

  • The workshop dataset — the Kang et al. PBMCs we call “ifnb”

  • Single-cell RNA-seq data of human peripheral blood mononuclear cells (PBMCs).

  • Two conditions: resting (CONTROL) and interferon-beta (IFN-β) stimulated (STIM) cells

  • Loaded via muscData::Kang18_8vs8()

  • Starts at ~24,500 annotated singlets across CTRL and STIM.

  • After the QC filters in Tutorial 01 you should retain ~23,000 cells.

Tip

Doublet detection in this tutorial runs per sample (CTRL then STIM separately) with scDblFinder.

  • The two per-sample objects are then merged back with merge() and JoinLayers() before the final filter.
  • Running doublet detection on the merged object would treat cross-sample cell pairs as doublets — they are not.

NOTE — percent.rb and when to use it

Tutorial 01 also computes percent.rb (ribosomal-protein fraction). It is not a standard QC filter; rather:

  • A very high percent.rb (>50%) in a cluster often signals a metabolically active or activated state, not a low-quality cell.
  • Monitor it per sample: if one sample has dramatically higher percent.rb, check whether library prep differed.
  • Do not filter on percent.rb by default — just note it.