kraken2 微生物物种分类

kraken 是微生物组分析进行物种分类的工具，目前已经是第二代 kraken2 了。kraken2 对比 kraken 重点优化了数据库创建速度和数据库大小，以及分类速度。

kraken 用 k-mer 方法对输入数据的每一条序列进行分类分析。将每一条序列分成多个 k-mers, 每个 k-mer 在分类数据库寻找 LCA (lowest common ancestor), 序列所有 k-mers 所在的分类及其祖先组成一个分类树————属于总分类树的子集，这个分类树每个节点的权重是序列 k-mers 分配到该分类的次数。分类树上每个 RTL (root-to-leaf) 路径得分是这些权重总分，得分最高的路径是分类路径。然后将该路径分类 left 标签分配给该序列。如果有多个路径得分相等，那么他们的共同祖先标签分配给该序列。
从上面原理可以看出，用 kraken2 分析时提供合适的分类数据库至关重要。比如你期望分类到物种（species）水平，那么 silva 的 16S 数据库是不合适的，因为 silva 数据库本身只分类到属（genus）。如果你提供的分类数据库根本不包含样品主要物种，那么分析肯定是错误的。
完成 kraken2 分析得到每个序列的分类，再用 bracken 对各分类进行丰度估计。

kraken 原理

分类数据库
分类数据库可以用 kraken2-build 命令创建，也可以从 Index zone by BenLangmead 下载现成的。默认下载的包含 50, 75, 100, 150, 200, 250, 300 k-mers 索引。

kraken2 提供现成的数据库下载

用 kraken2-build 命令创建时需要安装标记低复杂度区域的 dustmasker 工具，如果没有用 --no-masking 取消标记，否则工具会出错。创建标准的数据库直接用 --standard 命令，这将下载 NCBI 分类数据库，包含细菌、古细菌、病毒、人类和一些已知质粒载体。

kraken2-build --standard --threads 12 --no-masking --db ${db_dir}

普通数据库先下载 NCBI 分类数据，然后下载参考序列数据（archaea, bacteria, viral, human, fungi, ...）。这一步可以下载多个数据库到同一目录，最后一起创建成分类数据库。

kraken2-build --download-taxonomy --db ${db_dir}
kraken2-build --no-masking --download-library bacteria --db ${db_dir}
kraken2-build --build --threads 8 --db ${db_dir}

除 NCBI 分类数据库，也可以创建别的，比如 silva 16S 数据库，使用 --special 参数设置对应的 "TYPE". 比如 "silva/greengenes/rdp"，这些数据库是 kraken2 下载对应数据库，并处理为 kraken2 可用的格式。

kraken2-build --db $DBNAME --special TYPE

如果不是下载处理好的，需要自己 bracken-build 命令创建 bracken 定量数据库。

bracken-build -d kraken2db -t 12

分类与定量
kraken2 输入是需要分类的序列，可以是 fasta 格式也可以是 fastq 格式。

kraken2 --paired --gzip-compressed --use-names --threads 8 --output ${kraken2_dir}/${sample}.kraken \
--report ${kraken2_dir}/${sample}.kreport --db ${kraken2_db} \
${clean_dir}/${sample}_R1.fq.gz ${clean_dir}/${sample}_R2.fq.gz

设置 --confidence 给每个序列分类增加过滤阈值，根据 From defaults to databases: parameter and database choice dramatically impact the performance of metagenomic taxonomic classification tools 这个参数对结果有比较明显影响。
"confidence" 的得分计算：每个序列的每个分类得分为 其中 是该序列分配到某个 LCA 的 k-mer 数目， 是该序列不包含未知（ambiguous）碱基的 k-mer 数目。以下面第 5 条 LCA 记录为例，#561 是 #562 的父节点，#562 的分数是

#561 的分数是

如果设定的 --confidence 大于 16/21 那么 kraken2 会将序列分类到 #561; 如果大于 20/21 那么序列会是未分类。

kraken2 的输出每条序列分类结果，每一列内容是：

1. 第一列 "C/U" 分别表示该 reads 是否被分类
2. 第二列是序列 ID
3. 第三列是分类 ID, 如果 0 表示未分类，比如：Lactobacillus jensenii (taxid 109790)
4. 第四列是序列长度，双端测序的话用 "|" 分割，比如："100|98"
5. 第五列是序列 kmer 的 LCA 记录，比如 "562:13 561:4 A:31 0:1 562:3" 表示：
   • 前 13 kmer 分类到分类 ID #562
   • 接下去 4 kmer 分类到 #561
   • 接下去 31 kmer 包含 ambiguous 碱基
   • 下一个 kmer 不在分类数据库
   • 最后 3 kmer 分类到 #562
     注意如果是双端测序会出现 |:| 表示一个 reads 的结束和另一个 reads 的开始

用 kraken2 分类结果进行 bracken 定量分析。

bracken -d ${kraken2_db} -i ${kraken2_dir}/${sample}.kreport \
-o ${kraken2_dir}/${sample}.bracken -t 8

分析输出 *.bracken 后缀结果，未分类或低于分类层级的 reads 会被移除，剩下 reads 计算百分比。

name    taxonomy_id     taxonomy_lvl    kraken_assigned_reads   added_reads     new_est_reads   fraction_total_reads
Aureimonas sp. AU20     1349819 S       33      1       34      0.00000
Granulicatella adiacens 46124   S       8       0       8       0.00000
Bacillus velezensis     492670  S       35      23      58      0.00000
Clostridium manihotivorum       2320868 S       8       0       8       0.00000
Candidatus Fonsibacter ubiquis  1925548 S       11      0       11      0.00000

输出 *bracken_species.kreport 文件按分类学从大到小统计序列占比。至此完成了分类与定量，另外说一下 KrakenTools 提供了 kraken 分析相关的小工具，比如提取指定物种的序列。

参考资料
Wood, Derrick E., and Steven L. Salzberg. "Kraken: ultrafast metagenomic sequence classification using exact alignments." Genome biology 15.3 (2014): 1-12.
Home · DerrickWood/kraken2 Wiki · GitHub