功能注释后如何做富集分析

本文转载自《生信媛》，已获得授权。

本文是为了回答知识星球里的一个提问，他为了用 clusterProfiler做富集分析，打算构建一个 OrgDb，也就是物种数据库。

提问

我之前写过用Bioconductor对基因组注释，用 Bioconductor/AnnotationHub对模式植物的基因进行注释如何用bioconductor进行注释。昨天的推送，我讲过新物种的注释基本上都是基于同源相似性搜索数据库完成，最后得到的就是基因名和数据库中注释的对应关系。 Orgdb是 Bioconductor计划中其中一环，通过构建一个物种各个数据库注释条目和基因的对应关系数据库，方便在得到基因后对基因进行注释。

enrichGO的前三个参数 gene, OrgDb, keyType的目的是利用数据库将基因编号转换成GO号。 enrichKEGG的前三个参数 gene, organism, keyType 的目的也是为了基于物种名和基因编号直接爬取KEGG，将基因编号转换成KO号。

如果你只是为了做GO和KEGG富集分析，有必要构建物种数据库吗？我的答案是没有必要，因为不构建物种数据库也能够用 clusterProfiler做富集分析。

我相信Y叔一定提供了不通过 OrgDb，将转换基因编号为GO/KO编号，然后做富集分析的方法，所以我就去翻了Y叔为 clusterProfiler写的文档。于是我找到这一篇use clusterProfiler as an universal enrichment analysis tool, 这里面提到了一个通用的函数 enricher用于支持新注释物种.

一个优秀的R包一定会提供优质的帮助文档，何况clusterProfiler的引用已经破500了。单篇引用超过500了，来了解一下

核心参数两个 gene, TERM2GENE,前者表示的基因编号，后者是GO/KEGG条目和基因编号的对应关系

1. enricher(gene, pvalueCutoff = 0.05, pAdjustMethod = "BH", universe,

2.  minGSSize = 10, maxGSSize = 500, qvalueCutoff = 0.2, TERM2GENE,

3.  TERM2NAME = NA)

由于我只拿到了我的KEGG注释，GO注释还在运行中，这次就以KEGG富集分析作为例子。

我从KEGG上拿到的注释是下面这种情况，很明显，有些基因没有注释。这些没有注释的基因应该如何注释？Y叔的建议是不要，全部丢掉，原因去参考资料中找。

1. CAROC969890.1

2. CAROC969900.1 K12736

3. CAROC969910.1 K02943

4. CAROC969920.1 K13356

5. CAROC969930.1

6. CAROC969940.1

7. CAROC969950.1

8. CAROC969960.1

9. CAROC969970.1

10. CAROC969980.1

简单的grep就可以完成这个剔除工作， grep K query.ko > kegg.tsv,然后将 kegg.tsv导入到我们的R语言中

1. gene_ko <- read.table("C:/Users/DELL/Desktop/KEGG.tsv", header = FALSE,

2.  sep = " ")

然后我们随机抽样几个基因作为 gene输入，同时构建 TERM2GENE的输入

1. term2gene <- data.frame(TERM=gene_ko$V2, GENE=gene_ko$V1)

2. gene_sample <- sample(gene_ko$V1, 100)

3. enkegg <- enricher(gene_sample, TERM2GENE = term2gene, pAdjustMethod = "none")

我这里不用多重实验矫正的原因，因为我是随机抽的基因，很有可能是一个富集都找不到。。所以为了后续演示，就把矫正去掉了，真实情况下，你是要的。

参考资料

• why clusterProfiler fails