单细胞 RNA 测序（scRNA-seq）是解析细胞异质性的核心技术 —— 它能突破 bulk RNA-seq 的 “平均化” 局限，精准捕捉单个细胞的基因表达特征，揭示组织中隐藏的细胞亚群、细胞状态转换及罕见细胞类型。但 scRNA-seq 分析流程长、数据噪声高（如 dropout 效应、批次效应）、参数调节复杂，新手常陷入 “卡在哪一步” 的困境。

本文基于主流分析工具 Seurat（R 语言），结合 SingleR、Harmony 等辅助包，打造 “从原始数据到可视化报告” 的全流程实操手册。从环境搭建、数据质控到细胞分群注释，每一步都附带可直接复用的代码模板、参数选择逻辑和避坑指南，兼顾新手入门与进阶需求，帮你高效完成 scRNA-seq 分析并产出论文级结果。

一、分析核心流程与环境搭建

1. 核心分析流程

scRNA-seq 分析的核心逻辑是 “去噪→降维→聚类→注释→可视化”，具体流程如下：

plaintext

原始数据（10x Genomics等）→ 数据导入与Seurat对象构建 → 质控（过滤低质量细胞/基因）→ 数据预处理（归一化→高变基因筛选→标准化）→ 批次校正（多样本/多批次数据）→ 降维（PCA→UMAP/tSNE）→ 聚类分群 → 细胞类型注释（标记基因手动注释+SingleR自动注释）→ 可视化（分群图、标记基因图、功能富集图）→ 结果保存与导出

2. 核心工具与环境搭建

（1）必备 R 包

• 核心分析：Seurat（单细胞分析一站式工具）、SingleR（细胞类型自动注释）
• 数据处理：dplyr（数据清洗）、tidyr（格式转换）
• 可视化：ggplot2（基础绘图）、pheatmap（热图）、patchwork（图表拼接）
• 批次校正：Harmony（高效批次整合）、sctransform（标准化 + 批次校正一体）
• 功能富集：clusterProfiler（GO/KEGG 富集）

（2）安装代码（R/RStudio 中运行）

# 安装Bioconductor包管理器（若未安装）
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.18") # 适配R 4.3+版本

# 安装核心包
BiocManager::install(c("Seurat", "SingleR", "Harmony", "clusterProfiler", "org.Hs.eg.db", "org.Mm.eg.db"))
install.packages(c("dplyr", "tidyr", "ggplot2", "pheatmap", "patchwork", "viridis"))

# 加载所有包（后续分析直接复用）
library(Seurat)
library(SingleR)
library(Harmony)
library(clusterProfiler)
library(org.Hs.eg.db) # 人类注释包（小鼠用org.Mm.eg.db）
library(dplyr)
library(tidyr)
library(ggplot2)
library(pheatmap)
library(patchwork)
library(viridis)

# 论文级可视化主题（复用前文，统一风格）
theme_paper <- function() {
  theme_bw() +
    theme(
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold", family = "Arial"),
      axis.title = element_text(size = 10, family = "Arial"),
      axis.text = element_text(size = 8, family = "Arial"),
      legend.title = element_text(size = 9, family = "Arial"),
      legend.text = element_text(size = 8, family = "Arial"),
      panel.grid = element_blank()
    )
}

（3）环境验证

运行sessionInfo()查看包版本，确保 Seurat≥4.0、Harmony≥0.1.1，避免版本兼容问题。若安装失败，优先通过 conda 创建专属环境（conda create -n scRNA r=4.3 seurat harmony -c bioconda -c conda-forge）。

二、数据导入与 Seurat 对象构建

scRNA-seq 数据格式以 10x Genomics 为代表，输出 3 个核心文件：barcodes.tsv.gz（细胞条形码）、features.tsv.gz（基因名 / ID）、matrix.mtx.gz（表达矩阵）。若为其他平台（如 Smart-seq2），需先转换为 “基因 × 细胞” 的表达矩阵。

1. 10x Genomics 数据导入（最常用场景）

假设数据存储在data/filtered_feature_bc_matrix/目录下，直接用 Seurat 的Read10X函数导入：

# 数据导入（10x Genomics格式）
data_dir <- "data/filtered_feature_bc_matrix" # 替换为你的数据路径
counts <- Read10X(data.dir = data_dir) # 读取后为稀疏矩阵（行=基因，列=细胞）

# 构建Seurat对象（设置项目名，保存原始计数）
seurat_obj <- CreateSeuratObject(
  counts = counts,
  project = "scRNA_Sample", # 项目名，可自定义
  min.cells = 3, # 至少在3个细胞中表达的基因才保留（过滤低置信基因）
  min.features = 200 # 至少表达200个基因的细胞才保留（初步过滤空细胞）
)

# 查看Seurat对象结构
print(seurat_obj)
# 输出示例：An object of class Seurat 
# 31053 features across 2792 samples within 1 assay 
# Active assay: RNA (31053 features, 0 variable features)

2. 非 10x 数据导入（如表达矩阵）

若已有 “基因 × 细胞” 的 CSV/TSV 格式表达矩阵，直接读取并构建对象：

# 读取CSV格式表达矩阵（行=基因，列=细胞）
counts_df <- read.csv("data/expression_matrix.csv", row.names = 1)
counts_sparse <- as(as.matrix(counts_df), "dgCMatrix") # 转为稀疏矩阵（节省内存）

# 构建Seurat对象
seurat_obj <- CreateSeuratObject(
  counts = counts_sparse,
  project = "scRNA_Sample",
  min.cells = 3,
  min.features = 200
)

3. 关键预处理：添加细胞元数据

为后续质控和分群做准备，添加线粒体基因比例、核糖体基因比例等元数据：

# 计算线粒体基因比例（线粒体基因以"MT-"开头，人类；小鼠为"Mt-"）
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

# 计算核糖体基因比例（核糖体基因以"RPS"或"RPL"开头）
seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# 查看元数据前5行（包含nFeature_RNA=基因数、nCount_RNA=测序深度、percent.mt=线粒体比例）
head(seurat_obj@meta.data)

避坑指南：小鼠线粒体基因是 “Mt-”（大写 M + 小写 t），人类是 “MT-”（全大写），若物种不符会导致线粒体比例计算为 0，需提前核对基因名格式。

三、质控：过滤低质量细胞与基因

scRNA-seq 的噪声主要来自：空细胞（仅含游离 RNA）、死细胞（线粒体基因比例高）、低质量细胞（测序深度低、基因数少）。质控的核心是 “剔除异常值”，保留高质量细胞用于后续分析。

1. 质控指标与阈值选择

（1）细胞水平质控指标

指标	含义	推荐阈值（人类组织通用）	调整原则
nFeature_RNA	细胞表达的基因数	200 ~ 6000	上皮细胞基因数偏高，免疫细胞偏低
nCount_RNA	细胞测序深度（总 reads 数）	500 ~ 50000	与基因数正相关，比例约 10:1
percent.mt	线粒体基因表达比例	< 5%	肿瘤细胞可放宽至 < 10%
percent.ribo	核糖体基因表达比例	10% ~ 40%	比例过低可能是低质量细胞

（2）基因水平质控指标

• 仅保留在≥3 个细胞中表达的基因（CreateSeuratObject已设置min.cells=3）
• 剔除在所有细胞中均不表达的基因（自动过滤）

2. 质控可视化与过滤

先用小提琴图和散点图查看质控指标分布，再基于阈值过滤：

# 1. 可视化质控指标（小提琴图）
vln_plot <- VlnPlot(
  seurat_obj,
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.ribo"),
  ncol = 2, # 2列布局
  pt.size = 0.1 # 点大小（避免重叠）
) + theme_paper()
print(vln_plot)
ggsave("plots/quality_control_vln.pdf", vln_plot, width = 10, height = 8, dpi = 300)

# 2. 可视化基因数与测序深度的相关性（散点图）
feature_count_plot <- FeatureScatter(
  seurat_obj,
  feature1 = "nCount_RNA",
  feature2 = "nFeature_RNA",
  group.by = "orig.ident",
  pt.size = 0.3
) + theme_paper()
print(feature_count_plot)
ggsave("plots/feature_count_cor.pdf", feature_count_plot, width = 8, height = 6, dpi = 300)

# 3. 基于阈值过滤低质量细胞
seurat_obj_filtered <- subset(
  seurat_obj,
  subset = nFeature_RNA > 200 & nFeature_RNA < 6000 &
    percent.mt < 5 & percent.ribo > 10
)

# 查看过滤前后细胞数变化
cat("过滤前细胞数：", ncol(seurat_obj), "\n")
cat("过滤后细胞数：", ncol(seurat_obj_filtered), "\n")
# 示例输出：过滤前细胞数：2792 → 过滤后细胞数：2518（保留约90%，合理范围）

3. 质控避坑关键点

• 阈值不可 “一刀切”：如脑组织线粒体比例可放宽至 < 8%，血液免疫细胞基因数可低至 100+；
• 若过滤后细胞数 < 50%，需检查测序质量（如是否存在文库构建问题）；
• 避免过度过滤：如基因数阈值设为 5000 可能剔除高基因数的特异性细胞（如干细胞）。

四、数据预处理：归一化、标准化与高变基因筛选

质控后的原始表达矩阵存在 “测序深度差异”“dropout 效应” 等问题，需通过预处理转化为适合降维聚类的标准化数据。

1. 归一化（Normalization）

目的：消除不同细胞测序深度的差异，使细胞间基因表达可比。Seurat 默认用LogNormalize（log2 转换），适用于大多数场景：

# 归一化（默认LogNormalize：表达量=log2(Counts/文库大小×10000 + 1)）
seurat_obj_filtered <- NormalizeData(
  object = seurat_obj_filtered,
  normalization.method = "LogNormalize",
  scale.factor = 10000 # 缩放因子（每细胞总表达量标准化到10000）
)

# 查看归一化后的数据（前5行前5列）
head(GetAssayData(seurat_obj_filtered, assay = "RNA", slot = "data")[1:5, 1:5])

替代方案：若数据 dropout 严重（如低测序深度），可用sctransform标准化（同时处理异方差和 dropout）：

# sctransform标准化（需加载sctransform包）
# install.packages("sctransform")
library(sctransform)
seurat_obj_filtered <- SCTransform(
  object = seurat_obj_filtered,
  vars.to.regress = "percent.mt", # 回归线粒体比例（去除噪声）
  verbose = FALSE
)
# 后续分析需将assay切换为"SCT"（而非默认"RNA"）

2. 高变基因筛选（Variable Features）

目的：筛选 “在部分细胞中高表达、部分细胞中低表达” 的基因 —— 这些基因是细胞异质性的核心来源（如细胞类型标记基因），忽略低变异基因（如管家基因）可减少噪声。

# 筛选高变基因（默认筛选2000个，可调整nfeatures参数）
seurat_obj_filtered <- FindVariableFeatures(
  object = seurat_obj_filtered,
  selection.method = "vst", # 方差稳定转换（适合单细胞数据）
  nfeatures = 2000
)

# 可视化前20个高变基因
top10_var_genes <- head(VariableFeatures(seurat_obj_filtered), 10)
var_genes_plot <- VariableFeaturePlot(seurat_obj_filtered) +
  LabelPoints(
    plot = .,
    points = top10_var_genes,
    repel = TRUE, # 避免标签重叠
    size = 3
  ) + theme_paper()
print(var_genes_plot)
ggsave("plots/variable_features.pdf", var_genes_plot, width = 10, height = 6, dpi = 300)

3. 标准化（Scaling）

目的：将基因表达量标准化为均值 = 0、方差 = 1 的 Z-score，使不同基因表达范围可比，避免高表达基因主导降维结果。

# 标准化（默认对高变基因进行标准化，可通过features参数指定所有基因）
seurat_obj_filtered <- ScaleData(
  object = seurat_obj_filtered,
  features = VariableFeatures(object = seurat_obj_filtered), # 仅标准化高变基因
  vars.to.regress = "percent.mt", # 回归线粒体比例（进一步去除技术噪声）
  verbose = FALSE
)

# 查看标准化后的数据（Z-score）
head(GetAssayData(seurat_obj_filtered, assay = "RNA", slot = "scale.data")[1:5, 1:5])

关键参数说明：vars.to.regress可回归技术变量（如线粒体比例、测序深度）或生物变量（如细胞周期），若存在明显批次效应，可在此处加入批次变量（如vars.to.regress = c("percent.mt", "batch")）。

五、批次校正（多样本 / 多批次数据必备）

若分析数据来自多个样本、多个测序批次，会存在 “批次效应”（如不同批次的测序条件差异导致细胞聚类按批次分组，而非细胞类型），需通过批次校正消除。

1. 批次效应检测

先通过 PCA 和 UMAP 查看批次分布，判断是否需要校正：

# 假设元数据中存在"batch"列（记录每个细胞的批次信息）
# 若未添加，先手动添加：seurat_obj_filtered$batch <- c(rep("batch1", 1000), rep("batch2", 1518))

# 先做PCA（基于高变基因）
seurat_obj_filtered <- RunPCA(
  object = seurat_obj_filtered,
  features = VariableFeatures(object = seurat_obj_filtered),
  verbose = FALSE
)

# 可视化PCA的批次分布（若不同批次在PC1/PC2上明显分离，需校正）
pca_batch_plot <- DimPlot(
  seurat_obj_filtered,
  reduction = "pca",
  group.by = "batch",
  pt.size = 0.3
) + theme_paper()
print(pca_batch_plot)
ggsave("plots/pca_batch_effect.pdf", pca_batch_plot, width = 8, height = 6, dpi = 300)

2. 批次校正方法：Harmony（推荐）

Harmony 是目前最常用的批次校正工具，速度快、效果好，可直接集成到 Seurat 流程：

# 安装并加载Harmony（若未安装）
# BiocManager::install("Harmony")
library(Harmony)

# 运行Harmony批次校正（基于PCA结果）
seurat_obj_filtered <- RunHarmony(
  object = seurat_obj_filtered,
  group.by.vars = "batch", # 批次变量名
  assay.use = "RNA",
  reduction.use = "pca",
  dims.use = 1:30, # 使用前30个PCA维度
  verbose = FALSE
)

# 查看校正后的PCA结果（批次混合度提升）
harmony_pca_plot <- DimPlot(
  seurat_obj_filtered,
  reduction = "harmony",
  group.by = "batch",
  pt.size = 0.3
) + theme_paper()
print(harmony_pca_plot)
ggsave("plots/harmony_corrected_pca.pdf", harmony_pca_plot, width = 8, height = 6, dpi = 300)

替代方案：Seurat 内置的 CCA 方法（RunCCA），适合批次差异较小的数据，代码如下：

seurat_obj_filtered <- RunCCA(
  object = seurat_obj_filtered,
  group.by = "batch",
  features = VariableFeatures(object = seurat_obj_filtered),
  verbose = FALSE
)

3. 批次校正避坑

• 校正后需验证：若 UMAP 图中批次仍明显分离，需检查批次变量是否正确、是否存在样本质量差异；
• 避免过度校正：若不同批次存在真实生物差异（如不同处理组），不可用批次校正消除；
• 单批次数据可跳过此步骤，节省计算时间。

六、降维聚类：从高维数据到细胞分群

scRNA-seq 表达矩阵是 “万级基因 × 千级细胞” 的高维数据，需通过降维（PCA→UMAP/tSNE）将高维数据映射到 2D/3D 空间，再基于细胞间表达相似性聚类，得到细胞亚群。

1. 主成分分析（PCA）：线性降维

PCA 将高维基因表达数据压缩到少数主成分（PC），每个 PC 代表一个 “基因表达模式”，前 30 个 PC 基本能保留核心生物信息：

# 若已做Harmony校正，用harmony结果做后续分析；否则用pca
reduction_use <- ifelse("harmony" %in% Reductions(seurat_obj_filtered), "harmony", "pca")

# 可视化PCA的Elbow Plot（判断最佳PC数）
elbow_plot <- ElbowPlot(seurat_obj_filtered, ndims = 50) + theme_paper()
print(elbow_plot)
ggsave("plots/elbow_plot.pdf", elbow_plot, width = 8, height = 6, dpi = 300)

Elbow Plot 解读：曲线拐点处的 PC 数为最佳维度（如拐点在 15，选择前 15 个 PC），此时后续 PC 的方差贡献很小，继续增加 PC 数不会提升聚类效果。

2. UMAP/tSNE：非线性降维（可视化核心）

PCA 是线性降维，难以捕捉细胞间的非线性关系；UMAP（推荐）和 tSNE 能更好地保留局部聚类结构，适合可视化：

# 运行UMAP（基于最佳PC数，假设为1:15）
seurat_obj_filtered <- RunUMAP(
  object = seurat_obj_filtered,
  reduction = reduction_use,
  dims = 1:15, # 最佳PC数，需根据Elbow Plot调整
  verbose = FALSE
)

# 运行tSNE（可选，UMAP效果更稳定）
seurat_obj_filtered <- RunTSNE(
  object = seurat_obj_filtered,
  reduction = reduction_use,
  dims = 1:15,
  verbose = FALSE
)

# 可视化UMAP结果（按默认聚类分组，此时聚类尚未定义，暂按orig.ident分组）
umap_plot <- DimPlot(
  seurat_obj_filtered,
  reduction = "umap",
  group.by = "orig.ident",
  pt.size = 0.3,
  label = FALSE
) + theme_paper()
print(umap_plot)
ggsave("plots/umap_unclustered.pdf", umap_plot, width = 8, height = 6, dpi = 300)

3. 细胞聚类：基于表达相似性分群

Seurat 用 “共享最近邻（SNN）” 算法聚类，核心参数是resolution（分辨率），决定聚类颗粒度：

# 聚类分析（resolution=0.8，适合2000-3000个细胞）
seurat_obj_filtered <- FindNeighbors(
  object = seurat_obj_filtered,
  reduction = reduction_use,
  dims = 1:15,
  verbose = FALSE
) %>%
  FindClusters(
    resolution = 0.8 # 关键参数：分辨率越高，聚类数越多（推荐范围0.6-1.2）
  )

# 查看聚类结果（细胞数分布）
table(seurat_obj_filtered$seurat_clusters)
# 输出示例：
#  0   1   2   3   4   5   6   7   8 
# 320 289 256 231 210 198 187 176 151 

# 可视化UMAP聚类图（按聚类分群着色）
umap_clustered_plot <- DimPlot(
  seurat_obj_filtered,
  reduction = "umap",
  group.by = "seurat_clusters",
  pt.size = 0.3,
  label = TRUE, # 显示聚类编号
  label.size = 4
) + theme_paper()
print(umap_clustered_plot)
ggsave("plots/umap_clustered.pdf", umap_clustered_plot, width = 8, height = 6, dpi = 300)

分辨率选择技巧：

• 细胞数 < 1000：resolution=0.4-0.6；
• 细胞数 1000-5000：resolution=0.6-1.2；
• 细胞数 > 5000：resolution=1.0-1.5；
• 若聚类过细（如 20 + 群），降低分辨率；若聚类过粗（如 < 5 群），提高分辨率。

七、细胞类型注释：从聚类到生物学意义

聚类得到的 “数字编号群”（如 Cluster 0、Cluster 1）需要赋予生物学意义 —— 即注释为已知细胞类型（如 T 细胞、巨噬细胞、上皮细胞），核心是 “基于标记基因验证”。

1. 第一步：识别每个聚类的差异标记基因

用FindAllMarkers筛选每个聚类相对于其他聚类的高表达基因（差异标记基因）：

# 筛选每个聚类的标记基因（log2FC>1，padj<0.05）
cluster_markers <- FindAllMarkers(
  object = seurat_obj_filtered,
  only.pos = TRUE, # 仅保留上调基因
  min.pct = 0.25, # 至少在25%的细胞中表达
  logfc.threshold = 1, # 最小log2折叠变化
  test.use = "wilcox" # 检验方法（单细胞常用Wilcoxon秩和检验）
)

# 查看每个聚类的Top5标记基因
top_markers <- cluster_markers %>%
  group_by(cluster) %>%
  slice_head(n = 5) %>%
  ungroup()
print(top_markers[, c("cluster", "gene", "avg_log2FC", "padj")])

# 保存所有标记基因到CSV
write.csv(cluster_markers, "results/cluster_markers.csv", row.names = FALSE)

2. 第二步：基于标记基因手动注释

根据已知的细胞类型标记基因（如下表），匹配每个聚类的 Top 标记基因，手动注释细胞类型：

细胞类型	人类标记基因	小鼠标记基因
T 细胞	CD3D、CD3E、TRAC	Cd3d、Cd3e、Trac
B 细胞	MS4A1（CD20）、CD79A、IGHM	Ms4a1、Cd79a、Ighm
巨噬细胞	CD68、CSF1R、AIF1	Cd68、Csf1r、Aif1
上皮细胞	EPCAM、KRT19、KRT8	EpCAM、Krt19、Krt8
内皮细胞	CD31（PECAM1）、VWF、CLDN5	Cd31、Vwf、Cldn5
成纤维细胞	COL1A1、COL3A1、DCN	Col1a1、Col3a1、Dcn

手动注释代码：

# 创建细胞类型注释字典（根据你的标记基因结果调整）
cell_type_annot <- c(
  "0" = "T细胞",
  "1" = "巨噬细胞",
  "2" = "上皮细胞",
  "3" = "B细胞",
  "4" = "内皮细胞",
  "5" = "成纤维细胞",
  "6" = "NK细胞",
  "7" = "单核细胞",
  "8" = "树突状细胞"
)

# 给Seurat对象添加细胞类型注释
seurat_obj_filtered$cell_type <- unname(cell_type_annot[as.character(seurat_obj_filtered$seurat_clusters)])

# 可视化注释后的UMAP图
umap_celltype_plot <- DimPlot(
  seurat_obj_filtered,
  reduction = "umap",
  group.by = "cell_type",
  pt.size = 0.3,
  label = TRUE,
  label.size = 3,
  repel = TRUE # 避免标签重叠
) + theme_paper()
print(umap_celltype_plot)
ggsave("plots/umap_celltype_annotated.pdf", umap_celltype_plot, width = 10, height = 8, dpi = 300)

3. 第三步：基于参考数据集自动注释（SingleR）

手动注释依赖标记基因知识，容易出错，用 SingleR 基于公共参考数据集自动注释，交叉验证结果：

# 加载人类参考数据集（HumanPrimaryCellAtlasData，含多种原代细胞）
ref_data <- HumanPrimaryCellAtlasData()

# 提取Seurat对象的表达矩阵（标准化后的数据）
expr_data <- GetAssayData(seurat_obj_filtered, assay = "RNA", slot = "scale.data")

# 运行SingleR自动注释（按聚类注释，速度更快）
singleR_annot <- SingleR(
  test = expr_data,
  ref = ref_data,
  labels = ref_data$label.main, # 参考数据集的细胞类型标签
  clusters = seurat_obj_filtered$seurat_clusters, # 按聚类注释
  assay.type.test = "logcounts",
  assay.type.ref = "logcounts"
)

# 查看自动注释结果（每个聚类的预测细胞类型）
print(singleR_annot$labels)

# 将自动注释结果添加到Seurat对象
seurat_obj_filtered$cell_type_singleR <- unname(singleR_annot$labels[as.character(seurat_obj_filtered$seurat_clusters)])

# 对比手动注释和自动注释结果
annot_compare <- table(seurat_obj_filtered$cell_type, seurat_obj_filtered$cell_type_singleR)
print(annot_compare)

注释避坑：

• 若两种方法结果差异大，重新检查标记基因（如是否混淆人类 / 小鼠基因）；
• 罕见细胞类型可能无匹配参考数据，需依赖手动注释；
• 避免过度依赖自动注释，参考数据集可能不包含你的研究组织特异性细胞。

八、核心可视化：从分群到功能验证

可视化的核心是 “清晰展示生物学结论”，以下是 scRNA-seq 分析必做的 6 类图表，覆盖分群、标记基因、细胞比例等核心信息。

1. 细胞分群可视化（UMAP/tSNE）

除了基础分群图，可按样本、批次、细胞类型分别可视化，展示数据分布：

# 按样本+细胞类型联合可视化
umap_sample_celltype <- DimPlot(
  seurat_obj_filtered,
  reduction = "umap",
  group.by = c("orig.ident", "cell_type"),
  pt.size = 0.3,
  split.by = "orig.ident" # 按样本拆分图
) + theme_paper()
ggsave("plots/umap_sample_celltype.pdf", umap_sample_celltype, width = 12, height = 8, dpi = 300)

2. 标记基因表达可视化（小提琴图 + 气泡图）

验证细胞类型注释的准确性，展示标记基因在不同细胞群的表达：

# 选择关键标记基因（覆盖主要细胞类型）
marker_genes <- c("CD3D", "MS4A1", "CD68", "EPCAM", "PECAM1", "COL1A1")

# 1. 小提琴图（展示单个基因在各细胞群的表达分布）
vln_marker_plot <- VlnPlot(
  seurat_obj_filtered,
  features = marker_genes,
  group.by = "cell_type",
  pt.size = 0.1,
  ncol = 2
) + theme_paper()
ggsave("plots/marker_genes_vln.pdf", vln_marker_plot, width = 12, height = 10, dpi = 300)

# 2. 气泡图（展示多个基因在各细胞群的表达水平和比例）
dot_marker_plot <- DotPlot(
  seurat_obj_filtered,
  features = marker_genes,
  group.by = "cell_type",
  dot.scale = 8
) +
  theme_paper() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # x轴基因名旋转45度
ggsave("plots/marker_genes_dot.pdf", dot_marker_plot, width = 10, height = 8, dpi = 300)

3. 细胞类型比例可视化（堆叠柱状图）

比较不同样本 / 处理组的细胞类型比例差异：

# 计算每个样本的细胞类型比例
cell_type_ratio <- seurat_obj_filtered@meta.data %>%
  group_by(orig.ident, cell_type) %>%
  summarise(count = n(), .groups = "drop_last") %>%
  mutate(ratio = count / sum(count) * 100)

# 绘制堆叠柱状图
ratio_bar_plot <- ggplot(
  cell_type_ratio,
  aes(x = orig.ident, y = ratio, fill = cell_type)
) +
  geom_col(position = "stack", width = 0.7) +
  labs(x = "Sample", y = "Cell Type Ratio (%)", fill = "Cell Type") +
  theme_paper() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("plots/cell_type_ratio.pdf", ratio_bar_plot, width = 10, height = 6, dpi = 300)

4. 标记基因热图（Top 标记基因）

展示每个细胞群的 Top5 标记基因，直观反映聚类特异性：

# 提取每个聚类的Top5标记基因
top5_markers_per_cluster <- cluster_markers %>%
  group_by(cluster) %>%
  slice_head(n = 5) %>%
  pull(gene) %>%
  unique()

# 绘制热图（按细胞类型分组，行=基因，列=细胞）
heatmap_markers <- DoHeatmap(
  seurat_obj_filtered,
  features = top5_markers_per_cluster,
  group.by = "cell_type",
  size = 3, # 基因名字体大小
  angle = 0, # 基因名水平显示
  raster = FALSE # 非光栅化，放大不失真
) +
  scale_fill_viridis(option = "plasma") +
  theme_paper() +
  theme(legend.position = "right")
ggsave("plots/markers_heatmap.pdf", heatmap_markers, width = 12, height = 10, dpi = 300)

5. 功能富集可视化（GO/KEGG）

对关键细胞群（如差异细胞类型）做功能富集，揭示生物学功能：

# 以巨噬细胞（cell_type="巨噬细胞"）为例，提取高表达基因
macrophage_genes <- cluster_markers %>%
  filter(cluster == 1) %>% # 假设巨噬细胞是Cluster 1
  pull(gene)

# GO生物学过程富集分析
go_enrich <- enrichGO(
  gene = macrophage_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP", # 生物学过程
  pAdjustMethod = "fdr",
  qvalueCutoff = 0.05
)

# 可视化Top10富集通路（柱状图）
go_bar_plot <- barplot(
  go_enrich,
  showCategory = 10,
  title = "Macrophage GO Biological Process Enrichment"
) + theme_paper()
ggsave("plots/macrophage_go_enrich.pdf", go_bar_plot, width = 10, height = 6, dpi = 300)

九、结果保存与导出

分析完成后，需保存 Seurat 对象（便于后续重分析）、关键结果表格和可视化图表，形成完整的分析报告。

1. 保存 Seurat 对象

# 保存完整的Seurat对象（包含所有分析结果）
saveRDS(seurat_obj_filtered, "results/seurat_final.rds")

# 加载已保存的对象（后续分析直接复用）
# seurat_obj <- readRDS("results/seurat_final.rds")

2. 导出关键结果表格

# 导出细胞元数据（包含细胞类型、聚类、批次等信息）
write.csv(seurat_obj_filtered@meta.data, "results/cell_metadata.csv", row.names = TRUE)

# 导出差异标记基因
write.csv(cluster_markers, "results/cluster_markers.csv", row.names = FALSE)

# 导出细胞类型比例
write.csv(cell_type_ratio, "results/cell_type_ratio.csv", row.names = FALSE)

3. 导出表达矩阵（用于后续分析）

# 导出原始计数矩阵
counts_raw <- GetAssayData(seurat_obj_filtered, assay = "RNA", slot = "counts")
write.csv(as.matrix(counts_raw), "results/expression_raw.csv", row.names = TRUE)

# 导出归一化后的表达矩阵
counts_normalized <- GetAssayData(seurat_obj_filtered, assay = "RNA", slot = "data")
write.csv(as.matrix(counts_normalized), "results/expression_normalized.csv", row.names = TRUE)

十、完整代码模板（直接复用）

将上述步骤整合为完整代码，只需修改数据路径、质控阈值、标记基因等参数即可运行：

# 单细胞RNA-seq完整分析代码模板
# 作者：生信实操手册
# 适用场景：10x Genomics人类scRNA-seq数据（可调整为小鼠）

# -------------------------- 1. 加载包 --------------------------
library(Seurat)
library(SingleR)
library(Harmony)
library(clusterProfiler)
library(org.Hs.eg.db)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(patchwork)

theme_paper <- function() {
  theme_bw() +
    theme(
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold", family = "Arial"),
      axis.title = element_text(size = 10, family = "Arial"),
      axis.text = element_text(size = 8, family = "Arial"),
      legend.title = element_text(size = 9, family = "Arial"),
      legend.text = element_text(size = 8, family = "Arial"),
      panel.grid = element_blank()
    )
}

# -------------------------- 2. 数据导入与对象构建 --------------------------
data_dir <- "data/filtered_feature_bc_matrix" # 替换为你的数据路径
counts <- Read10X(data.dir = data_dir)
seurat_obj <- CreateSeuratObject(
  counts = counts,
  project = "scRNA_Sample",
  min.cells = 3,
  min.features = 200
)

# 添加元数据
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")
seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# -------------------------- 3. 质控 --------------------------
seurat_obj_filtered <- subset(
  seurat_obj,
  subset = nFeature_RNA > 200 & nFeature_RNA < 6000 &
    percent.mt < 5 & percent.ribo > 10
)

# -------------------------- 4. 数据预处理 --------------------------
seurat_obj_filtered <- NormalizeData(seurat_obj_filtered)
seurat_obj_filtered <- FindVariableFeatures(seurat_obj_filtered, nfeatures = 2000)
seurat_obj_filtered <- ScaleData(seurat_obj_filtered, vars.to.regress = "percent.mt")

# -------------------------- 5. 批次校正（按需开启） --------------------------
# seurat_obj_filtered$batch <- c(rep("batch1", 1000), rep("batch2", 1518)) # 手动添加批次
# seurat_obj_filtered <- RunPCA(seurat_obj_filtered, verbose = FALSE)
# seurat_obj_filtered <- RunHarmony(seurat_obj_filtered, group.by.vars = "batch", verbose = FALSE)
reduction_use <- ifelse("harmony" %in% Reductions(seurat_obj_filtered), "harmony", "pca")

# -------------------------- 6. 降维聚类 --------------------------
seurat_obj_filtered <- RunPCA(seurat_obj_filtered, verbose = FALSE)
seurat_obj_filtered <- RunUMAP(seurat_obj_filtered, reduction = reduction_use, dims = 1:15)
seurat_obj_filtered <- FindNeighbors(seurat_obj_filtered, reduction = reduction_use, dims = 1:15) %>%
  FindClusters(resolution = 0.8)

# -------------------------- 7. 细胞类型注释 --------------------------
# 手动注释（替换为你的标记基因对应的聚类）
cell_type_annot <- c(
  "0" = "T细胞", "1" = "巨噬细胞", "2" = "上皮细胞",
  "3" = "B细胞", "4" = "内皮细胞", "5" = "成纤维细胞",
  "6" = "NK细胞", "7" = "单核细胞", "8" = "树突状细胞"
)
seurat_obj_filtered$cell_type <- unname(cell_type_annot[as.character(seurat_obj_filtered$seurat_clusters)])

# -------------------------- 8. 可视化 --------------------------
# UMAP分群图
umap_celltype <- DimPlot(seurat_obj_filtered, reduction = "umap", group.by = "cell_type", label = TRUE) + theme_paper()
ggsave("plots/umap_celltype.pdf", umap_celltype, width = 10, height = 8, dpi = 300)

# 标记基因气泡图
marker_genes <- c("CD3D", "MS4A1", "CD68", "EPCAM", "PECAM1", "COL1A1")
dot_plot <- DotPlot(seurat_obj_filtered, features = marker_genes, group.by = "cell_type") + theme_paper()
ggsave("plots/marker_dot.pdf", dot_plot, width = 10, height = 8, dpi = 300)

# -------------------------- 9. 结果保存 --------------------------
saveRDS(seurat_obj_filtered, "results/seurat_final.rds")
write.csv(seurat_obj_filtered@meta.data, "results/cell_metadata.csv", row.names = TRUE)

十一、常见问题排查与避坑总结

1. 聚类模糊、UMAP 图重叠严重：

   • 解决方案：提高测序深度、调整质控阈值（如降低线粒体比例）、增加高变基因数量（nfeatures=3000）、调整 UMAP 的dims参数。

2. 细胞类型注释不准确：

   • 解决方案：交叉验证手动注释和 SingleR 结果、补充组织特异性标记基因、更换参考数据集（如小鼠用 MouseRNAseqData）。

3. 批次效应无法消除：

   • 解决方案：检查批次变量是否正确、增加 PCA 维度（dims=1:20）、使用 sctransform 标准化后再做 Harmony 校正。

4. 运行内存不足：

   • 解决方案：过滤低质量细胞（保留核心细胞）、使用稀疏矩阵（避免转为普通矩阵）、分批次处理（如先分析单个样本）。

总结

scRNA-seq 分析的核心是 “去噪→提特征→聚类型→释意义”，本文基于 Seurat 构建的全流程，覆盖了从原始数据到可视化报告的每一个关键步骤，且所有代码可直接复用。新手可先按模板运行，再根据自身数据（如小鼠、肿瘤组织）调整参数；进阶用户可扩展分析（如细胞通讯分析、拟时序分析）。

可视化是 scRNA-seq 结果的 “门面”，务必保证图表符合论文规范（矢量图、配色友好、标注清晰）。若能熟练掌握本文流程，可应对 80% 以上的 scRNA-seq 分析场景，为后续的机制研究和论文发表打下基础。