【生信实战】最详细的非模式物种的GO/KEGG富集分析全攻略

有时候不得不感慨，时间已经来到了2024年的尾声，但是对于许多非模式的物种来说，各种工具、数据集、参考等等依然还在寒冬。今天带大家手动的构建属于小众物种的富集分析全攻略。

1. 简介

相信每一个做过生信的人来说，对于GO/KEGG等富集分析并不陌生。有越来越多的工具和云平台可以提供这些分析。然而，对于许多对于许多农口或者微生物研究人员来说，由于低质量的参考基因组，和无人构建的数据集都使得这一“简单”的需求莫名的走了许多弯路。
其实，简单来说，对于任何一个物种都需要对应的基因->条目Term这样的关系表，有了这样的关系就可以构建我们的抽奖模型“超分布几何检验”。只不过对于大多数的模式物种来说，有大量的前辈已经把这些整理好了；前人栽树，后人乘凉即可。所以，理论上当我们没有现成的数据库时，就需要我们构建自己的物种的注释数据库，就可以轻松的完成今后的富集分析。

2.注释库构建

注释库的构建我们使用eggNOG，eggNOG（evolutionary genealogy of genes: Non-supervised Orthologous Groups）是一种用于基因组功能注释的数据库和工具，广泛应用于基因组学和生物信息学研究。目前已经更新到了eggNOG v5版本，他可以一次性的完整GO/COG/KEGG/PFAM等一系列的注释。

2.1 准备工作

1. 全转录本/基因的核酸或蛋白序列： 通常为fasta格式，里面包含了我们物种的所有基因/转录本。当然如果是组装程度极低或者甚至无参考的，可能是一些contig/scaffold。
2.eggNOG-maaper：该工具可以在github上下载https://github.com/eggnogdb/eggnog-mapper/releases/
3.eggNOG数据库：这部份建议手动下载，目前最新版本的地址如下：

http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog_proteins.dmnd.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/mmseqs.tar.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog.db.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/eggnog.taxa.tar.gz
http://eggnog6.embl.de/download/emapperdb-5.0.2/pfam.tar.gz

2.2 比对数据库

安装好eggNOG-maaper，和对应的python依赖后将我们下载的数据库解压到eggNOG-mapper的data目录下，结构大致如下：

data/
├── eggnog.db
├── eggnog_proteins.dmnd
├── eggnog.taxa.db
├── eggnog.taxa.db.traverse.pkl
├── mmseqs
│   ├── mmseqs.db
│   ├── mmseqs.db.dbtype
│   ├── mmseqs.db_h
│   ├── mmseqs.db_h.dbtype
│   ├── mmseqs.db_h.index
│   ├── mmseqs.db.index
│   ├── mmseqs.db.lookup
│   └── mmseqs.db.source
└── pfam
    ├── Pfam-A.clans.tsv.gz
    ├── Pfam-A.hmm
    ├── Pfam-A.hmm.h3f
    ├── Pfam-A.hmm.h3i
    ├── Pfam-A.hmm.h3m
    ├── Pfam-A.hmm.h3m.ssi
    ├── Pfam-A.hmm.h3p
    └── Pfam-A.hmm.idmap

如果我们的输入是蛋白序列则使用nohup emapper.py -i proteins.fa -o myeggnog --cpu 12 --dbmem -m diamond &
如果我们的是核酸序列则需要先翻译，此时使用命令nohup emapper.py -i gene.fa --itype genome -o myeggnog --translate --cpu 12 --dbmem -m diamond &
具体的其他参数可以通过-h进行查看，cpu数可以根据实际调整。

2.3 整理结果

当emapper比对完毕后会得到 myeggnog.emapper.annotations文件，该文件就是我们后面用于构建我们自己GO/KEGG数据库的关键，我们需要对结果进行进一步的整理。

library(tidyverse)
egg <- read_tsv('myeggnog.emapper.annotations', comment = '##') %>% mutate( `#query` = str_replace( `#query`, "_//d+$", "")) %>% group_by(`#query`) %>% top_n( n=1, wt = score)
egg_df <- egg %>% select(GID = `#query`, GOs, KEGG_ko) %>% filter(GOs != '-', KEGG_ko != '-') %>% separate_rows(GOs, sep = ",")  %>% separate_rows(KEGG_ko, sep = ",") 
egg_df <- egg_df %>% mutate(KEGG_ko = str_replace_all(KEGG_ko,'ko:','')) %>% rename(KO = KEGG_ko) %>% distinct()
library(KEGGREST)
kegg_pathways <- keggList("pathway", "ko") # "ko" 表示 KO 专用的 pathway
kegg_ko_pathways <- keggLink("ko", "pathway")
kegg_db <- tibble(Pathway = names(kegg_pathways), Name = kegg_pathways)
kegg_db <- tibble(Ko = sub("ko:","", kegg_ko_pathways), Pathway = sub("path:", "", names(kegg_ko_pathways))) %>% filter(str_detect(Pathway, '^ko'))  %>% inner_join(kegg_db, by = 'Pathway')
gene_info <- egg %>% select(GID = `#query`, GENENAME = seed_ortholog) 
gene2go <- egg_df %>% select(GID, GOs) %>% mutate(EVIDENCE = 'IEA') %>% distinct() %>% rename(GO = GOs)
koterms <- egg_df %>% select(GID, KO) %>% distinct()
gene2pathway <- kegg_db %>% inner_join(koterms, by = c('Ko' = 'KO')) %>% select(GID, Pathway) %>% distinct()
library(AnnotationForge)
makeOrgPackage(gene_info=gene_info, go=gene2go, ko=koterms,  pathway=gene2pathway, version="1.0.0", maintainer='Terry Xizzy ', author='Terry Xizzy ',outputDir=".", tax_id="12345", genus="Genus", species="species",goTable="go")
kegg_df <- koterms %>% inner_join(kegg_db, by = c('KO' = 'Ko')) %>% write_tsv('MyKEGG.xls')

在实际运行中只需要修改自己的输入文件名即可。这里的结果是基于5.0.2构建，后续可能出现列名变化的情况，需要根据实际进行调整。
代码运行完毕后会在目录下生成用于GO富集的Org会根据genus和species生成，这里就是org.Gspecies.eg.db，使用install.packages("org.Gspecies.eg.db",repos = NULL, type = "source")安装我们的包，后续分析时候只需要参考正常的clusterProfiler方法，详细的方法可以参考我以前的文章，核心步骤如下：

library(org.Gspecies.eg.db))
Db <- org.Gspecies.eg.db
enrichGO(gene_list, OrgDb = Db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 1,qvalueCutoff = 1,keyType = 'SYMBOL')

这里有一点需要注意的是，我们构建数据使用的GID应该和我们输入的gene_list的使用同一种。
而KEGG的分析则需要我们手动使用 enricher进行，关键代码如下：

library(tidyverse)
library(clusterProfiler)
kegg_df <- read_tsv('MyKEGG.xls')
#这里用于测试我们取前100个GID作为gene_list进行分析
gene_list <- kegg_df$GID[1:100]
kegg_res <- enricher(gene_list, TERM2GENE = select(kegg_df, Pathway, GID),  TERM2NAME = select(kegg_df, Pathway, Name), pvalueCutoff = 1,qvalueCutoff = 1)
kegg_res %>% as_tibble

输出如下：

kegg_res %>% as_tibble
# A tibble: 134 × 9
   ID      Description GeneRatio BgRatio   pvalue p.adjust   qvalue geneID Count
                                   
 1 ko01120 Microbial … 70/81     162/46… 3.75e-97 5.02e-95 3.04e-95 unass…    70
 2 ko00010 Glycolysis… 46/81     46/4637 1.58e-88 1.06e-86 6.41e-87 unass…    46
 3 ko01110 Biosynthes… 75/81     441/46… 1.05e-71 4.69e-70 2.84e-70 unass…    75
 4 ko01200 Carbon met… 51/81     113/46… 9.27e-67 3.10e-65 1.88e-65 unass…    51

至此，我们就可以对任意的物种轻松的进行GO/KEGG分析，一起动手试试吧！