最开始的思路是先构建OrgDb，然后使用enrichGO和enrichKEGG函数做分析。但是最后遇到报错 Error in get_GO_data(OrgDb, ont, keyType) : keytype is not supported...。因为这两个函数都需要指定一个keyType参数。如果物种是人可以根据基因名的类型来指定，比如symbol或者entrezid等。但是自己构建的OrgDB如何指定这个参数还不太清楚。后来发现不构建orgdb也可以做GO或者KEGG的富集分析，可以使用enricher()函数。下面记录利用eggnog-mapper软件注释结果做GO和KEGG富集分析做GO和KEGG富集分析的过程。

这里我使用 Schizosaccharomyces pombe 这个物种的蛋白数据做例子，搜了一下拉丁名好像是裂殖酵母。

第一步是使用eggnog-mapper做功能注释

conda activate emapper
python emapper.py -i orgdb_example/GCF_000002945.1_ASM294v2_protein.faa --output orgdb_example/out -m diamond --cpu 8

将注释结果下载到本地，手动删除前三行带井号的行，第四行开头的井号去掉，文件末尾带井号的行去掉。

利用R语言将注释结果整理成enricher函数需要的输入格式

GO富集

library(stringr)
library(dplyr)
eggegg[egg==""]gterms %
  select(query_name, GOs) %>%
  na.omit()
gene2go                       gene = character())
for (row in 1:nrow(gterms)) {
  gene_terms   gene_id   tmp                     term = gene_terms)
  gene2go }
head(gene2go)

> head(gene2go)
# A tibble: 6 x 2
  gene           term

1 NP_001018179.1 GO:0003674
2 NP_001018179.1 GO:0003824
3 NP_001018179.1 GO:0004418
4 NP_001018179.1 GO:0005575
5 NP_001018179.1 GO:0005622
6 NP_001018179.1 GO:0005623

获得一个两列的数据框，有了这个数据框就可以做GO富集分析了

在 https://www.jianshu.com/p/9c9e97167377 这篇文章里的评论区有人提到上面用到的for循环代码效率比较低，他提供的代码是

gene_ids eggnog_lines_with_go eggnog_lines_with_go
eggnog_annoations_go gene_to_go                                    times = sapply(eggnog_annoations_go, length)),
                         term = unlist(eggnog_annoations_go))
head(gene_to_go)

> head(gene_to_go)
            gene       term
1 NP_001018179.1 GO:0003674
2 NP_001018179.1 GO:0003824
3 NP_001018179.1 GO:0004418
4 NP_001018179.1 GO:0005575
5 NP_001018179.1 GO:0005622
6 NP_001018179.1 GO:0005623

用这个代码替换for循环，确实快好多。

接下来可以做GO富集分析了

首先准备一个基因列表，我这里选取gene2go中的前40个基因作为测试 还需要为TERM2GENE=参数准备一个数据框，第一列是term，第二列是基因ID，只需要把gene2go的列调换顺序就可以了。

library(clusterProfiler)
gene_listterm2genedf             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             TERM2GENE = term2gene)
head(df)
barplot(df)
dotplot(df)

y轴的标签通常是GO term (就是文字的那个)而不是GO id。clusterProfiler包同样提供了函数对ID和term互相转换。go2term()go2ont()

dfhead(df)
dim(df)
df1dim(df1)
head(df1)
df$termdf2dim(df2)
head(df2)
df$Onthead(df)
df3%
  select(c("term","Ont","pvalue"))
head(df3)
library(ggplot2)
ggplot(df3,aes(x=term,y=-log10(pvalue)))+
  geom_col(aes(fill=Ont))+
  coord_flip()+labs(x="")+
  theme_bw()

这里遇到一个问题：数据框如何分组排序？目前想到一个比较麻烦的办法是将每组数据弄成一个单独的数据框，排好序后再合并。

KEGG富集

library(stringr)
library(dplyr)
library(clusterProfiler)
eggegg[egg==""]gene2ko %
  dplyr::select(GID = query_name, Ko = KEGG_ko) %>%
  na.omit()
head(gene2ko)
head(gene2go)
gene2ko[,2]head(gene2ko)
#kegg_info.RData这个文件里有pathway2name这个对象
load(file = "kegg_info.RData")
pathway2gene % left_join(ko2pathway, by = "Ko") %>%
  dplyr::select(pathway=Pathway,gene=GID) %>%
  na.omit()
head(pathway2gene)
pathway2name
df             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             TERM2GENE = pathway2gene,
             TERM2NAME = pathway2name)
dotplot(df)
barplot(df)

以上最开始的输入文件是eggnog-mapper软件本地版注释结果，如果用在线版获得的注释结果，下载的结果好像没有表头，需要自己对应好要选择的列。

如果大家需要以上提到的任何文件，都可以给我的公众号留言。

欢迎大家关注我的公众号

小明的数据分析笔记本

中国加油！武汉加油！