2024-06-16

超详细的R中进行GSEA分析的流程！一看就会！

小花生信果 2024-05-27 19:00:09

转自公众号：生信花
http://mp.weixin.qq.com/s?__biz=MzkzODUxNjcwMA==&mid=2247483861&idx=1&sn=8ba2b9d488e296355533f99127860627

在日常研究中，有时需要用到GSEA分析来识别实验组和对照组之间的差异调节通路。GSEA（基因集富集分析）是用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献。它的基本思想是使用预定义的基因集（通常来自功能注释或先前实验的结果），将基因按照在两个分组中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。相较于GO和KEGG这种更加依赖差异基因的分析（实则是对一部分基因的分析，忽略差异不显著的基因)，GSEA则是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集，故能将那些GO/KEGG富集分信息中容易遗漏掉的差异表达不显著却有着重要生物学意义的基因包含在内，从而为分析提供更多有意义的数据信息。

定制生信分析

云服务器租赁

（加微信备注99领取使用）

GSEA官网提供了一个本地分析工具，但是这个工具在使用的时候准备输入文件比较麻烦，而且出图也是一个通路一张图，那有没有一种更好的方法呢？小果今天就带着大家来学习一下用R语言开展GSEA分析吧！

第一步：输入文件准备

这里我们可以直接使用limma差异分析后的结果文件，我们只需要前两列内容就可以了。如果不知道怎么使用limma进行差异分析的可以看一下小果以前的推文。

第二步：配置所需的R包

 所需R包均可以用BiocManager包安装，如果没有BiocManager包，可以用以下命令安装。

if (!requireNamespace("BiocManager", quietly = TRUE))    install.packages("BiocManager")    #安装BiocManager包
BiocManager::install("clusterProfiler")         #安装分析所需R包library(ReactomePA)    #加载所需R包library(tidyverse)library(data.table)library(org.Hs.eg.db)library(clusterProfiler)library(biomaRt)library(enrichplot)

第三步：导入数据并进行预处理

在进行GSEA分析之前我们还需要对输入文件进行一个预处理，主要就是将文件中第一列SYMBOL格式的ID换成ENTREZ格式的ID，同时将ENTREZID与logFC结合，并根据logFC的值进行降序排列，最后将结果整理成GSEA分析所需的格式，就可以用于后续分析了。下面跟着小果一起来操作一下吧！

setwd("G:/job/job1/GSEA/")   #进入存放数据的工作目录genelist_input <- fread(file="GSEA.txt", header = T, sep='t', data.table = F)   #导入数据genename <- as.character(genelist_input[,1])   #提取第一列基因名gene_map <- select(org.Hs.eg.db, keys=genename, keytype="SYMBOL", columns=c("ENTREZID"))  #将SYMBOL格式的ID换成ENTREZ格式的ID。

（注意：这一步运行结果正常情况就是1:1 mapping，如图所示。但是也有会出现1：many mapping的情况，如图示，这种情况可能会影响后续分析结果，同时有些分析也无法进行，因此我们可以进行去重复操作，再继续后续分析。）

non_duplicates_idx <- which(duplicated(gene_map$SYMBOL) == FALSE)  gene_map <- gene_map[non_duplicates_idx, ]   #去除重复值colnames(gene_map)[1]<-"Gene"  #输出mapping结果

#将ENTREZID与logFC对应起来，并根据logFC的值降序排列，最终生成结果如图所示。temp<-inner_join(gene_map,genelist_input,by = "Gene")temp<-temp[,-1]temp<-na.omit(temp)temp$logFC<-sort(temp$logFC,decreasing = T)

#最后我们将文件内容整理成GSEA分析所需的格式（如图所示），我们就可以开始GSEA分析了！geneList = temp[,2]names(geneList) = as.character(temp[,1]) geneList

第四步：进行GSEA分析，并保存结果。

Go_gseresult <- gseGO(geneList, 'org.Hs.eg.db', keyType = "ENTREZID", ont="BP", pvalueCutoff=0.05)   #使用GSEA进行GO富集分析（'org.Hs.eg.db'：对应物种的数据库；ont：选择输出条目，可选“BP,MF,CC或者ALL”，pvalueCutoff：设置P的阈值）KEGG_gseresult <- gseKEGG(geneList, pvalueCutoff=0.05)  #使用GSEA进行KEGG富集分析    #保存富集分析结果go_results<-as.data.frame(Go_gseresult)kegg_results<-as.data.frame(KEGG_gseresult)write.csv (go_results, file ="Go_gseresult.csv")write.csv (kegg_results, file ="KEGG_gseresult.csv")

第五步：结果可视化（以GO富集结果为例）

这里我们使用的是gseaplot2包进行结果可视化分析。它既可以输出指定单基因的图片，也可以输出包含多条基因的图片。

gseaplot2(Go_gseresult, 1:3, title = "Specific GO Biological Process in T2D group", pvalue_table = FALSE)  #1:3：这表示在图上显示前3条富集结果，也可以根据自己分析需要指定输出某一条结果；Go_gseresult：GO富集分析结果；title：加上标题；pvalue_table：是否在图上显示P值列表。