Seurat v5 加权最近邻分析 (WNN) 教程¶
基于 Seurat v5.5.1 官方教程，使用 PBMC 3k RNA 数据模拟 CITE-seq 多模态场景。
从 RNA 表达矩阵模拟 ADT 蛋白数据，完整演示 WNN 多模态分析流程。
目录¶
- 1. 背景与原理
- 2. 模拟 ADT 蛋白数据
- 3. 独立预处理每个模态
- 4. 构建多模态邻居图 (WNN)
- 5. WNN UMAP 与聚类
- 6. 单模态 UMAP 对比
- 7. Marker 基因与蛋白可视化
- 8. 模态权重分析
- 9. ADT 蛋白表达 DotPlot
- 10. 多模态联合可视化
- 11. 总结
- 附录: 完整脚本
1. 背景与原理¶
1.1 什么是 WNN？¶
加权最近邻（Weighted Nearest Neighbor, WNN）是 Seurat v5 引入的多模态分析方法，用于同时分析来自不同模态的数据（如 RNA + 蛋白）。
传统方法将不同模态的数据简单拼接后分析，但不同模态的信息量不同。WNN 的核心思想是：为每个细胞自动学习不同模态的权重，根据每个细胞中哪个模态信息更丰富来自适应加权。
1.2 WNN 分析流程¶
模态1 (RNA) -> 预处理 -> PCA
|
-> FindMultiModalNeighbors -> WNN 图 -> UMAP + 聚类
|
模态2 (ADT) -> 预处理 -> CLR + PCA
1.3 CITE-seq 数据特点¶
CITE-seq 同时测量 RNA 表达和细胞表面蛋白（ADT），两种模态需要不同的预处理：
- RNA：标准 LogNormalize + HVG + ScaleData + PCA
- ADT：CLR（Centered Log Ratio）归一化 + PCA
2. 模拟 ADT 蛋白数据¶
由于官方 bmcite 数据集下载超时，我们从 PBMC 3k 的 RNA 表达矩阵模拟 ADT 蛋白信号：
# 选择 25 个关键 marker 基因模拟对应 ADT 蛋白
adt_genes <- c("CD3D", "CD3E", "CD4", "CD8A", "MS4A1", "CD14", "FCGR3A", ...)
# 从 RNA counts 模拟 ADT 信号
for (i in seq_along(adt_genes)) {
gene <- adt_genes[i]
rna_expr <- as.numeric(rna_counts[gene, ])
# ADT 信号 = RNA 表达 * 随机系数 + 基线噪声
adt_signal <- rna_expr * runif(ncol(pbmc), 0.8, 1.2) +
rnorm(ncol(pbmc), mean = 1, sd = 0.5)
adt_matrix[i, ] <- round(pmax(adt_signal, 0) * 10)
}
# 创建 ADT assay
pbmc[["ADT"]] <- CreateAssayObject(counts = adt_matrix)
注意：真实 CITE-seq 的 ADT 数据来自抗体捕获，与 RNA 是独立的测量。这里使用 RNA 模拟仅用于演示 WNN 分析流程。
3. 独立预处理每个模态¶
3.1 RNA 预处理（标准流程）¶
DefaultAssay(pbmc) <- "RNA"
pbmc <- NormalizeData(pbmc, verbose = FALSE)
pbmc <- FindVariableFeatures(pbmc, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, verbose = FALSE)
3.2 ADT 预处理（CLR 归一化）¶
DefaultAssay(pbmc) <- "ADT"
VariableFeatures(pbmc) <- rownames(pbmc[["ADT"]])
# CLR 归一化是 CITE-seq ADT 数据的标准方法
pbmc <- NormalizeData(pbmc, normalization.method = "CLR", margin = 2, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, reduction.name = "apca", verbose = FALSE)
CLR（Centered Log Ratio） 是处理组成性数据（如蛋白占比）的标准归一化方法，相比 LogNormalize 更适合 ADT 数据。
4. 构建多模态邻居图 (WNN)¶
pbmc <- FindMultiModalNeighbors(
pbmc,
reduction.list = list("pca", "apca"),
dims.list = list(1:30, 1:min(18, nrow(pbmc[["ADT"]]) - 1)),
modality.weight.name = "RNA.weight"
)
关键参数：
- reduction.list：各模态的降维结果
- dims.list：各模态使用的维度数
- modality.weight.name：权重存储的列名
FindMultiModalNeighbors 的输出：
- weighted.nn：加权邻居图
- wsnn：WNN SNN 图（用于聚类）
- RNA.weight：每个细胞的 RNA 模态权重（0-1）
5. WNN UMAP 与聚类¶
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn", reduction.name = "wnn.umap",
reduction.key = "wnnUMAP_", verbose = FALSE)
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3,
resolution = 0.5, verbose = FALSE)
WNN UMAP：聚类 vs 已知细胞类型¶
WNN UMAP 中，聚类结果与已知细胞类型高度一致，说明多模态信息有效提升了聚类准确性。
6. 单模态 UMAP 对比¶
6.1 仅 RNA vs 仅 ADT¶
6.2 三模态对比：RNA / ADT / WNN¶
关键发现：
- 仅 RNA 的 UMAP 已经能较好区分主要细胞类型
- 仅 ADT 的 UMAP 聚类效果较差（因为模拟数据只有 25 个蛋白特征）
- WNN UMAP 综合了两个模态的信息，聚类效果最佳
- 在真实 CITE-seq 数据中，ADT 模态通常对 T 细胞亚型区分更有优势
7. Marker 基因与蛋白可视化¶
7.1 RNA Marker 基因¶
DefaultAssay(pbmc) <- "RNA"
FeaturePlot(pbmc, features = c("CD3D", "MS4A1", "CD14", "GNLY"),
reduction = "wnn.umap", ncol = 2, max.cutoff = 3)
7.2 ADT 蛋白 Marker¶
DefaultAssay(pbmc) <- "ADT"
# 注意：Seurat 会自动将下划线 _ 替换为短横线 -
FeaturePlot(pbmc, features = c("adt-CD3D", "adt-MS4A1", "adt-CD14", "adt-GNLY"),
reduction = "wnn.umap", ncol = 2, max.cutoff = 2)
重要：Seurat 会自动将特征名中的下划线
_替换为短横线-，因此引用 ADT 特征时需用adt-CD3D而非adt_CD3D。
8. 模态权重分析¶
8.1 RNA 权重小提琴图¶
WNN 为每个细胞计算了一个 RNA 模态权重（0-1），权重越高表示该细胞中 RNA 信息越重要：
VlnPlot(pbmc, features = "RNA.weight", group.by = "celltype", sort = TRUE)
8.2 RNA 权重在 UMAP 上着色¶
FeaturePlot(pbmc, features = "RNA.weight", reduction = "wnn.umap")
权重解读：
- RNA 权重接近 0.5 表示两个模态贡献均衡
- 权重偏向 RNA（>0.5）的细胞，RNA 信息更关键
- 权重偏向 ADT（<0.5）的细胞，蛋白信息更关键
- 在真实数据中，T 细胞亚型通常更依赖 ADT 权重来区分
9. ADT 蛋白表达 DotPlot¶
DefaultAssay(pbmc) <- "ADT"
DotPlot(pbmc, features = adt_markers_subset, group.by = "celltype") + coord_flip()
DotPlot 同时展示了各细胞类型中 ADT 蛋白的表达比例（点大小）和表达水平（颜色深浅），是验证细胞类型注释的利器。
10. 多模态联合可视化¶
10.1 CD3D：RNA vs ADT¶
同一基因/蛋白在 RNA 和 ADT 模态中的表达对比：
10.2 CD14：RNA vs ADT¶
RNA vs ADT 对比：
- RNA 表达（红色）更稀疏，有明显的零值
- ADT 信号（绿色）更连续，背景噪声更高
- 两者的高表达区域一致，验证了模拟数据的合理性
11. 总结¶
核心流程¶
RNA 数据 -> NormalizeData -> FindVariableFeatures -> ScaleData -> PCA
|
-> FindMultiModalNeighbors -> WNN -> UMAP + 聚类
|
ADT 数据 -> NormalizeData(CLR) -> ScaleData -> PCA (apca)
关键知识点¶
- WNN 自动学习模态权重：每个细胞都有独立的 RNA/ADT 权重，自适应地依赖信息量更大的模态
- CLR 归一化：ADT 数据的标准归一化方法，适用于组成性数据
- 多 Assay 管理：Seurat 对象可以同时包含 RNA 和 ADT 两个 assay，通过
DefaultAssay切换 - 特征名转换：Seurat 自动将
_替换为-，引用 ADT 特征时需注意 - 模拟数据局限：由于 ADT 是从 RNA 模拟的，两个模态高度相关，WNN 优势不如真实 CITE-seq 数据明显
附录: 完整脚本¶
wnn_vignette.R¶
# ============================================================
# Seurat v5 阶段4-B：加权最近邻分析 (WNN)
# ============================================================
# 基于官方教程：https://satijalab.org/seurat/articles/weighted_nearest_neighbor_analysis
# 由于 bmcite 数据集下载超时，使用 PBMC 3k RNA 数据模拟 CITE-seq 多模态场景
# 从 RNA 表达矩阵模拟 ADT（抗体衍生标签）蛋白数据，
# 演示完整的 WNN 多模态分析流程
# ============================================================
library(Seurat)
library(ggplot2)
library(patchwork)
library(Matrix)
# 设置工作目录
setwd("E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial")
output_dir <- "output/wnn"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)
cat("=== Seurat v5 阶段4-B：加权最近邻分析 (WNN) ===\n\n")
# 保存图片函数
save_plot <- function(filename, plot_obj, width = 10, height = 8) {
filepath <- file.path(output_dir, filename)
ggsave(filepath, plot_obj, width = width, height = height, dpi = 150, bg = "white")
cat(sprintf(" 已保存: %s\n", filename))
}
# ============================================================
# 1. 加载 PBMC 数据
# ============================================================
cat("## 1. 加载 PBMC 数据并准备多模态对象\n\n")
pbmc <- readRDS("output/pbmc3k_final.rds")
pbmc$celltype <- Idents(pbmc)
cat(sprintf("加载完成：共 %d 个细胞，%d 个基因\n", ncol(pbmc), nrow(pbmc)))
cat(sprintf("细胞类型：%s\n", paste(levels(pbmc$celltype), collapse = ", ")))
# ============================================================
# 2. 模拟 ADT（蛋白）数据
# ============================================================
cat("\n## 2. 模拟 ADT 蛋白数据\n\n")
# 提取 RNA counts 矩阵
rna_counts <- GetAssayData(pbmc, assay = "RNA", layer = "counts")
# 选择关键 marker 基因来模拟对应的 ADT 蛋白
# 模拟 CITE-seq 中常见的 25 个抗体面板
adt_genes <- c(
"CD3D", "CD3E", # T cell markers
"CD4", # CD4 T
"CD8A", # CD8 T
"MS4A1", # B cell (CD19/CD20)
"CD14", # Monocyte
"FCGR3A", # CD16 Mono
"NCAM1", # NK (CD56)
"GNLY", "NKG7", # NK markers
"IL7R", # CD127
"CCR7", # Naive T
"SELL", # CD62L
"ITGAX", # CD11c (DC)
"PPBP", # Platelet
"CD34", # Progenitor
"CD27", # Memory
"CD28", # Costimulation
"PDCD1", # PD-1
"ICOS", # Activation
"HLA-DRA", # MHC-II
"CD68", # Macrophage
"ITGAM", # CD11b
"CDK6", # Cell cycle
"CD69" # Early activation
)
# 只保留存在的基因
adt_genes <- adt_genes[adt_genes %in% rownames(rna_counts)]
cat(sprintf("模拟 ADT 面板：共 %d 个蛋白标记\n", length(adt_genes)))
# 从 RNA 表达模拟 ADT 信号
# ADT 信号通常比 RNA 更连续、噪声更低
set.seed(42)
adt_matrix <- matrix(0, nrow = length(adt_genes), ncol = ncol(pbmc))
rownames(adt_matrix) <- paste0("adt_", adt_genes)
colnames(adt_matrix) <- colnames(pbmc)
for (i in seq_along(adt_genes)) {
gene <- adt_genes[i]
# 取 RNA 表达作为基础
rna_expr <- as.numeric(rna_counts[gene, ])
# 模拟 ADT：加入蛋白层面噪声 + 基线信号
# ADT 信号通常更连续（对数正态分布）
adt_signal <- rna_expr * runif(ncol(pbmc), 0.8, 1.2) + # 表达量乘以随机系数
rnorm(ncol(pbmc), mean = 1, sd = 0.5) # 加入基线噪声
adt_signal <- pmax(adt_signal, 0) # 确保非负
adt_matrix[i, ] <- round(adt_signal * 10) # 放大信号
}
adt_matrix <- Matrix(adt_matrix, sparse = TRUE)
# 创建 ADT assay 并添加到 Seurat 对象
pbmc[["ADT"]] <- CreateAssayObject(counts = adt_matrix)
cat("ADT assay 已创建并添加到 Seurat 对象\n")
cat(sprintf("ADT 维度：%d 蛋白 x %d 细胞\n", nrow(pbmc[["ADT"]]), ncol(pbmc[["ADT"]])))
# ============================================================
# 3. 独立预处理每个模态
# ============================================================
cat("\n## 3. 独立预处理 RNA 和 ADT 模态\n\n")
# --- RNA 预处理 ---
cat("预处理 RNA...\n")
DefaultAssay(pbmc) <- "RNA"
pbmc <- NormalizeData(pbmc, verbose = FALSE)
pbmc <- FindVariableFeatures(pbmc, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, verbose = FALSE)
cat("RNA 预处理完成\n")
# --- ADT 预处理 ---
# 使用 CLR（Centered Log Ratio）归一化，这是 CITE-seq ADT 数据的标准方法
cat("预处理 ADT（CLR 归一化）...\n")
DefaultAssay(pbmc) <- "ADT"
# 使用所有 ADT 特征进行降维
VariableFeatures(pbmc) <- rownames(pbmc[["ADT"]])
pbmc <- NormalizeData(pbmc, normalization.method = "CLR", margin = 2, verbose = FALSE)
pbmc <- ScaleData(pbmc, verbose = FALSE)
pbmc <- RunPCA(pbmc, reduction.name = "apca", verbose = FALSE)
cat("ADT 预处理完成\n")
# ============================================================
# 4. 构建多模态邻居图 (FindMultiModalNeighbors)
# ============================================================
cat("\n## 4. 构建多模态邻居图 (WNN)\n\n")
# 计算加权最近邻
# RNA 使用 30 个 PC，ADT 使用 18 个 PC（因为蛋白面板较小）
cat("正在计算 WNN...\n")
pbmc <- FindMultiModalNeighbors(
pbmc,
reduction.list = list("pca", "apca"),
dims.list = list(1:30, 1:min(18, nrow(pbmc[["ADT"]]) - 1)),
modality.weight.name = "RNA.weight"
)
cat("WNN 计算完成\n")
# ============================================================
# 5. 基于 WNN 的 UMAP 和聚类
# ============================================================
cat("\n## 5. WNN UMAP 与聚类\n\n")
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn",
reduction.name = "wnn.umap",
reduction.key = "wnnUMAP_", verbose = FALSE)
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3,
resolution = 0.5, verbose = FALSE)
cat("WNN 聚类完成\n")
# WNN UMAP（聚类 vs 细胞类型）
p1 <- DimPlot(pbmc, reduction = "wnn.umap", label = TRUE, repel = TRUE,
label.size = 3) + NoLegend() + ggtitle("WNN UMAP（聚类）")
p2 <- DimPlot(pbmc, reduction = "wnn.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 3) + NoLegend() +
ggtitle("WNN UMAP（已知细胞类型）")
save_plot("01_wnn_umap.png", p1 + p2, width = 16, height = 7)
# ============================================================
# 6. 单模态 UMAP 对比
# ============================================================
cat("\n## 6. 单模态 UMAP 对比\n\n")
# RNA UMAP
pbmc <- RunUMAP(pbmc, reduction = "pca", dims = 1:30, assay = "RNA",
reduction.name = "rna.umap", reduction.key = "rnaUMAP_",
verbose = FALSE)
# ADT UMAP
pbmc <- RunUMAP(pbmc, reduction = "apca", dims = 1:min(18, nrow(pbmc[["ADT"]]) - 1),
assay = "ADT",
reduction.name = "adt.umap", reduction.key = "adtUMAP_",
verbose = FALSE)
p3 <- DimPlot(pbmc, reduction = "rna.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend() +
ggtitle("仅 RNA UMAP")
p4 <- DimPlot(pbmc, reduction = "adt.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend() +
ggtitle("仅 ADT UMAP")
save_plot("02_rna_vs_adt_umap.png", p3 + p4, width = 16, height = 7)
# 三模态对比
p5 <- DimPlot(pbmc, reduction = "rna.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 2.5) + ggtitle("RNA")
p6 <- DimPlot(pbmc, reduction = "adt.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 2.5) + ggtitle("ADT")
p7 <- DimPlot(pbmc, reduction = "wnn.umap", group.by = "celltype",
label = TRUE, repel = TRUE, label.size = 2.5) + ggtitle("WNN")
save_plot("03_three_modality_comparison.png",
(p5 + p6 + p7) & NoLegend() & theme(plot.title = element_text(hjust = 0.5)),
width = 21, height = 7)
# ============================================================
# 7. Marker 基因/蛋白可视化
# ============================================================
cat("\n## 7. Marker 基因与蛋白可视化\n\n")
DefaultAssay(pbmc) <- "RNA"
# RNA marker 基因
p8 <- FeaturePlot(pbmc, features = c("CD3D", "MS4A1", "CD14", "GNLY"),
reduction = "wnn.umap", ncol = 2,
max.cutoff = 3, cols = c("lightgrey", "red"))
save_plot("04_rna_markers_wnn.png", p8, width = 12, height = 10)
# ADT 蛋白 marker
DefaultAssay(pbmc) <- "ADT"
p9 <- FeaturePlot(pbmc, features = c("adt-CD3D", "adt-MS4A1", "adt-CD14", "adt-GNLY"),
reduction = "wnn.umap", ncol = 2,
max.cutoff = 2, cols = c("lightgrey", "darkgreen"))
save_plot("05_adt_markers_wnn.png", p9, width = 12, height = 10)
cat("Marker 可视化完成\n")
# ============================================================
# 8. 模态权重分析
# ============================================================
cat("\n## 8. 模态权重分析\n\n")
# RNA 权重分布
p10 <- VlnPlot(pbmc, features = "RNA.weight", group.by = "celltype",
sort = TRUE, pt.size = 0.1) +
NoLegend() +
labs(title = "各细胞类型 RNA 模态权重",
x = "细胞类型", y = "RNA 权重") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("06_rna_weight_violin.png", p10, width = 10, height = 7)
# RNA 权重在 WNN UMAP 上着色
p11 <- FeaturePlot(pbmc, features = "RNA.weight", reduction = "wnn.umap",
cols = c("#3B528B", "#21918C", "#5EC962", "#FDE725")) +
ggtitle("RNA 模态权重分布")
save_plot("07_rna_weight_umap.png", p11, width = 10, height = 8)
cat(sprintf("平均 RNA 权重: %.3f\n", mean(pbmc$RNA.weight)))
cat(sprintf("RNA 权重范围: [%.3f, %.3f]\n", min(pbmc$RNA.weight), max(pbmc$RNA.weight)))
# 各细胞类型平均 RNA 权重
cat("\n各细胞类型平均 RNA 权重:\n")
weight_by_type <- tapply(pbmc$RNA.weight, pbmc$celltype, mean)
weight_by_type <- sort(weight_by_type, decreasing = TRUE)
for (ct in names(weight_by_type)) {
cat(sprintf(" %s: %.3f\n", ct, weight_by_type[ct]))
}
# ============================================================
# 9. ADT 蛋白表达热图
# ============================================================
cat("\n## 9. ADT 蛋白表达热图\n\n")
DefaultAssay(pbmc) <- "ADT"
adt_markers_subset <- c("adt-CD3D", "adt-CD8A", "adt-MS4A1", "adt-CD14",
"adt-FCGR3A", "adt-GNLY", "adt-IL7R", "adt-NCAM1",
"adt-ITGAX", "adt-PPBP")
adt_markers_subset <- adt_markers_subset[adt_markers_subset %in% rownames(pbmc[["ADT"]])]
p12 <- DotPlot(pbmc, features = adt_markers_subset, group.by = "celltype") +
coord_flip() +
labs(title = "ADT 蛋白表达 DotPlot",
x = "蛋白标记", y = "细胞类型") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("08_adt_dotplot.png", p12, width = 10, height = 8)
# ============================================================
# 10. 多模态联合可视化
# ============================================================
cat("\n## 10. 多模态联合可视化\n\n")
# RNA vs ADT 同一基因对比
DefaultAssay(pbmc) <- "RNA"
p13 <- FeaturePlot(pbmc, features = "CD3D", reduction = "wnn.umap",
max.cutoff = 3, cols = c("lightgrey", "red")) +
ggtitle("RNA: CD3D")
DefaultAssay(pbmc) <- "ADT"
p14 <- FeaturePlot(pbmc, features = "adt-CD3D", reduction = "wnn.umap",
max.cutoff = 2, cols = c("lightgrey", "darkgreen")) +
ggtitle("ADT: CD3D")
save_plot("09_rna_vs_adt_cd3d.png", p13 | p14, width = 16, height = 7)
DefaultAssay(pbmc) <- "RNA"
p15 <- FeaturePlot(pbmc, features = "CD14", reduction = "wnn.umap",
max.cutoff = 3, cols = c("lightgrey", "red")) +
ggtitle("RNA: CD14")
DefaultAssay(pbmc) <- "ADT"
p16 <- FeaturePlot(pbmc, features = "adt-CD14", reduction = "wnn.umap",
max.cutoff = 2, cols = c("lightgrey", "darkgreen")) +
ggtitle("ADT: CD14")
save_plot("10_rna_vs_adt_cd14.png", p15 | p16, width = 16, height = 7)
# ============================================================
# 11. 保存结果
# ============================================================
cat("\n## 11. 保存结果\n\n")
# 恢复默认 assay
DefaultAssay(pbmc) <- "RNA"
saveRDS(pbmc, file.path(output_dir, "pbmc_wnn.rds"))
cat("已保存: pbmc_wnn.rds\n")
# 保存模态权重
write.csv(data.frame(
cell = colnames(pbmc),
celltype = pbmc$celltype,
RNA.weight = pbmc$RNA.weight
), file.path(output_dir, "modality_weights.csv"), row.names = FALSE)
cat("已保存: modality_weights.csv\n")
cat("\n=== 阶段4-B 加权最近邻分析 (WNN) 完成！===\n")
cat(sprintf("共生成 10 张图 + 1 个CSV + 1 个RDS文件\n"))
本教程基于 Seurat v5.5.1 官方文档整理
R 4.6.1 + Seurat 5.5.1 + Windows 环境
阶段4-B: 加权最近邻分析（WNN / CLR 归一化 / 模态权重 / 多模态可视化）