Seurat — R语言单细胞RNA-seq分析核心工具包¶

# 在R中安装Seurat v5（最新版v5.4.0）
install.packages("Seurat")          # 从CRAN安装

# 如果CRAN版本较旧，从GitHub安装最新版
install.packages("remotes")
remotes::install_github("satijalab/seurat", ref="seurat5")

# 安装常用配套包
install.packages(c(
    "ggplot2",    # 可视化
    "dplyr",      # 数据处理
    "patchwork"   # 图片拼合
))

# 加载库
library(Seurat)
library(ggplot2)
library(dplyr)

# 验证版本
packageVersion("Seurat")   # 应显示 5.5.0 或更新

# 读入CellRanger/STARsolo输出的10X格式矩阵
counts <- Read10X(
    data.dir = "sample/outs/filtered_feature_bc_matrix/"  # 矩阵目录
)

# 创建Seurat对象（最低细胞数和基因数过滤）
seurat_obj <- CreateSeuratObject(
    counts = counts,          # 计数矩阵
    project = "PBMC_3p",      # 项目名称
    min.cells = 3,            # 至少出现在3个细胞中的基因才保留
    min.features = 200        # 至少检测到200个基因的细胞才保留
)

# 查看基本信息
seurat_obj                    # 输出：细胞数、基因数等

# 计算线粒体基因百分比（质控关键指标）
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(
    seurat_obj,
    pattern = "^MT-"      # 人类线粒体基因以MT-开头（小鼠是^mt-）
)

# 可视化QC指标分布
VlnPlot(
    seurat_obj,
    features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),  # 三个关键指标
    ncol = 3
)

# 过滤低质量细胞（阈值根据数据调整）
seurat_obj <- subset(
    seurat_obj,
    subset = nFeature_RNA > 200 &   # 最少检测到200个基因
             nFeature_RNA < 5000 &  # 最多5000个基因（去除双细胞）
             percent.mt < 20        # 线粒体基因比例<20%（去除死细胞）
)

# 标准化（每细胞总表达量归一化到10000，再log转化）
seurat_obj <- NormalizeData(seurat_obj)

# 筛选高变基因（最具信息量的基因，用于后续分析）
seurat_obj <- FindVariableFeatures(
    seurat_obj,
    selection.method = "vst",  # VST方法（推荐）
    nfeatures = 2000           # 选取2000个高变基因
)

# 数据缩放（中心化，消除基因表达量差异）
seurat_obj <- ScaleData(seurat_obj)

# PCA降维
seurat_obj <- RunPCA(
    seurat_obj,
    features = VariableFeatures(seurat_obj)  # 使用高变基因做PCA
)

# 查看PCA结果，选择合适的PC数
ElbowPlot(seurat_obj, ndims = 50)  # 肘部图确定PC数（一般15-30个）

# 构建邻居图（基于PCA空间中的KNN）
seurat_obj <- FindNeighbors(
    seurat_obj,
    dims = 1:20    # 使用前20个PC（根据ElbowPlot调整）
)

# Leiden/Louvain聚类（resolution控制细胞类型粒度）
seurat_obj <- FindClusters(
    seurat_obj,
    resolution = 0.5    # 分辨率：0.1-0.3少聚类，0.5-1.0多聚类
)

# UMAP降维可视化
seurat_obj <- RunUMAP(
    seurat_obj,
    dims = 1:20    # 与FindNeighbors使用相同PC数
)

# 绘制UMAP图
DimPlot(
    seurat_obj,
    reduction = "umap",  # 使用UMAP坐标
    label = TRUE         # 在图上显示聚类编号
)

# 完整单细胞分析流程（PBMC数据）

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

# 1. 读取数据
counts <- Read10X("pbmc_cellranger_output/filtered_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts, project="PBMC", min.cells=3, min.features=200)

# 2. QC过滤
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern="^MT-")
pbmc <- subset(pbmc, subset=nFeature_RNA>200 & nFeature_RNA<4500 & percent.mt<15)
cat("过滤后细胞数:", ncol(pbmc), "\n")

# 3. 预处理
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, nfeatures=2000)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)

# 4. 聚类和可视化
pbmc <- FindNeighbors(pbmc, dims=1:15)
pbmc <- FindClusters(pbmc, resolution=0.5)
pbmc <- RunUMAP(pbmc, dims=1:15)

# 5. 可视化
p1 <- DimPlot(pbmc, reduction="umap", label=TRUE) +   # UMAP聚类图
      ggtitle("PBMC Clusters")
ggsave("pbmc_umap.pdf", p1, width=8, height=6)

# 6. 找每个聚类的marker基因
pbmc_markers <- FindAllMarkers(
    pbmc,
    only.pos = TRUE,          # 仅找上调基因（marker）
    min.pct = 0.25,           # 至少25%的细胞表达
    logfc.threshold = 0.25    # 至少0.25 log2倍数变化
)

# 查看每个聚类的top5 marker基因
pbmc_markers %>%
    group_by(cluster) %>%
    top_n(n=5, wt=avg_log2FC)

# 7. 基于marker基因注释细胞类型（以PBMC为例）
# 已知PBMC cell type markers:
# CD3D, CD4 -> CD4+ T cells
# CD3D, CD8A -> CD8+ T cells
# MS4A1 -> B cells
# CD14, LYZ -> Monocytes
# GNLY, NKG7 -> NK cells

FeaturePlot(pbmc,
    features=c("CD3D","CD4","CD8A","MS4A1","CD14","GNLY"),
    ncol=3)                   # 可视化各marker基因表达

# 8. 重命名聚类（根据marker注释结果）
new_cluster_ids <- c("CD4 T","CD14 Mono","CD4 T","B","CD8 T",
                     "NK","CD14 Mono","DC","Platelet")
names(new_cluster_ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new_cluster_ids)

# 9. 保存结果
saveRDS(pbmc, "pbmc_analyzed.rds")   # 保存Seurat对象