Chapter 5 — Multi-Sample Integration

Author

Single Cell RNA-seq Workshop

Note

Where this chapter sits. Companion to Lecture 05 and Tutorial 05. Builds on Chapter 2 (PCA + clustering) and Chapter 3 (manual annotation, used here for diagnostic checks).

5.1 The problem integration solves

Suppose you have two samples — say, peripheral blood from a control donor and from a patient — and you want a unified clustering and annotation. You run the Chapter 1 → Chapter 2 workflow on the joint object. You compute a UMAP. You color it by sample.

If the colors interleave within every cluster, you don’t need integration: the technical and biological structures are aligned, and downstream analyses can proceed.

If the samples form parallel structures — every cell type appearing twice, once per sample — you have a batch effect. PC1 of the joint dataset is now tracking which sample the cell came from rather than what cell type it is. Clusters at the standard resolution will be sample-specific rather than cell-type-specific. Pseudobulk DE between samples becomes non-comparable because you have no shared definition of “the same cell type” across them.

Integration is the family of methods that warp the embedding so the same cell type from different samples ends up in the same cluster, while preserving any real biological differences between samples (e.g. an activated state present in one sample but not another).

Note🔑 Key concept — what integration actually aligns

Integration does not align raw expression values; it aligns the geometric positions of cells in PC space. A batch effect in scRNA-seq is not a simple additive shift (as in bulk RNA-seq) but a complex non-linear warp: the same cell type measured in two libraries can sit in two distinct PC-space clouds, with different variances, different inter-cluster distances, and often different numbers of recovered clusters. Integration methods identify which cells from different batches are most likely the same type — either by finding mutual nearest-neighbor pairs across batches (Seurat anchors, MNN) or by penalizing cluster centroids that are batch-unbalanced (Harmony) — and apply a correction that moves those cells toward each other. The correction is applied in PC space (Harmony) or in a corrected expression assay (Seurat anchors), but in both cases the downstream biological signal — stored in the original RNA assay — is untouched. This is why you always switch back to DefaultAssay(seu) <- "RNA" before any differential expression or module scoring.

5.2 The diagnostic: when do you need integration?

Run the unintegrated PCA + UMAP from Chapter 2 and inspect three things:

  1. UMAP colored by sample. Do colors interleave or form islands?
  2. UMAP colored by author cell-type label (if available — the workshop’s ifnb ships them in seurat_annotations). Does each cell type form one region or multiple (one per sample)?
  3. table(seu$cluster, seu$sample) — what fraction of each cluster is each sample?

If samples are interleaved within most clusters and a balanced fraction shows up everywhere, you do not need integration. If samples form per-cluster islands and many clusters are >80% one sample, you do.

The workshop’s ifnb shows the second pattern dramatically: PC1 is essentially the IFN-β response axis. CTRL and STIM cells split every cell type into two sister clusters. This is the canonical “integration is needed” case.

5.3 The conceptual taxonomy of integration methods

The methods that have been developed for batch correction span a wide range of conceptual approaches, from simple linear regression to deep generative models. Table 1 places the five main families in context; the two the workshop implements in depth are Harmony and Seurat anchors.

Table 1: The five conceptual families of single-cell batch-correction methods. There is no universally best method: large benchmarks find that the right choice depends on dataset complexity and the trade-off between batch removal and biological-variation conservation1,2.
Family Examples Mechanism When to use
Linear / regression ComBat, sc.pp.regress_out Fit batch as a covariate; subtract its linear effect Simple cases; rarely sufficient on its own for scRNA-seq
Mutual nearest neighbors MNN3, fastMNN, Seurat anchors4 Find pairs of cells that are nearest neighbors across batches; learn a correction from those anchors Two or a few batches with similar cell-type compositions
Soft clustering in PC space Harmony5 Iteratively re-cluster cells with a batch-mixing penalty; shift centroids toward batch balance Most-used in 2026 — fast, scales to atlas sizes, works well for 2–50 batches
Variational autoencoders scVI6, scANVI7 Learn a latent representation conditioned on batch; the decoder can be ignored, the latent used like PCs Atlas-scale; especially with strong batch heterogeneity
Graph-based BBKNN, Conos Modify the kNN graph to balance cross-batch edges Fast; works well on already-PCA’d data

The workshop’s Tutorial 05 walks Harmony and Seurat anchors because they are the two dominant in-R methods. scVI requires Python; BBKNN is mainly used in Scanpy. Harmony, Seurat, scVI and scANVI all rank near the top under different conditions in the benchmarks, with no single winner across all datasets1,2.

5.4 Harmony in detail

Harmony5 operates in PC space — after you’ve run RunPCA on the joint object. It iterates four steps until convergence:

  1. Soft-cluster cells with a batch-diversity penalty: each cell is assigned soft membership in \(K\) clusters; the objective rewards clusters whose membership is balanced across batches
  2. Compute centroids: a global centroid per cluster, plus per-batch centroids per cluster
  3. Compute correction factors: the difference between each per-batch centroid and the global centroid
  4. Shift each cell by a soft-cluster-weighted combination of those correction factors

The output is a harmony reduction that is a drop-in replacement for pca in FindNeighbors and RunUMAP. The original PCs are preserved separately. The soft-clustering approach is what distinguishes Harmony from hard-assignment methods: because each cell belongs to multiple clusters with fractional membership, the correction is smooth and continuous — a cell on the boundary between two populations is shifted by a weighted average of both populations’ correction factors, rather than being abruptly assigned to one population’s correction. This continuity makes Harmony robust to rare populations that straddle cluster boundaries in the pre-correction PC space.

Harmony’s strengths:

  • Fast — minutes for atlas-scale (~M cells) data
  • Conservative — preserves rare populations well
  • Single-tunable: theta (default 2) controls how aggressively to mix batches; rarely needs adjustment

Pre-processing for Harmony: before running Harmony, the typical workflow is: (1) split by sample, (2) normalize and select 2,000 HVGs per sample, (3) identify features that are variable across most samples with SelectIntegrationFeatures(object.list = ...) — this avoids letting one dominant sample dictate the feature set, (4) merge the per-sample objects back with merge() and call JoinLayers() (Seurat v5 merge() leaves per-sample split layers; JoinLayers() collapses them so AddModuleScore and FindAllMarkers work on the RNA assay later), (5) run PCA on the integration features, then (6) call RunHarmony().

Harmony’s weaknesses:

  • Less effective on very heterogeneous batches (cell types absent in one but present in another); works best when most cell types are shared
Note⚙️ Key parameter — RunHarmony theta

RunHarmony(seu, group.by.vars = "stim", theta = 2). Default: theta = 2. Controls how aggressively Harmony penalizes batch imbalance in its soft-clustering objective. Higher theta (e.g. 4–6) forces more aggressive mixing across batches — useful for strong, hard-to-correct batch effects but risks erasing real biological variation between conditions. Lower theta (e.g. 1) is less aggressive and is appropriate when batches are nearly comparable already or when biological differences between samples are strong and should not be over-corrected. For most multi-sample scRNA-seq studies the default of 2 is a sensible starting point; the ISG-score diagnostic in §5.6 is the practical test of whether you’ve pushed too hard.

Note⚙️ Key parameter — SelectIntegrationFeatures nfeatures

SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 2000). Default: 2000. Selects genes that are among the top HVGs in the most samples, producing a feature set that is robust across batches rather than dominated by any single sample. Increasing nfeatures (e.g. 3000–5000) includes more genes and can recover rarer cell-type markers, at the cost of including more batch-specific noise genes. For most two-sample designs (like ifnb), 2000 is sufficient; for highly heterogeneous multi-tissue or multi-protocol atlases, raising to 3000–4000 is common practice. These features are then used to scale the merged object and run PCA before Harmony or IntegrateLayers.

5.5 Seurat anchors in detail

The Seurat anchor workflow4,8 works in two stages:

5.5.1 Find anchors

For each pair of samples:

  1. Project both into a shared low-rank space using canonical correlation analysis (CCA) + L2 normalization
  2. For each cell, find its mutual nearest neighbors across samples — pairs \((i_A, j_B)\) where \(i_A\) is among \(j_B\)’s nearest neighbors and vice versa
  3. Score each anchor by how consistent its local neighborhood is — anchors whose neighborhoods are mutually inhabited are scored high; isolated anchors (likely doublet-driven or technical) are scored low

The result is a list of anchored cell pairs and their confidence scores.

5.5.2 Apply the correction

Each cell in sample \(B\) is moved toward its anchored partner in sample \(A\), weighted by anchor score. The amount of movement is regularized so cells far from any anchor barely move.

The corrected expression lives in seu@assays$integrated@data. PCA on this corrected assay produces the integrated reduction used downstream.

Seurat anchors’ strengths:

  • Most rigorous of the alignment methods — anchor scoring filters bad correspondences
  • Preserves rare populations more conservatively than Harmony
  • Built-in to Seurat, no extra packages

Seurat anchors’ weaknesses:

  • Slow — for atlas data, easily an hour+; Harmony is 10× faster
  • Complex — multiple steps each with parameters; harder to tune

5.5b Harmony vs. Seurat anchors — a direct comparison

Because the workshop implements both methods, it is worth pausing to contrast them explicitly before working through their shared failure modes. Table 2 places them side by side on the dimensions that most affect the analysis decision.

Table 2: Direct comparison of Harmony and Seurat anchors on the dimensions most relevant to choosing between them. For most PBMC-scale studies (10k–200k cells), both work well; Harmony is the practical default for speed; Seurat anchors can provide more conservative correction when cell-type compositions are very unequal across batches.
Dimension Harmony Seurat anchors
Speed Fast — minutes for millions of cells Slow — hours for large datasets
Operating space PC space (corrects coordinates) Gene expression space (corrects an assay)
Output slot reductions$harmony (drop-in for pca) assays$integrated (separate assay)
Memory footprint Low High (stores corrected expression matrix)
Rare cell conservation Good Good to excellent (anchor scoring filters bad pairs)
Key tuning parameter theta (default 2; higher = more aggressive mixing) k.filter, k.anchor, dims
Downstream DE Always use RNA assay Always use RNA assay — integrated assay is not valid for DE
Atlas scale (>500k cells) Yes Not recommended
Note🔑 Key concept — the over-correction vs. under-correction trade-off

Integration methods operate on a fundamental tension: every algorithm that moves cells from different samples together risks also moving cells with different biology together. A Harmony run with theta = 10 will produce a perfectly sample-mixed UMAP, but at the cost of flattening the IFN-β response, the tumor microenvironment gradient, or any other biology that correlates with the sample label. The core diagnostic is therefore not “is the UMAP sample-mixed?” but “are real biological signals preserved within the mixed clusters?” The ISG-score gap test (§5.6) is the canonical check for the ifnb dataset: if STIM cells still score higher on interferon-stimulated genes than CTRL cells of the same cluster type, the biology is intact. If the gap is zero, you’ve over-corrected. No integration metric should ever be interpreted in isolation from a biology-preservation check.

5.6 The over-correction problem

The cardinal failure mode of integration is erasing real biological differences by being too aggressive. If your two samples are biologically different (CTRL vs IFN-β–stimulated), and integration spatially mixes them too completely, you also erase the IFN response — which is the actual biology of the experiment.

The right outcome: cells from both samples interleave on the UMAP but their gene-level expression still differs in the way the perturbation predicts. CTRL and STIM cells of the same cluster sit in the same UMAP region but the STIM cells still express more ISG15, IFI6, MX1, etc.

This is testable. The workshop’s diagnostic is the ISG score gap:

  1. Score each cell on a panel of interferon-stimulated genes using AddModuleScore (run on the original RNA assay, not the integrated assay — the integrated assay’s expression values aren’t valid for DE)
  2. For each cell type, compute mean ISG score in CTRL vs STIM
  3. Confirm STIM > CTRL within every cell type

If the gap is gone, integration over-corrected. Push less hard (lower Harmony theta, or use Seurat anchors with default settings). If the gap is preserved, integration was successful.

5.7 The under-correction problem

The opposite failure: clusters remain sample-specific after integration. There are two interpretations and they require different responses:

Two interpretations are possible and they require different responses. The first is genuinely sample-specific biology: a cell state that exists in one sample but not the other — for example, an interferon-activated T-cell subpopulation present only in the stimulated condition — should remain sample-specific after integration, because that’s real biological information, not a technical artifact. Integration that erases it is over-correcting, not under-correcting. The second is a residual batch effect: a cell type that is clearly present in both samples but still splits into two clusters after integration. Here the method has not fully aligned the shared population, and pushing integration harder (higher Harmony theta, or more stringent anchor filtering in Seurat) is warranted.

Distinguishing them requires biological knowledge and careful inspection of cluster markers. If the sample-specific cluster’s markers match an activation state (cytokines, ISGs) rather than a lineage (CD4, CD79A), it’s likely real biology. If the markers are the same lineage markers as a sample-mixed cluster, it’s under-correction.

5.8 The three diagnostics, in order

For every integrated dataset, run all three:

  1. UMAP by sample — colors should interleave within each cluster (after integration). If a cluster is >80% one sample, suspect either over- or under-correction
  2. UMAP by biology (cell-type labels, or a known signature score) — each cell type / state should occupy one region (not multiple, as in unintegrated data)
  3. Known biological signal preserved — score a known biology-associated module per cell, confirm it still differs between conditions within clusters as the experiment predicts

Methods sections of integration papers virtually always report (1). Many skip (2). Most skip (3). All three are equally important.

5.9 Pseudoreplicate caveat

Some users of integration methods are tempted to treat the integrated assay as a fully-corrected expression matrix and run downstream analyses (DE, marker tests) on it. This is wrong.

The integrated assay holds expression values that are valid for geometry — PCA, neighbor graph, UMAP, clustering. They are not valid for differential expression because the correction can dampen, distort, or invert real biological variation in a way that biases per-gene tests.

The discipline:

  • Use DefaultAssay(seu) <- "integrated" before RunPCA, FindNeighbors, FindClusters, RunUMAP
  • Switch back to DefaultAssay(seu) <- "RNA" after clustering, before FindMarkers, FindAllMarkers, marker scoring, or any DE work

This is one of the most-violated rules in scRNA-seq analysis. Always explicit.

Over-correction false positive. A subtle trap specific to the over-correction diagnostic: if you score cells on a biological signal (e.g. an ISG module score) after forgetting to switch back to the RNA assay, you will score on the integrated (batch-corrected) expression values. Since integration explicitly reduces the STIM–CTRL difference in the corrected space, the ISG score gap will appear to be zero or near-zero — falsely suggesting over-correction. The fix is DefaultAssay(seu) <- "RNA" immediately before AddModuleScore. The same applies to any FeaturePlot used to visualize a gene’s expression post-integration: always verify the active assay before drawing.

5.10 Common errors

Table 3 collects the errors most frequently encountered during integration and the integrated-object diagnostics. The single most important habit to develop is checking DefaultAssay before any expression-level operation.

Table 3: Common errors in multi-sample integration, their root causes, and fixes. The DefaultAssay confusion (integrated vs. RNA) is the single most common bug; always be explicit about which assay is active before each analysis step.
You see What’s wrong What to do
Harmony errors “no pca reduction” You skipped RunPCA Run PCA first; Harmony works on PCs
ISG score gap is zero after integration You scored on the integrated assay Set DefaultAssay(seu) <- "RNA" before AddModuleScore
Clusters after integration look identical to unintegrated FindNeighbors is using pca not harmony FindNeighbors(seu, reduction = "harmony", dims = 1:30)
IntegrateData returns a 5-dim object that won’t slot into RunPCA Seurat v5 / v4 API mismatch In v5, anchor integration produces a layered Assay5; consult the Seurat v5 vignette
FindAllMarkers returns empty results after integration DefaultAssay is still "integrated" DefaultAssay(seu) <- "RNA" before marker tests
JoinLayers errors on a non-merged object Object was never merge()-d JoinLayers is only needed after merge() creates split layers

5.11 Where this material is also discussed

5.12 Going further

The two method papers the tutorial uses are Harmony5 and Seurat anchors4,8; the deep-learning alternatives are scVI6 and its semi-supervised extension scANVI7, which also transfers labels. To choose among them, read the two benchmarks — of 14 batch-correction methods1 and of atlas-level integration2, the latter introducing the now-standard scoring that separates batch removal from biological conservation (metrics such as LISI and kBET). Mutual-nearest-neighbour correction3 is the conceptual ancestor of the anchor approach. The full curated list is on Key Papers, Reviews & Benchmarks.

References

1.
Tran, H. T. N. et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biology 21, 12 (2020).
2.
Luecken, M. D. et al. Benchmarking atlas-level data integration in single-cell genomics. Nature Methods 19, 41–50 (2022).
3.
Haghverdi, L., Lun, A. T. L., Morgan, M. D. & Marioni, J. C. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nature Biotechnology 36, 421–427 (2018).
4.
Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177, 1888–1902 (2019).
5.
Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nature Methods 16, 1289–1296 (2019).
6.
Lopez, R., Regier, J., Cole, M. B., Jordan, M. I. & Yosef, N. Deep generative modeling for single-cell transcriptomics. Nature Methods 15, 1053–1058 (2018).
7.
Xu, C. et al. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models. Molecular Systems Biology 17, e9620 (2021).
8.
Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587 (2021).
Back to top