生信云服务器实战：单细胞转录组差异表达与功能富集分析全流程指南

本文将系统介绍如何利用 Seurat 包对单细胞转录组数据进行差异表达分析，并结合 clusterProfiler 进行功能富集分析，帮助研究者深入理解细胞类型特异性的基因表达变化。

一、单细胞与常规转录组差异分析的核心区别

单细胞转录组（scRNA-seq）与传统的混合样本转录组（bulk RNA-seq）在差异分析方面存在显著差异：

• 分析粒度：传统转录组分析的是细胞群体的平均表达水平，可能掩盖个体细胞的异质性；而单细胞转录组能够揭示单个细胞的表达差异，捕捉细胞间的多样性。

• 数据特性：单细胞数据通常存在高比例的零值（dropout），数据稀疏且噪声较大，需采用专门的方法进行处理。

• 分析流程：传统方法如 DESeq2 和 edgeR 直接比较组间差异，而单细胞分析需先进行细胞聚类，再在特定细胞类型中进行差异分析，常用方法包括 Wilcoxon 秩和检验、MAST 等。

• 批次效应处理：单细胞数据的批次效应更为复杂，常用 Harmony、Seurat 的整合功能等方法进行校正。
  
  

二、差异表达分析实操流程（基于 Seurat）

以下以 Seurat 处理的胰腺单细胞数据集（panc8）为例，比较 “celseq2” 与 “smartseq2” 两种测序技术在 acinar 细胞类型中的差异表达。

1. 环境准备与数据加载

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

# 加载 Seurat 对象
seurat_obj <- readRDS("path/to/your_seurat_object.rds")

# 查看各组中细胞类型的分布
table(seurat_obj$cell_type, seurat_obj$group)

2. 设置分析参数

# 指定感兴趣的细胞类型
cell_types_of_interest <- c("acinar")

# 定义分组变量及比较组
group_var <- "group"
ident.1 <- "celseq2"
ident.2 <- "smartseq2"

# 设置差异分析参数
min.pct <- 0.1
logfc.threshold <- 0.25
test.use <- "wilcox"

3. 提取特定细胞类型并进行差异分析

# 提取 acinar 细胞
acinar_cells <- subset(seurat_obj, cell_type == "acinar")

# 进行差异表达分析
diff_acinar <- FindMarkers(acinar_cells,
                           group.by = group_var,
                           ident.1 = ident.1,
                           ident.2 = ident.2,
                           min.pct = min.pct,
                           logfc.threshold = logfc.threshold,
                           test.use = test.use)

# 保存结果
write.csv(diff_acinar, "DEG_acinar_celseq2_vs_smartseq2.csv", row.names = FALSE)

4. 绘制火山图

# 筛选显著差异基因
sig_genes <- diff_acinar[diff_acinar$p_val_adj < 0.05 & abs(diff_acinar$avg_log2FC) >= 0.5, ]

# 绘制火山图
p <- EnhancedVolcano(diff_acinar,
                     lab = rownames(diff_acinar),
                     x = "avg_log2FC",
                     y = "p_val_adj",
                     pCutoff = 0.05,
                     FCcutoff = 0.5,
                     pointSize = 2.0,
                     labSize = 3.0,
                     title = "DEGs in Acinar Cells: celseq2 vs smartseq2",
                     subtitle = "Wilcoxon Test",
                     caption = paste("Total DEGs:", nrow(sig_genes)),
                     legendPosition = "right")
print(p)
ggsave("volcano_acinar_celseq2_vs_smartseq2.pdf", p, width = 10, height = 8)

三、功能富集分析（GO 和 KEGG）

在获得显著差异基因后，可进一步进行功能富集分析，以揭示这些基因在生物过程和通路中的潜在作用。

library(clusterProfiler)
library(org.Hs.eg.db)

# 提取显著基因的符号
sig_gene_symbols <- rownames(sig_genes)

# 转换为 Entrez ID
sig_entrez <- bitr(sig_gene_symbols,
                   fromType = "SYMBOL",
                   toType = "ENTREZID",
                   OrgDb = org.Hs.eg.db)

# GO 富集分析（生物过程）
go_enrich <- enrichGO(gene = sig_entrez$ENTREZID,
                      OrgDb = org.Hs.eg.db,
                      ont = "BP",
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.2)

# KEGG 通路富集分析
kegg_enrich <- enrichKEGG(gene = sig_entrez$ENTREZID,
                          organism = "hsa",
                          pAdjustMethod = "BH",
                          pvalueCutoff = 0.05)

四、总结

通过上述流程，研究者可以系统地对单细胞转录组数据进行差异表达分析和功能富集分析，深入挖掘特定细胞类型在不同条件下的基因表达变化及其生物学意义。