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

经常收到粉丝提问，明明是按照我课程视频操作，也是按照我的代码在处理他自己的数据，但是做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)

很明显是这位同学没有去看函数的帮助文档，其实我推荐过 为R包写一本书（像Y叔致敬），然后这个同学的统计学也烂的一塌糊涂！

而且不能理解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')

文末友情宣传

强烈建议你推荐我们生信技能树给身边的博士后以及年轻生物学PI，帮助他们多一点数据认知，让科研更上一个台阶：

• 生信爆款入门-全球听（买一得五）（第5期）（可能是最后一期）你的生物信息学入门课
• (必看！)数据挖掘第3期（两天变三周，实力加量），医学生/临床医师首选技能提高课
• 生信技能树的2019年终总结  ，你的生物信息学成长宝藏
• 2020学习主旋律，B站74小时免费教学视频为你领路，还等什么，看啊！！！

推荐阅读