Archives

序列功能注释神器:eggNOG-mapper,KEGG/COG/KOG/GO/BiGG 一网打尽

一、引言

对于测序基因组进行 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

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算法:

image.png-462.8kB

图一、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 安装和使用:

安装比较简单,直接下载解压后即可使用, 可参考 InstallationBasic 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

Comments are closed.