bio-workflows-multiome-pipeline

name	bio-workflows-multiome-pipeline
description	End-to-end multiome workflow for joint scRNA-seq + scATAC-seq analysis. Covers data loading, separate modality processing, and WNN integration with Seurat/Signac. Use when analyzing joint scRNA+scATAC data.
tool_type	r
primary_tool	Seurat
workflow	true
depends_on	["single-cell/data-io","single-cell/preprocessing","single-cell/clustering","single-cell/multimodal-integration","single-cell/scatac-analysis"]
qc_checkpoints	[{"after_loading":"Both modalities detected per cell"},{"after_rna_qc":"RNA quality filters passed"},{"after_atac_qc":"TSS enrichment >2, nucleosome signal <4"},{"after_wnn":"Joint embedding separates cell types"}]

Multiome Pipeline

Complete workflow for 10X Multiome (joint scRNA + scATAC) analysis using Seurat and Signac.

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)

# Load RNA
rna_counts <- Read10X_h5('filtered_feature_bc_matrix.h5')
# For multiome, this returns a list with 'Gene Expression' and 'Peaks'

# Create Seurat object with RNA
seurat_obj <- CreateSeuratObject(
    counts = rna_counts$`Gene Expression`,
    assay = 'RNA'
)

# Load ATAC
atac_counts <- rna_counts$Peaks
# Or from fragments file
frags <- CreateFragmentObject('atac_fragments.tsv.gz', cells = colnames(seurat_obj))

# Create ChromatinAssay
atac_assay <- CreateChromatinAssay(
    counts = atac_counts,
    sep = c(':', '-'),
    fragments = frags,
    annotation = GetGRangesFromEnsDbensdb  EnsDb.Hsapiens.v86


seurat_obj  atac_assay

Step 2: RNA Quality Control and Processing

# QC metrics
seurat_obj[['percent.mt']] <- PercentageFeatureSet(seurat_obj, pattern = '^MT-')

# Filter
seurat_obj <- subset(seurat_obj,
    nCount_RNA > 1000 &
    nCount_RNA < 25000 &
    percent.mt < 20
)

# Normalize RNA
seurat_obj <- SCTransform(seurat_obj, assay = 'RNA', verbose = FALSE)

# PCA
seurat_obj <- RunPCA(seurat_obj, assay = 'SCT', verbose = FALSE)

Step 3: ATAC Quality Control and Processing

# ATAC QC metrics
DefaultAssay(seurat_obj) <- 'ATAC'

seurat_obj <- NucleosomeSignal(seurat_obj)
seurat_obj <- TSSEnrichment(seurat_obj)

# Visualize
VlnPlot(seurat_obj, features = c('nCount_ATAC', 'TSS.enrichment', 'nucleosome_signal'),
        pt.size = 0, ncol = 3)

# Filter ATAC
seurat_obj <- subset(seurat_obj,
    nCount_ATAC > 1000 &
    nCount_ATAC < 100000 &
    TSS.enrichment > 2 &
    nucleosome_signal < 4
)

# Normalize ATAC (TF-IDF + SVD = LSI)
seurat_obj <- RunTFIDF(seurat_obj)
seurat_obj <- FindTopFeatures(seurat_obj, min.cutoff  
seurat_obj  RunSVDseurat_obj


DepthCorseurat_obj

Step 4: Weighted Nearest Neighbors (WNN)

# Build WNN graph using both modalities
seurat_obj <- FindMultiModalNeighbors(
    seurat_obj,
    reduction.list = list('pca', 'lsi'),
    dims.list = list(1:30, 2:30),  # Skip LSI component 1 if depth-correlated
    modality.weight.name = 'RNA.weight'
)

# UMAP on WNN graph
seurat_obj <- RunUMAP(seurat_obj, nn.name = 'weighted.nn',
                       reduction.name = 'wnn.umap', reduction.key = 'wnnUMAP_')

# Cluster on WNN
seurat_obj <- FindClusters(seurat_obj, graph.name = 'wsnn',
                            algorithm = 3, resolution = 0.5, verbose  

Step 5: Visualization and Markers

# Compare modality-specific and joint embeddings
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

# Modality weights per cell
VlnPlot(seurat_obj, features = 'RNA.weight', group.by = 'seurat_clusters', pt.size = 0)

# Find markers (RNA)
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

# Link peaks to genes
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')  # Example genes
)

# Visualize links
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)

# Configuration
data_dir <- 'multiome_output'
output_dir <- 'multiome_results'
dir.create(output_dir, showWarnings = FALSE)

# === Load Data ===
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')
seurat_obj[['ATAC']] <- CreateChromatinAssay(
    counts = counts$Peaks
    sep   
    fragments  frags
    annotation  GetGRangesFromEnsDbensdb  EnsDb.Hsapiens.v86

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   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 

Related Skills

• single-cell/data-io - Loading 10X data
• single-cell/preprocessing - QC and normalization
• single-cell/multimodal-integration - WNN details
• single-cell/scatac-analysis - ATAC-specific processing