Bio Single Cell Multimodal Integration

Multimodal Integration

Analyze multi-modal single-cell data where multiple measurements are made per cell.

Common Modalities

CITE-seq Analysis (Seurat)

Load Data

library(Seurat)

# Read 10X data with antibody capture
data <- Read10X('filtered_feature_bc_matrix/')

# Separate RNA and ADT
rna_counts <- data$`Gene Expression`
adt_counts <- data$`Antibody Capture`

# Create Seurat object with both assays
obj <- CreateSeuratObject(counts = rna_counts, assay = 'RNA')
obj[['ADT']] <- CreateAssayObject(counts = adt_counts)

QC and Normalization

# RNA QC (standard)
obj <- PercentageFeatureSet(obj, pattern = '^MT-', col.name = 'percent.mt')
obj <- subset(obj, nFeature_RNA > 200 & percent.mt < 20)

# Normalize RNA
obj <- NormalizeData(obj, assay = 'RNA')
obj <- FindVariableFeatures(obj, assay = 'RNA')
obj <- ScaleData(obj, assay = 'RNA')

# Normalize ADT (CLR normalization)
obj <- NormalizeData(obj, assay = 'ADT', normalization.method = 'CLR', margin = 2)
obj <- ScaleData(obj, assay = 'ADT')

Weighted Nearest Neighbor (WNN) Clustering

# Dimensionality reduction for each modality
obj <- RunPCA(obj, assay = 'RNA', reduction.name = 'pca')
obj <- RunPCA(obj, assay = 'ADT', reduction.name = 'apca',
              features = rownames(obj[['ADT']]))

# WNN graph combining both modalities
obj <- FindMultiModalNeighbors(obj,
    reduction.list = list('pca', 'apca'),
    dims.list = list(1:30, 1:18))

# Cluster on WNN graph
obj <- FindClusters(obj, graph.name = 'wsnn', resolution = 0.5)

# UMAP on WNN
obj <- RunUMAP(obj, nn.name = 'weighted.nn', reduction.name = 'wnn.umap')

Visualize

# UMAP colored by cluster
DimPlot(obj, reduction = 'wnn.umap', label = TRUE)

# ADT expression on UMAP
FeaturePlot(obj, features = c('adt_CD3', 'adt_CD19', 'adt_CD14'),
            reduction = 'wnn.umap')

# Compare modality weights
VlnPlot(obj, features = 'RNA.weight', group.by = 'seurat_clusters')

10X Multiome (RNA + ATAC)

Load Data

library(Seurat)
library(Signac)

# Read RNA counts
rna_counts <- Read10X_h5('filtered_feature_bc_matrix.h5')$`Gene Expression`

# Read ATAC fragments
atac_counts <- Read10X_h5('filtered_feature_bc_matrix.h5')$Peaks
fragments <- CreateFragmentObject('atac_fragments.tsv.gz')

# Create multiome object
obj <- CreateSeuratObject(counts = rna_counts, assay = 'RNA')
obj[['ATAC']] <- CreateChromatinAssay(counts = atac_counts, fragments = fragments,
                                       genome = 'hg38', min.cells = 5)

Process ATAC

# ATAC QC
obj <- NucleosomeSignal(obj)
obj <- TSSEnrichment(obj)

# ATAC normalization
obj <- RunTFIDF(obj, assay = 'ATAC')
obj <- FindTopFeatures(obj, assay = 'ATAC', min.cutoff = 'q0')
obj <- RunSVD(obj, assay = 'ATAC')

Joint Analysis

# RNA processing
DefaultAssay(obj) <- 'RNA'
obj <- NormalizeData(obj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()

# WNN integration
obj <- FindMultiModalNeighbors(obj, reduction.list = list('pca', 'lsi'),
                                dims.list = list(1:30, 2:30))
obj <- RunUMAP(obj, nn.name = 'weighted.nn', reduction.name = 'wnn.umap')
obj <- FindClusters(obj, graph.name = 'wsnn')

Scanpy/MuData (Python)

CITE-seq with MuData

import scanpy as sc
import muon as mu
from muon import prot as pt

# Load multimodal data
mdata = mu.read_10x_h5('filtered_feature_bc_matrix.h5')

# Access modalities
rna = mdata.mod['rna']
prot = mdata.mod['prot']

# Process RNA
sc.pp.filter_cells(rna, min_genes=200)
sc.pp.normalize_total(rna, target_sum=1e4)
sc.pp.log1p(rna)
sc.pp.highly_variable_genes(rna)
sc.tl.pca(rna)

# Process protein (CLR normalization)
pt.pp.clr(prot)

# Multi-omics factor analysis
mu.tl.mofa(mdata, n_factors=20)

# Joint UMAP
mu.tl.umap(mdata)
mu.pl.umap(mdata, color=['rna:leiden', 'prot:CD3'])

Integration Metrics

Modality Weights

# Check how much each modality contributes per cell
weights <- obj@reductions$wnn@misc$weights

# Average weight by cluster
aggregate(weights, by = list(obj$seurat_clusters), mean)

Correlation Between Modalities

import numpy as np

# Correlate RNA and protein for same genes/proteins
common = set(rna.var_names) & set(prot.var_names)
for gene in common:
    rna_expr = rna[:, gene].X.toarray().flatten()
    prot_expr = prot[:, gene].X.toarray().flatten()
    corr = np.corrcoef(rna_expr, prot_expr)[0, 1]
    print(f'{gene}: r={corr:.3f}')

Marker Discovery

Multi-Modal Markers

# Find markers using both modalities
DefaultAssay(obj) <- 'RNA'
rna_markers <- FindAllMarkers(obj, only.pos = TRUE)

DefaultAssay(obj) <- 'ADT'
adt_markers <- FindAllMarkers(obj, only.pos = TRUE)

# Combine
all_markers <- rbind(
    transform(rna_markers, modality = 'RNA'),
    transform(adt_markers, modality = 'ADT')
)

Related Skills

• single-cell/data-io - Loading single-cell data
• single-cell/clustering - Clustering methods
• single-cell/markers-annotation - Cell type annotation
• chip-seq/peak-calling - For ATAC peak calling