目录
- 背景
- 工作环境加载
- 数据前处理
背景
在进行差异分析过程中,尤其是食品中一些功能因子的添加可能影响并不在几个基因上,而是对于一系列基因有着一定的影响,这些影响可能集中在某一些通路上面。因此,联合分析某一些通路上,或者某些特定组合上的基因变化,对于寻找差异不是特别明显,但是却有实际作用的基因特别重要。
工作环境加载
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")