单细胞转录组：差异基因分析和富集分析
一、单细胞与普通转录组差异基因分析的区别
单细胞转录组(scRNA-seq)与普通转录组(bulk RNA-seq)在差异基因分析方面存在本质区别：
1. 分析单位不同：
普通转录组：以组织或细胞混合物为单位，结果反映的是平均表达水平
单细胞转录组：以单个细胞为单位，能够揭示细胞间的异质性
2. 数据特性不同：
普通转录组：信号强度高，噪音相对较低
单细胞转录组：存在大量零值(dropout)，数据稀疏性高，技术噪音大
3. 分析方法不同：
普通转录组：DESeq2、edgeR等方法直接比较组间差异
单细胞转录组：需要先分群，再进行细胞类型特异性分析，可使用MAST、Wilcoxon等特殊方法
4. 解析能力不同：
普通转录组：无法区分细胞类型特异性变化
单细胞转录组：可以揭示特定细胞类型的变化，避免少数细胞类型的信号被稀释
5. 批次效应处理不同：
单细胞数据批次效应更为复杂，需要特殊的整合方法(如Harmony、Seurat整合)
二、实践操作：基于Seurat进行单细胞转录组差异基因分析
以下是针对已完成Seurat标准流程并注释好细胞类型的数据，进行不同组间同一细胞类型差异基因分析的完整流程：
我使用的数据是Seurat处理胰腺单细胞转录组数据（panc8数据集），把”celseq2”和“smartseq2”两种技术之间的差异作为分组。所以seurat_obj$group在本次处理过程就是seurat_obj$tech。
1. 准备工作：加载包和数据
# 加载必要的R包
library(Seurat)
library(dplyr)
library(ggplot2)
library(EnhancedVolcano)# 加载处理好的Seurat对象
seurat_obj <- readRDS("path/to/your_seurat_object.rds")# 查看数据结构
print(table(seurat_obj$cell_type, seurat_obj$group))
这一步加载了必要的R包和Seurat对象。table函数帮助我们检查每个细胞类型在不同组中的细胞数量，确保有足够的细胞进行比较。
2. 设置差异分析参数
# 设置感兴趣的细胞类型
cell_types_of_interest <- c("acinar") # 这里仅分析acinar，可扩展到其他类型# 设置分组变量和比较组
group_var <- "group" # 元数据中的分组变量名
ident.1 <- "treatment" # 处理组名称
ident.2 <- "control" # 对照组名称# 设置差异分析参数
min.pct <- 0.1 # 至少在一个组中有10%的细胞表达
logfc.threshold <- 0.25 # 对数倍变化阈值
test.use <- "wilcox" # 使用Wilcoxon秩和检验
这一步设置了分析参数，包括感兴趣的细胞类型、分组变量和比较组，以及差异分析的参数。
3.处理选中的细胞类型
# 创建结果存储列表
deg_results <- list()#DEGs
Fibroblast <- subset(seurat_obj, celltype=="acinar")diff_Fibroblast <- FindMarkers(Fibroblast,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_Fibroblast, "DEG_acinar_celseq2_vs_smartseq2.csv", row.names = FALSE)
这一步根据前面设置感兴趣的细胞类型，提取子集，进行差异基因分析，并将结果存储在列表中。
4.绘制火山图
sig_genes <- diff_Fibroblast[diff_Fibroblast$p_val_adj < 0.05 & abs(diff_Fibroblast$avg_log2FC) >= 0.5,]
p <- EnhancedVolcano(diff_Fibroblast,lab = rownames(diff_Fibroblast),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:", length(sig_genes)),legendPosition = "right",
#selectLab = sig_genes # 只标记显著基因
)
print(p)
ggsave("volcano_acinar_celseq2_vs_smartseq2.pdf", p, width = 10, height = 8)
这一步为指定的细胞类型绘制火山图，用于可视化差异基因的分布。火山图x轴表示对数倍变化，y轴表示调整后的p值。
5.筛选显著差异基因并进行功能富集分析
# 加载包
library(clusterProfiler)
library(org.Hs.eg.db)#筛选显著差异基因
#sig_genes <- diff_Fibroblast[diff_Fibroblast$p_val_adj < 0.05 & abs(diff_Fibroblast$avg_log2FC) >= 0.5,]# 从 sig_genes 中提取基因符号（行名）
sig_gene_symbols <- rownames(sig_genes)# 将基因符号转换为 Entrez ID
sig_entrez <- bitr(sig_gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)# 查看转换结果
head(sig_entrez)
这一步筛选显著差异基因，因为在绘制火山图之前我已经执行过这一步，这里我注释掉了。可以查看一下sig_genes的情况，提取之后还有5863个基因。
富集分析：
# 4. 进行 GO 富集分析（Biological Process）
go_enrich <- enrichGO(gene = sig_entrez$ENTREZID,OrgDb = org.Hs.eg.db,ont = "BP", # Biological Process，可改为 "MF" 或 "CC"pAdjustMethod = "BH", # Benjamini-Hochberg 校正pvalueCutoff = 0.05,qvalueCutoff = 0.2)# 查看 GO 富集结果前几行
head(go_enrich)#进行 KEGG 通路富集分析
kegg_enrich <- enrichKEGG(gene = sig_entrez$ENTREZID,organism = "hsa", # 人类pAdjustMethod = "BH",pvalueCutoff = 0.05,qvalueCutoff = 0.2)# 查看 KEGG 富集结果前几行
head(kegg_enrich)# 保存富集分析结果为 CSV 文件
write.csv(as.data.frame(go_enrich), "go_enrichment.csv", row.names = FALSE)
write.csv(as.data.frame(kegg_enrich), "kegg_enrichment.csv", row.names = FALSE)# 可视化（可选）
# GO 富集结果条形图
barplot(go_enrich, showCategory = 10) # 显示前 10 个显著类别
ggsave("go_enrichment_barplot.pdf", width = 10, height = 8)# KEGG 富集结果点图
dotplot(kegg_enrich, showCategory = 10) # 显示前 10 个显著通路
ggsave("kegg_enrichment_dotplot.pdf", width = 10, height = 8)
注意事项
确保Seurat对象中的
cell_type和group变量名与自己的数据一致，必要时调整代码中的变量名根据数据特点调整参数，如
min.pct、logfc.threshold和test.use单细胞数据中，每个细胞类型的细胞数量对结果可靠性有重要影响，建议在分析前检查细胞数量
考虑使用特定于单细胞数据的差异分析方法，如MAST
对于稀有细胞类型，可能需要调整参数或使用伪时间分析（pseudobulk）方法
通过以上步骤，可以系统地分析不同组间特定细胞类型的差异基因表达，并通过火山图直观展示结果。这种细胞类型特异性的分析是单细胞转录组技术的独特优势，可以揭示在普通转录组分析中可能被掩盖的细胞类型特异性变化。当然，你也可以选择分析不同类群（亚群）间的差异基因表达。后面富集分析参照普通转录组的做法。