Kraken和Kraken2（Kraken的迭代版本）是微生物组分析中最常用的软件，与其他同功能的软件相比，在速度快的前提下，准确性也很高。那么，Kraken究竟是如何做到的呢？Kraken2和Kraken相比有哪些改进？今天将在这篇文章中先对Kraken进行详细介绍。（实际上是前阵子组会被“安排”讲解Kraken和Kraken2，我就决定偷懒顺便整理发布到简书啦。）

特别说明
Kraken软件已经出了非常久，所以网上有不少相关资料，因此关于Kraken的解读基本参考网上的资料+文献本身。资料来源：http://jackywu.site/technology/kraken-code-analysis/
（如有侵权删）

从作者说起

这个实验室十分钟情于Genome Biology，Kraken和Kraken2都发表该杂志上。Kraken发表于2014年，作者和通讯正是上图所圈选出来的2位。

可能有人已经知道这位通讯作者了，没错是个大佬了。

Steven Lloyd Salzberg

搜了一下该实验室之前发表的论文，可以看到生信领域常用的软件Bowtie2、TopHat、TopHat2、FLASH和HISAT等都出自该实验室。

Publications

当然了，正如之前所说，Kraken这个软件也是在不断更新之中。

作者Derrick E Wood在2019年发布了Kraken的新版本Kraken2。（嗯，看了作者的照片，果然是学习改变人生）

那么其实在2014-2019年之间，该实验室也有其他人对Kraken进行过迭代更新，在2019年发表了KrakenUniq，没错依然是在Genome Biology。

KrakenUniq与Kraken比较主要是对使用了独特的k-mer counts（也就是对k-mer进行了优化），所以在速度和准确性上得到了一定的提升。不过，这次我们并不会深入讲解KrakenUniq，有兴趣的小伙伴可以自己去读相关的paper。

这次讲解主要想要回答5个问题，也就是参考的资料里面提出来的5个问题：
• 为什么Kraken的分析速度那么快?
• 为什么Kraken的数据库有几百G那么大?
• 为什么Kraken建库的速度非常慢?
• 为什么Kraken数据库的载入速度非常慢?
• Kraken的数据库能否拆分使得其能够分布式运行?

The Kraken sequence classification algorithm

首先来看一下Kraken的基本算法。简单来讲，使用Kraken软件有2步：准备（建库）+鉴定。

准备(建库)
• 建立k-mer（Box1）对应的taxon数据库
• 将数据库和索引文件映射到内存

实际上建库的工作只需要在第一次运行该软件时进行即可，再次使用的时候，因为已经做好了准备工作，所以只需要直接对序列进行鉴定即可。

鉴定
• 将待鉴定序列切成k-mer
• 将k-mer比对到数据库上获得其LCA_taxon（Box2）以及比对上的次数。
• 将上述数据构建成Classification tree ，然后计算每条root-to-leaf上的所有权重和，最大者即为该条序列的分类树。

举个例子，如上图，如果输入的query sequence切割成k-mer后，与LCA进行mapping，最终发现可以map到的k-mer有16条（即标记为紫色、蓝色、橘黄色和红色的k-mer），然后对各个节点进行统计。发现紫色节点有1个k-mer，蓝色的有10个，橘黄色的4个，红色的1个，那么最终就有2条路径。
而紫色-蓝色-橘黄色的这条路径总分为15，而紫色-黑色-红色的总分为2，因此，前者是得分更高的路径。所以这条序列就会被认为是橘黄色节点对应的物种。

Box1：什么是k-mer
k-mer指的是将一条read,连续切割,挨个碱基划动得到的一序列长度为K的核苷酸序列。
比如，以下这条read为例：
ATCGTTGCTTAATGACGTCAGTCGAATGCGATGACGTGACTGACTG
如果是k-mer=13的话
ATCGTTGCTTAAT
TCGTTGCTTAATG
CGTTGCTTAATGA
GTTGCTTAATGAC
……
对基因组进行k-mer分析，可以为我们提供一些信息：
1.基因组大小
2.基因组杂合度
3.基因组重复片段大小

Box2：什么是LCA？
LCA的全称是Lowest Common Ancestor，中文译为最近公共祖先，是指在一个树或者有向无环图中同时拥有x和y作为后代的最深的节点。
例子:
在右图中，x与y的最近公共祖先被标记为深绿色，其他公共祖 先被标记为浅绿色。

计算最近公共祖先和根节点的长度往往是有用的。比如，为了计算树中两个节点x和y之间的距离，可以使用以下方法:分别计算由x到根节点和y到根节点的距离，两者之和减去最近公共祖先到根节点的距 离的两倍即可得到x到y的距离。

建库过程

具体地，从上图中我们可以看到建库过程一共有6步，其中Step4在目前的版本中已经不需要进行了。
根据不同步骤所花费的时间可以发现，建库耗时主要集中在Step3 sort set和Step6 set LCA values。

那么这是为什么呢？我们一步一步来看看建库究竟干了什么。

Step0：Download Database

Standard Kraken Database:
NCBI taxonomic information, the complete genomes in RefSeq for the bacterial, archaeal, and viral domains.
所以可以看到标准的库下载的是细菌、古菌以及病毒的RefSeq数据。
但是实际上我们知道，就人体微生物组数据而言，其实真菌也是很重要的组成部分，因此我们可以自主添加真菌数据库。

Custom Database:

#If you need to modify the taxonomy, edits can be made to the names.dmp and nodes.dmp files in this directory
kraken-build --download-taxonomy --db $DBNAME
kraken-build --download-library bacteria --db $DBNAME
kraken-build --add-to-library chr1.fa --db $DBNAME

Step1：Create k-mer set: Jellyfish -> database.jdb

步骤2就是利用Jllyfish软件切割k-mer，生成database.jdb文件, 文件内容是 “k-mer: count”。

Jllyfish是CBCB（Center for Bioinformatics and Computational Biology）的Guillaume Marçais 和 Carl Kingsford 研发的一款计数 DNA 的 k-mers 的软件。该软件运用 Hash 表来存储数据，同时能多线程运行，速度快，内存消耗小。该软件只能运行在64位的Linux系统下。其文章于2011年发表在杂志 Bioinformatics上。

Step2: reduce database, optional and skipped
这部分顾名思义，就是对k-mer数据库进行了一个优化，缩减大小。

Step3：Sort set: database.kdb + databse.idx
这就是我们刚才说到的特别慢的一步。这一步是干啥的呢？具体可以分为两步：
Step3.1: 对database.jdb进行排序
Step3.2: 生成索引文件的
由于第一步排序使用了快速排序算法（Box3），因此就特别慢。

Box3：什么是快速排序算法？
快速排序算法是在起泡排序的基础上进行改进的一种算法，其实现的基本思想是：通过一次排序将整个无序表分成相互独立的两部分，其中一部分中的数据都比另一部分中包含的数据的值小，然后继续沿用此方法分别对两部分进行同样的操作，直到每一个小部分不可再分，所得到的整个序列就成为了有序序列。
例如，对无序表{49，38，65，97，76，13，27，49}进行快速排序，大致过程为：
首先从表中选取一个记录的关键字作为分割点（称为“枢轴”或者支点，一般选择第一个关键字），例如选取 49；
将表格中大于 49 的放置于 49 的右侧，小于 49 的放置于 49 的左侧，假设完成后的无序表为：{27，38，13，49，65，97，76，49}；
以 49 为支点，将整个无序表分割成了两个部分，分别为{27，38，13}和{65，97，76，49}，继续采用此种方法分别对两个子表进行排序；
前部分子表以 27 为支点，排序后的子表为{13，27，38}，此部分已经有序；后部分子表以 65 为支点，排序后的子表为{49，65，97，76}；
此时前半部分子表中的数据已完成排序；后部分子表继续以 65为支点，将其分割为{49}和{97，76}，前者不需排序，后者排序后的结果为{76，97}；
通过以上几步的排序，最后由子表{13，27，38}、{49}、{49}、{65}、{76，97}构成有序表：{13，27，38，49，49，65，76，97}；
本部分来源：http://data.biancheng.net/view/71.html

那么第二步就是生成索引文件。具体的我们可以看下面这张图。

这里Kraken用到了一个叫做“minimizer”的玩意。所谓的minimizer就是，由于k-mer的切割方式所以导致临近的k-mer实际是很相似的，所以就可以用一个minimizer来表征一组k-mer （我是这么理解的，如有误欢迎指正）。

具体怎么取minimizer，我没有仔细研究，大家可以看下面这篇文章：

当然minimizer的长度就自然会影响到所需要的minimizer的个数以及所需要的存储空间。Kraken默认的是15-bp，但是也可以修改。有些人由于设备的限制可能会建立一个MiniKraken，那么在MiniKraken中作者就采用了13-bp，以保证最后建立的数据库大小在4GB以内。

通过设置minimizer可以发现，在鉴定的时候对于一个k-mer的搜索范围会大大缩小，从而减少运行时间。

Step4：现在已经木有啦
Step5: Sequence ID to taxon map
将SeqID和TxaonID进行匹配。

Step6：Set LCA Value
那么我们重点来讲一下Step6。可以看一下文章中关于设定LCA Value的描述。

其实在这一步就是要将k-mer构建成Taxonomy Tree。并给每一个节点分配一个LCA Value。

这个LCA Value用taxonomic ID number确定的，然后database.kdb中k-mer: count会变为:k-mer:lca_taxon_id。

具体地就是在之前的步骤完成后，会遍历一遍你下载的数据库。然后对于某一物种的序列的所有k-mer，其LCA Value都会设置成这个序列所对应的物种taxonomic ID number。

那么，这个时候就会有一个问题。如果物种A和物种B共有某一段序列，那么这个序列的LCA值要怎么确定呢？

这时候，LCA Value就会依据两个物种的taxonomic ID number进行计算，而这也就是如何形成Taxonomy Tree的过程。

举个粗暴的例子：
比如菌种a1和菌种a2都属于菌属A，并且共有一部分序列，那么假设序列S是a1和a2共有的。
一开始先遍历到了a1，所以序列S的k-mer的LCA Value是a1的taxonomic ID number，然后这个时候轮到了a2，又出现了这段序列S，这时候这些k-mer的LCA Value就会根据已有的LCA Value和a2的taxonomic ID number进行计算，就会变成它们对应的属A的taxonomic ID number。而这个节点就会从最底层往上走一个分类水平，变成了属。
那么如果a1和a2不是共同属于一个属A，而是共同属于科A，那么对应这个共有序列S的所有k-mer的LCA Value就会变成科A的taxonomic ID number，这个节点也就更接近树的根节点了。

这就是建库的全部过程。

鉴定

鉴定的命令十分简单。

kraken --db $DBNAME seqs.fa

当然，为了加速鉴定速度，kraken也提供了preload指令将数据库和索引预先加载到内存里去，具体使用的方法是mmap（Box4）

kraken --preload --db $DBNAME seqs.fa

Box4：什么是mmap？
mmap是一种内存映射文件的方法，即将一个文件或者其它对象映射到进程的地址空间，实现文件磁盘地址和进程虚拟地址空间中一段虚拟地址的一一对映关系。实现这样的映射关系后，进程就可以采用指针的方式读写操作这一段内存，而系统会自动回写页面到对应的文件磁盘上，即完成了对文件的操作而不必再调用read,write等系统调用函数。相反，内核空间对这段区域的修改也直接反映用户空间，从而可以实现不同进程间的文件共享。
常规文件操作需要从磁盘到页缓存再到用户主存的两次数据拷贝。而mmap操控文件，只需要从磁盘到用户主存的一次数据拷贝过程。
（Source: https://www.cnblogs.com/huxiao-tee/p/4660352.html）

那么最后我们来回到最开始的5个问题：
• 为什么Kraken的分析速度那么快? 使用mmap映射+索引 (C++ &Perl)
• 为什么Kraken的数据库有几百G那么大? 原始database大+k-mer大小
• 为什么Kraken建库的速度非常慢? 数据库qsort排序+创建index+遍历设置LCA value
• 为什么Kraken数据库的载入速度非常慢? 数据库大
• Kraken的数据库能否拆分使得其能够分布式运行? 原代码不能，但是可以DIY?

好了，今天就讲到这里啦，下次分享一下Kraken2的原理～