Biostack.org

一、引言

对于测序基因组进行 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的数据库主要更新历史：

1. 从 1997 年 第一个公布版本,7个完整基因组，720个COG分类， 包含原核基因组和单细胞真核基因组（酵母），2003 年和2014 年进行了版本升级，最后只保留了细菌和古菌，包含了711个基因组以及4,631个COG分类， 26个功能分类。
2. 2013 年构建真核分支COG（KOG, Eukaryotic orthologous groups）；
3. 2007 年构建古菌分支COG（arCOG， Archaeal Clusters of Orthologous Genes），2012 年和2014 年arCOG进一步升级，arCOG比较适合用于古菌基因组注释；
4. 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 大致通过以下几个步骤完成：

1. 选择基因组
2. all-against-all SW序列比对，形成SIMAP库
3. 物种分类水平选择
4. 对已选择分类水平构建直系同源簇，图无监督聚类算法
5. 根据物种系统分类，建立直系同源簇的等级关系， 根据算法对OG进行合并和分拆使得不同分类水平的OG保持一致。
6. 对直系同源簇进行功能注释， 包括了KEG/GO/COG/SMART等
7. 对直系同源簇构建系统进化关系，确立直系同源和旁系同源关系， 并构建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