Bio Applied Cite Seq Integration

CITE-seq and Multiome Data Integration

When to Use

• Normalizing antibody-derived tag (ADT) counts from a CITE-seq experiment before clustering
• Building a joint RNA+protein embedding (WNN) where surface markers should outweigh noisy mRNA for some cell types
• Integrating 10x Genomics Multiome (paired snRNA-seq + snATAC-seq from the same nucleus)
• Deciding between CLR and DSB normalization based on whether empty-droplet ("background") data is available
• Reproducing or extending an Azimuth/Seurat WNN PBMC reference (e.g. Hao et al. 2021, GSE164378)

Version Compatibility

Python: muon >=0.1.6, scanpy >=1.10, anndata >=0.10, Python >=3.10. R: Seurat >=5.0 (FindMultiModalNeighbors for WNN), SeuratData for reference objects.

Prerequisites

• pip install muon scanpy anndata (Python) or install.packages('Seurat') (R, v5+)
• Filtered and raw (raw_feature_bc_matrix.h5) count matrices from Cell Ranger if using DSB
• Familiarity with bio-applied-single-cell-scanpy (basic AnnData preprocessing) and bio-applied-scatac-chromatin (for the ATAC half of Multiome)

ADT Normalization

Goal: convert raw antibody-derived tag (ADT) counts into values comparable across proteins and cells, correcting for per-protein background that ordinary RNA normalization does not model.

Approach: ADT counts differ fundamentally from RNA — high background (~100-500 counts even in negative cells), essentially no true zeros, and a bimodal (negative/positive) distribution per protein. Use CLR when no empty-droplet data exists; use DSB when it does, since DSB explicitly subtracts a per-protein background estimated from empty droplets.

import numpy as np


def clr_normalize(X: np.ndarray) -> np.ndarray:
    """Centered log-ratio normalization for ADT counts.

    X: cells x proteins matrix of raw ADT counts.
    Matches Seurat's NormalizeData(assay='ADT', method='CLR', margin=2)
    (margin=2 = normalize across proteins within each cell).
    """
    X_ps = X + 0.5  # pseudocount avoids log(0)
    geo_mean = np.exp(np.log(X_ps).mean(axis=1, keepdims=True))
    return np.log(X_ps / geo_mean)


def dsb_normalize(X: np.ndarray, background: np.ndarray) -> np.ndarray:
    """Simplified DSB (Denoised and Scaled by Background) normalization.

    X: cells x proteins raw ADT counts (real cells).
    background: empty-droplet x proteins raw ADT counts (ambient noise only),
    pulled from the *raw* (unfiltered) Cell Ranger matrix.
    Values > ~3-4 indicate a clearly positive cell.
    """
    bg_mean = np.log1p(background).mean(axis=0)
    bg_std = np.log1p(background).std(axis=0) + 1e-6
    return (np.log1p(X) - bg_mean) / bg_std


if __name__ == "__main__":
    rng = np.random.default_rng(0)
    counts = rng.poisson(50, size=(200, 6)).astype(float)
    empty = rng.poisson(80, size=(500, 6)).astype(float)

    clr = clr_normalize(counts)
    dsb = dsb_normalize(counts, empty)
    assert clr.shape == counts.shape and np.isfinite(clr).all()
    assert dsb.shape == counts.shape and np.isfinite(dsb).all()
    print("CLR range:", clr.min().round(2), clr.max().round(2))
    print("DSB positive-cell fraction (>3):", (dsb > 3).mean().round(3))

Caveat: CLR assumes a roughly balanced antibody panel — if 80% of antibodies target one lineage, the per-cell geometric mean is skewed and CLR values become miscalibrated.

WNN Integration (Hao et al. 2021)

Goal: combine RNA and protein (or RNA and ATAC) into one neighbor graph where each cell gets its own per-modality weight, rather than one fixed weight for the whole dataset.

Approach: weighted-nearest-neighbor (WNN) learns, per cell, how well each modality's own neighbors predict the other modality — a clean, stable signal (e.g. CD4 surface protein) gets upweighted, while a noisy or low-abundance one (e.g. CD4 mRNA) gets downweighted for that cell. In muon, mu.pp.neighbors implements this automatically once each modality has its own PCA/neighbors computed.

import scanpy as sc
import muon as mu

# adata_rna, adata_prot: same cell barcodes, RNA and ADT AnnData objects
mdata = mu.MuData({"rna": adata_rna, "prot": adata_prot})

# RNA: standard log-normalize + PCA
sc.pp.normalize_total(mdata["rna"], target_sum=1e4)
sc.pp.log1p(mdata["rna"])
sc.pp.pca(mdata["rna"], n_comps=30)
sc.pp.neighbors(mdata["rna"])

# Protein: CLR (via muon's own implementation) + PCA
mu.prot.pp.clr(mdata["prot"])
sc.pp.pca(mdata["prot"], n_comps=15)
sc.pp.neighbors(mdata["prot"])

# WNN: per-cell modality weighting across both graphs
mu.pp.neighbors(mdata, key_added="wnn")
sc.tl.umap(mdata, neighbors_key="wnn")
sc.tl.leiden(mdata, neighbors_key="wnn", key_added="wnn_clusters")

10x Multiome (RNA + ATAC)

Goal: jointly analyze RNA and chromatin accessibility measured in the same nucleus, enabling direct peak-to-gene regulatory inference (unlike CITE-seq, which pairs RNA with protein).

Approach: treat ATAC like a second modality feeding into the same WNN machinery, but preprocess it with TF-IDF + LSI instead of PCA, and drop the first LSI component (it typically just tracks sequencing depth).

mdata = mu.read_10x_h5("filtered_feature_bc_matrix.h5")

# RNA: standard preprocessing
sc.pp.normalize_total(mdata["rna"], target_sum=1e4)
sc.pp.log1p(mdata["rna"])
sc.pp.pca(mdata["rna"])

# ATAC: TF-IDF + LSI (component 1 correlates with sequencing depth, drop it downstream)
mu.atac.pp.tfidf(mdata["atac"], scale_factor=1e4)
sc.tl.pca(mdata["atac"])

mu.pp.neighbors(mdata, key_added="wnn")
sc.tl.umap(mdata, neighbors_key="wnn")

R: Seurat FindMultiModalNeighbors

Goal: run the reference WNN implementation from Hao et al. 2021, as used to build the Azimuth PBMC atlas.

Approach: SCTransform the RNA assay, CLR-normalize the ADT assay across cells (margin = 2), run separate PCAs, then let FindMultiModalNeighbors learn per-cell modality weights.

library(Seurat)

# obj has 'RNA' and 'ADT' assays with matching cell barcodes
obj <- SCTransform(obj, assay = "RNA", verbose = FALSE)
obj <- RunPCA(obj, reduction.name = "pca", npcs = 30, verbose = FALSE)

DefaultAssay(obj) <- "ADT"
obj <- NormalizeData(obj, assay = "ADT", normalization.method = "CLR", margin = 2)
obj <- ScaleData(obj, assay = "ADT")
obj <- RunPCA(obj, reduction.name = "apca", npcs = 18, verbose = FALSE)

obj <- FindMultiModalNeighbors(
  obj,
  reduction.list = list("pca", "apca"),
  dims.list = list(1:30, 1:18)
)
obj <- RunUMAP(obj, nn.name = "weighted.nn")
obj <- FindClusters(obj, graph.name = "wsnn", algorithm = 3)

Pitfalls

• DSB needs the raw (unfiltered) matrix: filtered_feature_bc_matrix.h5 alone has no empty droplets — you must load raw_feature_bc_matrix.h5 to get background counts for DSB.
• CLR margin direction matters: Seurat's margin=2 (across proteins per cell) is standard for ADT; margin=1 (across cells per protein) is for RNA-like CLR and gives different, usually wrong, results here.
• Isotype controls: include isotype/IgG control antibodies in the panel and inspect them — a positive isotype signal in a cell indicates nonspecific binding, not a true positive.
• Multiome is nuclei, not cells: cytoplasmic/mitochondrial RNA QC thresholds from whole-cell scRNA-seq do not transfer directly; nuclear RNA content is lower and different filtering cutoffs apply.
• Same barcodes, same order: MuData/WNN assumes RNA and protein (or ATAC) AnnDatas share identical cell barcodes — mismatched or reordered barcodes silently pair the wrong cells.
• Batch effects across antibody lots: ADT panels are especially sensitive to lot-to-lot titration differences; check for batch confounding before trusting WNN clusters as purely biological.

See Also

• bio-applied-single-cell-scanpy — base AnnData/scanpy preprocessing before CITE-seq normalization
• bio-applied-scatac-chromatin — ATAC-specific QC and peak calling for the Multiome ATAC half
• bio-applied-cell-type-annotation — annotating clusters from the WNN embedding
• bio-applied-data-harmonization — batch correction when combining multiple CITE-seq/Multiome runs