Multi-Omics: CITE-seq & WNN Integration

CITE-seq measures RNA and antibody-derived tags (ADTs) in the same cells. This template follows the two assays from matched count matrices through separate preprocessing, a Seurat weighted-nearest-neighbor (WNN) graph, and checks that help you interpret modality weights without treating them as biological truth.

Learning objectives & prerequisites

By the end of this lesson you will be able to:

• Prepare a Seurat CITE-seq object containing RNA and antibody-derived tag (ADT) assays.
• Normalize and preprocess ADT data using CLR (centered log ratio) and run modality-specific dimensionality reduction.
• Run Seurat's Weighted Nearest Neighbor (WNN) workflow to integrate RNA and protein measurements, generate a joint embedding, and inspect per-cell modality weights.
• Interpret modality contributions and evaluate when protein or RNA drives neighborhood structure.

Prerequisites:

• Familiarity with single-cell RNA-seq and R. The function names below follow Seurat v4.3 and the current Seurat reference; check the installed version if an argument or object accessor differs.
• Installed software: R and Seurat. ggplot2 is needed only for plotting. The code uses the Seurat v4.3 WNN API and current documented function names; it is not a record of a tested software installation.
• Input files: two RDS matrices, with features in rows and cell barcodes in columns. The RNA matrix has genes; the ADT matrix has antibody features. They must use the same barcode convention. This example intersects their barcodes before object construction and reports how many cells are retained. It does not load an already-built Seurat object.

Why analyze RNA and ADT together?

CITE-seq pairs transcript counts with DNA-tagged antibody counts from the same cell. RNA and ADT assays have different feature sets and measurement properties, so concatenating their features directly can make one modality dominate a distance calculation. Seurat WNN instead learns cell-specific modality weights from modality-specific reductions and builds a shared neighbor graph.

WNN is an unsupervised way to combine local neighborhood information. A weight describes the relative utility of a modality for a cell’s neighborhood in this analysis; it is not a direct measure of assay quality, biological importance, or truth. Check each assay and the experimental context independently.

Build the Seurat object from matched count matrices

Replace the two input paths with your RDS files. The matrices should have unique feature and barcode names, with cells in columns. The code keeps only the shared barcodes and aligns ADT columns to the RNA object before adding the second assay. Check this retained-cell count against your intended filtering before proceeding.

Files you should have

• RNA count matrix (genes × cells), saved as an RDS matrix
• ADT count matrix (ADT features × cells), saved as an RDS matrix
• Optional: feature annotations and a barcode list for auditing the matched inputs.

No matrices are bundled. For a worked source-to-object example, follow the official Seurat multimodal vignette (CITE-seq CBMC RNA and ADT matrices); its CSV inputs must be converted/prepared as RDS objects for this script. Inputs may be base R or Matrix sparse matrices, with nonnegative finite raw counts, unique feature names and unique cell barcodes. Do not use normalized/scaled values as counts. Work from a project directory root so relative data/ and results/ paths resolve. Check prerequisites and create directories:

CBMC species checkpoint: the named GEO RNA matrix uses HUMAN_ and MOUSE_ prefixes. For this human-focused example, the script follows Seurat's CollapseSpeciesExpressionMatrix preparation before object creation: it removes the human prefix and retains 100 highly expressed mouse control features. This is not a joint analysis of both species. Do not strip prefixes blindly in a different study; confirm the species and feature annotation first. Without preparation, Seurat changes underscores to hyphens and a plain ^MT- or CD4 lookup misses the human features.

R --version
Rscript -e 'stopifnot(requireNamespace("Seurat", quietly=TRUE), requireNamespace("ggplot2", quietly=TRUE)); cat("Seurat", as.character(packageVersion("Seurat")), "\n")'
mkdir -p data results

For the vignette's downloadable CSV matrices (features in rows, barcodes in columns), convert them from this project directory with data/ already created; substitute the downloaded filenames. Keep the same cell-barcode names in both files. This small conversion checks matching cells and writes the RDS inputs expected below:

rna <- as.matrix(read.csv("RNA_umi.csv", row.names = 1, check.names = FALSE))
adt <- as.matrix(read.csv("ADT_umi.csv", row.names = 1, check.names = FALSE))
stopifnot(length(intersect(colnames(rna), colnames(adt))) > 0)
saveRDS(rna, "data/rna_counts.rds")
saveRDS(adt, "data/adt_counts.rds")

For large data, prefer a sparse-matrix-preserving reader/conversion rather than forcing dense matrices. Place/produce inputs at data/rna_counts.rds and data/adt_counts.rds. Save the following R code as cite_seq_wnn.R in this directory and invoke Rscript cite_seq_wnn.R. Stop if files or preflight checks fail; do not continue with a partial or empty barcode intersection.

R script: build a Seurat CITE-seq object and preprocess

# Load libraries
library(Seurat)
library(ggplot2)

# Read matrices with features in rows and cells in columns
rna_mat <- readRDS("data/rna_counts.rds")
adt_mat <- readRDS("data/adt_counts.rds")
stopifnot(is.matrix(rna_mat) || inherits(rna_mat, "Matrix"),
          is.matrix(adt_mat) || inherits(adt_mat, "Matrix"))
stopifnot(!is.null(rownames(rna_mat)), !is.null(colnames(rna_mat)),
          !is.null(rownames(adt_mat)), !is.null(colnames(adt_mat)))
stopifnot(!anyDuplicated(rownames(rna_mat)), !anyDuplicated(colnames(rna_mat)),
          !anyDuplicated(rownames(adt_mat)), !anyDuplicated(colnames(adt_mat)))
stopifnot(all(is.finite(rna_mat)), all(rna_mat >= 0),
          all(is.finite(adt_mat)), all(adt_mat >= 0))

# Keep shared cells in the same order in both matrices
shared_cells <- intersect(colnames(rna_mat), colnames(adt_mat))
if (length(shared_cells) == 0) stop("RNA and ADT have no matching cell barcodes")
rna_mat <- rna_mat[, shared_cells, drop = FALSE]
adt_mat <- adt_mat[, shared_cells, drop = FALSE]
message(length(shared_cells), " shared barcodes retained")

# Prepare the named human CBMC source before Seurat normalizes feature names.
# This route assumes HUMAN_ genes are of interest and MOUSE_ genes are controls.
if (any(startsWith(rownames(rna_mat), "HUMAN_")) &&
    any(startsWith(rownames(rna_mat), "MOUSE_"))) {
  rna_mat <- CollapseSpeciesExpressionMatrix(
    rna_mat, prefix = "HUMAN_", controls = "MOUSE_", ncontrols = 100
  )
  if (anyDuplicated(rownames(rna_mat))) stop("Species preparation produced duplicate feature names")
  message("Prepared human CBMC gene symbols; retained 100 mouse control features")
}

# Add ADT counts after aligning to the cells retained in the RNA object
seurat_obj <- CreateSeuratObject(counts = rna_mat, assay = "RNA", project = "CITE_project")
adt_mat <- adt_mat[, colnames(seurat_obj), drop = FALSE]
seurat_obj[["ADT"]] <- CreateAssayObject(counts = adt_mat)

dir.create("results", showWarnings = FALSE, recursive = TRUE)
saveRDS(seurat_obj, file = "results/seurat_initial.rds")

Check the input checkpoint:

• In-memory Seurat object with assays "RNA" and "ADT". The file results/seurat_initial.rds is a reproducible checkpoint.

Preprocess RNA (standard): filtering, normalization, PCA

# Basic QC: percent mitochondrial, per-cell counts, per-cell features
# Human-style symbols commonly use MT-; mouse-style symbols commonly use mt-.
# Inspect rownames and verify percent.mt is nonzero/plausible before filtering.
mito_pattern <- "^MT-"  # use "^mt-" for mouse-style symbols; adapt to reference
mito_features <- grep(mito_pattern, rownames(seurat_obj), value = TRUE)
if (length(mito_features) == 0) stop("No mitochondrial features matched; inspect names and set mito_pattern")
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, features = mito_features)
if (all(seurat_obj$percent.mt == 0)) stop("percent.mt is all zero; verify species-specific feature names")
summary(seurat_obj$percent.mt)

# Illustrative thresholds only: inspect QC distributions and adapt to tissue, species, chemistry, and study design
seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nCount_RNA > 500 & percent.mt < 20)

# Normalize and find variable features
DefaultAssay(seurat_obj) <- "RNA"
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 1e4)
seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 3000)
seurat_obj <- ScaleData(seurat_obj, features = VariableFeatures(seurat_obj))
seurat_obj <- RunPCA(seurat_obj, features = VariableFeatures(seurat_obj), reduction.name = "pca")

# Save
saveRDS(seurat_obj, file = "results/seurat_rna_pca.rds")

Produces:

• PCA embedding stored in seurat_obj[["pca"]].
• Results checkpoint results/seurat_rna_pca.rds.

Preprocess ADT: CLR normalization and protein PCA

# Switch to ADT assay
DefaultAssay(seurat_obj) <- "ADT"

# CLR normalization. margin = 2 performs per-cell centering across proteins (common default for CITE-seq)
# Note: margin = 1 normalizes across cells for each protein; choose based on experimental design.
seurat_obj <- NormalizeData(seurat_obj, normalization.method = "CLR", margin = 2)

# Scale and run PCA on all ADT features (usually few proteins; do not select HVGs)
adt_features <- rownames(seurat_obj[["ADT"]])
seurat_obj <- ScaleData(seurat_obj, features = adt_features)
adt_npcs <- min(18L, length(adt_features) - 1L, ncol(seurat_obj) - 1L)
if (adt_npcs < 2L) stop("Too few ADT features or cells for this PCA/WNN example")
seurat_obj <- RunPCA(seurat_obj, features = adt_features, npcs = adt_npcs, reduction.name = "apca")

# Save
saveRDS(seurat_obj, file = "results/seurat_adt_pca.rds")

Produces:

• ADT PCA stored in seurat_obj[["apca"]].
• Checkpoint results/seurat_adt_pca.rds.

Run WNN integration, embedding, and clustering

Use only PCs that actually exist in each reduction. The named CBMC panel has 13 ADT features, not enough for 18 protein PCs. The caps below are software limits, not a biological selection rule: inspect the reductions and choose fewer informative PCs when appropriate.

# Ensure RNA is default for the function that requires "pca" as an entry
DefaultAssay(seurat_obj) <- "RNA"

# Bound each requested range by the dimensions actually produced.
rna_available <- ncol(Embeddings(seurat_obj, reduction = "pca"))
adt_available <- ncol(Embeddings(seurat_obj, reduction = "apca"))
if (rna_available < 2L || adt_available < 2L) stop("Need at least two PCs per modality for this example")
rna_dims <- seq_len(min(30L, rna_available))
adt_dims <- seq_len(min(18L, adt_available))
message("WNN PCs: RNA=", length(rna_dims), "; ADT=", length(adt_dims))

# Find multi-modal neighbors with the inspected, available dimensions.
seurat_obj <- FindMultiModalNeighbors(
  seurat_obj,
  reduction.list = list("pca", "apca"),
  dims.list = list(rna_dims, adt_dims),
  modality.weight.name = c("RNA.weight", "ADT.weight")
)

# Build a UMAP on the weighted neighbor graph and cluster using the WNN graph
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)

# Save final integrated object
# Record exactly the clustering result produced by this run as active identities.
cluster_cols <- grep("^wsnn_res\\.", colnames(seurat_obj[[]]), value = TRUE)
if (length(cluster_cols) == 0) stop("No wsnn resolution column found; inspect metadata after FindClusters")
cluster_col <- tail(cluster_cols, 1)
Idents(seurat_obj) <- seurat_obj[[cluster_col, drop = TRUE]]
stopifnot(nlevels(Idents(seurat_obj)) > 1)
print(table(Idents(seurat_obj)))
saveRDS(seurat_obj, file = "results/seurat_wnn_integrated.rds")

Produces:

• Weighted neighbor graph stored internally (graph name "wsnn"), a UMAP reduction "wnn.umap", per-cell modality weights in metadata columns RNA.weight and ADT.weight, as named in this call.
• Checkpoint results/seurat_wnn_integrated.rds. FindClusters may create a newly named or suffixed metadata column on reruns; explicitly choose the intended run's cluster column as active identities before plotting.

Output inspection and interpretation

After running the pipeline, inspect the following elements. Do not treat any single plot or statistic as definitive; use them together to form a consistent interpretation.

Make concrete inspection plots

Run after integration; this writes five PNGs under results/. It assumes CD4 occurs in both assays; replace with features actually present if needed.

stopifnot(all(c("wnn.umap") %in% Reductions(seurat_obj)),
          all(c("RNA.weight", "ADT.weight") %in% colnames(seurat_obj[[]])))
print(summary(seurat_obj$RNA.weight)); print(summary(seurat_obj$ADT.weight))
stopifnot("CD4" %in% rownames(seurat_obj[["RNA"]]), "CD4" %in% rownames(seurat_obj[["ADT"]]))
p_wnn <- DimPlot(seurat_obj, reduction = "wnn.umap", group.by = "ident")
ggplot2::ggsave("results/wnn_umap.png", p_wnn, width = 7, height = 5, bg = "white")
p_rna_weight <- FeaturePlot(seurat_obj, features = "RNA.weight", reduction = "wnn.umap")
ggplot2::ggsave("results/rna_weight.png", p_rna_weight, width = 7, height = 5, bg = "white")
p_adt_weight <- FeaturePlot(seurat_obj, features = "ADT.weight", reduction = "wnn.umap")
ggplot2::ggsave("results/adt_weight.png", p_adt_weight, width = 7, height = 5, bg = "white")
DefaultAssay(seurat_obj) <- "RNA"
ggplot2::ggsave("results/cd4_rna.png", FeaturePlot(seurat_obj, features = "CD4", reduction = "wnn.umap"), width = 7, height = 5, bg = "white")
DefaultAssay(seurat_obj) <- "ADT"
ggplot2::ggsave("results/cd4_adt.png", FeaturePlot(seurat_obj, features = "CD4", reduction = "wnn.umap"), width = 7, height = 5, bg = "white")
DefaultAssay(seurat_obj) <- "RNA"

The UMAP should show points colored by active WNN cluster identity; weight plots show per-cell relative modality contributions, and the paired CD4 plots display separate assay measurements on the same coordinates. These checks establish that plots were generated, not that clusters are cell types, a weight is quality, or RNA and protein should agree.

Key items to inspect

• UMAP embedding: confirm plotted active identities correspond to the intended clustering run and graph.
• Per-cell modality weights: RNA.weight and ADT.weight are metadata columns named in this call. Inspect summaries and spatial distributions; a smooth gradient is not biological validation.
• Feature concordance: paired CD4 plots explicitly switch DefaultAssay between RNA and ADT on the same embedding. Compare only measured features; assay signals need not agree exactly.
• Cluster marker identification: Use FindMarkers on the integrated cluster labels but consider modality. Identify markers in RNA and ADT separately to verify.

What to look for and what it means

• When a cell has a higher ADT weight, that says ADT contributes more to its local neighbor definition under the chosen preprocessing and dimensions. It does not by itself explain why RNA and protein differ or establish a biological mechanism.
• Cells with high RNA weight: indicates neighborhoods where transcriptome structure is more informative than protein measurements (for example, transcriptionally distinct states without clear surface-protein differences).
• When both modality weights contribute, both reductions inform the local neighborhood. This does not demonstrate that RNA and protein agree on a cell identity; inspect the underlying measurements and markers.

What the weights and clusters do not establish

• WNN output does not establish causality: high protein weight does not imply protein measurement is biologically ground truth; it reflects local signal-to-noise relationships.
• Cluster labels inferred from WNN should not be considered validated cell types without verifying canonical markers, metadata, and experimental context.
• Discrepancies between RNA and ADT may arise from technical issues (e.g., ambient antibody background, poor antibody specificity) rather than biology.

Troubleshooting: symptoms, likely causes, and fixes

Below are common issues and concrete steps to diagnose/fix them.

1. ADT signal looks uniformly high or dominated by a few proteins

• Symptom: After CLR normalization, some ADT features dominate UMAP or PCA.
• Likely causes: Strongly expressed proteins, technical antibody overloading, or poor CLR choice (margin)
• Fixes:
  • Plot raw ADT counts by protein to identify outliers: boxplot/violin of counts per ADT.
  • Consider margin = 1 (per-protein across cells) if comparing the same protein across cell types is the primary goal; margin = 2 (default) standardizes each cell's protein profile.
  • Remove problematic antibody features (after careful documentation) if an antibody is clearly non-specific or saturated.

2. Few or no informative ADT PCs

• Symptom: APCA (ADT PCA) variances are flat and first PCs explain little variance.
• Likely causes: Low ADT complexity (few proteins) or poor signal-to-noise.
• Fixes:
  • Use more ADT features if available; include isotype controls to assess background.
  • Inspect per-protein mean and variance; if most proteins are near background, consider re-evaluating experimental protocol or antibody panel.
  • If you have very few proteins (e.g., <10), be cautious with PCA-based dimensionality reduction; you may use all proteins directly in neighbor calculations instead of many PCs.

3. WNN UMAP shows RNA-dominated structure or vice versa

• Symptom: The WNN UMAP resembles the RNA-only UMAP.
• Possible causes include the selected dimensions, weak or noisy ADT signal, or the UMAP projection. Similar-looking UMAPs alone do not prove that WNN ignored ADT.
• Fixes:
  • Re-evaluate dims.list to ensure ADT dims capture useful signal (try fewer RNA PCs or more ADT PCs if justified).
  • The call below names the two metadata columns `RNA.weight` and `ADT.weight`. Inspect both distributions and verify the selected reduction dimensions. Do not infer a defective assay from extreme weights alone.

4. Mismatched barcodes between RNA and ADT

• Symptom: errors when adding ADT assay or neighbor finding produces fewer cells than expected.
• Likely causes: barcode naming conventions differ (prefixes/suffixes) or missing barcodes due to filtering.
• Fixes:
  • Harmonize barcodes prior to CreateAssayObject. Use `intersect(colnames(rna_mat), colnames(adt_mat))`, then subset both matrices in that same order before object creation. Stop if no barcodes match; investigate prefixes, suffixes, and filtering rather than silently accepting an empty or unexpectedly small intersection.
  • Document any barcode filtering steps and save filtered matrices.

Limitations, version caveats, and responsible-use notes

Limitations to keep in mind:

• WNN is an unsupervised integration method. It reflects relative local informativeness but does not correct modality-specific artifacts such as antibody cross-reactivity or ambient RNA contamination.
• Normalization choices (CLR margin, RNA normalization method) substantially affect results. Carefully document settings, and, where possible, test alternatives (e.g., SCTransform for RNA).
• Results are sensitive to dim selection (dims.list). Inspect the available dimensions and their variance patterns before setting `dims.list`; the example values (30 RNA PCs and 18 ADT PCs) are placeholders, not validated choices for every dataset. Avoid requesting PCs that do not exist.
• Seurat APIs can change between versions. This workflow follows the official Seurat v4.3 example and current reference for these calls; check your installed version and save session information with `saveRDS(sessionInfo(), file = "results/session_info.rds")`.

Responsible-use notes:

• Do not treat clustering or automated labels as clinical diagnoses. Validate cell labels with orthogonal data or experiments when possible.
• Document all preprocessing and filtering steps. Share intermediate checkpoints and seed settings to enable reproducibility.

Concise next-step checklist

1. Confirm input matrices and barcode matching; save filtered inputs.
2. Run RNA and ADT preprocessing pipelines separately and save checkpoints (RDS files).
3. Run FindMultiModalNeighbors with a few sensible dims.list choices; inspect modality-weight distributions.
4. Visualize WNN UMAP and compare to RNA-only and ADT-only embeddings.
5. Identify cluster markers in both modalities; validate with known biology and additional controls.
6. Save final Seurat object and session info for provenance.

References

• Seurat Weighted Nearest Neighbor analysis vignette: official Seurat documentation and examples.
• Stoeckius et al., 2017: Simultaneous epitope and transcriptome measurement in single cells (Nature Biotechnology)
• Seurat project site: package documentation and release notes.
• Seurat NormalizeData reference: CLR margin argument and normalization methods.