20260912 121059 Seurat V5 Scrna Seq 与 Scatac Seq 整合教程

Seurat v5 scRNA-seq 与 scATAC-seq 整合教程 ¶

基于 Seurat v5.5.1 官方教程，使用 PBMC 3k RNA 数据模拟 ATAC 基因活性矩阵。
演示跨模态标签转移和共嵌入（Co-embedding）的核心流程。

目录 ¶

1. 背景与原理
2. RNA 参考数据准备
3. 模拟 ATAC 基因活性数据
4. ATAC 预处理
5. 查找跨模态锚点
6. 标签转移：RNA 到 ATAC
7. 评估预测准确性
8. 共嵌入（Co-embedding）
9. 总结
附录: 完整脚本

1. 背景与原理 ¶

1.1 为什么要整合 scRNA-seq 与 scATAC-seq？¶

单细胞 ATAC-seq 测量染色质可及性，但无法直接识别细胞类型（因为没有基因表达信息）。整合思路：

用已注释的 scRNA-seq 数据作为参考
通过基因活性（Gene Activity）将 ATAC 的峰转化为基因级别的特征
利用跨模态锚点将 RNA 的细胞类型注释转移到 ATAC 细胞

1.2 核心流程 ¶

scRNA-seq（已注释）-> 预处理 -> PCA
                                |
                    FindTransferAnchors(reduction = "cca")
                                |
scATAC-seq -> 基因活性 -> 预处理 -> PCA
                                |
                    TransferData -> ATAC 细胞类型预测
                                |
                    Co-embedding -> RNA + ATAC 共同 UMAP

1.3 与 WNN 的区别 ¶

指标	WNN（阶段4-B）	ATAC 整合（本教程）
数据类型	同一细胞同时测 RNA + ADT	不同细胞分别测 RNA 和 ATAC
模态关系	配对（paired）	非配对（unpaired）
核心方法	FindMultiModalNeighbors	FindTransferAnchors(CCA)
目标	联合聚类	标签转移 + 共嵌入

2. RNA 参考数据准备 ¶

pbmc.rna <- readRDS("output/pbmc3k_final.rds")
pbmc.rna$seurat_annotations <- Idents(pbmc.rna)

# 标准预处理
pbmc.rna <- NormalizeData(pbmc.rna, verbose = FALSE)
pbmc.rna <- FindVariableFeatures(pbmc.rna, verbose = FALSE)
pbmc.rna <- ScaleData(pbmc.rna, verbose = FALSE)
pbmc.rna <- RunPCA(pbmc.rna, verbose = FALSE)
pbmc.rna <- RunUMAP(pbmc.rna, dims = 1:30, verbose = FALSE)

RNA 参考集 UMAP（已知注释）¶

RNA参考集UMAP

3. 模拟 ATAC 基因活性数据 ¶

由于 Signac 包和 pbmcMultiome 数据集下载超时，我们从 RNA 数据模拟 ATAC 基因活性：

# ATAC 基因活性 = 基因体 + 2kb 上游的染色质可及性
# 特点：比 RNA 更稀疏、有更多 dropout、但与 RNA 相关

for (i in seq_along(hv_genes)) {
  gene <- hv_genes[i]
  rna_expr <- as.numeric(rna_counts[gene, ])
  base_signal <- rna_expr * runif(n_cells, 0.3, 0.8)
  dropout_mask <- runif(n_cells) > 0.4  # 40% dropout
  base_signal <- base_signal * dropout_mask
  base_signal <- base_signal + rnorm(n_cells, mean = 0.5, sd = 0.3)
  gene_activity[i, ] <- pmax(round(base_signal * 5), 0)
}

pbmc.atac <- CreateSeuratObject(counts = gene_activity, project = "ATAC")

ATAC vs RNA 数据差异：
- ATAC 基因活性稀疏度更高（约 40% dropout）
- ATAC 信号与 RNA 正相关但噪声更大
- 真实 ATAC 分析需要 Signac 包处理峰（peak）级别数据

4. ATAC 预处理 ¶

# 使用标准 Seurat 流程模拟 ATAC 预处理
# 真实场景中应使用 TF-IDF + LSI（Signac 包）
pbmc.atac <- NormalizeData(pbmc.atac, verbose = FALSE)
pbmc.atac <- FindVariableFeatures(pbmc.atac, verbose = FALSE)
pbmc.atac <- ScaleData(pbmc.atac, verbose = FALSE)
pbmc.atac <- RunPCA(pbmc.atac, verbose = FALSE)
pbmc.atac <- RunUMAP(pbmc.atac, dims = 1:30, reduction.name = "umap.atac")

注意：真实 scATAC-seq 分析应使用 TF-IDF + LSI（Latent Semantic Indexing）代替 LogNormalize + PCA。本教程使用标准 Seurat 流程模拟，核心整合步骤相同。

ATAC UMAP（未注释）¶

ATAC UMAP

RNA vs ATAC 对比 ¶

RNA vs ATAC

ATAC 模拟数据的聚类结构不如 RNA 清晰，这是因为 ATAC 基因活性信号更稀疏、噪声更大。

5. 查找跨模态锚点 ¶

# 使用 CCA 降维来查找跨模态锚点
transfer.anchors <- FindTransferAnchors(
  reference = pbmc.rna,
  query = pbmc.atac,
  features = VariableFeatures(pbmc.rna),
  reference.assay = "RNA",
  query.assay = "RNA",
  reduction = "cca",    # 跨模态使用 CCA
  dims = 1:30
)

CCA（Canonical Correlation Analysis） 在跨模态锚点查找中的作用：
- 寻找两个模态之间的最大相关方向
- 在这些方向上匹配相似的细胞
- 适合 RNA vs ATAC 这种数据类型差异大的场景

6. 标签转移：RNA 到 ATAC ¶

celltype.predictions <- TransferData(
  anchorset = transfer.anchors,
  refdata = pbmc.rna$seurat_annotations,
  weight.reduction = pbmc.atac[["pca"]],
  dims = 1:30
)
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)

预测标签 vs 真实标签 ¶

预测vs真实

由于模拟数据保留了真实标签，可以直接对比预测结果。准确率为 73.9%，低于阶段4-A 的 92.4%，这是因为 ATAC 基因活性信号比 RNA 更稀疏，跨模态转移难度更大。

7. 评估预测准确性 ¶

7.1 混淆矩阵 ¶

混淆矩阵

混淆矩阵显示：
- B 细胞和 CD14+ Mono 预测较好
- T 细胞亚型（Naive CD4 / Memory CD4 / CD8）容易混淆
- 这与 ATAC 数据中 T 细胞亚型基因活性差异较小有关

7.2 预测分数分布 ¶

预测分数

正确预测的细胞平均置信度更高，错误预测的细胞置信度较低，说明置信度可以作为预测可靠性的指标。

8. 共嵌入（Co-embedding）¶

共嵌入将 RNA 和 ATAC 细胞放在同一个 UMAP 空间中，便于联合可视化：

# 用锚点将 RNA 表达数据插补到 ATAC 细胞中
imputation <- TransferData(
  anchorset = transfer.anchors,
  refdata = refdata,
  weight.reduction = pbmc.atac[["pca"]],
  dims = 1:30
)
pbmc.atac[["RNA"]] <- imputation

# 合并两个数据集
coembed <- merge(x = pbmc.rna, y = pbmc.atac)
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use)
coembed <- RunUMAP(coembed, dims = 1:30)

8.1 按模态着色 ¶

共嵌入-模态

RNA 和 ATAC 细胞在共嵌入空间中充分混合，说明跨模态整合成功。

8.2 按细胞类型着色 ¶

共嵌入-细胞类型

不同细胞类型的 RNA 和 ATAC 细胞在同一区域聚集，验证了标签转移的合理性。

8.3 按模态分面 ¶

共嵌入-分面

9. 总结 ¶

核心流程 ¶

scRNA-seq（已注释）                  scATAC-seq（未注释）
      |                                    |
  预处理 + PCA                      基因活性 -> 预处理 + PCA
      |                                    |
      +---------- FindTransferAnchors(CCA) --------+
                        |
                  TransferData -> ATAC 细胞类型预测
                        |
                  Co-embedding -> 联合 UMAP 可视化

关键知识点 ¶

基因活性是连接 ATAC 和 RNA 的桥梁：将 ATAC 峰映射到基因级别
CCA 降维适合跨模态锚点查找，能处理数据类型差异大的场景
标签转移准确率（73.9%）低于同模态映射（92.4%），因为 ATAC 信号更稀疏
共嵌入将不同模态的细胞放在同一空间，便于联合分析和可视化
真实 ATAC 分析需要 Signac 包进行峰级别处理和 TF-IDF + LSI 降维

附录: 完整脚本 ¶

atac_integration_vignette.R ¶

# ============================================================
# Seurat v5 阶段4-C：scRNA-seq 与 scATAC-seq 整合
# ============================================================
# 基于官方教程：https://satijalab.org/seurat/articles/seurat5_atacseq_integration_vignette
# 由于 Signac 包和 pbmcMultiome 数据集下载超时，
# 使用 PBMC 3k RNA 数据模拟 ATAC 基因活性矩阵，
# 演示跨模态标签转移和共嵌入的核心流程
# ============================================================

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

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

cat("=== Seurat v5 阶段4-C：scRNA-seq 与 scATAC-seq 整合 ===\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. 加载并准备 RNA 参考数据
# ============================================================
cat("## 1. 加载 RNA 参考数据\n\n")

pbmc.rna <- readRDS("output/pbmc3k_final.rds")
pbmc.rna$seurat_annotations <- Idents(pbmc.rna)
cat(sprintf("RNA 参考集：%d 细胞，%d 基因\n", ncol(pbmc.rna), nrow(pbmc.rna)))
cat(sprintf("细胞类型：%s\n", paste(levels(pbmc.rna$seurat_annotations), collapse = ", ")))

# RNA 标准预处理
pbmc.rna <- NormalizeData(pbmc.rna, verbose = FALSE)
pbmc.rna <- FindVariableFeatures(pbmc.rna, verbose = FALSE)
pbmc.rna <- ScaleData(pbmc.rna, verbose = FALSE)
pbmc.rna <- RunPCA(pbmc.rna, verbose = FALSE)
pbmc.rna <- RunUMAP(pbmc.rna, dims = 1:30, verbose = FALSE)
cat("RNA 预处理完成\n")

# RNA UMAP
p1 <- DimPlot(pbmc.rna, group.by = "seurat_annotations", label = TRUE, repel = TRUE) + 
  NoLegend() + ggtitle("RNA 参考集（已知注释）")
save_plot("01_rna_reference_umap.png", p1)

# ============================================================
# 2. 模拟 ATAC 基因活性数据
# ============================================================
cat("\n## 2. 模拟 ATAC 基因活性数据\n\n")

# 提取 RNA counts 矩阵
rna_counts <- GetAssayData(pbmc.rna, assay = "RNA", layer = "counts")

# 选择高变基因作为基因活性基础
hv_genes <- VariableFeatures(pbmc.rna)
# 只取前 2000 个高变基因来模拟
hv_genes <- hv_genes[1:min(2000, length(hv_genes))]
hv_genes <- hv_genes[hv_genes %in% rownames(rna_counts)]
cat(sprintf("使用 %d 个高变基因模拟 ATAC 基因活性\n", length(hv_genes)))

# 模拟 ATAC 基因活性矩阵
# ATAC 基因活性 = 基因体 + 2kb 上游的染色质可及性
# 特点：比 RNA 更稀疏、有更多 dropout、但与 RNA 相关
set.seed(123)
n_cells <- ncol(pbmc.rna)
n_genes <- length(hv_genes)

# 从 RNA 表达模拟 ATAC 基因活性
# 添加 dropout（50% 的零值）和噪声
gene_activity <- matrix(0, nrow = n_genes, ncol = n_cells)
rownames(gene_activity) <- hv_genes
colnames(gene_activity) <- colnames(pbmc.rna)

for (i in seq_along(hv_genes)) {
  gene <- hv_genes[i]
  rna_expr <- as.numeric(rna_counts[gene, ])
  # 模拟 ATAC：与 RNA 相关但更稀疏
  # 1. 基础信号 = RNA 表达 * 随机系数
  base_signal <- rna_expr * runif(n_cells, 0.3, 0.8)
  # 2. 加入 dropout（随机置零）
  dropout_mask <- runif(n_cells) > 0.4  # 40% dropout
  base_signal <- base_signal * dropout_mask
  # 3. 加入随机噪声
  base_signal <- base_signal + rnorm(n_cells, mean = 0.5, sd = 0.3)
  # 4. 确保非负
  gene_activity[i, ] <- pmax(round(base_signal * 5), 0)
}

gene_activity <- Matrix(gene_activity, sparse = TRUE)

# 创建 ATAC Seurat 对象（使用基因活性矩阵）
pbmc.atac <- CreateSeuratObject(counts = gene_activity, project = "ATAC")
# 添加真实标签（用于验证预测准确性）
pbmc.atac$seurat_annotations <- pbmc.rna$seurat_annotations[match(colnames(pbmc.atac), colnames(pbmc.rna))]

cat(sprintf("ATAC 模拟数据：%d 细胞，%d 基因活性特征\n", ncol(pbmc.atac), nrow(pbmc.atac)))
cat(sprintf("稀疏度：%.1f%%\n", round(sum(gene_activity == 0) / length(gene_activity) * 100, 1)))

# ============================================================
# 3. ATAC 预处理（TF-IDF + LSI 模拟）
# ============================================================
cat("\n## 3. ATAC 预处理\n\n")

# 由于没有 Signac，使用标准 Seurat 流程模拟 ATAC 预处理
# TF-IDF 近似：先 LogNormalize 再做 PCA（模拟 LSI）
DefaultAssay(pbmc.atac) <- "RNA"
pbmc.atac <- NormalizeData(pbmc.atac, verbose = FALSE)
pbmc.atac <- FindVariableFeatures(pbmc.atac, verbose = FALSE)
pbmc.atac <- ScaleData(pbmc.atac, verbose = FALSE)
pbmc.atac <- RunPCA(pbmc.atac, verbose = FALSE)

# 使用 PCA 近似 LSI
# 在真实场景中，这里使用 RunSVD（LSI），但 PCA 也能演示流程
pbmc.atac <- RunUMAP(pbmc.atac, dims = 1:30, reduction.name = "umap.atac",
                     reduction.key = "atacUMAP_", verbose = FALSE)
cat("ATAC 预处理完成\n")

# ATAC UMAP（未注释）
p2 <- DimPlot(pbmc.atac, reduction = "umap.atac", group.by = "orig.ident") + 
  NoLegend() + ggtitle("ATAC 模拟数据 UMAP（未注释）")
save_plot("02_atac_umap_unannotated.png", p2)

# RNA vs ATAC 对比
p3 <- DimPlot(pbmc.rna, group.by = "seurat_annotations", label = TRUE, repel = TRUE) + 
  NoLegend() + ggtitle("RNA")
p4 <- DimPlot(pbmc.atac, reduction = "umap.atac", group.by = "seurat_annotations", 
              label = TRUE, repel = TRUE) + NoLegend() + ggtitle("ATAC（真实标签）")
save_plot("03_rna_vs_atac.png", p3 + p4, width = 16, height = 7)

# ============================================================
# 4. 查找跨模态锚点
# ============================================================
cat("\n## 4. 查找跨模态锚点\n\n")

# 使用 CCA 降维来查找跨模态锚点
# 将 ATAC 的基因活性作为 query，RNA 表达作为 reference
cat("正在查找跨模态锚点（CCA）...\n")
transfer.anchors <- FindTransferAnchors(
  reference = pbmc.rna,
  query = pbmc.atac,
  features = VariableFeatures(pbmc.rna),
  reference.assay = "RNA",
  query.assay = "RNA",
  reduction = "cca",
  dims = 1:30
)
cat("锚点查找完成\n")

# ============================================================
# 5. 标签转移：从 RNA 注释到 ATAC
# ============================================================
cat("\n## 5. 标签转移：RNA → ATAC\n\n")

# 转移细胞类型标签
celltype.predictions <- TransferData(
  anchorset = transfer.anchors,
  refdata = pbmc.rna$seurat_annotations,
  weight.reduction = pbmc.atac[["pca"]],
  dims = 1:30
)
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)
cat("标签转移完成\n")

# 预测结果
cat(sprintf("\n预测细胞类型分布:\n"))
print(table(pbmc.atac$predicted.id))

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

pbmc.atac$annotation_correct <- pbmc.atac$predicted.id == pbmc.atac$seurat_annotations
correct <- sum(pbmc.atac$annotation_correct)
total <- ncol(pbmc.atac)
cat(sprintf("正确预测：%d/%d (%.1f%%)\n", correct, total, correct / total * 100))

# 预测 vs 真实标签 UMAP 对比
p5 <- DimPlot(pbmc.atac, reduction = "umap.atac", group.by = "predicted.id",
              label = TRUE, repel = TRUE) + NoLegend() + ggtitle("ATAC 预测标签")
p6 <- DimPlot(pbmc.atac, reduction = "umap.atac", group.by = "seurat_annotations",
              label = TRUE, repel = TRUE) + NoLegend() + ggtitle("ATAC 真实标签")
save_plot("04_predicted_vs_truth.png", p5 + p6, width = 16, height = 7)

# 混淆矩阵
conf_matrix <- as.data.frame(table(pbmc.atac$seurat_annotations, pbmc.atac$predicted.id))
conf_matrix <- conf_matrix[conf_matrix$Freq > 0, ]
colnames(conf_matrix) <- c("True", "Predicted", "Count")
# 归一化
conf_matrix$Freq <- conf_matrix$Count / ave(conf_matrix$Count, conf_matrix$True, FUN = sum)

p7 <- ggplot(conf_matrix, aes(x = Predicted, y = True, fill = Freq)) +
  geom_tile() +
  scale_fill_gradient(name = "比例", low = "#ffffc8", high = "#7d0025") +
  labs(title = "混淆矩阵：RNA 注释 vs ATAC 预测",
       x = "预测标签（ATAC）", y = "真实标签（RNA）") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot("05_confusion_matrix.png", p7, width = 10, height = 9)

# 预测分数分布
p8 <- ggplot(pbmc.atac@meta.data, aes(x = prediction.score.max, fill = annotation_correct)) +
  geom_density(alpha = 0.5) +
  labs(title = "预测分数分布",
       x = "预测分数", y = "密度", fill = "注释正确") +
  theme_bw()
save_plot("06_prediction_scores.png", p8, width = 10, height = 7)

cat(sprintf("平均预测分数：%.3f\n", mean(pbmc.atac$prediction.score.max)))
cat(sprintf("正确预测平均分数：%.3f\n", 
    mean(pbmc.atac$prediction.score.max[pbmc.atac$annotation_correct])))
cat(sprintf("错误预测平均分数：%.3f\n", 
    mean(pbmc.atac$prediction.score.max[!pbmc.atac$annotation_correct])))

# ============================================================
# 7. 共嵌入 (Co-embedding)
# ============================================================
cat("\n## 7. 共嵌入 RNA + ATAC\n\n")

# 用锚点将 RNA 表达数据插补到 ATAC 细胞中
genes.use <- VariableFeatures(pbmc.rna)
refdata <- GetAssayData(pbmc.rna, assay = "RNA", layer = "data")[genes.use, ]

cat("正在插补 RNA 表达到 ATAC 细胞...\n")
imputation <- TransferData(
  anchorset = transfer.anchors,
  refdata = refdata,
  weight.reduction = pbmc.atac[["pca"]],
  dims = 1:30
)
pbmc.atac[["RNA"]] <- imputation
cat("RNA 插补完成\n")

# 合并两个数据集
# 给每个对象添加来源标签
pbmc.rna$modality <- "RNA"
pbmc.atac$modality <- "ATAC"

coembed <- merge(x = pbmc.rna, y = pbmc.atac)
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE, verbose = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30, verbose = FALSE)
cat("共嵌入完成\n")

# 共嵌入 UMAP
p9 <- DimPlot(coembed, group.by = "modality", 
              cols = c("RNA" = "#E41A1C", "ATAC" = "#377EB8")) +
  ggtitle("共嵌入：RNA + ATAC（按模态着色）")
save_plot("07_coembed_modality.png", p9, width = 10, height = 8)

p10 <- DimPlot(coembed, group.by = "seurat_annotations", label = TRUE, 
               repel = TRUE, label.size = 3) + NoLegend() +
  ggtitle("共嵌入：RNA + ATAC（按细胞类型着色）")
save_plot("08_coembed_celltype.png", p10, width = 10, height = 8)

# 按模态拆分
p11 <- DimPlot(coembed, group.by = "modality", split.by = "modality",
               cols = c("RNA" = "#E41A1C", "ATAC" = "#377EB8")) +
  ggtitle("共嵌入：按模态分面")
save_plot("09_coembed_split.png", p11, width = 16, height = 7)

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

saveRDS(pbmc.atac, file.path(output_dir, "atac_annotated.rds"))
saveRDS(coembed, file.path(output_dir, "coembed.rds"))
cat("已保存: atac_annotated.rds, coembed.rds\n")

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

cat("\n=== 阶段4-C scRNA-seq 与 scATAC-seq 整合 完成！===\n")
cat(sprintf("共生成 9 张图 + 1 个CSV + 2 个RDS文件\n"))
cat(sprintf("标签转移准确率：%.1f%%\n", correct / total * 100))

本教程基于 Seurat v5.5.1 官方文档整理
R 4.6.1 + Seurat 5.5.1 + Windows 环境
阶段4-C: scRNA-seq 与 scATAC-seq 整合（跨模态锚点 / 标签转移准确率 73.9% / 共嵌入）