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

简介kraken是基于k-mer精确比对，并采用最LCA投票结果快速宏基因组DNA序列进行物种注释的软件。图. Kraken2分类基本原理该文章于2014年发表于Genome Biolo...

文章共4,915字 · 阅读需要大约17分钟

一键AI生成摘要，助你高效阅读

问答

刘永鑫Adam

36501人浏览 · 2021-01-04 07:00:00

刘永鑫Adam · 2021-01-04 07:00:00 发布

简介

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

图. Kraken2分类基本原理

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

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

Cell：20种宏基因组学物种分类工具大比拼

主页：
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标准输出格式

五列表

C/U代表分类classified或非分类unclassifed
序列ID
物种注释
比序列注释的区域，如98|94代表左端98bp，右端94bp比对至数据库
LCA比对结果，如”562:13 561:4”代表13 k-mer比对至物种#562，4 k-mer比对至#561物种

报告输出格式

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

百分比
count
count最优
(U)nclassified, (R)oot, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. “G2”代表位于属一种间
NCBI物种ID
科学物种名

自定义数据库

默认的数据库修改：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

Wood, Derrick E., and Steven L. Salzberg. “Kraken: ultrafast metagenomic sequence classification using exact alignments.” Genome biology 15.3 (2014): R46.
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
三种NCBI常见数据库 https://www.sohu.com/a/211785974_99908466

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。