20260912 111444 Seurat V5 数据整合分析教程

# 读取原始 10X 数据
pbmc_data <- Read10X(data.dir = data_dir)
pbmc <- CreateSeuratObject(counts = pbmc_data, project = "pbmc3k",
                            min.cells = 3, min.features = 200)

# QC 过滤
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
# 过滤后细胞数: 2638

# 随机分为两个批次
set.seed(42)
pbmc$batch <- sample(c("batch1", "batch2"), size = ncol(pbmc), replace = TRUE)
# batch1: 1348 个细胞 | batch2: 1290 个细胞

# 为 batch2 添加技术偏差：
# 1) 乘法噪声（0.7-1.3 的基因特异噪声因子，模拟不同捕获效率）
# 2) Dropout 500 个基因的 30% 细胞（模拟检测率差异）

# 按 batch 分割 RNA assay 的 counts layer
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)
# 分割后 pbmc 对象包含 counts.batch1 和 counts.batch2 两个 layer

pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)
pbmc <- FindNeighbors(pbmc, dims = 1:30, reduction = "pca")
pbmc <- FindClusters(pbmc, resolution = 0.8, cluster.name = "unintegrated_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

pbmc <- IntegrateLayers(
    object = pbmc,
    method = CCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.cca",
    verbose = FALSE
)

# 整合后必须合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# 基于整合后的降维结果聚类
pbmc <- FindNeighbors(pbmc, reduction = "integrated.cca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "cca_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.cca", reduction.name = "umap.cca")

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 5, order_by = avg_log2FC)

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
                     "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
pbmc$cell_type <- Idents(pbmc)

FeaturePlot(pbmc, features = c("IL7R", "CD14", "MS4A1", "CD8A", "GNLY", "FCGR3A"),
            reduction = "umap.cca", ncol = 3)

# 注意：FindConservedMarkers 需要 layers 处于合并状态
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

Idents(pbmc) <- "cca_clusters"
# 选择最大的聚类进行分析
conserved_markers <- FindConservedMarkers(
    pbmc,
    ident.1 = target_cluster,  # 聚类1
    grouping.var = "batch",     # 按批次分组
    verbose = FALSE
)

VlnPlot(pbmc, features = top_conserved, split.by = "batch", pt.size = 0, ncol = 3)

DotPlot(pbmc, features = markers_to_plot, split.by = "batch") + RotatedAxis()

# 创建 celltype.batch 组合变量
pbmc$celltype.batch <- paste(pbmc$cell_type, pbmc$batch, sep = "_")
Idents(pbmc) <- "celltype.batch"

# 比较第一个共有细胞类型的两个批次
de_genes <- FindMarkers(pbmc, ident.1 = "Cluster0_batch2", ident.2 = "Cluster0_batch1")

VlnPlot(pbmc, features = top_de, split.by = "batch", pt.size = 0, ncol = 3)

# RPCA 需要 layers 处于分割状态
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)

pbmc <- IntegrateLayers(
    object = pbmc,
    method = RPCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.rpca",
    verbose = FALSE
)

# 合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

原始数据 → QC 过滤 → 模拟/识别批次
                           ↓
                    split(pbmc[["RNA"]], f = pbmc$batch)
                           ↓
              ┌─────────────┴──────────────┐
              ↓                           ↓
        未整合分析                    CCA/RPCA 整合
        (展示批次效应)              IntegrateLayers()
              ↓                           ↓
        UMAP（批次分离）           JoinLayers() + 聚类
                                      ↓
                              UMAP（批次混合）
                                      ↓
                    ┌─────────────────┴──────────────────┐
                    ↓                                    ↓
           FindConservedMarkers                    FindMarkers
           (跨批次保守 marker)                   (跨批次差异基因)

###############################################################################
# Seurat v5 阶段3：数据整合分析（Integration）
# 基于官方教程:
#   - https://satijalab.org/seurat/articles/integration_introduction
#   - https://satijalab.org/seurat/articles/seurat5_integration
#
# 由于 ifnb 数据集（394MB）从 Seurat 官方服务器下载过慢，
# 本脚本使用已有的 PBMC 3k 数据，通过人为分割模拟两个"批次"，
# 完整演示 Seurat v5 整合流程的所有核心步骤：
#   1. 模拟批次效应
#   2. 未整合分析（展示批次效应）
#   3. CCA 整合
#   4. 整合前后对比
#   5. 保守 Marker 基因（FindConservedMarkers）
#   6. 跨条件差异表达（FindMarkers）
#   7. 整合方法对比（CCA vs RPCA）
###############################################################################

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

# === 路径配置 ===
base_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial"
data_dir <- file.path(base_dir, "data/filtered_gene_bc_matrices/hg19")
output_dir <- file.path(base_dir, "output")
integ_dir <- file.path(output_dir, "integration")
dir.create(integ_dir, showWarnings = FALSE, recursive = TRUE)

# ============================================================================
# 1. 加载数据 & 模拟批次效应
# ============================================================================
cat("\n========== 1. 加载数据 & 模拟批次 ==========\n")

# 读取原始 10X 数据
pbmc_data <- Read10X(data.dir = data_dir)

# 创建 Seurat 对象
pbmc <- CreateSeuratObject(counts = pbmc_data, project = "pbmc3k",
                            min.cells = 3, min.features = 200)

# 计算线粒体百分比
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# QC 过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("过滤后细胞数:", ncol(pbmc), "\n")

# --- 模拟批次效应 ---
# 将细胞随机分为两个"批次"
set.seed(42)
cell_ids <- colnames(pbmc)
batch_labels <- sample(c("batch1", "batch2"), size = length(cell_ids), replace = TRUE)
names(batch_labels) <- cell_ids
pbmc$batch <- batch_labels

# 为 batch2 人为添加技术偏差：
# 1) 对 batch2 的部分基因表达量添加随机噪声（模拟批次效应）
# 2) 降低 batch2 中某些基因的检测率
cat("正在模拟批次效应...\n")

# 获取 counts 矩阵
counts_mat <- GetAssayData(pbmc, layer = "counts")

# 找到 batch2 的细胞
batch2_cells <- which(pbmc$batch == "batch2")
batch1_cells <- which(pbmc$batch == "batch1")
cat("batch1:", length(batch1_cells), "个细胞\n")
cat("batch2:", length(batch2_cells), "个细胞\n")

# 对 batch2 添加乘法噪声（模拟不同的捕获效率）
set.seed(123)
noise_factor <- 0.7 + 0.6 * runif(nrow(counts_mat))  # 0.7-1.3 的基因特异噪声
noise_matrix <- matrix(noise_factor, nrow = nrow(counts_mat), ncol = length(batch2_cells))
# 对 batch2 的 counts 乘以噪声因子（向下取整）
counts_batch2 <- ceiling(as.matrix(counts_mat[, batch2_cells]) * noise_matrix)
counts_mat[, batch2_cells] <- as(counts_batch2, "dgCMatrix")

# 为 batch2 随机 dropout 部分基因（模拟检测率差异）
set.seed(456)
dropout_genes <- sample(rownames(pbmc), size = 500)  # 随机选 500 个基因
for (gene in dropout_genes) {
    dropout_cells <- sample(batch2_cells, size = floor(length(batch2_cells) * 0.3))
    counts_mat[gene, dropout_cells] <- 0
}

# 写回修改后的 counts
pbmc[["RNA"]]$counts <- counts_mat
cat("批次效应模拟完成\n")

# ============================================================================
# 2. 按批次分割 Layers（Seurat v5 核心特性）
# ============================================================================
cat("\n========== 2. 分割 Layers ==========\n")

# Seurat v5 核心特性：在一个对象中通过 layers 管理多个批次
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)
cat("分割后 layers:\n")
print(pbmc)

# ============================================================================
# 3. 未整合分析（展示批次效应）
# ============================================================================
cat("\n========== 3. 未整合分析 ==========\n")

# 标准流程（对每个 layer 独立归一化和找高变基因）
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)

# 聚类
pbmc <- FindNeighbors(pbmc, dims = 1:30, reduction = "pca")
pbmc <- FindClusters(pbmc, resolution = 0.8, cluster.name = "unintegrated_clusters")

# UMAP（未整合）
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")

# 3.1 未整合 UMAP - 按批次着色
p_unintegrated_batch <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "batch") +
    ggtitle("未整合 - 按批次着色")
ggsave(file.path(integ_dir, "01_unintegrated_batch.png"), plot = p_unintegrated_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 01_unintegrated_batch.png\n")

# 3.2 未整合 UMAP - 按聚类着色
p_unintegrated_cluster <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "unintegrated_clusters") +
    ggtitle("未整合 - 按聚类着色")
ggsave(file.path(integ_dir, "02_unintegrated_cluster.png"), plot = p_unintegrated_cluster, width = 8, height = 7, dpi = 150)
cat("已保存: 02_unintegrated_cluster.png\n")

# 3.3 未整合 - split.by 批次
p_unintegrated_split <- DimPlot(pbmc, reduction = "umap.unintegrated", split.by = "batch") +
    ggtitle("未整合 - 按批次分割")
ggsave(file.path(integ_dir, "03_unintegrated_split.png"), plot = p_unintegrated_split, width = 14, height = 6, dpi = 150)
cat("已保存: 03_unintegrated_split.png\n")

# ============================================================================
# 4. CCA 整合
# ============================================================================
cat("\n========== 4. CCA 整合 ==========\n")
cat("正在运行 CCA 整合（可能需要几分钟）...\n")

# IntegrateLayers 一行代码完成整合
pbmc <- IntegrateLayers(
    object = pbmc,
    method = CCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.cca",
    verbose = FALSE
)
cat("CCA 整合完成\n")

# 重新合并 layers（整合后必须执行）
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])
cat("Layers 已重新合并\n")

# 基于整合后的降维结果聚类
pbmc <- FindNeighbors(pbmc, reduction = "integrated.cca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "cca_clusters")

# 整合后 UMAP
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.cca", reduction.name = "umap.cca")

# 4.1 整合后 UMAP - 按批次着色
p_integrated_batch <- DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合后 - 按批次着色")
ggsave(file.path(integ_dir, "04_integrated_batch.png"), plot = p_integrated_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 04_integrated_batch.png\n")

# 4.2 整合后 UMAP - 按聚类着色
p_integrated_cluster <- DimPlot(pbmc, reduction = "umap.cca", group.by = "cca_clusters", label = TRUE) +
    ggtitle("CCA 整合后 - 按聚类着色")
ggsave(file.path(integ_dir, "05_integrated_cluster.png"), plot = p_integrated_cluster, width = 8, height = 7, dpi = 150)
cat("已保存: 05_integrated_cluster.png\n")

# 4.3 整合后 - split.by 批次
p_integrated_split <- DimPlot(pbmc, reduction = "umap.cca", split.by = "batch") +
    ggtitle("CCA 整合后 - 按批次分割")
ggsave(file.path(integ_dir, "06_integrated_split.png"), plot = p_integrated_split, width = 14, height = 6, dpi = 150)
cat("已保存: 06_integrated_split.png\n")

# ============================================================================
# 5. 整合前后对比图
# ============================================================================
cat("\n========== 5. 整合前后对比 ==========\n")

# 5.1 并排对比 - 按批次着色
p_compare_batch <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "batch") +
    ggtitle("未整合") +
DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合后")
ggsave(file.path(integ_dir, "07_compare_batch.png"), plot = p_compare_batch, width = 16, height = 6, dpi = 150)
cat("已保存: 07_compare_batch.png\n")

# 5.2 并排对比 - 按聚类着色
p_compare_cluster <- DimPlot(pbmc, reduction = "umap.unintegrated", group.by = "unintegrated_clusters") +
    ggtitle("未整合聚类") +
DimPlot(pbmc, reduction = "umap.cca", group.by = "cca_clusters", label = TRUE) +
    ggtitle("整合后聚类")
ggsave(file.path(integ_dir, "08_compare_cluster.png"), plot = p_compare_cluster, width = 16, height = 6, dpi = 150)
cat("已保存: 08_compare_cluster.png\n")

# ============================================================================
# 6. 细胞类型注释
# ============================================================================
cat("\n========== 6. 细胞类型注释 ==========\n")

# 使用经典 marker 进行注释（与阶段1一致）
# 先寻找每个聚类的 marker 基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top_markers <- pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 5, order_by = avg_log2FC)
cat("各聚类 top 5 marker:\n")
print(top_markers)

# 保存 marker 表
write.csv(pbmc.markers, file.path(integ_dir, "integration_markers.csv"), row.names = FALSE)
cat("已保存: integration_markers.csv\n")

# 基于阶段1的知识进行注释
# 查看各聚类的 marker 来确定细胞类型
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B",
                      "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
# 根据实际聚类数调整
n_clusters <- length(unique(pbmc$cca_clusters))
cat("聚类数量:", n_clusters, "\n")

# 如果聚类数与注释数不匹配，使用通用注释
if (n_clusters == length(new.cluster.ids)) {
    names(new.cluster.ids) <- levels(pbmc)
    pbmc <- RenameIdents(pbmc, new.cluster.ids)
} else {
    cat("聚类数与预设注释数不一致，使用通用注释\n")
    # 使用数字标签
    new.cluster.ids <- paste0("Cluster", 0:(n_clusters - 1))
    names(new.cluster.ids) <- levels(pbmc)
    pbmc <- RenameIdents(pbmc, new.cluster.ids)
}
pbmc$cell_type <- Idents(pbmc)

# 6.1 注释后 UMAP
p_annotated <- DimPlot(pbmc, reduction = "umap.cca", label = TRUE, pt.size = 0.5) + NoLegend() +
    ggtitle("整合后 - 细胞类型注释")
ggsave(file.path(integ_dir, "09_annotated_umap.png"), plot = p_annotated, width = 8, height = 7, dpi = 150)
cat("已保存: 09_annotated_umap.png\n")

# 6.2 Marker FeaturePlot
p_markers_fp <- FeaturePlot(pbmc, features = c("IL7R", "CD14", "MS4A1", "CD8A", "GNLY", "FCGR3A"),
                             reduction = "umap.cca", ncol = 3)
ggsave(file.path(integ_dir, "10_marker_featureplot.png"), plot = p_markers_fp, width = 14, height = 10, dpi = 150)
cat("已保存: 10_marker_featureplot.png\n")

# ============================================================================
# 7. 保守 Marker 基因（FindConservedMarkers）
# ============================================================================
cat("\n========== 7. 保守 Marker 基因 ==========\n")

# FindConservedMarkers：在两个批次中分别寻找 marker，然后用 meta-analysis 合并 p 值
# 这可以找到不受批次效应影响的稳定 marker 基因

# 注意：FindConservedMarkers 需要 layers 处于合并状态
# （它在内部按 grouping.var 分组分析，不需要预先 split）
# 确保 layers 已合并
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# 以某个聚类为例（选择最大的聚类）
Idents(pbmc) <- "cca_clusters"
target_cluster <- as.character(which.max(table(pbmc$cca_clusters)))
cat("分析聚类", target_cluster, "的保守 marker\n")

conserved_markers <- FindConservedMarkers(
    pbmc,
    ident.1 = target_cluster,
    grouping.var = "batch",
    verbose = FALSE
)

cat("保守 marker 前 15 个:\n")
print(head(conserved_markers, 15))

# 保存保守 marker
write.csv(conserved_markers, file.path(integ_dir, "conserved_markers.csv"), row.names = TRUE)
cat("已保存: conserved_markers.csv\n")

# 7.1 保守 marker 小提琴图（按批次分割）
top_conserved <- rownames(head(conserved_markers, 6))
p_conserved_vln <- VlnPlot(pbmc, features = top_conserved, split.by = "batch", pt.size = 0, ncol = 3) +
    ggtitle("保守 Marker - 按批次分割")
ggsave(file.path(integ_dir, "11_conserved_markers_violin.png"), plot = p_conserved_vln, width = 14, height = 10, dpi = 150)
cat("已保存: 11_conserved_markers_violin.png\n")

# 7.2 保守 marker DotPlot（按批次分割）
# layers 已合并，直接使用
markers_to_plot <- rownames(head(conserved_markers, 12))
p_conserved_dot <- DotPlot(pbmc, features = markers_to_plot, split.by = "batch") + RotatedAxis() +
    ggtitle("保守 Marker DotPlot - 按批次对比")
ggsave(file.path(integ_dir, "12_conserved_markers_dot.png"), plot = p_conserved_dot, width = 12, height = 6, dpi = 150)
cat("已保存: 12_conserved_markers_dot.png\n")

# ============================================================================
# 8. 跨批次差异表达分析
# ============================================================================
cat("\n========== 8. 跨批次差异表达 ==========\n")

# 创建 celltype.batch 组合变量
pbmc$celltype.batch <- paste(pbmc$cell_type, pbmc$batch, sep = "_")
Idents(pbmc) <- "celltype.batch"

# 选择一个细胞类型，比较两个批次的差异
# 尝试找到 CD14+ Mono 或第一个细胞类型
unique_idents <- unique(pbmc$celltype.batch)
cat("可用的 celltype.batch 组合:\n")
print(table(pbmc$celltype.batch))

# 选择第一个细胞类型进行跨批次比较
# 找到同时存在于两个批次的细胞类型
batch1_types <- unique(pbmc$cell_type[pbmc$batch == "batch1"])
batch2_types <- unique(pbmc$cell_type[pbmc$batch == "batch2"])
common_types <- intersect(batch1_types, batch2_types)
cat("两个批次共有的细胞类型:\n")
print(common_types)

if (length(common_types) > 0) {
    target_type <- common_types[1]
    ident1 <- paste0(target_type, "_batch2")
    ident2 <- paste0(target_type, "_batch1")

    if (ident1 %in% unique_idents && ident2 %in% unique_idents) {
        de_genes <- FindMarkers(pbmc, ident.1 = ident1, ident.2 = ident2, verbose = FALSE)
        cat("\n", target_type, "的跨批次差异基因（前 15 个）:\n")
        print(head(de_genes, 15))

        write.csv(de_genes, file.path(integ_dir, "batch_de_genes.csv"), row.names = TRUE)
        cat("已保存: batch_de_genes.csv\n")

        # 8.1 差异基因小提琴图
        top_de <- rownames(head(de_genes, 6))
        if (length(top_de) > 0) {
            p_de_vln <- VlnPlot(pbmc, features = top_de, split.by = "batch", pt.size = 0, ncol = 3) +
                ggtitle(paste0("跨批次差异基因 - ", target_type))
            ggsave(file.path(integ_dir, "13_batch_de_violin.png"), plot = p_de_vln, width = 14, height = 10, dpi = 150)
            cat("已保存: 13_batch_de_violin.png\n")
        }
    }
}

# ============================================================================
# 9. RPCA 整合（更快速、更保守的整合方法）
# ============================================================================
cat("\n========== 9. RPCA 整合 ==========\n")

# RPCA 整合需要按批次分割 layers
pbmc[["RNA"]] <- split(pbmc[["RNA"]], f = pbmc$batch)

# RPCA 整合（比 CCA 更快、更保守）
cat("正在运行 RPCA 整合...\n")
pbmc <- IntegrateLayers(
    object = pbmc,
    method = RPCAIntegration,
    orig.reduction = "pca",
    new.reduction = "integrated.rpca",
    verbose = FALSE
)
cat("RPCA 整合完成\n")

# 合并 layers
pbmc[["RNA"]] <- JoinLayers(pbmc[["RNA"]])

# RPCA 聚类和 UMAP
pbmc <- FindNeighbors(pbmc, reduction = "integrated.rpca", dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5, cluster.name = "rpca_clusters")
pbmc <- RunUMAP(pbmc, dims = 1:30, reduction = "integrated.rpca", reduction.name = "umap.rpca")

# 9.1 RPCA 整合后 UMAP - 按批次着色
p_rpca_batch <- DimPlot(pbmc, reduction = "umap.rpca", group.by = "batch") +
    ggtitle("RPCA 整合后 - 按批次着色")
ggsave(file.path(integ_dir, "14_rpca_batch.png"), plot = p_rpca_batch, width = 8, height = 7, dpi = 150)
cat("已保存: 14_rpca_batch.png\n")

# 9.2 CCA vs RPCA 对比
p_cca_rpca <- DimPlot(pbmc, reduction = "umap.cca", group.by = "batch") +
    ggtitle("CCA 整合") +
DimPlot(pbmc, reduction = "umap.rpca", group.by = "batch") +
    ggtitle("RPCA 整合")
ggsave(file.path(integ_dir, "15_cca_vs_rpca.png"), plot = p_cca_rpca, width = 16, height = 6, dpi = 150)
cat("已保存: 15_cca_vs_rpca.png\n")

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

saveRDS(pbmc, file.path(output_dir, "pbmc3k_integrated.rds"))
cat("已保存整合后对象: pbmc3k_integrated.rds\n")

# ============================================================================
# 总结
# ============================================================================
cat("\n")
cat("=======================================================\n")
cat("  阶段3 数据整合教程完成!\n")
cat("=======================================================\n")
cat("  输出目录:", integ_dir, "\n")
cat("  图片数量: 15 张\n")
cat("  覆盖内容:\n")
cat("    01-03: 未整合分析（展示批次效应）\n")
cat("    04-06: CCA 整合结果\n")
cat("    07-08: 整合前后对比\n")
cat("    09-10: 细胞类型注释\n")
cat("    11-12: 保守 Marker 基因\n")
cat("    13:    跨批次差异表达\n")
cat("    14-15: RPCA 整合与 CCA 对比\n")
cat("  核心知识点:\n")
cat("    - Seurat v5 layers 机制（split / Join）\n")
cat("    - IntegrateLayers 一行代码整合\n")
cat("    - CCA vs RPCA 两种整合方法\n")
cat("    - FindConservedMarkers 跨批次保守 marker\n")
cat("    - FindMarkers 跨批次差异表达\n")
cat("    - 整合前后对比可视化\n")
cat("=======================================================\n")