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

经常收到粉丝提问，明明是按照我课程视频操作，也是按照我的代码在处理他自己的数据，但是做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，帮助他们多一点数据认知，让科研更上一个台阶：
推荐阅读

赞 (0)

相关推荐