Kraken2：宏基因组快速物种注释神器

简介

kraken是基于k-mer精确比对，并采用最LCA投票结果快速宏基因组DNA序列进行物种注释的软件。

图. Kraken2分类基本原理

该文章于2014年发表于Genome Biology，目前引用过两千次[1]。详见《Kraken：使用精确比对的超快速宏基因组序列分类软件》

kraken1对内存要求大限制了部分用户使用，kraken2优化了这一问题，而且输出结果格式更友好，下游分析更方便，文章于2019年再次发表于Genome Biology，一年被引231次[2]。

主页：
https://github.com/DerrickWood/kraken2

帮助文档：
https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown

软件安装

kraken2程序主要包括kraken2 和 kraken2-build两条命令，还依赖一些小程序和脚本。所有程序将安装于同一目录。如果安装后你移动安装文件格，需要修改kraken2脚本中的$KRAKEN2_DIR对应新的位置。

推荐conda安装，一条命令搞定

conda install kraken2

另一种安装方法：可选git或wget下载，不推荐

cd ~/software
git clone git@github.com:DerrickWood/kraken2.git

# 或wget下载并解压
wget https://github.com/DerrickWood/kraken2/archive/master.zip
unzip master.zip

# 安装，设置你要安装的位置
KRAKEN2_DIR=~/soft/kraken2
./install_kraken2.sh $KRAKEN2_DIR

# 添加到环境变量
cp $KRAKEN2_DIR/kraken2{,-build,-inspect} $HOME/bin

版本和升级

记录软件版本：

kraken2 --version # Kraken version 2.0.9-beta

升级

conda update kraken2 -c bioconda

kraken2 --version # 2.1.1 Copyright 2013-2020

Kraken2数据库

Kraken2数据库至少包括3个文件

  • hash.k2d: 物种地图，即所有序列与物种的数据库

  • opts.k2d: 数据库构建的信息

  • taxo.k2d: 数据库的分类学信息

以上文件为快速读取，全为二进制文件。
如果仅使用kraken2，除以上三个文件外，其它的文件在空间有限下可以删除。如果要使用bracken时，仍需保留其他文件，用于构建索引等。

标准数据库构建

需要下载50G数据，过程将消耗超100GB的空间。

# 下载数据库，先设置存放位置
DBNAME=~/db/kraken2/201127
mkdir -p $DBNAME
# 可调线程数据，加速构建过程
kraken2-build --standard --threads 24 --db $DBNAME

此步下载数据>50GB，下载速度由网络决定。索引时间4小时33分，多线程最快35min完成

标准模式下只下载5种数据库：古菌archaea、细菌bacteria、人类human、载体UniVec_Core、病毒viral

个性数据库构建

更多数据库构建的参数，详见kraken2-build -h

使用方法：kraken2-build 任务类型 可选参数
Usage: kraken2-build [task option] [options]

以下任务选择其中之一
Task options (exactly one must be selected):
--download-taxonomy 下载物种信息Download NCBI taxonomic information
--download-library TYPE 下载数据库类型，包括14种Download partial library
(TYPE = one of "archaea", "bacteria", "plasmid",
"viral", "human", "fungi", "plant", "protozoa",
"nr", "nt", "env_nr", "env_nt", "UniVec",
"UniVec_Core")
--special TYPE 特殊数据库，包括3种Download and build a special database
(TYPE = one of "greengenes", "silva", "rdp")
--add-to-library FILE 添加文件到库中 Add FILE to library
--build 建索引Create DB from library
(requires taxonomy d/l'ed and at least one file
in library)
--clean 清理Remove unneeded files from a built database
--standard 下载标准库Download and build default database
--help 帮助Print this message
--version 版本信息Print version information

Options:
--db NAME 库位置Kraken 2 DB/library name (mandatory except for
--help/--version)
--threads # 线程数Number of threads (def: 1)
--kmer-len NUM K-mer长度 length in bp/aa (build task only;
def: 35 nt, 15 aa)
--minimizer-len NUM Minimizer length in bp/aa (build task only;
def: 31 nt, 15 aa)
--minimizer-spaces NUM Number of characters in minimizer that are
ignored in comparisons (build task only;
def: 6 nt, 0 aa)
--protein Build a protein database for translated search
--no-masking Used with --standard/--download-library/
--add-to-library to avoid masking low-complexity
sequences prior to building; masking requires
dustmasker or segmasker to be installed in PATH,
which some users might not have.
--max-db-size NUM Maximum number of bytes for Kraken 2 hash table;
if the estimator determines more would normally be
needed, the reference library will be downsampled
to fit. (Used with --build/--standard/--special)

kraken2的库类型包括：”archaea”, “bacteria”, “plasmid”, “viral”, “human”, “fungi”, “plant”, “protozoa”, “nr”, “nt”, “env_nr”, “env_nt”, “UniVec”, “UniVec_Core”

archaea bacteria plasmid viral human fungi plant protozoa nr nt env_nr env_nt UniVec

数据越多，需要的内存也多，请谨慎选择。

以下实例供参考

# 下载数据库，先设置存放位置
DBNAME=~/db/kraken2/201201
mkdir -p DBNAME
cd $DBNAME

# 下载物种注释，259
kraken2-build --download-taxonomy --threads 24 --db $DBNAME

# 下载非默认数据库中的真菌库 1.33G
kraken2-build --download-library fungi --threads 24 --db $DBNAME

# 批量下载，非标准数据库
for i in archaea bacteria plasmid viral human fungi plant protozoa nr nt env_nr env_nt UniVec; do
kraken2-build --download-library $i --threads 24 --db $DBNAME
done

# 确定的库建索引
kraken2-build --build --threads 24 --db $DBNAME

注：NCBI上的三个重要的数据库—NR/NT，Taxonomy和RefSeq。NR(Non-Redundant Protein Sequence Database)非冗余蛋白库，所有GenBank+EMBL+DDBJ+PDB中的非冗余蛋白序列，对于所有已知的或可能的编码序列，NR记录中都给出了相应的氨基酸序列（通过已知或可能的读码框推断而来）以及专门蛋白数据库中的序列号。NR库相当于一个以核酸序列为基础的交叉索引，将核酸数据和蛋白数据联系起来。NT(Nucleotide Sequence Database),核酸序列数据库，是NR库的子集。NCBI的分类数据库，包括大于7万余个物种的名字和种系，这些物种都至少在遗传数据库中有一条核酸或蛋白序列。其目的是为序列数据库建立一个一致的种系发生分类学。RefSeq(the reference sequence database,https://www.ncbi.nlm.nih.gov/refseq/ ).参考序列数据库，包含RefSeq_genomic(NCBI genomic reference sequences)，RefSeq_protein(NCBI protein reference sequences)和RefSeq transpans(NCBI transpans reference sequences)具有生物意义上的非冗余基因,转录本和蛋白质序列，是经过NCBI和其他组织校正的数据库，使用人类基因命名委员会定义的术语，并且包括了官方的基因符号和可选的符号。[3]

序列分类

kraken2 --db $DBNAME seqs.fa

其它参数解析：

  • 多线程：—threads NUM

  • 快速操作：—quick 检索第一次匹配即停止；—min-hits 多个匹配结果才确定；

  • 序列过滤：分类—classified-out、末分类的—unclassified-out结果输出到文件

  • 输出结果：为标准输出，可|重定向，也可—output写入文件

  • 输入文件：默认为fasta，可—fastq-input指定为fastq

  • 输入压缩文件：—gzip-compressed or —bzip2-compressed

  • 输入文件自动检测：默认为自动检测，你可以帮忙指定类型 —fasta-input, —fastq-input, —gzip-compressed, and/or —bzip2-compressed

  • 双端数据：2开始检测双端数据了，  kraken2 —paired —classified-out cseqs#.fq seqs_1.fq seqs_2.fq

  • --use-names 输出物种名(ID)

  • --report输出报告

  • --report-zero-counts 输出所有物种，方便样本合并比较(sort -k5,5n)

  • --use-mpa-style输出metaphlan2格式

输出文件格式

Kraken标准输出格式

五列表

  1. C/U代表分类classified或非分类unclassifed

  2. 序列ID

  3. 物种注释

  4. 比序列注释的区域，如98|94代表左端98bp，右端94bp比对至数据库

  5. LCA比对结果，如”562:13 561:4”代表13 k-mer比对至物种#562，4 k-mer比对至#561物种

报告输出格式

包括6列，方便整理下游分析。

  1. 百分比

  2. count

  3. count最优

  4. (U)nclassified, (R)oot, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. “G2”代表位于属一种间

  5. NCBI物种ID

  6. 科学物种名

自定义数据库

默认的数据库修改：names.dmp、nodes.dmp和*.accession2taxid。

—download-library 下载数据库

—add-to-library

—skip-maps 跳过某数据库

常用数据库

  • archaea: 古菌686M, RefSeq complete archaeal genomes/proteins

  • bacteria: 细菌47G, RefSeq complete bacterial genomes/proteins

  • plasmid: 质粒，RefSeq plasmid nucleotide/protein sequences

  • viral: 病毒262M, RefSeq complete viral genomes/proteins

  • human: 人3.1G, GRCh38 human genome/proteins

  • fungi: 真菌，RefSeq complete fungal genomes/proteins

  • plant: 植物，RefSeq complete plant genomes/proteins

  • protozoa: 原始动物，RefSeq complete protozoan genomes/proteins

  • nr: 非冗余蛋白库，NCBI non-redundant protein database

  • nt: 非冗余核酸库，NCBI non-redundant nucleotide database

  • env_nr: 非冗余环境蛋白，NCBI non-redundant protein database with sequences from large environmental sequencing projects

  • env_nt: 非冗余环境核酸，NCBI non-redundant nucleotide database with sequences from large environmental sequencing projects

  • UniVec: 常用污染序列，如载体、头、引物等NCBI-supplied database of vector, adapter, linker, and primer sequences that may be contaminating sequencing projects and/or assemblies

  • UniVec_Core:载体核心库2M，用于去除污染序列, A subset of UniVec chosen to minimize false positive hits to the vector database

标准模式下只下载5种数据库：古菌archaea、细菌bacteria、人类human、载体UniVec_Core、病毒viral

我们可以手动下载指定数据库

kraken2-build --download-library plasmid --db $DBNAME
# Downloading plasmid files from FTP...
kraken2-build --download-library fungi --db $DBNAME

下载后需要--build来建索引，nr和env_nr需要—protein参数, 而UniVec 和UniVec_Core不能用—protein选项

（可选）自行添加基因组，需要满足以下两点要求

  • fasta格式，可以多个文件

  • 必须包括NCBI物种ID，如>sequence16|kraken:taxid|32630  Adapter sequence

# 添加单个序列
kraken2-build --add-to-library chr1.fa --db $DBNAME

# 批量添加序列
for file in chr*.fa
do
kraken2-build --add-to-library $file --db $DBNAME
done

数据库确定后，建索引

kraken2-build --build --db $DBNAME

  • —kmer-len 和 minimizer-len，要求核酸至少31，蛋白至少15

  • kraken2-build —help查看更多参数

非NCBI数据库

支持核酸数据库Greengene、RDP、SILVA数据库制作。它们的物种注释不遵循NCBI的Taxonomy标准。

  • Greengenes (greengenes), 所有可用的16S

  • RDP (rdp), 细菌和古菌16S

  • SILVA (silva), 核糖体小亚基99%非冗余NR99序列集

以里以SILVA为例，因为它不仅包括16S，还包括18S，更适合宏基因组数据

kraken2-build --db $DBNAME --special silva

检查数据库内容(可选)

kraken2-inspect —db EXAMPLE_DB | head -5

更多需求，请阅读官方帮助文档：https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown

常见问题

找不到Perl模块

kraken2 --version

错误信息：Fcntl.c: loadable library and perl binaries are mismatched (got handshake key 0xdb00080, needed 0xde00080)

解决方法：手动设置PERL5LIB

# 查看软件安装目录
type kraken2
kraken2 is hashed (/mnt/bai/yongxin/miniconda2/envs/kraken2/bin/kraken2)

# 设置环境位置env
e=/mnt/bai/yongxin/miniconda2/envs/kraken2
PERL5LIB=${e}/lib/5.26.2:${e}/lib/5.26.2/x86_64-linux-thread-multi

找不到libiconv.so.2

方法1. 直接conda安装库

conda install libiconv

方法2. 查找安装位置，链接到env的lib目录

locate libiconv.so.2
ln -s 找到的位置 软件环境的lib

Reference

  1. Wood, Derrick E., and Steven L. Salzberg. “Kraken: ultrafast metagenomic sequence classification using exact alignments.” Genome biology 15.3 (2014): R46.

  2. Derrick E. Wood, Jennifer Lu & Ben Langmead. (2019). Improved metagenomic analysis with Kraken 2. Genome Biology 20, 257, doi: https://doi.org/10.1186/s13059-019-1891-0

  3. 三种NCBI常见数据库 https://www.sohu.com/a/211785974_99908466

赞 (0)

相关推荐