如果你的差异基因集(500个左右的基因)本身确实没有一些生物学功能的偏好性,非常有可能是在统计学显著阈值条件限定下,就是拿不到富集的kegg通路。这个时候,大家可以修改代码,比如:
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
# 需要自己差异分析筛选得到 gene_down 基因集 ,然后进行超几何分布检验
拿到通路后自己去调整阈值。
或者可以去MSigDB(Molecular Signatures Database)看看其它功能基因集是否可以富集到结果,MSigDB数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列(Collection),每个系列分别是:
有些时候,大家挑选到的差异基因集数量过少或者过多,又不知道如何调整阈值,得到的结果会稀奇古怪,不可理喻。也就是说,这个分析严重依赖于差异基因集的挑选,因人而异,这样不好。
在R里面的实现方式也是非常的简单:
# 需要自己对全部的基因进行一个排序,排序后的变量是geneList
library(clusterProfiler)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
tmp=kk_gse@result
gsea_kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
也可以对MSigDB数据库的全部基因集进行批量gsea分析,然后细致的去挑选结果。
# 选择gmt文件(MigDB中的全部基因集)
# 自己在 MigDB 官网下载 gmt 文件哦!
d='~/biosoft/MSigDB/symbols/'
gmts <- list.files(d,pattern = 'all')
gmts
#GSEA分析
library(GSEABase) # BiocManager::install('GSEABase')
## 下面使用lapply循环读取每个gmt文件,并且进行GSEA分析
## 如果存在之前分析后保存的结果文件,就不需要重复进行GSEA分析。
f='gsea_results.Rdata'
if(!file.exists(f)){
gsea_results <- lapply(gmts, function(gmtfile){
# gmtfile=gmts[2]
geneset <- read.gmt(file.path(d,gmtfile))
print(paste0('Now process the ',gmtfile))
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
head(egmt)
# gseaplot(egmt, geneSetID = rownames(egmt[1,]))
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
})
dev.off()
gsea_results_df <- do.call(rbind, gsea_results_list)
write.csv(gsea_results_df,file = 'gsea_results_df.csv')
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')
如果你差异基因没办法富集到通路,使用GSEA也许会有不一样的结果哦!
比如下面的这个通路,如果你使用上下调基因集去注释,通常呢你看不到它被富集到,因为top上下调基因都不是这个通路里面的,但是gsea分析就可以看到,它其实是显著下调的!