0
0

差异基因富集分析（R语言——GO&KEGG&GSEA）

接着上次的内容，上篇内容给大家分享了基因表达量怎么做分组差异分析，从而获得差异基因集，想了解的可以去看一下，这篇主要给大家分享一下得到显著差异基因集后怎么做一下通路富集。

1.准备差异基因集

我就直接把上次分享的拿到这边了。我们一般都把差异基因分为上调基因和下调基因分别做通路富集分析。下面上代码，可能包含我的一些个人习惯，勿怪。显著差异基因的筛选条件根据个人需求设置哈。

  1. ##载入所需R包
  2. library(readxl)
  3. library(DOSE)
  4. library(org.Hs.eg.db)
  5. library(topGO)
  6. library(pathview)
  7. library(ggplot2)
  8. library(GSEABase)
  9. library(limma) 
  10. library(clusterProfiler)
  11. library(enrichplot)
  13. ##edger
  14. edger_diff <- diff_gene_Group
  15. edger_diff_up <- rownames(edger_diff[which(edger_diff$logFC > 0.584962501),])
  16. edger_diff_down <- rownames(edger_diff[which(edger_diff$logFC < -0.584962501),])
  18. ##deseq2
  19. deseq2_diff <- diff_gene_Group2
  20. deseq2_diff_up <- rownames(deseq2_diff[which(deseq2_diff$log2FoldChange > 0.584962501),])
  21. deseq2_diff_down <- rownames(deseq2_diff[which(deseq2_diff$log2FoldChange < -0.584962501),])
  23. ##将差异基因集保存为一个list
  24. gene_diff_edger_deseq2 <- list()
  25. gene_diff_edger_deseq2[["edger_diff_up"]] <- edger_diff_up
  26. gene_diff_edger_deseq2[["edger_diff_down"]] <- edger_diff_down
  27. gene_diff_edger_deseq2[["deseq2_diff_up"]] <- deseq2_diff_up
  28. gene_diff_edger_deseq2[["deseq2_diff_down"]] <- deseq2_diff_down

2.进行通路富集分析

这里主要介绍普通的GO&KEGG&GSEA的简单富集。筛选显著富集通路的筛选条件也是根据自己的需求决定，一般是矫正后P值小于0.05。我这里是省事，写了各list循环。

  1. for (i in 1:length(gene_diff_edger_deseq2)){
  2.   keytypes(org.Hs.eg.db)
  3.   
  4.   entrezid_all = mapIds(x = org.Hs.eg.db,  
  5.                         keys = gene_diff_edger_deseq2[[i]], 
  6.                         keytype = "SYMBOL", #输入数据的类型
  7.                         column = "ENTREZID")#输出数据的类型
  8.   entrezid_all  = na.omit(entrezid_all)  #na省略entrezid_all中不是一一对应的数据情况
  9.   entrezid_all = data.frame(entrezid_all) 
  10.   
  11.   ##GO富集##
  12.   GO_enrich = enrichGO(gene = entrezid_all[,1],
  13.                        OrgDb = org.Hs.eg.db, 
  14.                        keyType = "ENTREZID", #输入数据的类型
  15.                        ont = "ALL", #可以输入CC/MF/BP/ALL
  16.                        #universe = 背景数据集我没用到它。
  17.                        pvalueCutoff = 1,qvalueCutoff = 1, #表示筛选的阈值，阈值设置太严格可导致筛选不到基因。可指定 1 以输出全部
  18.                        readable = T) #是否将基因ID映射到基因名称。
  19.   
  20.   GO_enrich_data  = data.frame(GO_enrich)
  21.   write.csv(GO_enrich_data,paste('GO_enrich_',names(gene_diff_edger_deseq2)[i], '.csv', sep = ""))
  22.   GO_enrich_data <- GO_enrich_data[which(GO_enrich_data$p.adjust < 0.05),]
  23.   write.csv(GO_enrich_data,paste('GO_enrich_',names(gene_diff_edger_deseq2)[i], '_filter.csv', sep = ""))
  24.   
  25.   
  26.   ###KEGG富集分析###
  27.   KEGG_enrich = enrichKEGG(gene = entrezid_all[,1], #即待富集的基因列表
  28.                            keyType = "kegg",
  29.                            pAdjustMethod = 'fdr',  #指定p值校正方法
  30.                            organism= "human",  #hsa，可根据你自己要研究的物种更改，可在https://www.kegg.jp/brite/br08611中寻找
  31.                            qvalueCutoff = 1, #指定 p 值阈值（可指定 1 以输出全部）
  32.                            pvalueCutoff=1) #指定 q 值阈值（可指定 1 以输出全部）
  33.   KEGG_enrich_data  = data.frame(KEGG_enrich)
  34.   write.csv(KEGG_enrich_data, paste('KEGG_enrich_',names(gene_diff_edger_deseq2)[i], '.csv', sep = ""))
  35.   KEGG_enrich_data <- KEGG_enrich_data[which(KEGG_enrich_data$p.adjust < 0.05),]
  36.   write.csv(KEGG_enrich_data, paste('KEGG_enrich_',names(gene_diff_edger_deseq2)[i], '_filter.csv', sep = ""))
  37. }

3.通路富集情况可视化

这里只介绍一种简单的气泡图，当然还有其他的自己去了解吧。

  1. ##GO&KEGG富集BP\CC\MF\KEGG分面绘图需要分开处理一下，富集结果里的ONTOLOGYL列修改
  2. GO_enrich_data_BP <- subset(GO_enrich_data, subset = GO_enrich_data$ONTOLOGY == "BP")
  3. GO_enrich_data_CC <- subset(GO_enrich_data, subset = GO_enrich_data$ONTOLOGY == "CC")
  4. GO_enrich_data_MF <- subset(GO_enrich_data, subset = GO_enrich_data$ONTOLOGY == "MF")
  6. ##提取GO富集BP\CC\MF的top5
  7. GO_enrich_data_filter <- rbind(GO_enrich_data_BP[1:5,], GO_enrich_data_CC[1:5,], GO_enrich_data_MF[1:5,])
  9. ##重新整合进富集结果
  10. GO_enrich@result <- GO_enrich_data_filter
  12. ##处理KEGG富集结果
  13. KEGG_enrich@result <- KEGG_enrich_data
  14. ncol(KEGG_enrich@result)
  15. KEGG_enrich@result$ONTOLOGY <- "KEGG"
  16. KEGG_enrich@result <- KEGG_enrich@result[,c(10,1:9)]
  18. ##整合GO KEGG富集结果
  19. ego_GO_KEGG <- GO_enrich
  20. ego_GO_KEGG@result <- rbind(ego_GO_KEGG@result, KEGG_enrich@result[1:5,])
  21. ego_GO_KEGG@result$ONTOLOGY <- factor(ego_GO_KEGG@result$ONTOLOGY, levels = c("BP", "CC", "MF","KEGG"))##规定分组顺序
  23. ##简单画图
  24. pdf("edger_diff_up_dotplot.pdf", width = 7, height = 7)
  25. dotplot(ego_GO_KEGG, split = "ONTOLOGY", title="UP-GO&KEGG", label_format = 60, color = "pvalue") + 
  26.   facet_grid(ONTOLOGY~., scale = "free_y")+
  27.   theme(plot.title = element_text(hjust = 0.5, size = 13, face = "bold"), axis.text.x = element_text(angle = 90, hjust = 1))
  28. dev.off()

4.气泡图如图所示

做了些处理，真实图片，左侧pathway是跟后面气泡一一对应的，当然还有其他可视化方式那就需要各位自己去探索了，谢谢！

5.GSEA富集分析

这里也是做一下简单的GSEA

  1. ##GSEA官方网站下载背景gmt文件并读入
  2. geneset <- list()
  3. geneset[["c2_cp"]] <- read.gmt("c2.cp.v2023.2.Hs.symbols.gmt")
  4. geneset[["c2_cp_kegg_legacy"]] <- read.gmt("c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")
  5. geneset[["c2_cp_kegg_medicus"]] <- read.gmt("c2.cp.kegg_medicus.v2023.2.Hs.symbols.gmt")
  6. geneset[["c2_cp_reactome"]] <- read.gmt("c2.cp.reactome.v2023.2.Hs.symbols.gmt")
  7. geneset[["c3_tft"]] <- read.gmt("c3.tft.v2023.2.Hs.symbols.gmt")
  8. geneset[["c4_cm"]] <- read.gmt("c4.cm.v2023.2.Hs.symbols.gmt")
  9. geneset[["c5_go_bp"]] <- read.gmt("c5.go.bp.v2023.2.Hs.symbols.gmt")
  10. geneset[["c5_go_cc"]] <- read.gmt("c5.go.cc.v2023.2.Hs.symbols.gmt")
  11. geneset[["c5_go_mf"]] <- read.gmt("c5.go.mf.v2023.2.Hs.symbols.gmt")
  12. geneset[["c6"]] <- read.gmt("c6.all.v2023.2.Hs.symbols.gmt")
  13. geneset[["c7"]] <- read.gmt("c7.all.v2023.2.Hs.symbols.gmt")
  15. ##进行GSEA富集分析，这里也是写了个循环
  16. gsea_results <- list()
  17. for (i in names(gene_diff)){
  18.   geneList <- gene_diff[[i]]$logFC
  19.   names(geneList) <- toupper(rownames(gene_diff[[i]]))
  20.   geneList <- sort(geneList,decreasing = T)
  21.   for (j in names(geneset)){
  22.     listnames <- paste(i,j,sep = "_")
  23.     gsea_results[[listnames]] <- GSEA(geneList = geneList,
  24.                          TERM2GENE = geneset[[j]],
  25.                          verbose = F,
  26.                          pvalueCutoff = 0.1,
  27.                          pAdjustMethod = "none",
  28.                          eps=0)
  29.   }
  30. }
  32. ##批量绘图，注意这里如果有空富集通路，会报错！
  33. for (j in 1:nrow(gsea_results[[i]]@result)) {
  34.     p <- gseaplot2(x=gsea_results[[i]],geneSetID=gsea_results[[i]]@result$ID[j], title = 
  35.     gsea_results[[i]]@result$ID[j]) 
  36.     
  37.     pdf(paste(paste(names(gsea_results)[i], gsea_results[[i]]@result$ID[j], sep = 
  38.     "_"),".pdf",sep = ""))
  39.     print(p)
  40.     dev.off()
  41.   }

6.GSEA富集最简单图形如下

分享到此结束了，希望对大家有所帮助。

标签： r语言 数据分析 数据挖掘
本文转载自: https://blog.csdn.net/m0_73145613/article/details/140145314
版权归原作者 R语言打工人 所有， 如有侵权，请联系我们删除。
登录后发布评论

“差异基因富集分析（R语言——GO&KEGG&GSEA）”的评论:

还没有评论
关于作者
...
overfit同步小助手
文章同步
相关阅读

网络安全法-网络运行安全

使用selenium/drissionpage时如何阻止chrome自动跳转http到https

docker desktop 里部署的Open WebUI 管理员密码忘记了的处理方法

在ubuntu20.04中搭建onsite比赛运行环境

利用开源的低代码表单设计器FcDesigner高效管理和渲染复杂表单结构

Kafka学习笔记

【前端】浏览器输入url到页面呈现发生了什么？

文章导航
没有找到文章内容分段。。。