单细胞转录组 | 差异基因分析和富集分析

# 加载必要的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))

# 设置感兴趣的细胞类型
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秩和检验

# 创建结果存储列表
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)

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)

# 加载包
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)

# 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)