20260912 121054 Seurat V5 参考数据集映射教程

# 加载已注释的 PBMC 数据
pbmc <- readRDS("output/pbmc3k_final.rds")
pbmc$seurat_annotations <- Idents(pbmc)

# 随机拆分
set.seed(42)
ref_idx <- sample(1:ncol(pbmc), size = floor(ncol(pbmc) * 0.7))
pancreas.ref <- subset(pbmc, cells = cell_ids[ref_idx])      # 1846 细胞
pancreas.query <- subset(pbmc, cells = cell_ids[-ref_idx])   # 792 细胞

pancreas.ref <- NormalizeData(pancreas.ref, verbose = FALSE)
pancreas.ref <- FindVariableFeatures(pancreas.ref, verbose = FALSE)
pancreas.ref <- ScaleData(pancreas.ref, verbose = FALSE)
pancreas.ref <- RunPCA(pancreas.ref, verbose = FALSE)
pancreas.ref <- FindNeighbors(pancreas.ref, dims = 1:30, verbose = FALSE)
pancreas.ref <- FindClusters(pancreas.ref, verbose = FALSE)
pancreas.ref <- RunUMAP(pancreas.ref, dims = 1:30, verbose = FALSE)

pancreas.anchors <- FindTransferAnchors(
  reference = pancreas.ref,
  query = pancreas.query,
  dims = 1:30,
  reference.reduction = "pca"
)

predictions <- TransferData(
  anchorset = pancreas.anchors,
  refdata = pancreas.ref$seurat_annotations,
  dims = 1:30
)
pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

markers <- c("MS4A1", "CD3D", "CD3E", "IL7R", "CD14", "LYZ", "GNLY", "NKG7", "FCGR3A", "PPBP")
VlnPlot(pancreas.query, features = markers, group.by = "predicted.id", ncol = 5, pt.size = 0)

# 参考集 UMAP 需要保存模型
pancreas.ref <- RunUMAP(pancreas.ref, dims = 1:30, reduction = "pca",
                        return.model = TRUE, verbose = FALSE)

# MapQuery 一步到位
pancreas.query <- MapQuery(
  anchorset = pancreas.anchors,
  reference = pancreas.ref,
  query = pancreas.query,
  refdata = list(celltype = "seurat_annotations"),
  reference.reduction = "pca",
  reduction.model = "umap"
)

ggplot(pancreas.query@meta.data, aes(x = predicted.id, y = prediction.score.max)) +
  geom_boxplot(fill = "#4C9BE8", alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.5)

参考集（已注释）-> 预处理 + PCA + UMAP（保存模型）
                           |
查询集（未注释）-> 预处理 -> FindTransferAnchors
                           |
                    TransferData -> 预测标签 + 置信度
                           |
                    MapQuery -> 投影到参考 UMAP 空间

# ============================================================
# Seurat v5 阶段4-A：参考数据集映射 (Mapping & Annotation)
# ============================================================
# 基于官方教程：https://satijalab.org/seurat/articles/integration_mapping
# 由于 panc8 数据集下载超时，使用已有的 PBMC 3k 数据模拟参考映射场景
# 将 PBMC 数据拆分为参考集（70%）和查询集（30%），
# 演示完整的标签转移 + UMAP 投影流程
# ============================================================

library(Seurat)
library(ggplot2)
library(patchwork)

# 设置工作目录
setwd("E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial")
output_dir <- "output/mapping"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

cat("=== Seurat v5 阶段4-A：参考数据集映射 ===\n\n")

# ============================================================
# 1. 加载已注释的 PBMC 数据
# ============================================================
cat("## 1. 加载已注释的 PBMC 数据\n\n")

# 加载阶段1完成的带有细胞类型注释的 Seurat 对象
pbmc <- readRDS("output/pbmc3k_final.rds")
# 将 active.ident（细胞类型注释）保存为 metadata 列
pbmc$seurat_annotations <- Idents(pbmc)
cat(sprintf("加载完成：共 %d 个细胞，%d 个基因\n", ncol(pbmc), nrow(pbmc)))
cat(sprintf("已知细胞类型：%s\n", paste(levels(pbmc$seurat_annotations), collapse = ", ")))

# 保存图片函数
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))
}

# ============================================================
# 2. 拆分为参考集和查询集
# ============================================================
cat("\n## 2. 拆分为参考集和查询集\n\n")

set.seed(42)
cell_ids <- colnames(pbmc)
n_cells <- length(cell_ids)
# 70% 参考集，30% 查询集
ref_idx <- sample(1:n_cells, size = floor(n_cells * 0.7))
ref_cells <- cell_ids[ref_idx]
query_cells <- cell_ids[-ref_idx]

cat(sprintf("参考集：%d 个细胞\n", length(ref_cells)))
cat(sprintf("查询集：%d 个细胞\n", length(query_cells)))

# 提取参考集和查询集
pancreas.ref <- subset(pbmc, cells = ref_cells)
pancreas.query <- subset(pbmc, cells = query_cells)

# 确保使用 seurat_annotations 作为 Idents
Idents(pancreas.ref) <- "seurat_annotations"
Idents(pancreas.query) <- "seurat_annotations"

# 记录参考集的真实细胞类型（用于后续验证）
cat(sprintf("\n参考集细胞类型分布:\n"))
print(table(pancreas.ref$seurat_annotations))

cat(sprintf("\n查询集真实细胞类型分布（用于验证预测准确性）:\n"))
print(table(pancreas.query$seurat_annotations))

# 保存查询集的真实标签用于后续比较
query_true_labels <- pancreas.query$seurat_annotations

# ============================================================
# 3. 构建参考集（预处理 + 聚类 + UMAP）
# ============================================================
cat("\n## 3. 构建参考集\n\n")

# 参考集预处理
pancreas.ref <- NormalizeData(pancreas.ref, verbose = FALSE)
pancreas.ref <- FindVariableFeatures(pancreas.ref, verbose = FALSE)
pancreas.ref <- ScaleData(pancreas.ref, verbose = FALSE)
pancreas.ref <- RunPCA(pancreas.ref, verbose = FALSE)
pancreas.ref <- FindNeighbors(pancreas.ref, dims = 1:30, verbose = FALSE)
pancreas.ref <- FindClusters(pancreas.ref, verbose = FALSE)
pancreas.ref <- RunUMAP(pancreas.ref, dims = 1:30, verbose = FALSE)

cat("参考集预处理完成\n")

# 参考集 UMAP（按已知细胞类型着色）
p1 <- DimPlot(pancreas.ref, reduction = "umap", group.by = "seurat_annotations",
              label = TRUE, label.size = 4, repel = TRUE) +
  ggtitle("参考集 UMAP（已知注释）")
save_plot("01_reference_umap.png", p1)

# ============================================================
# 4. 查找转移锚点 (FindTransferAnchors)
# ============================================================
cat("\n## 4. 查找转移锚点\n\n")

# 查询集预处理
pancreas.query <- NormalizeData(pancreas.query, verbose = FALSE)

# 查找锚点
cat("正在查找转移锚点...\n")
pancreas.anchors <- FindTransferAnchors(
  reference = pancreas.ref,
  query = pancreas.query,
  dims = 1:30,
  reference.reduction = "pca"
)
cat("锚点查找完成\n")

# ============================================================
# 5. 标签转移 (TransferData)
# ============================================================
cat("\n## 5. 标签转移 (TransferData)\n\n")

# 转移细胞类型标签
predictions <- TransferData(
  anchorset = pancreas.anchors,
  refdata = pancreas.ref$seurat_annotations,
  dims = 1:30
)
pancreas.query <- AddMetaData(pancreas.query, metadata = predictions)

cat("标签转移完成\n")
cat(sprintf("\n预测细胞类型分布:\n"))
print(table(pancreas.query$predicted.id))

# ============================================================
# 6. 评估预测准确性
# ============================================================
cat("\n## 6. 评估预测准确性\n\n")

# 比较预测标签和真实标签
pancreas.query$prediction.match <- pancreas.query$predicted.id == pancreas.query$seurat_annotations
match_result <- table(pancreas.query$prediction.match)
cat(sprintf("预测结果:\n"))
print(match_result)
cat(sprintf("\n预测准确率: %.1f%%\n", 
    round(match_result["TRUE"] / sum(match_result) * 100, 1)))

# 按细胞类型查看预测准确性
cat("\n各细胞类型预测准确性:\n")
acc_by_type <- sapply(levels(pancreas.query$seurat_annotations), function(ct) {
  idx <- pancreas.query$seurat_annotations == ct
  if (sum(idx) > 0) {
    correct <- sum(pancreas.query$predicted.id[idx] == ct)
    sprintf("%s: %d/%d (%.0f%%)", ct, correct, sum(idx), correct/sum(idx)*100)
  } else {
    sprintf("%s: 0/0", ct)
  }
})
print(acc_by_type)

# 预测准确性柱状图
acc_df <- data.frame(
  celltype = levels(pancreas.query$seurat_annotations),
  accuracy = sapply(levels(pancreas.query$seurat_annotations), function(ct) {
    idx <- pancreas.query$seurat_annotations == ct
    if (sum(idx) > 0) sum(pancreas.query$predicted.id[idx] == ct) / sum(idx) * 100
    else 0
  })
)
p2 <- ggplot(acc_df, aes(x = reorder(celltype, -accuracy), y = accuracy)) +
  geom_bar(stat = "identity", fill = "#4C9BE8") +
  geom_text(aes(label = sprintf("%.0f%%", accuracy)), vjust = -0.5, size = 3.5) +
  ylim(0, 105) +
  labs(title = "各细胞类型预测准确率", x = "细胞类型", y = "准确率 (%)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("02_prediction_accuracy.png", p2)

# 混淆矩阵热图
conf_matrix <- as.data.frame(table(pancreas.query$seurat_annotations, pancreas.query$predicted.id))
colnames(conf_matrix) <- c("True", "Predicted", "Count")
p3 <- ggplot(conf_matrix, aes(x = Predicted, y = True, fill = Count)) +
  geom_tile() +
  geom_text(aes(label = Count), size = 3) +
  scale_fill_gradient(low = "white", high = "#2166AC") +
  labs(title = "混淆矩阵：真实标签 vs 预测标签", x = "预测标签", y = "真实标签") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("03_confusion_matrix.png", p3, width = 10, height = 9)

# ============================================================
# 7. Marker 基因验证
# ============================================================
cat("\n## 7. Marker 基因验证\n\n")

# 用经典 Marker 基因验证预测
markers <- c("MS4A1", "CD3D", "CD3E", "IL7R", "CD14", "LYZ", "GNLY", "NKG7", "FCGR3A", "PPBP")
p4 <- VlnPlot(pancreas.query, features = markers, group.by = "predicted.id",
              ncol = 5, pt.size = 0) +
  plot_layout(ncol = 1) &
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8))
save_plot("04_marker_violin.png", p4, width = 15, height = 12)

cat("Marker 基因小提琴图已保存\n")

# ============================================================
# 8. UMAP 投影 (MapQuery)
# ============================================================
cat("\n## 8. UMAP 投影 (MapQuery)\n\n")

# 重新运行 UMAP 并保存模型
pancreas.ref <- RunUMAP(pancreas.ref, dims = 1:30, reduction = "pca",
                         return.model = TRUE, verbose = FALSE)

# 使用 MapQuery 一步完成标签转移 + 嵌入整合 + UMAP 投影
cat("正在执行 MapQuery...\n")
pancreas.query <- MapQuery(
  anchorset = pancreas.anchors,
  reference = pancreas.ref,
  query = pancreas.query,
  refdata = list(celltype = "seurat_annotations"),
  reference.reduction = "pca",
  reduction.model = "umap"
)
cat("MapQuery 完成\n")

# ============================================================
# 9. 参考集与查询集 UMAP 对比
# ============================================================
cat("\n## 9. 参考集与查询集 UMAP 对比\n\n")

# 参考集 UMAP
p5 <- DimPlot(pancreas.ref, reduction = "umap", group.by = "seurat_annotations",
              label = TRUE, label.size = 3, repel = TRUE) +
  NoLegend() +
  ggtitle("参考集注释")
# 查询集投影 UMAP
p6 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "predicted.celltype",
              label = TRUE, label.size = 3, repel = TRUE) +
  NoLegend() +
  ggtitle("查询集预测标签")
p_combo <- p5 + p6
save_plot("05_ref_vs_query_umap.png", p_combo, width = 16, height = 7)

# 查询集：真实标签 vs 预测标签对比
p7 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "seurat_annotations",
              label = TRUE, label.size = 3, repel = TRUE) +
  ggtitle("查询集真实标签")
p8 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "predicted.celltype",
              label = TRUE, label.size = 3, repel = TRUE) +
  ggtitle("查询集预测标签")
p_truth_vs_pred <- p7 + p8
save_plot("06_truth_vs_predicted.png", p_truth_vs_pred, width = 16, height = 7)

# ============================================================
# 10. 预测置信度分析
# ============================================================
cat("\n## 10. 预测置信度分析\n\n")

# 查看预测分数分布
p9 <- ggplot(pancreas.query@meta.data, aes(x = predicted.id, y = prediction.score.max)) +
  geom_boxplot(fill = "#4C9BE8", alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 0.5) +
  labs(title = "各细胞类型预测置信度", x = "预测细胞类型", y = "预测分数") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("07_prediction_scores.png", p9, width = 10, height = 7)

cat(sprintf("平均预测分数: %.3f\n", mean(pancreas.query$prediction.score.max)))

# ============================================================
# 11. 合并参考集和查询集可视化
# ============================================================
cat("\n## 11. 合并参考集和查询集可视化\n\n")

# 给参考集和查询集添加来源标签
pancreas.ref$dataset <- "Reference"
pancreas.query$dataset <- "Query"

# 合并 UMAP 可视化
p10 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "predicted.celltype",
               label = TRUE, label.size = 3, repel = TRUE) +
  labs(title = "查询集投影到参考 UMAP", color = "预测细胞类型") +
  theme(plot.title = element_text(size = 14))
save_plot("08_query_projected_umap.png", p10, width = 10, height = 8)

# 按数据来源分组着色
p11 <- DimPlot(pancreas.query, reduction = "ref.umap", group.by = "dataset",
               cols = c("Query" = "#E41A1C")) +
  labs(title = "查询集在参考 UMAP 空间中的分布", color = "数据集")
save_plot("09_query_on_ref_space.png", p11, width = 10, height = 8)

# ============================================================
# 12. 保存结果
# ============================================================
cat("\n## 12. 保存结果\n\n")

saveRDS(pancreas.ref, file.path(output_dir, "reference_dataset.rds"))
saveRDS(pancreas.query, file.path(output_dir, "query_mapped.rds"))
cat("已保存: reference_dataset.rds, query_mapped.rds\n")

# 保存预测结果
write.csv(pancreas.query@meta.data[, c("seurat_annotations", "predicted.id", 
                                        "prediction.score.max", "prediction.match")],
          file.path(output_dir, "mapping_predictions.csv"), row.names = TRUE)
cat("已保存: mapping_predictions.csv\n")

cat("\n=== 阶段4-A 参考数据集映射 完成！===\n")
cat(sprintf("共生成 9 张图 + 1 个CSV + 2 个RDS文件\n"))
cat(sprintf("预测准确率: %.1f%%\n", 
    round(match_result["TRUE"] / sum(match_result) * 100, 1)))