R语言——基于clusterProfiler的GSEV数据分析

背景

在进行差异分析过程中，尤其是食品中一些功能因子的添加可能影响并不在几个基因上，而是对于一系列基因有着一定的影响，这些影响可能集中在某一些通路上面。因此，联合分析某一些通路上，或者某些特定组合上的基因变化，对于寻找差异不是特别明显，但是却有实际作用的基因特别重要。

工作环境加载

library(clusterProfiler)
library(org.Hs.eg.db)
library(dplyr)

数据前处理

input <- read.table("GSEA_analysis_input.csv",sep=",",header=T,check.names=F)
colnames(input) <- c("ENSEMBL","logFC")
###################################
#############对编码进行转化#########
input_gene <- bitr(input$gene_id,
                fromType = "ENSEMBL",
                toType = c("ENSEMBL","ENTREZID"),
                OrgDb = org.Hs.eg.db)# 对编码进行转化。

head(gene.df)

转化之后，一个严重的问题，在于很多基因名称在转化过程中并没有对应，我们现在不知道原因，
但是为了下一步的数据处理，我们只可以将没有转化成功的予以丢弃，同时把转化成功的数据和对应的logFC值组成下一步处理需要的文献。

input_all = left_join(input_gene,input,by="ENSEMBL")

转化为需要的数据格式，这里要注意，

## 1: 提取logFC值，并储存在一个向量中
geneList = input_all$logFC # 按照logFC值对基因进行排序
## 2: 对geneList进行命名
names(geneList) = as.character(input_all[,2])
head(geneList)
## 3: 根据logFC值降序排列
geneList = sort(geneList, decreasing = TRUE)

进行基因富集分析

########## GO的GSEA富集分析：gseGO ##########
gsego <- gseGO(geneList     = geneList,#排序后的基因列表，一般根据logFC进行排序
              OrgDb        = org.Hs.eg.db,
              ont          = "CC",#可选择bp.MF,CC,ALL
              nPerm        = 1000,#置换检验的次数，默认1000，保持默认即可
              minGSSize    = 100,#最小基因集的基因数
              maxGSSize    = 500,#最大基因集的基因数
              pvalueCutoff = 0.05,#p值的阈值 这里对于Met数据要做调整，否者没有结果。
              verbose      = FALSE)#是否输出提示信息，默认为false
head(gsego)

# 保存GO的GSEA分析结果到文件
write.table(gsego,file="data/GSEA_GO_result.txt",sep="\t",
            quote=F,row.names = F)

KEGG的GSEA富集分析

########## KEGG的GSEA富集分析：gseKEGG ##########
gsekk <- gseKEGG(geneList     = geneList,
               organism     = 'hsa',
               nPerm        = 1000,
               minGSSize    = 120,
               pvalueCutoff = 0.05,
               verbose      = FALSE)
head(gsekk)

# 保存KEGG的GSEA分析结果到文件
write.table(gsekk,file="data/GSEA_KEGG_result.txt",sep="\t",
            quote=F,row.names = F)

读取制定的gmt文件

# 读取上面指定的gmt文件
all_msigdb <- read.gmt(file.path(msigdb_GMTs,msigdb))# 这里注意下载正确的文件，并且读入进来。

gsegmt <- GSEA(geneList, TERM2GENE=all_msigdb, verbose=FALSE)#TERM2GENE，数据集条目名称TERM与基因名的对应关系，一般是两列的数据框格式
head(gsegmt)

# 保存KEGG的GSEA分析结果到文件
write.table(gsegmt,file="data/GSEA_MSigDb_result.txt",sep="\t",
            quote=F,row.names = F)

# 将以上所有结果保存为R二进制格式，方便快速加载
save(geneList, gsego, gsekk, gsegmt, file = "data/gsea_result.Rda")