Chapter 1 — Preprocessing: From a Counts Matrix to a Filtered, Normalized Object

Note

Where this chapter sits. This is the long-form companion to Lecture 01 — Pre-processing and Tutorial 01 — QC & Preprocessing. Read this for the deeper “why”; the lecture is the slide-deck summary; the tutorial is the runnable code. Prerequisites: Appendix A — Single-Cell Biology Refresher, Appendix B — Statistical Foundations §§ 1–3 (Poisson/NB, mean–variance), P2 — R & RStudio (vectors, sparse matrices).

1.1 Why preprocessing exists

Single-cell RNA-seq starts with a counts matrix: rows are genes, columns are barcoded droplets, entries are integer UMI counts. This is the output that Cell Ranger or STARsolo produces from a raw 10x sequencing run. It is not an analysis-ready object. It contains, in roughly equal numbers, four categories of objects we need to separate:

1. Real cells, the population of interest, with reasonable RNA content
2. Empty droplets that captured only ambient RNA from the cell suspension
3. Doublets — droplets that captured two cells each, producing a chimeric expression profile
4. Damaged or dying cells, identifiable by elevated mitochondrial RNA

Preprocessing is the act of cleaning out (2)–(4) before downstream analysis. In a typical 10x run, perhaps 70–85% of barcodes that pass Cell Ranger’s first-pass filter are real, well-behaved cells. The remaining 15–30% are noise that — if left in — will drive PCs, clusters, and apparent cell types in your dataset. You will spend more analytical effort tracking them down later than you would spend filtering them out now.

This is not a minor housekeeping step. Preprocessing choices propagate through every downstream module: the genes you keep define the space PCA operates in, the cells you keep define the clusters you can find, and the normalization you choose sets the scale on which every later test of differential expression and abundance is computed. The single most consistent finding across community best-practice surveys is that QC and normalization decisions have outsized, hard-to-reverse effects on the final biology^1,2. The field’s honest reckoning with how much these computational choices shape conclusions is worth reading early^3,4.

This chapter walks each of these decisions: what the noise looks like, what it costs to leave it in, and how we recognize it in the data.

1.2 The standard Seurat-flavored workflow, named

Across the field there is broad agreement on the steps and on roughly the order. The Seurat eight-step workflow^5,6 — which the workshop tutorials use — is:

1. Load the counts matrix into a Seurat object
2. Compute QC metrics per cell: total UMI count, total detected gene count, mitochondrial fraction
3. (Optional) Ambient-RNA correction with SoupX
4. Doublet detection with scDblFinder
5. Filter cells using thresholds chosen from the QC distributions
6. Normalize counts to remove per-cell sequencing-depth differences
7. Find highly variable genes (HVGs) to focus downstream on biological signal
8. Scale the HVG matrix (z-score across cells) so each gene is given equal weight in PCA

The Scanpy⁷ and Bioconductor scran+scater workflows⁸ are essentially identical with different function names. The intuition behind each step is what matters; once you understand the why, swapping toolkits is mechanical. The order is not entirely arbitrary, and two ordering questions recur:

• Ambient correction and doublet detection come before filtering, because both estimate their models from the full set of barcodes (including borderline and empty droplets). Filter first and you throw away the information those models need.
• Normalization comes before HVG selection and scaling, because “variability” is only meaningful once per-cell depth has been removed — otherwise the most variable genes are simply the ones in the most deeply sequenced cells.

Tip

Read along with the tutorial. The Tutorial 01 chunks are organized in this same eight-step order. Each chunk’s Think about it prompt is meant to be answered with the conceptual material in this chapter.

1.3 Quality-control metrics — what they measure and why

Three numbers per cell carry most of the QC information, and a fourth (percent.rb) is often worth adding. The original systematic treatment of which features actually separate low-quality cells from good ones found that these library-complexity and mitochondrial metrics dominate⁹. Table 1 summarizes the four metrics and how to read their tails; the subsections that follow explain each in turn.

Table 1: Per-cell QC metrics, their typical healthy ranges for PBMCs, and how to read each distribution’s tails. Ranges are tissue- and chemistry-dependent — always set thresholds from your own per-sample distributions rather than from this table.

Metric	What it measures	Healthy range (PBMC)	A high value means	A low value means
nCount_RNA	Total UMIs per cell (sequencing depth)	~5,000–15,000	Very active cell or a doublet	Empty, lysed, or undersequenced droplet
nFeature_RNA	Unique genes detected	~1,500–4,000	Doublet (two cells’ transcriptomes)	Low-complexity or dying cell
percent.mt	Mitochondrial-UMI fraction	~2–10%	Dying/damaged cell (membrane leaked)	Expected; near-zero is normal for snRNA-seq
percent.rb	Ribosomal-protein fraction	tissue-dependent	High translational activity (state, not quality)	Translationally quiet cell state

1.3.1 nCount_RNA — total UMIs per cell

This is just the sum of the cell’s column in the counts matrix: how many distinct RNA molecules (UMIs, not reads — see Appendix A §4) were captured for this barcode. It serves as a rough proxy for “how deeply was this cell sequenced?”

The distribution across cells is typically log-normal-shaped, with a peak at 5,000–15,000 UMIs for a healthy 10x 3’ v3 run. Cells in the left tail are either undersequenced, lysed, or empty; cells in the right tail are either particularly transcriptionally active or are doublets (two cells contributing UMIs to one barcode).

Numerical signpost: a healthy cell has at least ~1,000 UMIs; below 500 is almost always noise. The exact lower threshold depends on the tissue (some cell types — neutrophils, plasma cells — are biologically transcriptionally lower than lymphocytes). The distinction between UMIs and reads matters here: UMI counting collapses PCR duplicates so that each counted molecule is (ideally) a real captured transcript, which is what makes nCount_RNA a meaningful biological quantity rather than a sequencing artefact. The UMI deduplication that produces these counts is itself a non-trivial step (Chapter 11).

Reading the UMI-vs-gene scatter. Plotting nCount_RNA on x against nFeature_RNA on y produces one of the most informative single QC plots. Healthy cells fall on a tight, near-linear trend. Points far below the trend (many UMIs but few unique genes) are low-complexity — over-amplifying a small repertoire; common in degraded libraries or red-blood-cell contamination. Points far above the trend in the upper-right corner have both high UMI counts and many unique genes — the classic doublet signature (two cells’ worth of RNA in one barcode). Coloring the scatter by percent.mt overlays the dying-cell signal: points that sit low on the trend and show high mitochondrial fraction are apoptotic cells, not droplets to keep.

1.3.2 nFeature_RNA — number of unique genes detected

The count of distinct genes with at least one UMI in the cell. Strongly correlated with nCount_RNA but plateaus at the saturation of the cell’s expressed transcriptome. A typical PBMC has 1,500–4,000 detected genes; a metabolically active cell or doublet pushes higher.

nFeature_RNA and nCount_RNA together let you see the cell’s complexity: a cell with many UMIs but few unique genes is over-amplifying a small set of highly-expressed transcripts (a bad sign — sometimes a hemoglobin-saturated red blood cell, sometimes a degraded library). Plot them against each other; healthy cells fall on a single linear-ish trend. The ratio log10(nFeature)/log10(nCount) — sometimes called library complexity or log10GenesPerUMI — is a useful derived metric: cells far below the cloud are low-complexity and frequently low quality.

1.3.3 percent.mt — mitochondrial fraction

Compute the fraction of UMIs that map to mitochondrial-encoded genes (in human, gene symbols starting with MT-; in mouse, mt-). The mitochondrion has its own genome with ~37 genes; when a cell’s plasma membrane breaks down (apoptosis, dissociation stress), cytoplasmic mRNAs leak out faster than mitochondrial ones (which are protected by the mitochondrial membrane), so dying cells show elevated percent.mt.

Healthy cells: 2–10% (tissue-dependent). Tumor cells often run higher (15–25%). Cells above 30% are nearly always damaged. The rule is plot the distribution per sample first, then choose a threshold from the data — there is no universal cutoff. Critically, several tissues legitimately run high mitochondrial content: cardiomyocytes, skeletal muscle, brown adipose, and proliferating tissue can sit at 15–30% in perfectly healthy cells. A fixed 5% gate that is correct for PBMCs will silently delete the biology in a heart sample. The mirror-image case is single-nucleus RNA-seq (snRNA-seq), where the cytoplasm (and most mitochondria) is removed during nuclear isolation: there, percent.mt should be near zero, and elevated mitochondrial reads instead flag cytoplasmic contamination of the nuclei prep.

1.3.4 percent.rb — ribosomal-protein fraction

Ribosomal-protein genes (RPS*/RPL* in human, Rps*/Rpl* in mouse) are among the most highly expressed transcripts in any cell, and their fraction varies systematically with cell state — actively translating cells (e.g. proliferating or recently activated lymphocytes) carry a higher ribosomal load. percent.rb is therefore less a quality flag than a biological-state covariate: a cluster distinguished only by percent.rb is often a translational-activity difference rather than a distinct cell type. It is computed exactly like percent.mt (the tutorial adds both) and is worth inspecting per cluster when an annotation looks suspicious.

For full statistical detail on why we use these metrics specifically and not, say, total reads or duplication rate, see Appendix B §1 and Appendix A §6.

1.3.5 Fixed thresholds vs. data-driven outlier detection

There are two philosophies for turning these metrics into a keep/discard decision. Fixed thresholds (nFeature_RNA > 200 & < 2500 & percent.mt < 5) are transparent and reproducible but require you to look at the data and justify each number. Data-driven outlier detection — scater::isOutlier() / scuttle::isOutlier(), which flags cells more than k median-absolute-deviations (MADs) from the per-sample median — adapts to each sample automatically and is the recommended default in the Bioconductor best-practice workflow^1,8. A common, defensible choice is “remove cells > 3 MADs above the median percent.mt, and > 3 MADs below the median nFeature/nCount,” computed per sample so a single low-quality sample doesn’t distort the cutoff for the others. Whichever you choose, compute it per sample and write the exact numbers into your methods (Appendix G).

1.4 The barcode-rank knee plot — empty droplets vs cells

A complementary diagnostic is the knee plot: sort barcodes by total UMI count descending and plot rank vs total UMIs on a log–log scale.

Real cells form a high plateau on the left. Empty droplets form a low plateau on the right. Between them, the curve drops sharply through a knee that separates the two. Cell Ranger’s empty-droplet calling uses this knee (plus an EmptyDrops likelihood test below it) to call cells. The filtered_feature_bc_matrix/ you start the tutorials with has had this filter applied.

The classic knee/inflection heuristic keeps only barcodes above the steep drop, but it discards real cells that happen to have low RNA content (small resting lymphocytes, for instance). EmptyDrops¹⁰ improves on this by testing each ambiguous barcode against a statistical model of the ambient profile: a barcode is called a cell if its expression deviates significantly from the soup, even if its total count is modest. This recovers low-RNA cell types that a pure knee cutoff would lose, and it is the reason modern Cell Ranger and DropletUtils::emptyDrops() outperform a naive threshold.

The plot is still worth looking at for your own data because pathological samples produce pathological knees:

• Two knees — usually two cell-size populations (e.g. a clean lymphocyte cell line + a contaminating macrophage line)
• No clear knee — the barcode-loading was too sparse, or the cell suspension was very heterogeneous in size
• A short plateau and a long tail — over-loaded chip; many doublets; the upper plateau is contaminated with multiplets

The knee plot is the single most informative diagnostic figure for an scRNA-seq run. Save it with your raw data. DropletUtils::barcodeRanks() (R) re-derives the knee and inflection points for you.

1.5 Ambient-RNA correction

When cells lyse during dissociation (and some always do), their cytoplasmic mRNA spills into the suspension. Every droplet — including good cells — captures a small amount of this ambient RNA. The contamination is in proportion to the ambient profile (mostly the most-expressed genes from the most-abundant cell type) and is often 5–15% of a cell’s measured counts; in hemorrhagic or high-death samples it can be much higher.

The practical effect: a cell that doesn’t express gene X actually shows non-zero X due to the ambient bath. If gene X is a highly-discriminative marker for some cell type, this can blur cell-type boundaries, falsely suggest cell–cell communication, and inflate marker p-values.

1.5.1 SoupX, DecontX, and CellBender — three design philosophies

There are three widely used correction tools, and they differ in what they assume:

• SoupX¹¹ estimates the ambient profile from the empty-droplet counts and subtracts a learned fraction (ρ, the contamination rate) from each cell. ρ is estimated automatically (autoEstCont()) by looking at the expression of genes that should be zero in a given cluster but show low non-zero counts because of the soup. SoupX is fast, interpretable, and the most common choice, but it needs preliminary clusters and both the raw and filtered matrices.
• DecontX¹² (in the celda package) models each cell’s counts as a mixture of its “own” expression distribution and a contamination distribution, estimating a per-cell contamination fraction by Bayesian inference. It does not require the raw matrix, which makes it convenient when only the filtered matrix is available.
• CellBender remove-background¹³ uses a deep generative (variational autoencoder) model of the full unfiltered matrix to jointly call cells and remove ambient background and chimeric (barcode-swapped) counts. It is the most powerful but also the most computationally demanding (GPU recommended) and operates on the raw .h5.

When to run any of them:

• Always, if your suspension was complex (tumor microenvironment, multiple tissues, high cell-death rate)
• Probably not necessary for clean immune cell populations (PBMCs from healthy donors usually have ρ < 0.05)
• Required if you intend to do cell–cell communication analysis (Chapter 14) — that field is acutely sensitive to ambient bleed, because a ligand “expressed” by the wrong cell type produces a phantom interaction

SoupX and CellBender need the raw Cell Ranger matrix because the soup profile is learned from empty droplets. The workshop’s ifnb dataset only ships the filtered matrix, so the SoupX block in Tutorial 01 is illustrative; run it for real on your own 10x output (which has both matrices), or on the Full Talapas run from raw FASTQs.

For the broader file-format context, see Appendix E — Gene-Expression Databases & File Formats.

1.6 Doublet detection

A doublet is one barcode’s worth of UMIs from two cells trapped in one droplet. The 10x chemistry produces ~0.8–1% multiplet rate per 1,000 cells loaded, so a 10k-cell run is roughly 8–10% doublets and a 20k-cell run ~16–20%¹⁴. Homotypic doublets (two cells of the same type) are mostly harmless — they look like a slightly bigger cell of that type. Heterotypic doublets (two different types) are the dangerous ones: a CD4 T + B-cell doublet expresses both CD3 and MS4A1 and can form a spurious “intermediate” cluster that gets mistaken for a real transitional state. The basic gene-count filter catches the extreme cases, but similar-sized heterotypic doublets pass it — you need a model.

1.6.1 How simulation-based detectors work

The dominant detectors all share one idea: simulate artificial doublets by averaging the profiles of random pairs of real cells, embed real and simulated cells together, and flag real cells that sit in neighbourhoods dense with simulated doublets. They differ mainly in implementation and ecosystem:

• scDblFinder¹⁵ (R/Bioconductor) simulates doublets, then trains a gradient-boosted-tree classifier on a set of features (library size, neighbourhood doublet density in PCA space, etc.). Fast and accurate; the workshop default.
• DoubletFinder¹⁶ (R) is the classic Seurat-ecosystem tool; it requires a parameter sweep (pK) and an estimate of the expected doublet rate.
• Scrublet¹⁷ (Python) is the standard in Scanpy workflows and uses the same simulate-and-score logic with a k-NN classifier.

A systematic benchmark of nine such methods found that simulation-based detectors clearly outperform simple count-threshold rules, with scDblFinder and a small number of others leading on both accuracy and speed¹⁸.

1.6.2 Practical notes and the per-sample rule

• Run detection per sample, never on a merged object. Cells from different 10x channels can never co-encapsulate, so an inter-sample “doublet” is a logical impossibility — running on merged data confuses the simulator and silently mislabels integration structure as doublets.
• Doublet removal is destructive. A real activated T cell or a megakaryocyte can look slightly doublet-like by sheer transcriptional output. The conservative approach is to remove flagged doublets after clustering, only when an entire cluster is doublet-enriched and its markers make no biological sense.
• Sanity-check the rate against loading. 5,000 cells loaded → ~5% expected; 20,000 → ~20%. A 25% call from a 5k-cell run is a red flag for over-loading or debris, not a real doublet rate.
• Set a seed. All three tools are stochastic; the tutorial sets set.seed(2026) immediately before the call.

Experimental alternatives to computational detection: if you pool samples with distinct genotypes (multiple donors) or species (a human–mouse “barnyard” mix), cross-genotype doublets can be detected directly by genetic demultiplexing (demuxlet, souporcell, Vireo) or by counting reads mapping to both genomes. This both demultiplexes the pool and measures the true multiplet rate empirically — the design used by the ifnb dataset’s parent study¹⁹.

For more on why doublets specifically (and not, say, all noise) deserve a dedicated detector, see Appendix A §11.

1.7 Filtering: thresholds chosen from the distribution, not the textbook

The big mistake at this step: copying thresholds from a tutorial without looking at your data. The right thresholds depend on the tissue, the chemistry, the dissociation protocol, and the run quality. Three rules:

1. Plot the per-sample distributions before choosing thresholds. A violin or histogram per sample, on log scale for nCount and nFeature. The left tail’s shoulder is a defensible cutoff; the bare minimum is whatever clearly cuts away the empty-droplet remnants.
2. Use median ± k·MAD as a starting point. The scater::isOutlier() function does this automatically. For percent.mt, the rule “anything more than 3 MADs above the per-sample median” is conservative but data-driven⁸.
3. Document your choices. Record the exact numerical thresholds in your methods section (Appendix G). “We filtered low-quality cells” is not a defensible methods sentence.

The workshop’s ifnb filter — nFeature_RNA > 200, nFeature_RNA < 2500, percent.mt < 5, scDblFinder.class == "singlet" — is reasonable for healthy PBMCs but would over-filter a tumor or BAL fluid sample. Treat any specific thresholds as starting points. A useful discipline is to record how many cells each criterion removes (the tutorial writes a per-step cell-count table): if a single threshold is deleting a quarter of your data, that is a signal to investigate the sample, not a result to accept silently.

1.8 Normalization

After filtering, every cell still has a different nCount_RNA — its sequencing depth. Without correction, a deeper-sequenced cell has higher absolute counts for every gene, regardless of biology. Normalization removes this technical variance. This is the step where the count nature of the data — and its mean–variance relationship — matters most; the statistical background is in Appendix B §§2–3.

Note🔑 Key concept — why depth normalization comes first

A cell’s total UMI count reflects how deeply that droplet happened to be sequenced, not how much RNA the cell biologically contained. Until that technical depth is divided out, the “most variable” genes are simply those in the most deeply sequenced cells, and PC1 tracks library size instead of biology. Normalization is the hinge of the whole workflow: every later step — HVG selection, scaling, PCA, clustering, and differential expression — assumes it has already happened, which is also why you must never normalize the same object twice.

1.8.1 Log-normalization (the default)

The standard transformation is log-normalize:

\[\tilde y_{ij} = \log\!\left(\frac{y_{ij}}{C_j}\cdot 10^4 + 1\right)\]

where \(y_{ij}\) is gene \(i\)’s count in cell \(j\), \(C_j\) is the cell’s total UMI count, and the \(+1\) keeps zeros mapping to zero. This is what Seurat’s NormalizeData() does; the result lives in seu@assays$RNA$data. The “counts-per-10,000” scale factor is conventional, not magic; the log compresses the long right tail of highly expressed genes and makes the data approximately homoskedastic enough for PCA. Post-normalization, every cell has a normalized counts vector with no strong dependence on its raw library size.

1.8.2 SCTransform — variance stabilization by regression

An alternative is SCTransform²⁰, which fits a per-gene regularized negative-binomial regression of raw counts on cell depth and uses the Pearson residuals as the normalized values. Because the negative binomial explicitly models the overdispersion of count data, the residuals have approximately constant variance across the mean-expression range — which makes downstream PCA and clustering more stable, especially for low-count cells. SCTransform is preferred for very low-depth datasets and is the default for spatial transcriptomics (Visium spots pool 1–10 cells with wider count variation; Chapter 16). It is a single call that replaces NormalizeData → FindVariableFeatures → ScaleData; running it and the three-step chain on the same object stacks transforms and breaks downstream assumptions.

1.8.3 scran pooling and the broader benchmark

For very heterogeneous data with many zeros, Bioconductor’s scran computes size factors by pooling cells before deconvolving per-cell factors, which is more robust than a simple library-size factor when many genes are undetected²¹. Which transformation is “best” has been benchmarked directly: a systematic comparison found that the simple log-normalization (sometimes called shifted-log or log1p) is a remarkably strong, hard-to-beat default, with the more elaborate methods offering modest, dataset-dependent gains²². The three approaches are compared in Table 2. The practical takeaway: consistency matters more than the choice — pick one, apply it identically across all samples, and record it. Do not normalize the same object twice.

Table 2: Normalization approaches in the Seurat/Bioconductor ecosystem. Benchmarks find the simple log-normalization a strong, hard-to-beat default²²; the elaborate methods help most for low-count or highly heterogeneous data.

Method	Underlying model	Best suited to	Replaces the 3-step chain?
Log-normalize (NormalizeData)	\(\log(\text{CP10k} + 1)\) scaling	General-purpose default; well-behaved 10x data	No — follow with FindVariableFeatures + ScaleData
SCTransform	Regularized negative-binomial; Pearson residuals	Low-depth data and spatial (Visium)	Yes — single call replaces all three steps
scran pooling	Pooled-then-deconvolved size factors	Very heterogeneous data with many zeros	No — produces size factors for normalization

1.9 Highly variable genes (HVGs)

Once normalized, every gene contributes some variance to the data. Most of it is uninformative: housekeeping genes, ribosomal proteins, ubiquitously-expressed genes whose variance reflects sequencing noise more than biology. We want to focus PCA and clustering on the genes that actually distinguish cell types.

The standard tool is variance-stabilized feature selection (FindVariableFeatures(method = "vst")):

1. Fit the relationship between each gene’s mean expression and its variance across cells (a loess fit on log-mean vs log-variance)
2. For each gene, compute a standardized variance = observed variance / fitted variance-given-mean
3. Take the top N (typically 2,000) by standardized variance — these are the HVGs

The shortlist almost always includes recognizable cell-type markers (CD14, MS4A1, GNLY, etc.) plus condition-induced genes if your dataset has a strong perturbation (interferon-stimulated genes in ifnb, hypoxia markers in tumor data). Seeing genes you expect near the top of the HVG list is a quick sanity check that QC and normalization worked.

Why 2,000? Empirically, somewhere between 1,500 and 3,000 is the sweet spot for most 10x datasets:

• Below ~1,000: rare cell types disappear because their distinguishing genes don’t make the cut
• Above ~5,000: the list becomes diluted with low-information genes; PCA performance plateaus or degrades

Going lower (e.g. 500 HVGs) is sometimes justified for atlas-scale data where speed matters more than fine resolution; going higher (3,000–5,000) helps for highly heterogeneous tissue with many subtle subtypes.

Note⚙️ Key parameter — nfeatures (number of HVGs)

FindVariableFeatures(nfeatures = 2000) sets how many genes enter PCA. Default: 2000. Too few (≲1,000) and rare cell types vanish because their distinguishing genes never make the list; too many (≳5,000) and low-information genes dilute the signal so PCA plateaus or degrades. Raise it (3,000–5,000) for highly heterogeneous tissue with many subtle subtypes; lower it (≈500) only for atlas-scale data where speed outweighs fine resolution. Whatever you choose, select HVGs after normalization and — for multi-sample work — in a batch-aware way with SelectIntegrationFeatures() (Chapter 5).

Design choices worth knowing:

• Exclude technical genes. Mitochondrial and ribosomal genes are highly variable for technical reasons; including them as HVGs lets technical variation drive your PCs. Many pipelines remove ^MT-/^RP[SL] patterns from the variable list before PCA.
• Cell cycle. If cycling is a confounder, score and (optionally) regress it out before HVG selection, otherwise cell-cycle genes dominate the variable list. Whether to regress cell cycle at all is a judgement call — in a study of proliferation you would keep it.
• HVGs and integration. When you integrate multiple samples (Chapter 5), the convention is to select HVGs in a way that is robust across batches (e.g. SelectIntegrationFeatures), so batch-specific noise genes don’t dominate.

1.10 Scaling

The last step before PCA: z-score each gene across cells so it has mean 0 and SD 1. Why? Because PCA finds directions of maximum variance. Without scaling, a few highly-expressed genes (whose variance is large simply because their absolute expression is large) would dominate PC1 — and we want PC1 to track biology, not raw expression magnitude.

By default ScaleData() scales only the ~2,000 HVGs, which is all PCA needs and keeps the scale.data slot small. If you later want a non-HVG gene on a heatmap, re-run ScaleData(seu, features = <those genes>).

Optional: regress out unwanted covariates. ScaleData(seu, vars.to.regress = c("percent.mt", "nCount_RNA")) removes the linear effect of those covariates from each gene before z-scoring. Useful when a covariate is visibly driving early PCs; usually unnecessary for clean data. Note that regressing a batch/condition variable here (e.g. vars.to.regress = "stim") is a crude form of batch correction — it removes a single linear effect uniformly across all cell types, whereas real batch effects are non-linear and cell-type-specific. Use a proper integration method instead (Chapter 5).

A note: SCTransform does scaling implicitly — its Pearson residuals are already on a comparable scale across genes — so if you used SCTransform, you don’t run ScaleData.

1.11 What you have at the end

After steps 1–8, the Seurat object holds:

• @assays$RNA$counts — raw integer UMI counts (the input)
• @assays$RNA$data — log-normalized values
• @assays$RNA$scale.data — z-scored HVG matrix (the input to PCA)
• VariableFeatures(seu) — vector of ~2,000 HVG names
• @meta.data — per-cell metadata (sample, QC metrics, doublet calls)

Typical cell counts for the workshop ifnb dataset. Starting from muscData::Kang18_8vs8() and keeping only annotated singlets yields ~24,500 cells across both conditions (~12,000 CTRL + ~12,000 STIM). After the QC filter (nFeature_RNA > 200 & < 2500, percent.mt < 5, scDblFinder.class == "singlet"), the retained count is approximately ~23,000 cells. These post-QC cells are what every downstream module in the workshop uses.

In Seurat v5 the assay is layered — after merge() you may see counts.CTRL/counts.STIM layers that must be recombined with JoinLayers() before single-matrix tools (e.g. FindMarkers, FindVariableFeatures, AddModuleScore) will work; the tutorial does this explicitly after the per-sample doublet-detection merge. If you skip JoinLayers(), downstream calls will either error with “data layers are not joined” or silently skip every gene and return an empty table. The next chapter (Chapter 2 — Dimensionality Reduction & Clustering) takes the scaled HVG matrix and turns it into clusters and a UMAP.

1.12 Common errors and what they mean

Table 3 collects the errors you are most likely to hit at this stage, what each one means, and the fix.

Table 3: Common preprocessing errors, their underlying cause, and the fix.

You see	What’s wrong	What to do
Read10X returns a list, not a matrix	The Cell Ranger output has multiple feature types (RNA + ADT + CMO)	Pull out the slice you want: counts$``Gene Expression``
Filter removes 80% of cells	Your thresholds are wrong, or QC distributions are bad	Plot per-sample distributions; recompute thresholds from the data
percent.mt is 0 for every cell	The mitochondrial gene pattern doesn’t match	Check rownames(seu) for MT- (human) or mt- (mouse)
scDblFinder flags 30% of cells	Either you over-loaded the chip, or you ran it on a merged object	Run per sample; if rate is genuinely 30% then your run is poor quality
HVG list is dominated by ribosomal genes	Your data had no strong biology, or QC was insufficient	Re-check QC; consider excluding ^RP[SL] from HVGs; choose a tissue with cleaner heterogeneity
FindMarkers errors after merge()	The v5 assay still has split layers	Run JoinLayers(seu) first

For a longer list of preprocessing pitfalls, see the Common errors callout at the top of Tutorial 01.

1.13 Where this material is also discussed

• Lecture: Lecture 01 — Pre-processing
• Tutorial: Tutorial 01 — QC & Preprocessing
• Biology background: Appendix A — Single-Cell Biology Refresher
• Stats background: Appendix B — Statistical Foundations, §§ 1–3
• File formats and what’s in filtered_feature_bc_matrix/: Appendix E
• R basics for the code chunks: P2 — R & RStudio
• Methods-section template covering this chapter’s choices: Appendix G

1.14 Going further

The foundational best-practice tutorials remain the single best next read: the Bioconductor “Orchestrating Single-Cell Analysis” workflow⁸, the original Seurat/Scanpy best-practices papers^1,7, and the modern, modality-spanning update². For the specific tools in this chapter, go to the primary methods: EmptyDrops¹⁰, SoupX¹¹, DecontX¹² and CellBender¹³ for ambient RNA; scDblFinder¹⁵, DoubletFinder¹⁶ and Scrublet¹⁷ for doublets, with the comparative benchmark to choose between them¹⁸; and sctransform²⁰, scran pooling²¹ and the transformation benchmark²² for normalization. The full curated list is on Key Papers, Reviews & Benchmarks.

References

1.

Luecken, M. D. & Theis, F. J. Current best practices in single-cell RNA-seq analysis: A tutorial. Molecular Systems Biology 15, e8746 (2019).

2.

Heumos, L. et al. Best practices for single-cell analysis across modalities. Nature Reviews Genetics 24, 550–572 (2023).

3.

Kharchenko, P. V. The triumphs and limitations of computational methods for scRNA-seq. Nature Methods 18, 723–732 (2021).

4.

Lähnemann, D. et al. Eleven grand challenges in single-cell data science. Genome Biology 21, 31 (2020).

5.

Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial reconstruction of single-cell gene expression data. Nature Biotechnology 33, 495–502 (2015).

6.

Hao, Y. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nature Biotechnology 42, 293–304 (2024).

7.

Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: Large-scale single-cell gene expression data analysis. Genome Biology 19, 15 (2018).

8.

Amezquita, R. A. et al. Orchestrating single-cell analysis with Bioconductor. Nature Methods 17, 137–145 (2020).

9.

Ilicic, T. et al. Classification of low quality cells from single-cell RNA-seq data. Genome Biology 17, 29 (2016).

10.

Lun, A. T. L. et al. EmptyDrops: Distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data. Genome Biology 20, 63 (2019).

11.

Young, M. D. & Behjati, S. SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data. GigaScience 9, giaa151 (2020).

12.

Yang, S. et al. Decontamination of ambient RNA in single-cell RNA-seq with DecontX. Genome Biology 21, 57 (2020).

13.

Fleming, S. J. et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nature Methods 20, 1323–1335 (2023).

14.

Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nature Communications 8, 14049 (2017).

15.

Germain, P.-L., Lun, A., Garcia Meixide, C., Macnair, W. & Robinson, M. D. Doublet identification in single-cell sequencing data using scDblFinder. F1000Research 10, 979 (2021).

16.

McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: Doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Systems 8, 329–337 (2019).

17.

Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: Computational identification of cell doublets in single-cell transcriptomic data. Cell Systems 8, 281–291 (2019).

18.

Xi, N. M. & Li, J. J. Benchmarking computational doublet-detection methods for single-cell RNA sequencing data. Cell Systems 12, 176–194 (2021).

19.

Kang, H. M. et al. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology 36, 89–94 (2018).

20.

Hafemeister, C. & Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biology 20, 296 (2019).

21.

Lun, A. T. L., Bach, K. & Marioni, J. C. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biology 17, 75 (2016).

22.

Ahlmann-Eltze, C. & Huber, W. Comparison of transformations for single-cell RNA-seq data. Nature Methods 20, 665–672 (2023).