Chapter 12: Single-Cell Multi-Omics: CITE-Seq, SHARE-Seq & Single-Cell ATAC-Seq

Chapter 12: Single-Cell Multi-Omics: CITE-Seq, SHARE-Seq & Single-Cell ATAC-Seq

While scRNA-seq captures gene transcription, cellular function is governed by multiple interconnected molecular layers: cell-surface protein abundance, epigenetic chromatin accessibility, and transcription factor binding. Single-cell multi-omics assays capture multiple molecular modalities from the exact same physical cell.

1. CITE-Seq Surface Protein Analysis in Seurat

Because protein tags (ADTs) have completely different background noise and sequencing dynamics than mRNA, ADT counts must be normalized using Centered Log-Ratio (CLR) normalization across cells, not standard library size scaling.

library(Seurat)

# Create multimodal Seurat object containing RNA and ADT assays
pbmc_multimodal <- CreateSeuratObject(counts = rna_counts)
pbmc_multimodal[["ADT"]] <- CreateAssayObject(counts = adt_counts)

# Normalize RNA with SCTransform
pbmc_multimodal <- SCTransform(pbmc_multimodal, assay = "RNA")

# Normalize ADT using Centered Log-Ratio (CLR) across cells (margin = 2)
DefaultAssay(pbmc_multimodal) <- "ADT"
pbmc_multimodal <- NormalizeData(pbmc_multimodal, normalization.method = "CLR", margin = 2)
pbmc_multimodal <- ScaleData(pbmc_multimodal)

# Multi-Modal Integration using Weighted Nearest Neighbor (WNN) Analysis
# WNN calculates cell-specific modality weights based on signal-to-noise ratio
pbmc_multimodal <- FindMultiModalNeighbors(
    pbmc_multimodal, 
    reduction.list = list("pca", "apca"), # PCA on RNA, APCA on ADT
    dims.list = list(1:30, 1:18)
)

# Run WNN UMAP
pbmc_multimodal <- RunUMAP(pbmc_multimodal, nn.name = "weighted.nn", reduction.name = "wnn.umap")
DimPlot(pbmc_multimodal, reduction = "wnn.umap", label = TRUE)

2. Single-Cell Chromatin Accessibility with Signac

Chromatin accessibility data is inherently binary (a locus is either cleaved by Tn5 or not). Analyzing scATAC-seq requires Term Frequency-Inverse Document Frequency (TF-IDF) normalization and Singular Value Decomposition (SVD), termed Latent Semantic Indexing (LSI).

library(Signac)
library(GenomicRanges)

# Load fragment file and peak matrix
chrom_assay <- CreateChromatinAssay(
    counts = peak_matrix, 
    sep = c(":", "-"), 
    fragments = "fragments.tsv.gz", 
    annotation = annotations_granges
)
pbmc_atac <- CreateSeuratObject(counts = chrom_assay, assay = "peaks")

# 1. TF-IDF Normalization
pbmc_atac <- RunTFIDF(pbmc_atac)

# 2. Select top accessible features and run SVD (LSI)
pbmc_atac <- FindTopFeatures(pbmc_atac, min.cutoff = 'q0')
pbmc_atac <- RunSVD(pbmc_atac) # Dim 1 typically correlates with sequencing depth

# 3. UMAP on LSI dimensions 2 to 30 (excluding depth-correlated Dim 1)
pbmc_atac <- RunUMAP(pbmc_atac, reduction = 'lsi', dims = 2:30)

3. Transcription Factor Motif Accessibility & Footprinting

By mapping accessible peaks against the JASPAR transcription factor motif database, tools like chromVAR calculate motif accessibility deviation scores across single cells, revealing which regulatory master factors (e.g. GATA, PU.1, NF-\(\kappa\)B) are actively binding DNA.

library(JASPAR2020)
library(TFBSTools)

# Query vertebrate PWM matrices from JASPAR
pfm <- getMatrixSet(x = JASPAR2020, opts = list(collection = "CORE", tax_group = "vertebrates", all_versions = FALSE))

# Add motif annotations to chromatin assay
pbmc_atac <- AddMotifs(object = pbmc_atac, genome = BSgenome.Hsapiens.UCSC.hg38, pfm = pfm)