571 Seurat v5 单细胞分析从零入门¶

Scanpy:  read → filter → normalize → log1p → HVG → scale → PCA → neighbors → UMAP → leiden → markers
Seurat:  Read10X → CreateSeuratObject → filter → Normalize → FindVariableFeatures → Scale → PCA → FindNeighbors → FindClusters → UMAP → FindAllMarkers

# v4 的写法（旧）
seurat_obj@assays$RNA@counts     # 原始计数
seurat_obj@assays$RNA@data       # 归一化数据

# v5 的写法（新）
LayerData(seurat_obj, layer = "counts")    # 原始计数
LayerData(seurat_obj, layer = "data")      # 归一化数据

# 安装 Seurat v5（当前最新版 5.5.0）
install.packages('Seurat')  # 从 CRAN 安装，会自动安装依赖

# 安装常用辅助包
install.packages('patchwork')    # 拼图工具，把多个 ggplot 拼在一起
install.packages('dplyr')        # 数据处理神器
install.packages('ggplot2')      # 画图核心包

# SingleR 在 Bioconductor 上，需要用 BiocManager 安装
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")  # 先安装 BiocManager

BiocManager::install("SingleR")            # 安装 SingleR（自动细胞类型注释）
BiocManager::install("celldex")            # 安装参考数据集包
BiocManager::install("SingleCellExperiment")  # SingleR 依赖的数据结构

library(Seurat)              # 加载 Seurat
packageVersion("Seurat")     # 查看版本号，应该是 5.5.0

# 如果你更习惯 conda
conda create -n seurat_env r-base=4.4 -y  # 创建 R 环境
conda activate seurat_env                   # 激活环境
conda install -c conda-forge r-seurat -y    # 安装 Seurat
# 注意：conda 版本可能不是最新，建议在 R 里用 install.packages

# ---- 加载所需的 R 包 ----
library(Seurat)       # 单细胞分析主包
library(dplyr)        # 数据处理
library(patchwork)    # 图片拼接
library(ggplot2)      # 画图

# ---- 下载 PBMC 3k 数据 ----
# 10x 官方提供的经典数据集，约 2700 个外周血单个核细胞
# 方法一：从网上下载（需要联网）
# 下载地址：https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
# 解压后得到 filtered_gene_bc_matrices/hg19/ 目录

# 方法二：用 SeuratData 包（更方便）
# install.packages("SeuratData")
# library(SeuratData)
# InstallData("pbmc3k")
# data("pbmc3k")

# ---- 读取 10x 数据 ----
# 假设数据解压在 data/ 目录下
pbmc.data <- Read10X(data.dir = "data/filtered_gene_bc_matrices/hg19/")
# Read10X 会自动读取 barcodes.tsv, genes.tsv, matrix.mtx 三个文件
# 返回一个稀疏矩阵（行=基因, 列=细胞）

# 查看数据维度
dim(pbmc.data)  # 应该是 32738 行（基因）× 2700 列（细胞）

# 看前几个基因和细胞
pbmc.data[1:5, 1:3]  # 前 5 个基因 × 前 3 个细胞

# ---- 创建 Seurat 对象 ----
pbmc <- CreateSeuratObject(
    counts = pbmc.data,     # 原始计数矩阵
    project = "pbmc3k",     # 项目名称（会显示在 metadata 里）
    min.cells = 3,          # 一个基因至少在 3 个细胞中表达才保留
    min.features = 200      # 一个细胞至少检测到 200 个基因才保留
)

# 查看 Seurat 对象
pbmc  # 会显示类似：An object of class Seurat, 13714 features across 2700 samples

# ---- 计算线粒体基因比例 ----
# 线粒体基因以 "MT-" 开头（人类），高比例说明细胞可能已经死亡
pbmc[["percent.mt"]] <- PercentageFeatureSet(
    pbmc,
    pattern = "^MT-"  # 正则表达式：以 MT- 开头的基因
)

# 查看 QC 指标
head(pbmc@meta.data, 5)  # 查看前 5 个细胞的元数据
# nCount_RNA   = 总 UMI 数（类似 Scanpy 的 total_counts）
# nFeature_RNA = 基因数（类似 Scanpy 的 n_genes_by_counts）
# percent.mt   = 线粒体比例

# ---- 画 QC 小提琴图 ----
VlnPlot(
    pbmc,
    features = c("nFeature_RNA",   # 基因数
                 "nCount_RNA",     # UMI 数
                 "percent.mt"),    # 线粒体比例
    ncol = 3                       # 三张图并排
)
ggsave("qc_violin.png", width = 12, height = 4)  # 保存图片

# ---- 画散点图看相关性 ----
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2  # 用 patchwork 拼图
ggsave("qc_scatter.png", width = 12, height = 5)

# ---- 过滤 ----
pbmc <- subset(
    pbmc,
    subset = nFeature_RNA > 200 &     # 基因数 > 200（去掉空液滴）
             nFeature_RNA < 2500 &     # 基因数 < 2500（去掉 doublet）
             percent.mt < 5            # 线粒体比例 < 5%（去掉死细胞）
)

cat("质控后:", ncol(pbmc), "个细胞,", nrow(pbmc), "个基因\n")
# 典型输出：质控后: 2638 个细胞, 13714 个基因

# ---- 标准归一化 ----
pbmc <- NormalizeData(
    pbmc,
    normalization.method = "LogNormalize",  # 方法：先除以总 UMI，再乘以 scale.factor，再取 log
    scale.factor = 10000                     # 缩放因子（和 Scanpy 的 target_sum=1e4 一样）
)
# 归一化后的数据存在 Layers 的 "data" 层

# ---- SCTransform v2 归一化（推荐） ----
# 同时完成归一化 + 高变基因选择 + 缩放，一步到位
# pbmc <- SCTransform(
#     pbmc,
#     vars.to.regress = "percent.mt",  # 回归掉线粒体比例的影响
#     verbose = FALSE                   # 不输出进度信息
# )
# 注意：如果用 SCTransform，可以跳过 Step 4 和 Step 5 的 FindVariableFeatures + ScaleData

# ---- 找高变基因（Highly Variable Features, HVFs）----
pbmc <- FindVariableFeatures(
    pbmc,
    selection.method = "vst",   # 使用 VST 方法（方差稳定化转换）
    nfeatures = 2000            # 选 2000 个高变基因（默认值）
)

# 查看 top 10 高变基因
top10 <- head(VariableFeatures(pbmc), 10)  # 获取前 10 个高变基因名
print(top10)

# 画高变基因图
plot1 <- VariableFeaturePlot(pbmc)                   # 散点图
plot2 <- LabelPoints(plot = plot1, points = top10)   # 标注 top10
plot2
ggsave("hvg.png", width = 10, height = 5)

# ---- 缩放：让每个基因的均值为 0，标准差为 1 ----
all.genes <- rownames(pbmc)     # 获取所有基因名
pbmc <- ScaleData(
    pbmc,
    features = all.genes,       # 对所有基因缩放（默认只对高变基因缩放）
    vars.to.regress = "percent.mt"  # 可选：回归掉线粒体比例
)
# 缩放后的数据存在 Layers 的 "scale.data" 层

# ---- PCA 降维 ----
pbmc <- RunPCA(
    pbmc,
    features = VariableFeatures(object = pbmc),  # 只用高变基因做 PCA
    npcs = 50                                     # 计算 50 个主成分
)

# 查看前 2 个 PC 的 top 基因
print(pbmc[["pca"]], dims = 1:2, nfeatures = 5)  # 每个 PC 显示 5 个基因

# 画 PCA 图
DimPlot(pbmc, reduction = "pca")  # 在 PC1-PC2 平面画细胞
ggsave("pca.png", width = 7, height = 5)

# 画肘部图（Elbow Plot）决定用多少个 PC
ElbowPlot(pbmc, ndims = 50)  # 显示前 50 个 PC 的标准差
ggsave("elbow.png", width = 7, height = 5)
# 一般在第 10-15 个 PC 附近出现拐点

# ---- 构建 KNN 邻居图 ----
pbmc <- FindNeighbors(
    pbmc,
    dims = 1:10  # 使用前 10 个 PC（根据肘部图决定）
)

# ---- 聚类（Leiden / Louvain）----
pbmc <- FindClusters(
    pbmc,
    resolution = 0.5,       # 分辨率：越大分的群越多
    algorithm = 4            # 算法：1=Louvain, 4=Leiden（推荐 Leiden）
)

# 查看聚类结果
head(Idents(pbmc))  # 显示前几个细胞的聚类标签
table(Idents(pbmc)) # 统计每个 cluster 的细胞数

# ---- UMAP 降维 ----
pbmc <- RunUMAP(
    pbmc,
    dims = 1:10  # 使用前 10 个 PC
)

# 画 UMAP 聚类图
DimPlot(
    pbmc,
    reduction = "umap",      # 使用 UMAP 坐标
    label = TRUE,             # 在每个 cluster 上显示编号
    pt.size = 0.5             # 点的大小
) + NoLegend()                # 去掉图例（因为标签已经在图上了）
ggsave("umap_clusters.png", width = 7, height = 6)

# ---- 找所有 cluster 的 marker 基因 ----
pbmc.markers <- FindAllMarkers(
    pbmc,
    only.pos = TRUE,          # 只看上调的基因（正的 log2FC）
    min.pct = 0.25,           # 至少在 25% 的细胞中表达
    logfc.threshold = 0.25,   # log2 fold change 阈值
    test.use = "wilcox"       # 使用 Wilcoxon 检验（和 Scanpy 一样）
)

# 查看每个 cluster 的 top 5 marker
top5 <- pbmc.markers %>%
    group_by(cluster) %>%          # 按 cluster 分组
    filter(p_val_adj < 0.05) %>%   # 只保留显著的（调整后 p < 0.05）
    slice_head(n = 5)              # 每组取前 5 个
print(top5, n = 40)                # 打印结果

# 画 marker 基因热图
DoHeatmap(
    pbmc,
    features = top5$gene    # 用 top5 marker 基因
) + NoLegend()
ggsave("heatmap_markers.png", width = 12, height = 8)

# ---- PBMC 常见细胞类型的 marker 基因 ----
# 先用 FeaturePlot 看关键 marker 在 UMAP 上的表达
FeaturePlot(
    pbmc,
    features = c("MS4A1",   # B 细胞
                 "CD14",    # CD14 单核细胞
                 "LYZ",     # 单核细胞
                 "IL7R",    # CD4 T 细胞
                 "CD8A",    # CD8 T 细胞
                 "GNLY",    # NK 细胞
                 "FCGR3A",  # CD16 单核细胞
                 "FCER1A",  # 树突状细胞
                 "PPBP"),   # 血小板
    ncol = 3                # 每行 3 图
)
ggsave("featureplot_markers.png", width = 15, height = 12)

# 用 DotPlot 看每个 cluster 的 marker 表达
DotPlot(
    pbmc,
    features = c("IL7R", "CD4", "CD8A", "CD8B", "MS4A1", "CD79A",
                 "GNLY", "NKG7", "CD14", "LYZ", "FCGR3A", "MS4A7",
                 "FCER1A", "CST3", "PPBP"),
    dot.scale = 8               # 点的最大大小
) + RotatedAxis()               # 横轴标签旋转
ggsave("dotplot_markers.png", width = 12, height = 5)

# ---- 手动注释：给每个 cluster 贴标签 ----
# 根据 FeaturePlot 和 DotPlot 结果判断（具体对应关系需要看你的图）
new.cluster.ids <- c(
    "Naive CD4 T",        # cluster 0
    "CD14+ Mono",         # cluster 1
    "Memory CD4 T",       # cluster 2
    "B",                  # cluster 3
    "CD8 T",              # cluster 4
    "FCGR3A+ Mono",       # cluster 5
    "NK",                 # cluster 6
    "DC",                 # cluster 7
    "Platelet"            # cluster 8
)
names(new.cluster.ids) <- levels(pbmc)     # 把标签和 cluster 编号对应
pbmc <- RenameIdents(pbmc, new.cluster.ids) # 重命名

# 画注释后的 UMAP
DimPlot(
    pbmc,
    reduction = "umap",
    label = TRUE,           # 显示标签
    pt.size = 0.5
) + NoLegend()
ggsave("umap_celltype.png", width = 8, height = 6)

# ---- 使用 SingleR 自动注释（需要先安装 SingleR 和 celldex）----
library(SingleR)                    # 加载 SingleR
library(celldex)                    # 加载参考数据集
library(SingleCellExperiment)       # 加载数据结构包

# 获取参考数据集（人类免疫细胞）
ref <- celldex::HumanPrimaryCellAtlasData()  # 下载人类细胞图谱参考数据

# 提取 Seurat 对象的表达矩阵，转为 SingleCellExperiment 格式
sce <- as.SingleCellExperiment(pbmc)  # Seurat 对象转 SCE 对象

# 运行 SingleR 自动注释
pred <- SingleR(
    test = sce,              # 你的数据
    ref = ref,               # 参考数据集
    labels = ref$label.main  # 使用主要标签（粗分类）
)

# 把预测结果加到 Seurat 对象的 metadata 里
pbmc$SingleR_labels <- pred$labels  # 添加注释列

# 画 SingleR 注释结果
DimPlot(
    pbmc,
    reduction = "umap",
    group.by = "SingleR_labels",  # 用 SingleR 标签着色
    label = TRUE,
    repel = TRUE,                  # 标签不重叠
    pt.size = 0.5
) + NoLegend()
ggsave("umap_singler.png", width = 10, height = 7)

# ---- 保存 Seurat 对象 ----
saveRDS(pbmc, file = "pbmc3k_analyzed.rds")   # 保存为 RDS 文件
# 下次加载：pbmc <- readRDS("pbmc3k_analyzed.rds")

# ---- 导出 metadata 为表格 ----
write.csv(
    pbmc@meta.data,                            # 细胞元数据
    file = "pbmc3k_metadata.csv",              # 保存路径
    row.names = TRUE                            # 保留 barcode 作为行名
)

Seurat 对象结构：
┌──────────────────────────────────────────────────┐
│  @assays                                          │
│  └── RNA (Assay5 对象)                            │
│      ├── layers$counts     → 原始计数矩阵         │  ← 类似 Scanpy 的 adata.layers['counts']
│      ├── layers$data       → 归一化后的数据        │  ← 类似 Scanpy 的 adata.X
│      └── layers$scale.data → 缩放后的数据          │
├──────────────────────────────────────────────────┤
│  @meta.data (DataFrame)                           │
│  每行一个细胞，包含：                               │
│  - nCount_RNA      总 UMI 数                       │  ← 类似 adata.obs['total_counts']
│  - nFeature_RNA    基因数                          │  ← 类似 adata.obs['n_genes_by_counts']
│  - percent.mt      线粒体比例                      │  ← 类似 adata.obs['pct_counts_mt']
│  - seurat_clusters 聚类标签                        │  ← 类似 adata.obs['leiden']
├──────────────────────────────────────────────────┤
│  @reductions                                      │
│  ├── pca   → PCA 降维结果                          │  ← 类似 adata.obsm['X_pca']
│  ├── umap  → UMAP 降维结果                         │  ← 类似 adata.obsm['X_umap']
│  └── tsne  → t-SNE 降维结果（如果计算了的话）       │  ← 类似 adata.obsm['X_tsne']
├──────────────────────────────────────────────────┤
│  @graphs                                          │
│  ├── RNA_nn    → KNN 图（最近邻）                   │  ← 类似 adata.obsp['distances']
│  └── RNA_snn   → SNN 图（共享最近邻）               │  ← 类似 adata.obsp['connectivities']
├──────────────────────────────────────────────────┤
│  @commands   → 记录你执行过的每一步操作（可复现）     │
│  @misc       → 杂物柜（自定义数据）                 │  ← 类似 adata.uns
└──────────────────────────────────────────────────┘

pbmc                           # 打印对象概况
dim(pbmc)                      # 维度：基因数 × 细胞数（注意和 Scanpy 反过来）
head(pbmc@meta.data)           # 查看细胞元数据
Idents(pbmc)                   # 查看当前的细胞标签（默认是聚类结果）
Layers(pbmc)                   # 查看有哪些数据层
Reductions(pbmc)               # 查看有哪些降维结果
colnames(pbmc)[1:5]            # 前 5 个细胞的 barcode
rownames(pbmc)[1:5]            # 前 5 个基因名

# 原因：Seurat v5 的 Assay5 数据结构中，某些操作还没生成对应的 layer
# 解决：确保在运行 NormalizeData 之后再使用需要 "data" 层的函数

# 如果是从旧版 Seurat 对象升级来的，可以试：
pbmc[["RNA"]] <- as(pbmc[["RNA"]], "Assay5")  # 转换为 v5 格式

# 原因：在 FindNeighbors 之前没有跑 RunPCA
# 解决：按顺序执行流程

# 正确顺序：
# NormalizeData → FindVariableFeatures → ScaleData → RunPCA → FindNeighbors

# 原因：系统缺少 C++ 编译工具或依赖库
# 解决：

# Ubuntu/Debian:
# sudo apt-get install libcurl4-openssl-dev libssl-dev libxml2-dev libhdf5-dev

# CentOS/RHEL:
# sudo yum install libcurl-devel openssl-devel libxml2-devel hdf5-devel

# macOS:
# xcode-select --install
# brew install hdf5

# 原因：SCTransform 会重新筛选基因，导致和之前的 Assay 基因不一致
# 解决：这通常是正常警告，不用担心
# 如果需要某些特定基因，可以在 SCTransform 之前就过滤好

# 或者使用标准归一化流程避免此问题：
# NormalizeData → FindVariableFeatures → ScaleData

# 原因：数据太大（几十万到百万细胞）
# 解决方案：

# 方法一：使用 BPCells 后端（磁盘上操作，不全加载到内存）
# install.packages("BPCells")
# library(BPCells)
# mat <- open_matrix_dir("path/to/bpcells_dir")
# pbmc <- CreateSeuratObject(counts = mat)

# 方法二：使用 Sketch 工作流（先抽样分析再映射回全量）
# pbmc <- SketchData(
#     object = pbmc,
#     ncells = 5000,           # 抽取 5000 个代表性细胞
#     method = "LeverageScore"
# )

# 读 → 建对象 → 质控 → 归一化 → 高变基因 → 缩放 → PCA → 邻居 → 聚类 → UMAP → Marker → 注释
Read10X → CreateSeuratObject → subset → NormalizeData → FindVariableFeatures → ScaleData → RunPCA → FindNeighbors → FindClusters → RunUMAP → FindAllMarkers → 手动/SingleR