Multi-Sample Integration Tutorial
Merge multiple scRNA-seq samples into a single dataset, remove batch effects with Harmony, and produce an integrated UMAP where biology — not sequencing batch — drives clustering.
What you'll learn
- Concatenate multiple H5AD files with Merge Data, preserving each cell's sample identity in obs["sample"].
- Run batch-aware HVG selection (
batch_key="sample") so feature selection isn't dominated by batch artifacts. - Apply Harmony correction, route Neighbors to
use_rep="X_pca_harmony", and verify the integrated UMAP shows cells of both samples interleaved within each cluster.
Prerequisites
- • Oreumics installed (install guide).
- • Familiarity with the standard scRNA workflow — complete the PBMC tutorial first.
- • Harmony tutorial (2-sample case, follows the shipped scRNA Harmony Integration template) recommended — this tutorial generalizes that pattern to N samples.
- • Sample data:
sampledata/scrna/pbmc3k_raw.h5ad(loaded twice, labelled per input in Merge Data).
When do you need multi-sample integration?
- • Multiple donors / conditions: control vs. disease, timepoint A vs. B, patient cohort studies
- • Multiple library preps: data produced on different days, different 10x chemistry versions, or different labs
- • Different sequencing runs: samples sequenced across multiple flow cells or instruments
If your data comes from a single sample / single run, skip integration and use the standard PBMC tutorial.
Dataset
Sample files (pre-packaged)
sampledata/scrna/pbmc3k_raw.h5adsampledata/scrna/pbmc3k_raw.h5ad ← loaded a second time as "sample_B"What this tutorial demonstrates
Two independent H5AD files are loaded, merged into a single AnnData with a sample batch column, preprocessed, and then batch-corrected with Harmony. The tutorial shows how to verify that batches are mixed correctly in the final UMAP.
Load and Merge Samples
Load each sample independently, then use Merge Data to concatenate them into a single AnnData. Merge Data adds a sample column to obs so you can identify which cell came from which sample — this is the batch_key you will use later in Harmony. The label for each input is entered on the Merge Data node itself (one card per connected input, prefilled from the source node's name) — Load H5AD has no sample-name parameter.
Pipeline
Load the first sample. Rename the node to "Sample A" (double-click its title) — Merge Data uses node names to prefill the per-input labels.
file_pathsampledata/scrna/pbmc3k_raw.h5adWhat you'll see
AnnData: 2,700 cells × 32,738 genes.
Load the second sample the same way and name the node "Sample B". Both outputs connect to the single Merge Data node.
file_pathsampledata/scrna/pbmc3k_raw.h5adWhat you'll see
AnnData: 2,700 cells × 32,738 genes.
Concatenate both AnnData objects into a single dataset. batch_key names the obs column that receives the per-input labels; the labels themselves are typed into the input cards (batch_labels_list). Cell barcodes are de-duplicated automatically, and the live pre-merge preview shows the resulting cell and gene counts before you run.
batch_keysamplecolumn name in the merged obsbatch_labels_listsample_A, sample_Bone card per inputjoininnerkeep only genes present in ALL samples (default: outer)What you'll see
5,400 cells × 32,738 genes. obs["sample"] identifies each cell's origin. Use join="outer" if samples have different gene sets and you want to keep all genes (missing values become 0).
Tip
Prefer join="inner" unless samples were aligned to different annotations. "outer" padding with zeros can inflate sparsity and distort downstream statistics.
Quality Control
QC is applied to the merged dataset. Because both samples have already been concatenated, the thresholds apply globally. For real multi-sample experiments, inspect QC distributions per sample before filtering — set color_by='sample' on a label-mode Quality Control node to get per-sample violin plots.
Pipeline
Compute per-cell QC metrics across all merged cells without removing anything, and draw one violin per sample. Human mitochondrial prefix is MT- (uppercase); for mouse data use mt- (lowercase).
modelabelmito_prefixMT-ribo_prefixRPS,RPLplot_typeviolinmetrics["n_genes_by_counts", "total_counts", "pct_counts_mt"]color_bysampleone violin per sampleTip
If one sample has a much lower n_genes median (< 500), it may be a low-quality library that warrants separate investigation or exclusion before merging. The QC Workspace shows MAD-based threshold recommendations for the pooled data.
Filter cells and genes using thresholds derived from the per-sample violin plots. These values work for typical human PBMC data.
modefiltermin_genes200max_pct_mito20filter_min_cellstruemin_cells3Standard Preprocessing
Preprocessing is identical to a single-sample workflow. Run these nodes in strict order — normalize → log1p → HVG → scale → PCA. Harmony is inserted after PCA.
Pipeline
Normalize each cell to the same total count. save_counts: true stores the raw integer counts in adata.layers["counts"] — needed if you want seurat_v3 HVG selection or a counts-based test later.
target_sum10000save_countstrueLog-transform the normalized counts. This stabilizes variance and makes the data approximately normally distributed — a prerequisite for PCA.
basenullnatural logSelect 2,000 highly variable genes. For multi-sample data, set batch_key so HVGs are chosen within each sample and then combined — this stops genes that merely differ between batches from being selected. flavor "seurat" runs on the log-normalized matrix; "seurat_v3" is an alternative that expects raw counts, so pair it with layer: counts.
n_top_genes2000flavorseuratbatch_keysampleper-sample HVG selectionTip
Providing batch_key selects genes that are variable inside samples, not just because of batch differences. This is the recommended approach for multi-sample datasets.
Standardize each gene to zero mean and unit variance. With hvg_only: true (default) X keeps only the HVGs while the full log-normalized gene set is saved to adata.raw. After scaling adata.X contains z-scores — never use it directly for expression visualization (the plot nodes read adata.raw or layers["normalized"] via layer="auto").
max_value10zero_centertruehvg_onlytrueRun PCA on scaled HVGs. At this stage the PCA embedding (X_pca) still contains batch effects — Harmony will correct it in the next section.
n_comps50use_highly_variabletrueTip
Add a PCA Visualization node (plot_type: pca_scatter, color_by: sample) after PCA. If the first PCs separate samples instead of cell types, batch effects are strong — Harmony correction is essential.
Harmony Batch Correction
Harmony operates on the PCA embedding and outputs a corrected version in obsm['X_pca_harmony']. The gene expression matrix (adata.X) is never modified — only the PCA coordinates are corrected. You must then pass use_rep='X_pca_harmony' to the Neighbors node, otherwise the batch correction has no effect downstream.
Critical: Harmony output key
Harmony outputs obsm['X_pca_harmony'] — it does NOT overwrite obsm['X_pca']. The most common mistake is forgetting to set use_rep='X_pca_harmony' in the Neighbors node. Without this, UMAP and clustering will use the uncorrected PCA and batch effects will remain.
Pipeline
Run Harmony integration on the PCA embedding. batch_key must match the obs column you assigned during Merge Data (e.g., "sample"; the node default is "batch"). Outputs obsm["X_pca_harmony"] — the batch-corrected PCA embedding.
batch_keysamplemust match obs column from Merge Datatheta2.0higher = more aggressive mixingmax_iter_harmony10What you'll see
obsm["X_pca_harmony"] added. This embedding should show cells from both samples interleaved, while uncorrected X_pca still separates them.
Build the kNN graph using the Harmony-corrected embedding. This is the most critical step — use_rep MUST be set to X_pca_harmony, not X_pca.
n_neighbors15use_repX_pca_harmony← CRITICAL: use Harmony outputCompute UMAP from the Harmony-corrected neighbor graph. The resulting embedding reflects biology (cell types), not batch of origin.
min_dist0.3spread1.0Run Leiden clustering on the Harmony-corrected neighbor graph.
methodleidenresolution0.8UMAP colored by cluster. Add a second Embedding Plot with color_by: sample to verify that batches are well-mixed within clusters — this is the integration quality check. (Cluster Visualization has no "umap" plot_type; it draws composition/proportions bars — its plot_type: proportions with groupby: sample gives a per-sample cluster breakdown.)
plot_typeumapcolor_byleidenthen a second node with "sample"What you'll see
A well-integrated UMAP: each cluster contains cells from both sample_A and sample_B. If one cluster is dominated by a single sample, try increasing theta in Harmony (up to 4) for stronger correction.
Tip
Good integration = clusters are not sample-specific. Because the demo loads the same file twice, mixing is trivially perfect — with real donors you should still expect some sample-specific clusters where biology differs.
Downstream Analysis
After integration, downstream analysis proceeds identically to a single-sample workflow. DEG Analysis identifies cluster-specific genes. The sample column can be used for condition-aware analyses (e.g., comparing cluster composition between groups).
Pipeline
Rank differentially expressed genes per cluster (formerly "Find Markers"). Works the same as in the single-sample workflow; leave n_genes empty to rank all genes. The Wilcoxon test is non-parametric and works well for integrated data.
groupbyleidenmethodwilcoxonn_genesnullall genesDotplot of the top N marker genes per cluster. Use this to annotate cell types and then use Edit Clusters to rename clusters.
plot_typerank_genes_groups_dotplotn_genes5Save the fully processed, integrated AnnData. The saved file contains the Harmony embedding, UMAP coordinates, cluster labels, and DEG results — ready for downstream comparative analyses.
output_pathresults/integrated_pbmc.h5adcompressiongzipIntegration Quality Checklist
Run Embedding Plot with color_by='sample' and check:
- ✓ Batch mixing: Each UMAP cluster contains cells from multiple samples — no cluster is dominated by a single batch
- ✓ Cell type preservation: Known cell type markers still show expected patterns (CD3D for T cells, LYZ for monocytes)
- ✓ Cluster coherence: Clusters visible in the Harmony-corrected UMAP are biologically interpretable, not random mixtures
If integration quality is poor:
- • Increase
thetain Harmony (2 → 3 → 4) for more aggressive mixing - • Increase
max_iter_harmonyto 20 if convergence is slow - • Re-run Highly Variable Genes with
batch_key='sample'(and, if you saved counts,flavor='seurat_v3'+layer='counts') for better feature selection - • Check if samples have very different total cell counts — extreme size imbalance can reduce integration quality
Common Pitfalls
Forgetting use_rep in Neighbors
The #1 mistake: running Harmony then leaving Neighbors at the default use_rep="X_pca". The batch correction only takes effect if Neighbors uses "X_pca_harmony".
Correcting when you should not
If "batch" differences actually reflect biology (e.g., comparing disease vs. control), over-correcting with high theta removes the signal you want to detect. Use lower theta or compare with uncorrected results.
Merging before per-sample QC
For real experiments, run a label-mode Quality Control node on each sample separately before merging. A single failed library can contaminate the merged dataset.
Using join="outer" without checking gene overlap
Outer join pads missing genes with zeros. For samples aligned to different genomes (e.g., hg38 vs mm10 in xenograft experiments), use inner join or pre-filter to the shared gene set.
Going Further
- Use the Harmony + Doublet tutorial for an explicit 2-batch walk-through with a doublet-screening front-end.
- Run cluster-level differential expression between samples (condition, treatment) using DEG Analysis with groupby="sample", or restrict the test to chosen clusters with its
groupsparameter. - Map cell-type labels from a reference dataset onto your integrated AnnData with the Cell Mapper node (the Reference Label Transfer template), or label cells directly with CellTypist Annotation — both bypass manual cluster-by-marker annotation.
Glossary
- batch_key
- An obs column that labels each cell's batch / sample / donor of origin. Set by Merge Data; consumed by Highly Variable Genes, Harmony, and condition-aware DEG.
- join (inner vs. outer)
- Merge Data's gene-set behavior. inner = keep only genes present in all samples (recommended). outer = keep all genes, padding missing values with 0 (the node default).
- layers["counts"]
- The raw integer counts that Normalize keeps aside when save_counts is on. seurat_v3 HVG selection and counts-based tests read from here; adata.raw (set by Scale) holds the log-normalized full gene set instead.
- obsm["X_pca_harmony"]
- Harmony's output — the batch-corrected PCA embedding. Pass to Neighbors via use_rep, otherwise correction is silently lost.
- Integration quality
- A subjective check — render the UMAP twice, once colored by leiden (clusters look biological?) and once by sample (cells from each batch interleaved within each cluster?).