使用SigDB(Molecular Signatures Database)基因集进行富集分析,包含8个系列
相较于KEGG,SigDB数据集包含的功能更多
对 MigDB中的全部基因集 做GSEA分析。
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
# 安装需要的包
# BiocManager::install("GSEABase")
# 导入kegg已经注释好的数据
load(file = 'anno_DEG.Rdata')
geneList=DEG$logFC
names(geneList)=DEG$symbol
geneList=sort(geneList,decreasing = T)
# 导入MigDB全部基因集
# 需要在 MigDB 官网下载 为gmt文件
# 下载好的路径
d='MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
# 导入8个序列的基因集名
#
library(GSEABase)
## 使用lapply循环读取每个gmt文件并进行GSEA分析
# 考验你计算机能力的时候到了
f='gsea_results.Rdata'
if(!file.exists(f)){
#相比较apply,lapply较多的用于list的循环操作
gsea_results <- lapply(gmts, function(gmtfile){
#读取gmt文件
geneset <- read.gmt(file.path(d,gmtfile))
# 打印目前进行的基因库
print(paste0('Now process the ',gmtfile))
# 进行GSEA分析
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
head(egmt)
# 返回计算的结果
return(egmt)
})
# 保存结果到本地文件
save(gsea_results,file = f)
}
# 如果上述操作完成之后
# 再次运算就可以直接导入
load(file = f)
#提取gsea结果,
gsea_results_list<- lapply(gsea_results, function(x){
cat(paste(dim(x@result)),'\n')
x@result
})
## 3 11
## 996 11
## 186 11
## 233 11
## 671 11
## 95 11
## 1591 11
## 27 11
# docall 函数能够对list使用dataframe结构的函数,下行为合并结果
gsea_results_df <- do.call(rbind, gsea_results_list)
# 选择有差异的基因集进行画图,第一个参数为基因集,第二个参数需要对应前面参数
# 具体可以查看gsea_results里面的行名
gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE')
gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6')
gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE')