Version Compatibility
Reference examples tested with: Cell Ranger ARC 2.2+, Seurat 5.1+, Signac 1.14+, EnsDb.Hsapiens.v86, BSgenome.Hsapiens.UCSC.hg38 1.4+, ggplot2 3.5+ (AMULET via the standalone java/python tool or scDblFinder's amulet() in R -- ArchR and snapATAC2 ship their OWN simulation-based doublet callers, not AMULET; MultiVI via scvi-tools if using the Python path)
Before using code patterns, verify installed versions match. If versions differ:
- R:
packageVersion('<pkg>') then ?function_name to verify parameters
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Note: cellranger-arc count (NOT cellranger-atac) emits the paired RNA+ATAC per-nucleus barcodes. Seurat FindClusters(algorithm=3) is SLM, not Leiden (1=Louvain, 2=Louvain-multilevel, 3=SLM, 4=Leiden). The atac_fragments.tsv.gz must be block-gzipped + tabix-indexed; the Tn5 +4/-5 offset is already applied by 10x -- do not re-shift. Confirm in-tool before quoting.
Multiome Pipeline
"Analyze my 10X Multiome data jointly" -> Orchestrate Cell Ranger ARC processing, Seurat/Signac scRNA+scATAC integration via WNN, chromatin accessibility peak calling, motif enrichment, and gene regulatory network inference.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
Made-once commitments
| Commitment | Consequence inherited downstream |
|---|
cellranger-arc run (NOT cellranger-atac) | Only ARC emits the joined RNA+ATAC per-nucleus barcodes; -atac gives ATAC-only barcodes and the join is impossible |
| Shared cell-barcode join | RNA and ATAC QC pass DIFFERENT barcodes; the analyzable set is their INTERSECTION. A namespace mismatch (-1 suffix, RNA vs ATAC whitelist) silently empties the join |
| Same genome build for GEX + ATAC (EnsDb + BSgenome) | Gene activity, LinkPeaks, and motif coordinates require identical build, else peak-to-gene linking is garbage |
| Consensus peak set (multi-sample) | Peaks are dataset-specific; merging samples on discordant peaks fabricates batch structure -- re-quantify against a unified peak set |
| Intron inclusion (GEX half) | Multiome is nuclei (mostly unspliced) -- introns are essential for the RNA UMI totals |
Pipeline orchestration: the joint-modality decisions that make or break the result
Multiome's defining feature is that RNA and ATAC are measured in the SAME nucleus, so the two assays share one barcode universe and must be reconciled, not analyzed independently. The orchestration decisions:
- The paired-cell anchor is the whole point: keep only barcodes that pass QC in BOTH modalities. RNA and ATAC QC use different metrics (RNA: gene count, mito %; ATAC: TSS enrichment, nucleosome signal, fragment count) and are computed per modality, but the surviving cell set is their intersection. cellranger-arc (not cellranger-atac) produces the paired barcodes; their universes differ. See single-cell/preprocessing and single-cell/scatac-analysis.
- Per-modality QC and doublet detection run BEFORE the joint embedding. ATAC doublets are missed by RNA-based callers (scDblFinder/Scrublet, see single-cell/doublet-detection) and need a fragment-based detector such as AMULET (see single-cell/scatac-analysis); resolving them after WNN lets fake intermediate states drive the joint clustering.
- Drop the depth-correlated LSI component before joint analysis. The first ATAC LSI/SVD component usually (not always) captures sequencing depth rather than biology; check with DepthCor and exclude whichever component correlates with depth (commonly #1, hence dims 2:30 in WNN). Forgetting this lets depth dominate the ATAC contribution to the joint graph.
- WNN vs a generative joint embedding is a real choice. Seurat/Signac WNN learns a per-cell modality weight on precomputed PCA + LSI and is the default when both modalities are well processed; MultiVI (scvi-tools) jointly models RNA + ATAC counts end-to-end and handles batch and mosaic (RNA-only or ATAC-only) cells better. WNN cannot integrate cells missing a modality. See single-cell/multimodal-integration; verify current best practice against installed docs.
- Embed (WNN) then cluster then annotate, and annotate from RNA primarily. Gene-activity scores derived from ATAC are an approximation of expression, so cell-type labels come from RNA markers; ATAC informs the regulatory state, not the identity call. See single-cell/clustering, single-cell/markers-annotation, single-cell/cell-annotation.
- Peak-to-gene linking is correlational, not causal. LinkPeaks correlates peak accessibility with gene expression across cells within a window; a link is a hypothesis to validate, not a proven enhancer-target pair. For genome-wide enhancer-gene mapping use ABC/ENCODE-rE2G. See atac-seq/co-accessibility and atac-seq/enhancer-gene-linking.
- Cross-condition questions still need pseudobulk and separate composition testing. Condition DE on either modality aggregates RAW counts per sample x cell type (cells-as-replicates is pseudoreplication, Squair 2021); proportion shifts between conditions are tested separately and can masquerade as DE. See single-cell/differential-abundance and differential-expression/deseq2-basics.
Workflow Overview
10X Multiome data
|
v
[1. Load Data] ---------> Read RNA + ATAC
|
v
[2. RNA Processing] ----> Standard scRNA workflow
|
v
[3. ATAC Processing] ---> Peak calling, LSI
|
v
[4. WNN Integration] ---> Weighted nearest neighbors
|
v
[5. Joint Analysis] ----> Clustering, markers
|
v
[6. Linked Features] ---> Gene-peak links
|
v
Integrated multiome object
Step 1: Load Multiome Data
library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
rna_counts <- Read10X_h5('filtered_feature_bc_matrix.h5')
seurat_obj <- CreateSeuratObject(
counts = rna_counts$`Gene Expression`,
assay = 'RNA'
)
atac_counts <- rna_counts$Peaks
frags <- CreateFragmentObject('atac_fragments.tsv.gz', cells = colnames(seurat_obj))
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- 'UCSC'
atac_assay <- CreateChromatinAssay
counts atac_counts
sep
fragments frags
annotation annotations
seurat_obj atac_assay
Step 2: RNA Quality Control and Processing
seurat_obj[['percent.mt']] <- PercentageFeatureSet(seurat_obj, pattern = '^MT-')
seurat_obj <- subset(seurat_obj,
nCount_RNA > 1000 &
nCount_RNA < 25000 &
percent.mt < 20
)
seurat_obj <- SCTransform(seurat_obj, assay = 'RNA', verbose = FALSE)
seurat_obj <- RunPCA(seurat_obj, assay = 'SCT', verbose = FALSE)
Step 3: ATAC Quality Control and Processing
DefaultAssay(seurat_obj) <- 'ATAC'
seurat_obj <- NucleosomeSignal(seurat_obj)
seurat_obj <- TSSEnrichment(seurat_obj)
VlnPlot(seurat_obj, features = c('nCount_ATAC', 'TSS.enrichment', 'nucleosome_signal'),
pt.size = 0, ncol = 3)
seurat_obj <- subset(seurat_obj,
nCount_ATAC > 1000 &
nCount_ATAC < 100000 &
TSS.enrichment > 2 &
nucleosome_signal < 4
)
seurat_obj <- RunTFIDF(seurat_obj)
seurat_obj <- FindTopFeatures(seurat_obj, min.cutoff
seurat_obj RunSVDseurat_obj
DepthCorseurat_obj
Step 3b: Doublet detection (per modality, BEFORE WNN)
Remove doublets before the joint embedding, or fake intermediate states drive the joint clustering. RNA-based callers (scDblFinder/Scrublet) MISS ATAC doublets -- ATAC needs a fragment-based caller (AMULET), run on the same nuclei. Detect per modality, drop the union of doublets, then build WNN. Mechanism: single-cell/doublet-detection (RNA) and single-cell/scatac-analysis (AMULET).
Step 4: Weighted Nearest Neighbors (WNN)
seurat_obj <- FindMultiModalNeighbors(
seurat_obj,
reduction.list = list('pca', 'lsi'),
dims.list = list(1:30, 2:30),
modality.weight.name = 'RNA.weight'
)
seurat_obj <- RunUMAP(seurat_obj, nn.name = 'weighted.nn',
reduction.name = 'wnn.umap', reduction.key = 'wnnUMAP_')
seurat_obj <- FindClusters(seurat_obj, graph.name = 'wsnn',
algorithm = 3, resolution = 0.5, verbose
Step 5: Visualization and Markers
p1 <- DimPlot(seurat_obj, reduction = 'pca', label = TRUE) + ggtitle('RNA PCA')
p2 <- DimPlot(seurat_obj, reduction = 'lsi', label = TRUE) + ggtitle('ATAC LSI')
p3 <- DimPlot(seurat_obj, reduction = 'wnn.umap', label = TRUE) + ggtitle('WNN UMAP')
p1 + p2 + p3
VlnPlot(seurat_obj, features = 'RNA.weight', group.by = 'seurat_clusters', pt.size = 0)
DefaultAssayseurat_obj
rna_markers FindAllMarkersseurat_obj only.pos min.pct
DefaultAssayseurat_obj
atac_markers FindAllMarkersseurat_obj only.pos min.pct
test.use latent.vars
Step 6: Gene-Peak Linkage
DefaultAssay(seurat_obj) <- 'ATAC'
seurat_obj <- RegionStats(seurat_obj, genome = BSgenome.Hsapiens.UCSC.hg38)
seurat_obj <- LinkPeaks(
seurat_obj,
peak.assay = 'ATAC',
expression.assay = 'SCT',
genes.use = c('CD8A', 'CD4', 'MS4A1', 'CD14')
)
CoveragePlot(seurat_obj, region = 'CD8A', features = 'CD8A',
expression.assay = 'SCT', extend.upstream = 10000, extend.downstream = 10000)
Complete Workflow Script
library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(BSgenome.Hsapiens.UCSC.hg38)
library(ggplot2)
data_dir <- 'multiome_output'
output_dir <- 'multiome_results'
dir.create(output_dir, showWarnings = FALSE)
cat('Loading data...\n')
counts <- Read10X_h5(file.path(data_dir, 'filtered_feature_bc_matrix.h5'))
frags <- file.path(data_dir, 'atac_fragments.tsv.gz')
seurat_obj <- CreateSeuratObject(counts = counts$`Gene Expression`, assay = 'RNA')
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <-
seurat_obj CreateChromatinAssay
counts countsPeaks
sep
fragments frags
annotation annotations
cat ncolseurat_obj
cat
seurat_obj PercentageFeatureSetseurat_obj pattern
seurat_obj subsetseurat_obj nCount_RNA nCount_RNA percent.mt
cat
DefaultAssayseurat_obj
seurat_obj NucleosomeSignalseurat_obj
seurat_obj TSSEnrichmentseurat_obj
seurat_obj subsetseurat_obj nCount_ATAC TSS.enrichment nucleosome_signal
cat ncolseurat_obj
cat
DefaultAssayseurat_obj
seurat_obj SCTransformseurat_obj verbose
seurat_obj RunPCAseurat_obj verbose
cat
DefaultAssayseurat_obj
seurat_obj RunTFIDFseurat_obj
seurat_obj FindTopFeaturesseurat_obj min.cutoff
seurat_obj RunSVDseurat_obj
cat
seurat_obj FindMultiModalNeighborsseurat_obj
reduction.list
dims.list
modality.weight.name
seurat_obj RunUMAPseurat_obj nn.name
reduction.name reduction.key
seurat_obj FindClustersseurat_obj graph.name algorithm resolution verbose
saveRDSseurat_obj file.pathoutput_dir
pdffile.pathoutput_dir width height
DimPlotseurat_obj reduction label
dev.off
cat output_dir
cat uniqueseurat_objseurat_clusters
Common Errors
| Symptom | Cause | Fix |
|---|
| Near-empty joint object / no cells after join | RNA vs ATAC barcode namespace mismatch (-1 suffix, different whitelist) | Reconcile barcodes; intersect on identical strings; confirm cellranger-ARC (not -atac) |
| ATAC depth dominates the joint graph | Kept the depth-correlated LSI component | DepthCor -> drop it (WNN dims.list 2:30 for ATAC) |
| Fake intermediate joint clusters | ATAC doublets not removed (RNA caller is blind to them) | AMULET fragment-based doublet call per modality before WNN |
| Cell types mislabeled | Annotated from ATAC gene-activity | Annotate identity from RNA markers; activity is a cluster-level proxy |
| Spurious batch across multi-sample multiome | Merged on discordant peak sets | Unify peaks and re-quantify |
| "Enhancer regulates gene" overclaim | Read LinkPeaks correlation as causal | Treat as a composition-confounded hypothesis; validate |
| Inflated cross-condition DE | Tested cells as replicates on either modality | Pseudobulk RAW per sample x cell-type (Squair 2021) |
References
- Hao Y, Hao S, Andersen-Nissen E, et al (2021) Integrated analysis of multimodal single-cell data. Cell 184:3573-3587.e29. DOI 10.1016/j.cell.2021.04.048. (WNN.)
- Ashuach T, Gabitto MI, Koodli RV, et al (2023) MultiVI: deep generative model for the integration of multimodal data. Nature Methods 20:1222-1231. DOI 10.1038/s41592-023-01909-9. (mosaic-capable joint RNA+ATAC alternative to WNN.)
- Squair JW, Gautier M, Kathe C, et al (2021) Confronting false discoveries in single-cell differential expression. Nature Communications 12:5692. DOI 10.1038/s41467-021-25960-2. (pseudobulk for cross-condition DE.)
Related Skills
- single-cell/data-io - Loading 10X, h5ad, RDS, and h5mu formats
- single-cell/preprocessing - Per-modality QC and normalization choice
- single-cell/doublet-detection - RNA-based and hashing doublet removal
- single-cell/clustering - Resolution sweep and cluster validation on the joint graph
- single-cell/markers-annotation - Marker discovery, manual labeling, and pseudobulk condition DE
- single-cell/cell-annotation - Automated reference-based label transfer from the RNA modality
- single-cell/differential-abundance - Test whether cell-type proportions shifted between conditions
- single-cell/multimodal-integration - WNN, totalVI/MultiVI, and MOFA joint-embedding details
- single-cell/scatac-analysis - ATAC-specific processing, LSI, gene activity, and AMULET fragment-based doublet detection
- differential-expression/deseq2-basics - Pseudobulk condition DE engine for aggregated counts
- atac-seq/single-cell-atac - Signac / ArchR / SnapATAC2 ecosystem decision; AMULET; cellranger-arc
- atac-seq/co-accessibility - Cicero / ArchR getCoAccessibility for cis-regulatory inference
- atac-seq/enhancer-gene-linking - ABC / ENCODE-rE2G for enhancer-gene mapping
- atac-seq/motif-deviation - chromVAR for per-cell TF motif activity
- atac-seq/footprinting - scprinter for sc footprinting
- workflows/grn-pipeline - Downstream: the paired object feeds SCENIC+ enhancer-GRN inference (Path B)
