Multimodal Integration

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

Common Modalities

Technology Modalities Package

CITE-seq RNA + surface proteins (ADT) Seurat

10X Multiome RNA + ATAC Seurat, Signac, ArchR

SHARE-seq RNA + ATAC Seurat, Signac

Spatial (Visium) RNA + spatial coordinates Seurat, Squidpy

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