Archives

驯化VirSorter: 预测metagenome contigs的prophage

鉴定噬菌体的工具有Phage_Finder, Prophinder , PHASTPhiSpy,不过要讲的是 VirSorter。VirSorter适合不完整的基因组,单细胞基因组,宏基因组。

VirSorter 运行时间很长问题,主要问题是 HMMER的问题,HMMER支持多线程不理性,即使设置多线程,实际执行的时候基本都是单线程,导致运行时间比较长。

那解决这个问题的方式就是:将线程强制变成进程,根据hmmsearch的特点,将库文件拆分成指定ncpu份, 单独提交可以达到并行目的, 然后将拆分后的结果合并为输出文件即可。

为解决这个问题 Biostack,实现了 hmmsearch-virsorter 做为HMMSEARCH任务提交的中间件,替换掉VirSorter的提交方式,可以顺利进行真实的并行任务提交。

$ hmmsearch-virsorter


Program: hmmsearch-virsorter: HMM based annotation.
Version: 0.0.1
Contact: ZHANG LEI <zhanglei@logicinformatics.com>

Usage:   hmmsearch-virsorter [options] <sequence> <tblout> <output>

Options: -c INT    CPU number, default: [40]
         -d STR    database location, default: [/biostack/database/pfam/Pfam-A.hmm]

现在利用40线程,一个典型的细菌基因组基本3分钟就可以完成前噬菌体鉴定。

生物信息数据工具封装:序列比对之hmmscan-pipe

一、前言

为什么要封装这个hmmscan-pipe呢, hmmscanHMMER(当前最新版本 3.1b2)程序包的子程序,工作模式是: 蛋白序列对谱序列库(hmmpress构建索引), 常用的氨基酸序列功能注释工具, 常用功能信息数据库有: PfamSuperfamilydbCANSMARTTIGRFAM 等。

hmmscan 当前版本可使用 --cpu 指定使用的线程数,但是一般不能有效利用多核心资源,所以最佳实践为: 序列拆分。
使用 fastx-utils partition 分割成指定线程数文件,使用进程对每个文件单独提交,这样就比较适合在集群模式下工作。
hmmscan-pipe依赖hmmscan-utils 支持标准输入流,主要解决两个问题:

  1. hmmscan-utils domtblout, 格式化hmmscan输出格式,使用制表符分隔,采用默认的过滤模式:如果比对片段 >80aa Evalue 阈值使用1e-5, <80aa Evalue 阈值使用1e-4
  2. hmmscan-utils resolve, 解析HMM匹配区域的交叠问题,去除交叠区域比较大的比对, 实现了文章 A fast and automated solution for accurately resolving protein domain architectures BMC 算法

hmmscan-utils 程序提供了两子命令程序:

Usage: hmmscan-utils <command> <arguments>
supports:

domtblout       <domtblout>
resolve         <domtblout>

分别解决了上述两个问题。

二、hmmscan-pipe介绍

小封装 hmmscan-pipe 命令行接口:

$hmmscan-pipe

Program: hmmscan-pipe: HMM based annotation.
Version: 0.0.1
Contact: ZHANG LEI <zhanglei@logicinformatics.com>

Usage:   hmmscan-pipe [options] <sequence> <project>

Options: -c INT    CPU number, default: [56]
         -e double set evalue cutoff, default: [1e-10]
         -d STR    database location, default: [/biostack/database/pfam/Pfam-A.hmm]

最简单的任务提交模式:

$hmmscan-pipe   A189.faa  Pfam

程序执行序列比对任务:采用了上述默认参数, 序列库也是默认参数:”/biostack/database/pfam/Pfam-A.hmm”

稍微复杂点的可以为:

$hmmscan-pipe -c 30 -e 1e-5  -d /biostack/database/pfam/Pfam-A.hmm A189.faa  Pfam

封装的好处?

1. 使用'fastx-utils partition' 对序列进行分拆,多进程执行,提高CPU利用率,可有效支持集群模式;
2. 可增加很多辅助操作,将生成的文件直接转换成Excel文档;
3. 按照自己常用需求调整默认参数,减少命令行提交复杂度;

三、最后

Biostack 提供定制小封装服务, 如有需求请联系微信号: biostack

NCBI Taxonomy 数据库更新,提供lineage、host信息

我们分析metagenome数据离不开使用NCBI的Taxonomy数据,NCBI Taxonomy 提供了一棵物种树,其实每个节点(Node)都分配了一个数字标识符,可以唯一描述一个系统分类信息。

NCBI Taxonomy 数据库提供了一个 taxdump.tar.gz, 并记录了节点的描述信息(names.dmp)以及树的上下游信息(nodes.dmp), 刚刚发布的更新版本提供了额外的lineage信息(rankedlineage.dmp) 以及 host 信息。

另外NCBI已经不再给Strain水平分配这种数字标识符,所以NCBI Taxonomy 提供了 typematerial.dmp 文件用于关联种和菌株(strain)的映射关系。

利用新的数据库我们可以很容易对一些短序列分类器进行注释, 常用的操作如下:

1、 格式化数据库,一般可以使用 tsv-utils

cut -f1,5  fullnamelineage.dmp | sed 's/ $//' >fullnamelineage.db
cut -f1,5  taxidlineage.dmp | sed 's/ $//' >taxidlineage.db
cut -f1,3  host.dmp  >host.db

2、 典型使用场景

下面以Kraken为例子,介绍如何格式化为有效信息, kraken的结果:

C   E00552:27:HJ2JYALXX:4:1101:5233:1801    435590  203 816:40 435590:21 A:31 435590:13 0:53 435590:10 0:5
U   E00552:27:HJ2JYALXX:4:1101:5781:1801    0   252 0:116 A:31 0:75

四列组成:

第一列: 序列分类标识符, C为分类,U为未分类
第二列: 序列ID
第三列: 命中的NCBI Taxonomy 数字ID,未命中为0;
第四列: Kmer匹配的信息;比如: 816:40,40个Kmer匹配Taxonomy ID 816;

现在可以生成一个关于序列分类的列表:

E00552:27:HJ2JYALXX:4:1110:28615:27257  12509   vertebrates,human   \
Viruses; dsDNA viruses, no RNA stage; Herpesvirales; Herpesviridae; Gammaherpesvirinae; Lymphocryptovirus; Human gammaherpesvirus 4;

通过执行下面命令问题就解决了:

grep -P ^"C" GY.txt \
| cut -f2,3 \
| tsv-utils definition  -c 2 fullnamelineage.db  - \
| tsv-utils definition  -c 2  host.db   -  \
>GY.taxonomy.tsv

毒力因子注释Protocol:VFDB数据库

一、毒力因子

毒力因子(Virulence factor), 详细介绍参见 维基百科 Virulence_factor 页面, 细菌、病毒、真菌等生成的分子,并产生毒力(主要有侵袭力和毒素等),包括:

1. 在宿主定殖 (colonization),黏附在宿主消化道、呼吸道、生殖道、尿道及眼结膜等处,以免被肠蠕动、黏液分泌、呼吸道纤毛运动等作用所清除
2. 免疫逃避,逃避宿主的免疫应答
3. 免疫抑制,抑制宿主的免疫反应
4. 进入和退出细胞
5. 从宿主获得营养

毒力因子可编码在可移动遗传元件(比如质粒、基因岛、噬菌体等)上并进行水平基因转移(传播),使无害细菌变成危险的病原菌,所以在鉴定毒力因子时一般会考虑: 基因岛、分泌蛋白等。

二、病原菌毒力因子数据库 VFDB

vfdb

毒力因子数据库VFDB 由中国医学科学院研发,被广泛应用于毒力因子基因鉴定。
VFDB收集了包括30个属( 74个病原菌)的细菌毒力基因序列信息。

summary

VFDB提供了对应的毒力基因核酸和蛋白质序列信息,因此鉴定毒力基因最简单的办法就是序列比对(BLAST),

2.1 数据库预处理

数据预处理需要以下几个步骤:

VFDB的元信息可以通过序列文件以及提供的描述文件获得:

>VFG000676(gb|AAD32411) (lef) anthrax toxin lethal factor precursor [Anthrax toxin (VF0142)] [Bacillus anthracis str. Sterne]`

1、格式化序列文件,只保留毒力基因编号 VFG000676 并获得 VFG000676 -> VF0142 映射关系;
2、格式化序列库;
3、预处理VFDB的数据库描述信息:

VF_Name         -> Acinetobactin
VF_FullName     -> -
Bacteria        -> Acinetobacter baumannii
Characteristics -> -
Structure       -> An iron-chelating molecule composed of equimolar quantities \
                   of 2,3-dihydroxybenzoic acid (DHBA), L-threonine, and N-hydroxyhistamine
Function        -> High-affinity catechol-hydroxamate siderophore competing with host cells for iron
Mechanism       -> -
Keyword         -> Iron uptake; Siderophore
VFID            -> VF0467

这里我们比较感兴趣应该是: VFID列和Keyword列, 需要格式化成:
VFID Keyword 映射关系;
4、序列比对
5、keywords注释

2.2 VFDB 注释protocol

综合以上我们使用命令行实现如下:

vfdb-fmt  VFDB_setA_pro.fas.gz   >VFDB_setA.pep
makeblastdb  -in VFDB_setA.pep  -dbtype prot  -out  VFDB_setA  -title  VFDB_setA
blast-pipe  -d VFDB_setA  -c 40  A189.faa  A189
blast-utils  hits  -i  30  A189/align/blast.tsv | blast-utils best_hsp -   >A189.vfdb.tsv
tail -n +3  VFs.txt  | tabtk cut -r -f9,8 VFs.txt >vfdb-keyword
tail -n +3  VFs.txt  | tabtk cut -r -f9,6 VFs.txt >vfdb-fuction
zgrep '>' VFDB_setA_pro.fas.gz | vfdb-tab | cut -f1,6 >vfdb-map
tsv-utils  definition  -c 2  -t "VF" vfdb-map A189.vfdb.tsv \
| tsv-utils  definition  -c 3  -t "fuction" vfdb-fuction -  \
| tsv-utils  definition  -c 3 -t 'keyword'  vfdb-keyword  - >A189.vfdb-ann.tsv
tsv-utils tsv2xlsx A189.vfdb-ann.xlsx  VFDB:A189.vfdb-ann.tsv

如果数据库格式化好了其实只需三步:

blast-utils  hits  -i  30  A189/align/blast.tsv | blast-utils best_hsp -   >A189.vfdb.tsv
tsv-utils  definition  -c 2  -t "VF" vfdb-map A189.vfdb.tsv \
| tsv-utils  definition  -c 3  -t "fuction" vfdb-fuction -  \
| tsv-utils  definition  -c 3 -t 'keyword'  vfdb-keyword  - >A189.vfdb-ann.tsv
tsv-utils tsv2xlsx A189.vfdb-ann.xlsx  VFDB:A189.vfdb-ann.tsv

生物信息数据工具封装:序列比对之blast-pipe

一、前言

为了生物学家更加容易使用命令行模式的生物信息工具,在数据分析流程集成工具水平下我们设计了 “小封装” 模式,即封装一些常用的几个工作步骤,尽量使用优化后的默认参数,整个小封装依赖实现的 “tool-utils”, 比如 blast-utilsfastx-utilstsv-utilssam-utils, 这些小程序一般都是 C语言 实现的高性能应用。

二、blast-pipe介绍

小封装 blast-pipe 目的就是解决序列比对BLAST的提交任务, 命令行接口:

$blast-pipe

Program: blast-pipe: blast submit and parse protocol.
Version: 0.0.1
Contact: ZHANG LEI <zhanglei@logicinformatics.com>

Usage:   blast-pipe [options] <sequence> <project>

Options: -t STR    blast type. blastx|blastp|blastn, default [blastp],
                   for special task, can do like this: 'blastn -task megablast'
         -c INT    CPU number, default: [56]
         -e double set evalue cutoff, default: [1e-10]
         -i double set identity cutoff for filter, default: [0]
         -f STR    set outfmt parameter, default:
                   ['6 qseqid stitle pident length mismatch gapopen qstart qend sstart send evalue bitscore']
         -b double set bit score cutoff for filter, default: [60]
         -m INT    set max_target_seqs parameter, default: [1]
         -d STR    swiss-prot database location, default: [/biostack/database/swiss_prot/swiss_prot]

最简单的任务提交模式:

$blast-pipe   A189.faa  Swiss-prot

程序执行序列比对任务:采用了上述默认参数, 序列库也是默认参数:”/biostack/database/swiss_prot/swiss_prot”

稍微复杂点的可以为:

$blast-pipe -c 10 -m 10 -d /biostack/database/vfdb/vfdb A189.faa vfdb

为什么要对序列比对也进行封装呢?

1. 使用'fastx-utils partition' 对序列进行分拆,多进程执行,提高CPU利用率,可有效支持集群模式;
2. 可增加很多辅助操作,将生成的文件直接转换成Excel文档;
3. 按照自己常用需求调整默认参数,减少命令行提交复杂度;

三、最后

Biostack 提供定制小封装服务, 如有需求请联系微信号: biostack

序列功能注释神器: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

otutable-utils之计算OTU表汇总信息-otutable_summary和otutab_stats

本文介绍otutable-utils系列的一个小工具:otutab_summary (Github 链接 ), 就是计算OTU表中每个样本的OTU数和tag数,支持未注释的OTU表。

1. 命令行接口

$ otutable_summary
Usage:  otutable_summary <otutab>
version: 0.0.1

2. 命令行接口

otutable_summary OTU_table.txt  >OTU_table.summary.txt

3. 其它实现

3.1 otutab_stats

otutab_stats 是 Usearch 的命令行工具,接口也比较简单

usearch  -otutab_stats  OTU_table.txt   -output     OTU_table.report.txt

该程序对整个OTU表中样本进行统计分析, 示例结果如下:

   1624812  Reads (1.6M)
        27  Samples
      3565  OTUs

     96255  Counts
     44360  Count  =0  (46.1%)
     14975  Count  =1  (15.6%)
     14774  Count >=10 (15.3%)

       504  OTUs found in all samples (14.1%)
       739  OTUs found in 90% of samples (20.7%)
      1808  OTUs found in 50% of samples (50.7%)

Sample sizes: min 48273, lo 59502, med 60943, mean 60178.2, hi 61976, max 62871

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有.

Last Upate: 2017-11-29 10:18 PM

扩增子序列分析之物种分类: 注释到种水平,也许可以!

一、引言

DADA2 升级至 1.6 版本, 并正式支持 Species assignment, 使用了精确匹配模式判别是否可以分类到种水平,当然种以上分类水平还是要使用RDP算法或者Sintax算法。

DADA2

DADA2 Species assignment 逻辑很简单:

  1. 预处理16S库文件(去除未分类的序列或者分类不明确的序列, 保留属和种水平分类信息), DADA2 提供了rdp_species_assignment_16.fa.gzsilva_species_assignment_v128.fa.gz 两种库。
  2. 序列精确比对,DADA2使用了ShortRead进行序列比对;
  3. 判断,如果所有有效匹配(可设正向和反向选项)在 属+种 分类都一致,那么是可以直接分配该种分类等级, DADA2 提供了允许分配多个的选项,不推荐使用这个模式。

代码段:

assignSpecies <- function(seqs, refFasta, allowMultiple=FALSE, tryRC=FALSE, verbose=FALSE) {
  # Define number of multiple species to return
  if(is.logical(allowMultiple)) {
    if(allowMultiple) keep <- Inf
    else keep <- 1
  } else {
    keep <- as.integer(allowMultiple)
  }
  # Get character vector of sequences
  seqs <- getSequences(seqs)
  # Read in the reference fasta
  refsr <- readFasta(refFasta)
  ids <- as(id(refsr), "character")
  # Crude format check
  if(!length(unlist(strsplit(ids[[1]], "\\s"))) >= 3) {
    if(length(unlist(gregexpr(";", ids[[1]]))) >= 3) {
      stop("Incorrect reference file format for assignSpecies (this looks like a file formatted for assignTaxonomy).")
    } else {
      stop("Incorrect reference file format for assignSpecies.")
    }
  }
  genus <- sapply(strsplit(ids, "\\s"), `[`, 2)
  species <- sapply(strsplit(ids, "\\s"), `[`, 3)
  # Identify the exact hits
  hits <- vector("list", length(seqs))
  lens <- nchar(seqs)
  for(len in unique(lens)) { # Requires all same length sequences
    seqdict <- PDict(seqs[lens==len])
    vhit <- (vcountPDict(seqdict, sread(refsr))>0)
    if(tryRC) vhit <- vhit | (vcountPDict(seqdict, reverseComplement(sread(refsr)))>0)
    hits[lens==len] <- lapply(seq(nrow(vhit)), function(x) vhit[x,])
  }
  # Get genus species return strings
  rval <- cbind(unlist(sapply(hits, mapHits, refs=genus, keep=1)),
                unlist(sapply(hits, mapHits, refs=species, keep=keep)))
  colnames(rval) <- c("Genus", "Species")
  rownames(rval) <- seqs
  if(verbose) cat(sum(!is.na(rval[,"Species"])), "out of", length(seqs), "were assigned to the species level.\n")
  rval
}

该思路可以迁移其它扩增子数据分析平台,比如QIIME2uParse,为了更方面命令行模式下执行该分析,我们重新实现了一个独立版本: species_assign

二、 species_assign程序

2.1 准备数据库

数据库推荐使用silva_species_assignment_v128, 通过下面命令获取序列和分类映射表;

wget https://zenodo.org/record/824551/files/silva_species_assignment_v128.fa.gz
gunzip silva_species_assignment_v128.fa.gz
grep ">" silva_species_assignment_v128.fa | sed 's/>//' | cut -f1,2,3 -d ' ' >silva_species_assignment_v128.map

构建Bowtie索引库:

bowtie-build  example_species_assignment.fa  example_species_assignment
2.2 命令行接口

species_assign 提供了和dada2 assignSpecies 一致的可选项, 即

  1. 正链还是双链匹配模式
  2. 匹配多个物种, 一般不会选择这个参数。

命令行接口:

$ species_assign
Usage: species_assign [options] <map> <bowtie>
Options:
  -b  enable both strand. Default:[plus]
  -m  enable multiple species assignment. Default:[not allowed]
  -v  print version number

可选项:

-b: both strands.
-m: multiple species match.
2.3 程序示例

结合bowtie我们可以轻松完成dada2的物种注释功能:

bowtie  -f -v 0 -a silva_species_assignment_v128 example_seqs.fa 2>/dev/null |\
   species_assign  silva_species_assignment_v128.map -

解读:

bowtie 参数:

-f: 输入格式为fasta;
-v: end-to-end 比对模式, 允许0个错配。
-a:输出所有匹配;
2>/dev/null:将log信息丢掉;

species_assign 程序在 taxon-utils 仓库。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017/11/28 23:41

Linux 命令行小工具之tldr:更好的 man?

一、引言

Linux的伟大成功一个重要的原因就是有各种各样的小工具,比如ls, dutopdu等,可以完成Linux操作系统的各种文件操作、目录操作以及系统管理等,但是,以下几个因素也会成为我们习惯命令行操作的拦路虎:

  1. 可是使用的命令行工具太多了。
  2. 命令行接口(选项)太多了。

当安装这些小工具不是问题的时候,使用这些小工具就成为最大的问题了,首当其冲就是了解命令行的使用说明,最直接的肯定是 man 工具, 如下图,我们可以在命令行下输入man man就可以其实获得关于man的使用说明。

man

图一: man使用说明

当然很多程序, 可以通过直接执行命令行工具本身, 比如samtools 或者 bedtools --help 获取帮助信息。

不过今天介绍是一个小工具tldr()

二、TLDR介绍

tldr-pages 采用一种社区驱动模式, Github仓库地址: https://github.com/tldr-pages/tldr, Web 客户端: https://tldr.ostera.io/

Linux命令行工具运行命令如下:

图二: tldr使用模式

相对于man页面,tldr具有更多的例子。

三、安装

tldr 提供了各种语言的客户端: 下面列出几种比较常用的安装模式:

Node.js 客户端(需要先安装 Node.js) :

npm install -g tldr

Python 客户端

pip install tldr

GO 客户端:

wget https://github.com/pranavraja/tldr/releases/download/v1/tldr_linux_amd64.tar.gz
tar xzvf  tldr_linux_amd64.tar.gz
mv  tldr_linux_amd64 tldr-1.0
export PATH=$PATH:$PWD/tldr-1.0

四、使用

命令行接口(以Node.js客户端为例):

$ tldr

  Usage: tldr command [options]
  Simplified and community-driven man pages
  Options:
    -h, --help               output usage information
    -V, --version            output the version number
    -l, --list               List all commands for the chosen platform in the cache
    -a, --list-all           List all commands in the cache
    -1, --single-column      List single command per line (use with options -l or -a)
    -r, --random             Show a random command
    -e, --random-example     Show a random example
    -f, --render [file]      Render a specific markdown [file]
    -m, --markdown           Output in markdown format
    -o, --os [type]          Override the operating system [linux, osx, sunos]
    --linux                  Override the operating system with Linux
    --osx                    Override the operating system with OSX
    --sunos                  Override the operating system with SunOS
    -t, --theme [theme]      Color theme (simple, base16, ocean)
    -s, --search [keywords]  Search pages using keywords
    -u, --update             Update the local cache
    -c, --clear-cache        Clear the local cache

Examples

    $ tldr tar
    $ tldr du --os=linux
    $ tldr --search "create symbolic link to file"
    $ tldr --list
    $ tldr --list-all
    $ tldr --random
    $ tldr --random-example

To control the cache

    $ tldr --update
    $ tldr --clear-cache

To render a local file (for testing)

    $ tldr --render /path/to/file.md

下面介绍几种功能性可选项:

-l: 列出缓存中指定操作系统的所有命令;
-a: 列出缓存中所有命令;
-m: 输出格式为markdown 格式;
--linux/--osx/--sunos:重载操作系统类型
-t: 颜色主题模式;
-u: 升级本地缓存;
-c:清除本地缓存;

示例:

$ tldr tar

  tar

  Archiving utility.
  Often combined with a compression method, such as gzip or bzip.

  - Create an archive from files:
    tar cf target.tar file1 file2 file3

  - Create a gzipped archive:
    tar czf target.tar.gz file1 file2 file3

  - Extract an archive in a target folder:
    tar xf source.tar -C folder

  - Extract a gzipped archive in the current directory:
    tar xzf source.tar.gz

  - Extract a bzipped archive in the current directory:
    tar xjf source.tar.bz2

  - Create a compressed archive, using archive suffix to determine the compression program:
    tar caf target.tar.xz file1 file2 file3

  - List the contents of a tar file:
    tar tvf source.tar

很方便的一个小工具,如果忘记一个命令具体使用,请使用:tldr command

所支持的不同平台的命令行使用用例,请下载:

tldr-book

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017/11/27 21:34

taxon-utils之NCBI Taxonomy 标识符lineage注释: taxon-translation

一、引言:

本文介绍一个小工具,taxon-translation, 使用场景:

需要将NCBI Taxonomy数据库的数字标识符转换成整个Lineage信息。

taxon-utils 仓库地址: https://github.com/jameslz/taxon-utils

二:应用

2.1 数据库格式化 taxon-db

taxon-db程序将NCBI的Taxonomy数据库的 node.dmpname.dmp 合并成一个文件。

执行方法:

taxon-db  nodes.dmp  names.dmp  >ncbi.map
gzip ncbi.map

ncbi.map 也可以从 taxon-utils 仓库下载, ncbi.map.gz(版本为2017-09-15)链接

2.1 数据库格式化 taxon-translate

taxon-translate 命令行接口比较简单:

$ taxon-translate
Usage: taxon-translate [options] <ncbi-map> <tab>
Options:
  -c INT     target col for taxon translate, default: 1
  -v         print version number

可以通过 -c 指定列进行翻译。

测试文件(test.tsv):

read_1  68      Lysobacter enzymogenes

第一列为序列名称, 第二列为NCBI Taxon ID, 第三列为科学命名。

对第二列进行翻译:

taxon-translate -c 2  ncbi.map.gz test.tsv

我们可以获得:

read_1  69      Lysobacter enzymogenes  \
d:Bacteria;p:Proteobacteria;c:Gammaproteobacteria;  \
o:Xanthomonadales;f:Xanthomonadaceae;g:Lysobacter;s:Lysobacter enzymogenes

很直接,很简单。

该程序每次都要将ncbi.map.gz 读入内存, 因此可以先将ncbi.map.gz对于对应的数据结构序列化,然后直接内存映射,这样可以对taxon-translate进一步优化。

本文材料为 BASE (Biostack Applied bioinformatic SEies ) 课程 Linux Command Line Tools for Life Scientists 材料, 版权归 上海逻捷信息科技有限公司 所有。

Last Update: 2017-11-22 5:56 PM