一、引言
对于测序基因组进行 KEGG（Kyoto Encyclopedia of Genes and Genomes）和COG（clusters of orthologous groups，对直系同源基因进行聚类）功能注释，基本成为基因组注释的标配内容， 特别是微生物基因组基因注释，其基因功能注释逻辑基础是 直系同源基因具有相同的功能，最经典的鉴定直系同源基因策略是 BBH（bi-directional best hit）策略，但是通常最直接的直系同源基因很难鉴定，而对同源基因进行聚类并定义一个簇会是更好的策略：每一个簇会包含直系同源基因（伴随物种形成事件出现）和旁系同源基因（伴随拷贝事件出现），每一簇共享同一个功能， KO(KEGG Orholog), COG, eggNOG 等都是基于聚类的方式定义簇，并对簇进行注释。
今天要讲的是eggNOG, eggNOG的出现要从COG说开，下面看看NCBI COG的数据库主要更新历史：
- 从 1997 年 第一个公布版本,7个完整基因组，720个COG分类， 包含原核基因组和单细胞真核基因组（酵母），2003 年和2014 年进行了版本升级，最后只保留了细菌和古菌，包含了711个基因组以及4,631个COG分类， 26个功能分类。
- 2013 年构建真核分支COG（KOG, Eukaryotic orthologous groups）；
- 2007 年构建古菌分支COG（arCOG， Archaeal Clusters of Orthologous Genes），2012 年和2014 年arCOG进一步升级，arCOG比较适合用于古菌基因组注释；
- 2011 年构建Phage分支COG（POG，phage orthologous groups），2013 年进行了升级；
由于计算资源需求，NCBI COG 构建了不同系统分类分支的COG簇，比如arCOG，KOG, POG等，推荐使用这些分支对新测序基因组进行注释，其实eggNOG 尤其是4.x版本也使用了clade特异的聚类模式。
从2003年至2014年NCBI COG一直未更新，EMBL EggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)继承了NCBI COG的衣钵，极大的扩展了基因组信息。当前版本为 4.5.1 版本， 把包含了2,031个基因组， 其中 352病毒基因组， 190k个直系同源家族。
eggNOG 数据库包含了丰富的注释信息，除了COG/KOG/NOG的分类和注释信息外，还包含了KEGG/GO/SMART/PFAM信息。
最新版本的EggNOG 还提供了自动化注释工具eggnog-mapper，可很方便的完成基因组的功能注释，注释信息可以关联COG/KOG/KEGG/GO/BiGG等。
eggNOG在不同的系统分类水平都进行了构建直系同源簇，因此当我们可以确定基因组系统分类时，可以选择系统进化上最近的分支进行注释，以减少计算量， 比如注释的真菌组可以选择fuNOG, 古菌基因组可选择arNOG, 每一个分类水平都提供了HMM谱， 可使用HMMER快速比对。
最新版本eggNOG构建流程如下图：
eggNOG 大致通过以下几个步骤完成：
- 选择基因组
- all-against-all SW序列比对，形成SIMAP库
- 物种分类水平选择
- 对已选择分类水平构建直系同源簇，图无监督聚类算法
- 根据物种系统分类，建立直系同源簇的等级关系， 根据算法对OG进行合并和分拆使得不同分类水平的OG保持一致。
- 对直系同源簇进行功能注释， 包括了KEG/GO/COG/SMART等
- 对直系同源簇构建系统进化关系，确立直系同源和旁系同源关系， 并构建HMM谱
了解eggNOG数据库的构建思路有助于我们理解eggnog-mapper是如何通过实现的。
下面我们着重介绍eggnog-mapper 工具， 看看如何使用eggNOG资源实现序列（核酸和氨基酸都可以）的功能注释。
二、eggNOG-mapper介绍
2.1 eggNOG-mapper算法：
图一、eggNOG-mapper工作流
通过A/B/C/D 四个步骤完成eggNOG注释和功能关联。
下面组条解读：
1. 数据库搜索，可以通过HMMER和DIAMOND获取种子序列，DIAMOND直接与序列库进行相似性搜索，通过 HMMER 与特定 level 的HMMER谱搜索获取最佳HMM,然后使用 phmmer 获取最佳种子序列。
2. 解析直系同源基因 Orthologs, 存在几种模式:
a. one-to-one, 没有发生duplication 事件;
b. one-to-many, 物种形成后发生的duplication 事件;
c. many-to-many;
3. 去除远源Orthologs，选择系统分类上比较近的Orthologs
4. 对Orholog类型进行处理，功能转移（注释）
2.2 eggNOG-mapper数据准备：
当前版本： 1.0.3
下载库文件：
emapperdb-4.5.1
hmmdb_levels/ 29-May-2016 08:05 -
OG_fasta.tar.gz 13-Oct-2017 09:18 6G
eggnog.db.gz 17-Oct-2017 12:00 4G
eggnog_proteins.dmnd.gz 13-Oct-2017 09:18 2G
og2level.tsv.gz
注： 如果使用的是最新的Diamond版本，需要重新对序列库文件进行格式化，按照下面方式进行操作：
gunzip eggnog_proteins.dmnd
diamond getseq --db eggnog_proteins.dmnd >eggnog_proteins.pep
diamond makedb --db eggnog_proteins.dmnd --in eggnog_proteins.pep
更多使用帮助请参考： Basic-Usage
2.3 安装和使用：
安装比较简单，直接下载解压后即可使用, 可参考 Installation 和 Basic Usage：
git clone https://github.com/jhcepas/eggnog-mapper
export PATH=$PWD/eggnog-mapper:$PATH
emapper.py -i test/polb.fa --output polb_bact -d bact
emapper.py -i test/p53.fa --output p53_maNOG -m diamond
emapper 的主要命令行接口：
$ emapper.py -h
usage: emapper.py [-h] [--guessdb] [--database] [--dbtype {hmmdb,seqdb}]
[--data_dir] [--qtype {hmm,seq}] [--tax_scope]
[--target_orthologs {one2one,many2one,one2many,many2many,all}]
[--excluded_taxa]
[--go_evidence {experimental,non-electronic}]
[--hmm_maxhits] [--hmm_evalue] [--hmm_score]
[--hmm_maxseqlen] [--hmm_qcov] [--Z] [--dmnd_db DMND_DB]
[--matrix {BLOSUM62,BLOSUM90,BLOSUM80,BLOSUM50,BLOSUM45,PAM250,PAM70,PAM30}]
[--gapopen GAPOPEN] [--gapextend GAPEXTEND]
[--seed_ortholog_evalue] [--seed_ortholog_score] [--output]
[--resume] [--override] [--no_refine] [--no_annot]
[--no_search] [--report_orthologs] [--scratch_dir]
[--output_dir] [--temp_dir] [--no_file_comments]
[--keep_mapping_files] [-m {hmmer,diamond}] [-i]
[--translate] [--servermode] [--usemem] [--cpu]
[--annotate_hits_table] [--version]
optional arguments:
-h, --help show this help message and exit
--version
Target HMM Database Options:
--guessdb guess eggnog db based on the provided taxid
--database , -d specify the target database for sequence searches.
Choose among: euk,bact,arch, host:port, or a local
hmmpressed database
--dbtype {hmmdb,seqdb}
--data_dir Directory to use for DATA_PATH.
--qtype {hmm,seq}
Annotation Options:
--tax_scope Fix the taxonomic scope used for annotation, so only
orthologs from a particular clade are used for
functional transfer. By default, this is automatically
adjusted for every query sequence.
--target_orthologs {one2one,many2one,one2many,many2many,all}
defines what type of orthologs should be used for
functional transfer
--excluded_taxa (for debugging and benchmark purposes)
--go_evidence {experimental,non-electronic}
Defines what type of GO terms should be used for
annotation:experimental = Use only terms inferred from
experimental evidencenon-electronic = Use only non-
electronically curated terms
HMM search_options:
--hmm_maxhits Max number of hits to report. Default=1
--hmm_evalue E-value threshold. Default=0.001
--hmm_score Bit score threshold. Default=20
--hmm_maxseqlen Ignore query sequences larger than `maxseqlen`.
Default=5000
--hmm_qcov min query coverage (from 0 to 1). Default=(disabled)
--Z Fixed database size used in phmmer/hmmscan (allows
comparing e-values among databases).
Default=40,000,000
diamond search_options:
--dmnd_db DMND_DB Path to DIAMOND-compatible database
--matrix {BLOSUM62,BLOSUM90,BLOSUM80,BLOSUM50,BLOSUM45,PAM250,PAM70,PAM30}
Scoring matrix
--gapopen GAPOPEN Gap open penalty
--gapextend GAPEXTEND
Gap extend penalty
Seed ortholog search option:
--seed_ortholog_evalue
Min E-value expected when searching for seed eggNOG
ortholog. Applies to phmmer/diamond searches. Queries
not having a significant seed orthologs will not be
annotated. Default=0.001
--seed_ortholog_score
Min bit score expected when searching for seed eggNOG
ortholog. Applies to phmmer/diamond searches. Queries
not having a significant seed orthologs will not be
annotated. Default=60
Output options:
--output , -o base name for output files
--resume Resumes a previous execution skipping reported hits in
the output file.
--override Overwrites output files if they exist.
--no_refine Skip hit refinement, reporting only HMM hits.
--no_annot Skip functional annotation, reporting only hits
--no_search Skip HMM search mapping. Use existing hits file
--report_orthologs The list of orthologs used for functional transferred
are dumped into a separate file
--scratch_dir Write output files in a temporary scratch dir, move
them to final the final output dir when finished.
Speed up large computations using network file
systems.
--output_dir Where output files should be written
--temp_dir Where temporary files are created. Better if this is a
local disk.
--no_file_comments No header lines nor stats are included in the output
files
--keep_mapping_files Do not delete temporary mapping files used for
annotation (i.e. HMMER and DIAMOND search outputs)
Execution options:
-m {hmmer,diamond} Default:hmmer
-i Input FASTA file containing query sequences
--translate Assume sequences are genes instead of proteins
--servermode Loads target database in memory and keeps running in
server mode, so another instance of eggnog-mapper can
connect to this sever. Auto turns on the --usemem flag
--usemem If a local hmmpressed database is provided as target
using --db, this flag will allocate the whole database
in memory using hmmpgmd. Database will be unloaded
after execution.
--cpu
--annotate_hits_table
Annotatate TSV formatted table of query->hits. 4
fields required: query, hit, evalue, score. Implies
--no_search and --no_refine.
主要可选参数解读：
-i: 输入文件
-o: 输出文件前缀
-m: 使用HMMER策略还是DIAMOND策略，默认使用HMMER，如果eggNOG存在很近的序列，HMMER\
鉴定的种子子直系同源序列的准确性会较低，选用DIAMOND比较普适，尤其是metagenome \
序列
--cpu：使用的线程数
--translate：如使用的核酸序列，选择HMMER策略时需要先翻译成氨基酸序列，DIAMOND模式 \
会自动选择blastx序列比对模式，如果输入的是微生物基因组核酸序列,推荐使\
用Prodigal预测基因，如果是拼装的转录本，可以考虑使用TRANSDECODE预测 \
ORF
--usemem： 将 emapper.db 读入内存，他然后执行emapper， 默认为磁盘数据库查找
--servermode：使用服务器模式，即将emapper.db调入内存，供其它进程使用，分配监听端口。
--output_dir：输出文件目录
--report_orthologs： 列出所有进行功能转移的直系同源基因
--no_refine: 只针对HMMER模式，需要结合 --no_annot 使用
--no_annot: 只汇总鉴定的最佳seed序列以及相应的E值和Bitscore值
--no_search:可直接基于--no_annot的结果进行后续功能注释
--target_orthologs: one2one,many2one,one2many,many2many,all可选。
--data_dir: 数据库目录
--tax_scope: 指定选择的直系同源基因的物种分类范围，默认为自动判断。
从实现上看，emapper使用SQLite数据库存储注释信息，序列比对（HMMER和Diamond）和序列注释是独立的两步，可以分开执行，注释部分可直接使用序列比对的中间结果。
了解这些参数，执行起来就比较简单，比如：
HMMER模式
emapper.py -i test/polb.fa --output polb_bact -d bact
diamond模式
emapper.py -i test/p53.fa --output p53_maNOG -m diamond
Server模式：
python emapper.py -d bact --cpu 10 --servermode
python emapper.py -d bact:localhost:51500 -i test/polb.fa -o polb
2.4 结果解读：
针对生成文件： polb.emapper.annotations：
结果解读：
第一列：查询序列名称；
第二列：eggNOG种子序列；
第三列：eggNOG种子序列 evalue；
第四列：eggNOG种子序列 bit score；
第五列：预测基因名称；
第六列：GO_terms, 预测的GO，分号分隔；
第七列：KEGG_KO: 预测的KO,分号分隔；
第八列：BiGG_Reactions: BiGG代谢反应预测，分号分隔；
第九列：eggNOG Taxonomic Scope信息；
第十列：匹配的OGs;
第十一列：best_OG|evalue|score: Best matching Orthologous Groups (only in HMM mode)
第十二列：COG功能分类；
第十三列：eggNOG功能描述；
注：种子可通过两种方式获取，一种基于HMMER HMM谱 相似性搜索，一种基于DIAMOND序列相似性搜索。
针对Metagenome, 首选使用DIAMOND于整个eggNOG库进行相似性搜索。
细节可参考：Results Interpretation
三、功能关联
程序见: https://github.com/jameslz/ann-utils
GO/KO/BiGG注释
我们根据上述表中的第六列/第七列/第八列 实现了对 KEGG/BiGG/GO 通用的注释提取程序，用于生成类似 goslim 输出结果。
$ emapper-ann
Usage: emapper-ann <text>
Options:
-d STR parse annotation, [KEGG, BiGG, GO], default: [KEGG]
-v print version number
示例：
emapper-ann -d KEGG nr.emapper.annotations
KOG/COG注释
有时候我们只关注COG或者KOG分类，这样我们就可以根据 第九列选取对应的KOG/COG匹配信息。
$ emapper-eggnog
Usage: emapper-eggnog [options] <emapper>
Options:
-c catalog to report [COG, KOG, XOG], XOG for [KC]OG; default: [COG]
-v print version number
示例：
emapper-eggnog -c COG nr.emapper.annotations
到目前为止， 我们可以根据eggnog-mapper完成KEGG/BiGG/COG/KOG/GO的功能注释。
本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料， 版权归 上海逻捷信息科技有限公司 所有。
Last Update: 2017-12-01 11:09 AM
