20260912 091245 Seurat V5 Pbmc 3K 单细胞分析完整教程

读取数据 -> QC过滤 -> 归一化 -> 高变基因 -> 标准化 -> PCA
-> 确定维度 -> 聚类 -> UMAP -> 寻找Marker -> 细胞类型注释

# set a CRAN mirror (TUNA Tsinghua)
local({r <- getOption("repos")
      r["CRAN"] <- "https://mirrors.tuna.tsinghua.edu.cn/CRAN"
      options(repos=r)})

install.packages("Seurat")

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

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

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

# 创建 Seurat 对象
# min.cells = 3: 至少在3个细胞中检测到的基因才保留
# min.features = 200: 每个细胞至少检测到200个基因才保留
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200)

        AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 ...
CD3D    4                .                10               ...
TCL1A   .                .                .                ...
MS4A1   .                6                .                ...

# 计算线粒体基因表达比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 过滤条件:
#   nFeature_RNA > 200   (去除空液滴)
#   nFeature_RNA < 2500  (去除 doublet)
#   percent.mt < 5       (去除线粒体污染严重的细胞)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

# LogNormalize: 每个细胞表达量除以总UMI，乘10000，再取log
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# 默认选择 2000 个高变基因，用于下游 PCA
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 对每个基因做 Z-score 标准化（均值=0, 方差=1）
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

# 对高变基因做 PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# ElbowPlot: 按方差解释比例排序各 PC，拐点处即为建议维度
ElbowPlot(pbmc)

# 基于 PCA 空间构建 KNN 图
pbmc <- FindNeighbors(pbmc, dims = 1:10)

# Louvain 算法聚类，resolution=0.5 适合约3K细胞
pbmc <- FindClusters(pbmc, resolution = 0.5)

Number of nodes: 2638
Number of edges: 95927
Number of communities: 9

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

# 为每个聚类寻找特异性标记基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)

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)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

Seurat 对象
|-- assays
|   +-- RNA
|       |-- layers
|       |   |-- counts    # 原始 UMI 计数矩阵
|       |   |-- data       # 归一化后的矩阵
|       |   +-- scale.data # 标准化后的矩阵
|       +-- VariableFeatures  # 高变基因列表
|-- reductions
|   |-- pca               # PCA 降维结果
|   +-- umap              # UMAP 降维结果
|-- meta.data             # 细胞元数据（nCount, nFeature, percent.mt 等）
+-- active.ident          # 当前聚类标签

###############################################################################
# Seurat v5 PBMC 3k 标准分析流程
# 基于官方教程: https://satijalab.org/seurat/articles/pbmc3k_tutorial
# 数据: 10X Genomics PBMC 2,700 个外周血单核细胞
###############################################################################
#
# 学习目标: 跑通单细胞 RNA-seq 标准分析全流程
# 流程链路: 读取数据 -> QC过滤 -> 归一化 -> 高变基因 -> 标准化 -> PCA
#           -> 确定维度 -> 聚类 -> UMAP -> 寻找marker -> 细胞类型注释
#
###############################################################################


# ============================================================================
# 0. 环境准备
# ============================================================================

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

# 设置工作目录（请根据实际路径修改）
data_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/data/filtered_gene_bc_matrices/hg19"
output_dir <- "E:/BaiduNetdiskDownload/生物信息学/pbmc3k_tutorial/output"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)


# ============================================================================
# 1. 读取数据 & 创建 Seurat 对象
# ============================================================================
# Read10X() 读取 10X CellRanger 输出目录（包含 barcodes.tsv, genes.tsv, matrix.mtx）
# 返回一个稀疏矩阵: 行=基因, 列=细胞, 值=UMI计数

cat("\n========== 1. 读取数据 ==========\n")
pbmc.data <- Read10X(data.dir = data_dir)

# 创建 Seurat 对象
# min.cells = 3:   至少在 3 个细胞中检测到的基因才保留
# min.features = 200: 每个细胞至少检测到 200 个基因才保留
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k",
                           min.cells = 3, min.features = 200)
print(pbmc)

# 查看计数矩阵的稀疏性（大多数值为 0，用稀疏矩阵存储节省内存）
cat("\n--- 稀疏矩阵示例（前30个细胞中3个基因）---\n")
print(pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30])


# ============================================================================
# 2. 质量控制（QC）& 细胞过滤
# ============================================================================
# 常用 QC 指标:
#   - nFeature_RNA: 每个细胞检测到的独特基因数（过低=空液滴/低质量，过高=doublet）
#   - nCount_RNA:   每个细胞检测到的分子总数（UMI总数）
#   - percent.mt:   线粒体基因占比（过高=低质量/凋亡细胞）

cat("\n========== 2. 质量控制 ==========\n")

# 计算线粒体基因表达比例（基因名以 MT- 开头）
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 查看 QC 统计
cat("\n--- QC 指标（前5个细胞）---\n")
print(head(pbmc@meta.data, 5))

# 可视化 QC 指标（小提琴图）
p_qc_vln <- VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
ggsave(file.path(output_dir, "01_qc_violin.png"), plot = p_qc_vln, width = 12, height = 5, dpi = 150)
cat("已保存: 01_qc_violin.png\n")

# FeatureScatter: 查看指标间关系
p_scatter1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
p_scatter2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
p_qc_scatter <- p_scatter1 + p_scatter2
ggsave(file.path(output_dir, "02_qc_scatter.png"), plot = p_qc_scatter, width = 12, height = 5, dpi = 150)
cat("已保存: 02_qc_scatter.png\n")

# 过滤细胞:
#   nFeature_RNA > 200  （已有 min.features=200 过滤，这里二次确认）
#   nFeature_RNA < 2500 （去除 doublet）
#   percent.mt < 5      （去除线粒体污染严重的细胞）
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
cat("过滤后细胞数:", ncol(pbmc), "\n")


# ============================================================================
# 3. 数据归一化（Normalization）
# ============================================================================
# LogNormalize: 每个细胞的表达量除以该细胞总UMI，乘以10000，再取log
# 结果存储在 pbmc[["RNA"]]$data

cat("\n========== 3. 归一化 ==========\n")
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
cat("归一化完成\n")


# ============================================================================
# 4. 识别高变基因（Feature Selection）
# ============================================================================
# FindVariableFeatures: 识别在不同细胞间表达变异最大的基因
# 默认选择 2000 个高变基因，用于下游 PCA 降维

cat("\n========== 4. 高变基因识别 ==========\n")
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 查看前10个高变基因
top10 <- head(VariableFeatures(pbmc), 10)
cat("Top 10 高变基因:", paste(top10, collapse = ", "), "\n")

# 可视化高变基因
p_var <- VariableFeaturePlot(pbmc)
p_var_label <- LabelPoints(plot = p_var, points = top10, repel = TRUE)
ggsave(file.path(output_dir, "03_variable_features.png"), plot = p_var + p_var_label, width = 10, height = 7, dpi = 150)
cat("已保存: 03_variable_features.png\n")


# ============================================================================
# 5. 数据标准化（Scaling）
# ============================================================================
# ScaleData: 对每个基因做 Z-score 标准化（均值=0, 方差=1）
# 使下游 PCA 分析中各基因权重一致
# 结果存储在 pbmc[["RNA"]]$scale.data

cat("\n========== 5. 标准化 ==========\n")
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
cat("标准化完成（", length(all.genes), "个基因）\n")


# ============================================================================
# 6. 线性降维（PCA）
# ============================================================================
# 对标准化后的高变基因做 PCA
# 每个 PC 代表一组相关基因的「元特征」

cat("\n========== 6. PCA 降维 ==========\n")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

# 查看前5个 PC 的 top 基因
cat("\n--- 前5个主成分的标记基因---\n")
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

# 可视化 PCA
p_pca <- DimPlot(pbmc, reduction = "pca") + NoLegend()
ggsave(file.path(output_dir, "04_pca.png"), plot = p_pca, width = 7, height = 6, dpi = 150)
cat("已保存: 04_pca.png\n")

# PCA 负载图（前2个 PC）
p_pca_load <- VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
ggsave(file.path(output_dir, "05_pca_loadings.png"), plot = p_pca_load, width = 10, height = 8, dpi = 150)
cat("已保存: 05_pca_loadings.png\n")

# PCA 热图（前15个 PC）
p_pca_heat <- DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
ggsave(file.path(output_dir, "06_pca_heatmap.png"), plot = p_pca_heat, width = 12, height = 15, dpi = 150)
cat("已保存: 06_pca_heatmap.png\n")


# ============================================================================
# 7. 确定数据维度
# ============================================================================
# 选择用于下游分析的主成分数量
# ElbowPlot: 按方差解释比例排序各 PC，拐点处即为建议维度

cat("\n========== 7. 确定维度 ==========\n")
p_elbow <- ElbowPlot(pbmc)
ggsave(file.path(output_dir, "07_elbow_plot.png"), plot = p_elbow, width = 7, height = 5, dpi = 150)
cat("已保存: 07_elbow_plot.png\n")
cat("拐点约在 PC 9-10，本教程选取 10 个 PC\n")


# ============================================================================
# 8. 细胞聚类（Clustering）
# ============================================================================
# FindNeighbors: 基于 PCA 空间构建 KNN 图（K近邻图）
# FindClusters:  使用 Louvain 算法做模块化优化聚类
#   resolution: 聚类粒度，0.4-1.2 适合约 3K 细胞的数据集

cat("\n========== 8. 聚类 ==========\n")
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# 查看聚类结果
cat("\n--- 前5个细胞的聚类 ID ---\n")
print(head(Idents(pbmc), 5))
cat("聚类数量:", length(levels(pbmc)), "\n")


# ============================================================================
# 9. 非线性降维（UMAP）
# ============================================================================
# UMAP: 将高维数据投射到 2D 空间用于可视化
# 注意: UMAP 用于可视化探索，不应仅凭图形下生物学结论

cat("\n========== 9. UMAP ==========\n")
pbmc <- RunUMAP(pbmc, dims = 1:10)

p_umap <- DimPlot(pbmc, reduction = "umap")
ggsave(file.path(output_dir, "08_umap.png"), plot = p_umap, width = 8, height = 7, dpi = 150)
cat("已保存: 08_umap.png\n")

# 保存中间结果（可重新加载，无需重跑前面的步骤）
saveRDS(pbmc, file.path(output_dir, "pbmc_tutorial.rds"))
cat("已保存中间结果: pbmc_tutorial.rds\n")


# ============================================================================
# 10. 寻找差异表达基因（Marker Genes）
# ============================================================================
# FindAllMarkers: 为每个聚类寻找特异性标记基因（与其他所有细胞比较）
# presto 包加速 DE 分析（Seurat v5 自动调用）

cat("\n========== 10. Marker 基因 ==========\n")
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)

# 筛选 log2FC > 1 的显著 marker
pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) -> top_markers

cat("\n--- 各聚类 marker 数量 ---\n")
print(table(top_markers$cluster))

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

# 每个聚类取 top 10 基因做热图
top10_all <- pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup()

p_heatmap <- DoHeatmap(pbmc, features = top10_all$gene) + NoLegend()
ggsave(file.path(output_dir, "09_marker_heatmap.png"), plot = p_heatmap, width = 12, height = 15, dpi = 150)
cat("已保存: 09_marker_heatmap.png\n")

# FeaturePlot: 在 UMAP 上展示经典 marker 基因表达
p_feature <- FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14",
                                              "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
ggsave(file.path(output_dir, "10_feature_plot.png"), plot = p_feature, width = 14, height = 12, dpi = 150)
cat("已保存: 10_feature_plot.png\n")

# VlnPlot: 小提琴图展示特定 marker
p_vln <- VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
ggsave(file.path(output_dir, "11_violin_markers.png"), plot = p_vln, width = 8, height = 5, dpi = 150)
cat("已保存: 11_violin_markers.png\n")


# ============================================================================
# 11. 细胞类型注释
# ============================================================================
# 根据已知 canonical markers 将聚类映射到细胞类型
#
# 聚类  Marker 基因          细胞类型
# 0    IL7R, CCR7           Naive CD4+ T
# 1    CD14, LYZ            CD14+ Mono（单核细胞）
# 2    IL7R, S100A4         Memory CD4+ T
# 3    MS4A1                 B（B细胞）
# 4    CD8A                  CD8+ T
# 5    FCGR3A, MS4A7         FCGR3A+ Mono
# 6    GNLY, NKG7            NK（自然杀伤细胞）
# 7    FCER1A, CST3          DC（树突状细胞）
# 8    PPBP                  Platelet（血小板）

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

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)

# 带标签的 UMAP 图
p_umap_labeled <- DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ggsave(file.path(output_dir, "12_umap_labeled.png"), plot = p_umap_labeled, width = 8, height = 7, dpi = 150)
cat("已保存: 12_umap_labeled.png\n")

# 保存最终结果
saveRDS(pbmc, file.path(output_dir, "pbmc3k_final.rds"))
cat("已保存最终结果: pbmc3k_final.rds\n")


# ============================================================================
# 12. 总结
# ============================================================================

cat("\n")
cat("=======================================================\n")
cat("  PBMC 3k 标准分析流程完成!\n")
cat("=======================================================\n")
cat("  输入细胞数: 2,700\n")
cat("  过滤后细胞数:", ncol(pbmc), "\n")
cat("  聚类数:", length(levels(pbmc)), "\n")
cat("  细胞类型: Naive CD4 T / CD14+ Mono / Memory CD4 T / B /\n")
cat("           CD8 T / FCGR3A+ Mono / NK / DC / Platelet\n")
cat("  输出目录:", output_dir, "\n")
cat("  输出文件:\n")
cat("    01_qc_violin.png       - QC小提琴图\n")
cat("    02_qc_scatter.png      - QC散点图\n")
cat("    03_variable_features.png - 高变基因图\n")
cat("    04_pca.png             - PCA散点图\n")
cat("    05_pca_loadings.png    - PCA负载图\n")
cat("    06_pca_heatmap.png     - PCA热图\n")
cat("    07_elbow_plot.png      - 拐点图\n")
cat("    08_umap.png            - UMAP聚类图\n")
cat("    09_marker_heatmap.png  - Marker热图\n")
cat("    10_feature_plot.png    - Feature图\n")
cat("    11_violin_markers.png  - Marker小提琴图\n")
cat("    12_umap_labeled.png    - 注释后UMAP图\n")
cat("    all_markers.csv        - 全部Marker表\n")
cat("    pbmc_tutorial.rds      - 中间结果\n")
cat("    pbmc3k_final.rds       - 最终结果\n")
cat("=======================================================\n")