不要怀疑，你的基因就是没办法富集到统计学显著的通路

Posted on by

经常收到粉丝的提问，明明是按照我课程视频操作，也是按照我的代码在处理他自己的数据，但是做kegg数据库富集的时候，就是返回值为空。

代码如下:

#该包的gene需要识别ENTREZID，因此需要再次转换
library(org.Hs.eg.db)
gene.df<-bitr(gene,fromType="SYMBOL",
 toType=c("ENSEMBL","ENTREZID"),
 OrgDb=org.Hs.eg.db)
head(gene.df)
genelist=gene.df$ENTREZID
kegg <- enrichKEGG(genelist, organism = 'hsa', pvalueCutoff = 0.01)
head(kegg)
dim(kegg)
dotplot(kegg,showCategory=30)

有意思的是，我回复了其中一个粉丝，需要调整pvalueCutoff参数。

他给我的回答是：

请问这个pvalueCutoff应该怎么设置呢？我在之前筛选表达上调的基因的时候设定的cutoff值就是p.value<0.01，logFC>logFC_cutoff的基因呀。

但enrich的时候也将pvalueCutoff同样设置成0.01，但还是不行，试了0.001，也还是空白呢。

我想，这位同学应该是不明白P值是什么，居然把P值阈值越调越小了。其实应该是：

kk.down <- enrichKEGG(gene = gene_down,
 organism = 'hsa', 
 pvalueCutoff = 0.9,
 qvalueCutoff =0.9)

很明显是这位同学没有去看函数的帮助文档。

而且不能理解KEGG富集的原理就是超几何分布检验了，也就没办法接受为什么自己给定的基因集，在KEGG数据库里面，居然会无法富集到统计学显著的通路。

举个例子

中国有13亿人口，其中广东人是1个亿，现在发现某次比赛金牌获得者100名里面有13位广东人。你就需要计算这次比赛金牌获得者人群里面是否富集了广东人，超几何分布检验公式套进去，发现统计学P值很大，并不显著。

然后你不甘心，去继续检查是否富集了其它省的人，发现都不显著。

我勒个去，难道你分析错了吗？仅仅是因为这次比赛确实没有人口分布的偏向性啊！同样的道理，你的基因，也是有可能并没有kegg通路的偏向性的啊！

不过，你换一个方式，比如统计比赛金牌获得者，发现不同年龄的确有富集哦，比如100个都是青少年，但是青少年占中国人口肯定不是百分百啊！所以这个统计学显著的富集结论你就可以心安理得的下了，同理，kegg数据库富集不到，你可以换其它数据库哦，msigdb里面有很多哦。

pathway在我这里是基因集的别名，其中msigdb有着丰富的基因集，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

  • H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
  • C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
  • C2: curated gene sets：（专家）校验基因集合，基于通路、文献等；
  • C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分；
  • C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
  • C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）；
  • C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据；
  • C7: immunologic signatures: 免疫相关基因集合。

比如注释到GO数据库：

go <- enrichGO(gene_up, OrgDb = "org.Hs.eg.db", ont="all") 
library(ggplot2)
library(stringr)
barplot(go, split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free") 
barplot(go, split="ONTOLOGY",font.size =10)+ 
 facet_grid(ONTOLOGY~., scale="free") + 
 scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
 ggsave('gene_up_GO_all_barplot.png') 
go <- enrichGO(gene_down, OrgDb = "org.Hs.eg.db", ont="all") 
barplot(go, split="ONTOLOGY",font.size =10)+ 
 facet_grid(ONTOLOGY~., scale="free") + 
 scale_x_discrete(labels=function(x) str_wrap(x, width=50))+
 ggsave('gene_down_GO_all_barplot.png')
​

Comments are closed.