GSEA富集分析流程及解释（R代码）【待完善】

GSEA简介

• GSEA法基因功能富集分析原理详解
• 豆豆学习GSEA
• clusterprofiler|GSEA富集分析及可视化

数据准备

GSEA只需要gene id和log2FoldChange两列

library(clusterProfiler)
library(enrichplot)
library(org.Mm.eg.db) # mouse全基因组注释包，用于不同版本基因名的转换
# 这个是属于注释包，每个基因集可能对应的注释包不一样，要从基因集所在的平台找到对应的注释包，然后从bioconductor获取安装。
library(dplyr)
library(stringr)
library(msigdbr) # 提供多个物种的MSigDB数据

dif_genes <- read.csv("dif_genes.csv")
# 得到差异分析结果：all_different_genes
geneList <- dif_genes$avg_logFC
# 如果是Ensemble ID，并且如果还带着版本号，需要去除版本号，再进行基因ID转换，得到Entrez ID
names(geneList) <- dif_genes$gene_id
# 最后从大到小排序，得到一个字符串
geneList <- sort(geneList, decreasing = T) 

# 检查是否有重复基因名
genelist <- dif_genes$gene_id
genelist[duplicated(genelist)]

数据集准备

MSigDB数据库是GSEA的官方数据库

官网

misgdbr包

• Introduction to msigdbr (r package vignette)

msigdbr包可提供多个物种的MSigDB数据

# GSEA #mdigbd dataset
msigdbr_df <- msigdbr(species = "Mus musculus",category = "C2")
# msigdbr_df$gene_symbol:"SYMBOL";$entrez_gene:"ENTREZID"
# msigdbr_list <- split(x = msigdbr_df$entrez_gene, f = msigdbr_df$gs_name)
# msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, entrez_gene) %>% as.data.frame()
# msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, ensembl_gene) %>% as.data.frame()
msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()

用clusterProfiler实现GSEA

• 利用clusterProfiler进行富集分析

gsea_result <- GSEA(geneList = geneList,
                    minGSSize = 1,
                    maxGSSize = 1000,
                    pvalueCutoff = 1,
                    TERM2GENE = msigdbr_t2g)
# 将其中的基因名变成symbol ID
#gsea_result <- setReadable(gsea_result, OrgDb = org.Mm.eg.db)
# 还可以直接点击查看，只需要转成数据框
#gsea_result_df <- as.data.frame(gsea_result)
# 导出结果
write.csv(gsea_result@result, file = "gsea_result.csv")
# 对感兴趣的通路可视化
gseaplot2(gsea_result, geneSetID = c("PATHWAY_NAME"))

或进行GO分析

## gseGO
# 想看biological process(BP)，可以先不过滤p值
go_result <- gseGO(geneList = geneList,
                   ont = "BP", 
                   OrgDb = org.Mm.eg.db,
                   minGSSize = 10,
                   maxGSSize = 500,
                   pvalueCutoff = 1)
# 将其中的基因名变成symbol ID
go_result <- setReadable(go_result, OrgDb = org.Mm.eg.db)
# 还可以直接点击查看，只需要转成数据框
go_result_df <- as.data.frame(go_result)
# 导出结果
write.csv(go_result@result, file = "gseGO_result.csv")
# 对感兴趣的通路可视化
gseaplot2(go_result, geneSetID = c("GO:0032981"))

KEGG分析

##gseKEGG
kk <- gseKEGG(geneList = geneList,
              organism = 'hsa',
              nPerm = 1000,
              minGSSize = 10,
              maxGSSize = 500,
              pvalueCutoff = 1,
              verbose = FALSE)
gseaplot(kk, geneSetID = "hsa04110")

结果注释，可视化解读

enrichmentscore 或 NES 大于0表示上调，小于0表示下调
GSEA详细解释及结果解读

===================================================

完整代码

rm(list=ls())
setwd("C:\\Users\\DELL\\Desktop\\gsea")
memory.limit(102400)

library(clusterProfiler)
library(enrichplot)
library(org.Mm.eg.db) # mouse全基因组注释包，用于不同版本基因名的转换
library(dplyr)
library(stringr)
library(msigdbr) 
 
# GSEA #mdigbd dataset
msigdbr_df <- msigdbr(species = "Mus musculus", category = "C2")
# msigdbr_df$gene_symbol:"SYMBOL";$entrez_gene:"ENTREZID"
# msigdbr_list <- split(x = msigdbr_df$entrez_gene, f = msigdbr_df$gs_name)
# msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, entrez_gene) %>% as.data.frame()
# msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, ensembl_gene) %>% as.data.frame()
msigdbr_t2g <- msigdbr_df %>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()

files <- list.files(getwd(), pattern = "*.csv") 

# 多组数据，用for循环
for (i in 1:length(files)){
  file <- files[i]
  print(file)
  dif_genes <- read.csv(file)
  # 得到差异分析结果：all_different_genes
  geneList <- dif_genes$avg_logFC
  # 如果是Ensemble ID，并且如果还带着版本号，需要去除版本号，再进行基因ID转换，得到Entrez ID
  names(geneList) <- dif_genes$gene_id
  # 最后从大到小排序，得到一个字符串
  geneList <- sort(geneList, decreasing = T) 

  # 检查是否有重复基因名
  genelist <- dif_genes$gene_id
  genelist <- genelist[duplicated(genelist)]
  
  gsea_result <- GSEA(geneList = geneList,
                      minGSSize = 1,
                      maxGSSize = 1000,
                      pvalueCutoff = 1,
                      TERM2GENE = msigdbr_t2g)
  # 将其中的基因名变成symbol ID
  #gsea_result <- setReadable(gsea_result, OrgDb = org.Mm.eg.db)
  # 还可以直接点击查看，只需要转成数据框
  #gsea_result_df <- as.data.frame(gsea_result)
  # 导出结果
  write.csv(gsea_result@result, file = "gsea_result.csv")
}