Kraken和Kraken2（Kraken的迭代版本）是微生物組分析中最常用的軟件，與其他同功能的軟件相比，在速度快的前提下，準(zhǔn)確性也很高。那么，Kraken究竟是如何做到的呢？Kraken2和Kraken相比有哪些改進(jìn)？今天將在這篇文章中先對(duì)Kraken進(jìn)行詳細(xì)介紹。（實(shí)際上是前陣子組會(huì)被“安排”講解Kraken和Kraken2，我就決定偷懶順便整理發(fā)布到簡(jiǎn)書啦。）

特別說(shuō)明
Kraken軟件已經(jīng)出了非常久，所以網(wǎng)上有不少相關(guān)資料，因此關(guān)于Kraken的解讀基本參考網(wǎng)上的資料+文獻(xiàn)本身。資料來(lái)源：http://jackywu.site/technology/kraken-code-analysis/
（如有侵權(quán)刪）

從作者說(shuō)起

這個(gè)實(shí)驗(yàn)室十分鐘情于Genome Biology，Kraken和Kraken2都發(fā)表該雜志上。Kraken發(fā)表于2014年，作者和通訊正是上圖所圈選出來(lái)的2位。

可能有人已經(jīng)知道這位通訊作者了，沒錯(cuò)是個(gè)大佬了。

Steven Lloyd Salzberg

搜了一下該實(shí)驗(yàn)室之前發(fā)表的論文，可以看到生信領(lǐng)域常用的軟件Bowtie2、TopHat、TopHat2、FLASH和HISAT等都出自該實(shí)驗(yàn)室。

Publications

當(dāng)然了，正如之前所說(shuō)，Kraken這個(gè)軟件也是在不斷更新之中。

作者Derrick E Wood在2019年發(fā)布了Kraken的新版本Kraken2。（嗯，看了作者的照片，果然是學(xué)習(xí)改變?nèi)松?/p>

那么其實(shí)在2014-2019年之間，該實(shí)驗(yàn)室也有其他人對(duì)Kraken進(jìn)行過(guò)迭代更新，在2019年發(fā)表了KrakenUniq，沒錯(cuò)依然是在Genome Biology。

KrakenUniq與Kraken比較主要是對(duì)使用了獨(dú)特的k-mer counts（也就是對(duì)k-mer進(jìn)行了優(yōu)化），所以在速度和準(zhǔn)確性上得到了一定的提升。不過(guò)，這次我們并不會(huì)深入講解KrakenUniq，有興趣的小伙伴可以自己去讀相關(guān)的paper。

這次講解主要想要回答5個(gè)問題，也就是參考的資料里面提出來(lái)的5個(gè)問題：
? 為什么Kraken的分析速度那么快?
? 為什么Kraken的數(shù)據(jù)庫(kù)有幾百G那么大?
? 為什么Kraken建庫(kù)的速度非常慢?
? 為什么Kraken數(shù)據(jù)庫(kù)的載入速度非常慢?
? Kraken的數(shù)據(jù)庫(kù)能否拆分使得其能夠分布式運(yùn)行?

The Kraken sequence classification algorithm

首先來(lái)看一下Kraken的基本算法。簡(jiǎn)單來(lái)講，使用Kraken軟件有2步：準(zhǔn)備（建庫(kù)）+鑒定。

準(zhǔn)備(建庫(kù))
? 建立k-mer（Box1）對(duì)應(yīng)的taxon數(shù)據(jù)庫(kù)
? 將數(shù)據(jù)庫(kù)和索引文件映射到內(nèi)存

實(shí)際上建庫(kù)的工作只需要在第一次運(yùn)行該軟件時(shí)進(jìn)行即可，再次使用的時(shí)候，因?yàn)橐呀?jīng)做好了準(zhǔn)備工作，所以只需要直接對(duì)序列進(jìn)行鑒定即可。

鑒定
? 將待鑒定序列切成k-mer
? 將k-mer比對(duì)到數(shù)據(jù)庫(kù)上獲得其LCA_taxon（Box2）以及比對(duì)上的次數(shù)。
? 將上述數(shù)據(jù)構(gòu)建成Classification tree ，然后計(jì)算每條root-to-leaf上的所有權(quán)重和，最大者即為該條序列的分類樹。

舉個(gè)例子，如上圖，如果輸入的query sequence切割成k-mer后，與LCA進(jìn)行mapping，最終發(fā)現(xiàn)可以map到的k-mer有16條（即標(biāo)記為紫色、藍(lán)色、橘黃色和紅色的k-mer），然后對(duì)各個(gè)節(jié)點(diǎn)進(jìn)行統(tǒng)計(jì)。發(fā)現(xiàn)紫色節(jié)點(diǎn)有1個(gè)k-mer，藍(lán)色的有10個(gè)，橘黃色的4個(gè)，紅色的1個(gè)，那么最終就有2條路徑。
而紫色-藍(lán)色-橘黃色的這條路徑總分為15，而紫色-黑色-紅色的總分為2，因此，前者是得分更高的路徑。所以這條序列就會(huì)被認(rèn)為是橘黃色節(jié)點(diǎn)對(duì)應(yīng)的物種。

Box1：什么是k-mer
k-mer指的是將一條read,連續(xù)切割,挨個(gè)堿基劃動(dòng)得到的一序列長(zhǎng)度為K的核苷酸序列。
比如，以下這條read為例：
ATCGTTGCTTAATGACGTCAGTCGAATGCGATGACGTGACTGACTG
如果是k-mer=13的話
ATCGTTGCTTAAT
TCGTTGCTTAATG
CGTTGCTTAATGA
GTTGCTTAATGAC
……
對(duì)基因組進(jìn)行k-mer分析，可以為我們提供一些信息：
1.基因組大小
2.基因組雜合度
3.基因組重復(fù)片段大小

Box2：什么是LCA？
LCA的全稱是Lowest Common Ancestor，中文譯為最近公共祖先，是指在一個(gè)樹或者有向無(wú)環(huán)圖中同時(shí)擁有x和y作為后代的最深的節(jié)點(diǎn)。
例子:
在右圖中，x與y的最近公共祖先被標(biāo)記為深綠色，其他公共祖 先被標(biāo)記為淺綠色。

計(jì)算最近公共祖先和根節(jié)點(diǎn)的長(zhǎng)度往往是有用的。比如，為了計(jì)算樹中兩個(gè)節(jié)點(diǎn)x和y之間的距離，可以使用以下方法:分別計(jì)算由x到根節(jié)點(diǎn)和y到根節(jié)點(diǎn)的距離，兩者之和減去最近公共祖先到根節(jié)點(diǎn)的距 離的兩倍即可得到x到y(tǒng)的距離。

建庫(kù)過(guò)程

具體地，從上圖中我們可以看到建庫(kù)過(guò)程一共有6步，其中Step4在目前的版本中已經(jīng)不需要進(jìn)行了。
根據(jù)不同步驟所花費(fèi)的時(shí)間可以發(fā)現(xiàn)，建庫(kù)耗時(shí)主要集中在Step3 sort set和Step6 set LCA values。

那么這是為什么呢？我們一步一步來(lái)看看建庫(kù)究竟干了什么。

Step0：Download Database

Standard Kraken Database:
NCBI taxonomic information, the complete genomes in RefSeq for the bacterial, archaeal, and viral domains.
所以可以看到標(biāo)準(zhǔn)的庫(kù)下載的是細(xì)菌、古菌以及病毒的RefSeq數(shù)據(jù)。
但是實(shí)際上我們知道，就人體微生物組數(shù)據(jù)而言，其實(shí)真菌也是很重要的組成部分，因此我們可以自主添加真菌數(shù)據(jù)庫(kù)。

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文件, 文件內(nèi)容是 “k-mer: count”。

Jllyfish是CBCB（Center for Bioinformatics and Computational Biology）的Guillaume Mar?ais 和 Carl Kingsford 研發(fā)的一款計(jì)數(shù) DNA 的 k-mers 的軟件。該軟件運(yùn)用 Hash 表來(lái)存儲(chǔ)數(shù)據(jù)，同時(shí)能多線程運(yùn)行，速度快，內(nèi)存消耗小。該軟件只能運(yùn)行在64位的Linux系統(tǒng)下。其文章于2011年發(fā)表在雜志 Bioinformatics上。

Step2: reduce database, optional and skipped
這部分顧名思義，就是對(duì)k-mer數(shù)據(jù)庫(kù)進(jìn)行了一個(gè)優(yōu)化，縮減大小。

Step3：Sort set: database.kdb + databse.idx
這就是我們剛才說(shuō)到的特別慢的一步。這一步是干啥的呢？具體可以分為兩步：
Step3.1: 對(duì)database.jdb進(jìn)行排序
Step3.2: 生成索引文件的
由于第一步排序使用了快速排序算法（Box3），因此就特別慢。

Box3：什么是快速排序算法？
快速排序算法是在起泡排序的基礎(chǔ)上進(jìn)行改進(jìn)的一種算法，其實(shí)現(xiàn)的基本思想是：通過(guò)一次排序?qū)⒄麄€(gè)無(wú)序表分成相互獨(dú)立的兩部分，其中一部分中的數(shù)據(jù)都比另一部分中包含的數(shù)據(jù)的值小，然后繼續(xù)沿用此方法分別對(duì)兩部分進(jìn)行同樣的操作，直到每一個(gè)小部分不可再分，所得到的整個(gè)序列就成為了有序序列。
例如，對(duì)無(wú)序表{49，38，65，97，76，13，27，49}進(jìn)行快速排序，大致過(guò)程為：
首先從表中選取一個(gè)記錄的關(guān)鍵字作為分割點(diǎn)（稱為“樞軸”或者支點(diǎn)，一般選擇第一個(gè)關(guān)鍵字），例如選取 49；
將表格中大于 49 的放置于 49 的右側(cè)，小于 49 的放置于 49 的左側(cè)，假設(shè)完成后的無(wú)序表為：{27，38，13，49，65，97，76，49}；
以 49 為支點(diǎn)，將整個(gè)無(wú)序表分割成了兩個(gè)部分，分別為{27，38，13}和{65，97，76，49}，繼續(xù)采用此種方法分別對(duì)兩個(gè)子表進(jìn)行排序；
前部分子表以 27 為支點(diǎn)，排序后的子表為{13，27，38}，此部分已經(jīng)有序；后部分子表以 65 為支點(diǎn)，排序后的子表為{49，65，97，76}；
此時(shí)前半部分子表中的數(shù)據(jù)已完成排序；后部分子表繼續(xù)以 65為支點(diǎn)，將其分割為{49}和{97，76}，前者不需排序，后者排序后的結(jié)果為{76，97}；
通過(guò)以上幾步的排序，最后由子表{13，27，38}、{49}、{49}、{65}、{76，97}構(gòu)成有序表：{13，27，38，49，49，65，76，97}；
本部分來(lái)源：http://data.biancheng.net/view/71.html

那么第二步就是生成索引文件。具體的我們可以看下面這張圖。

這里Kraken用到了一個(gè)叫做“minimizer”的玩意。所謂的minimizer就是，由于k-mer的切割方式所以導(dǎo)致臨近的k-mer實(shí)際是很相似的，所以就可以用一個(gè)minimizer來(lái)表征一組k-mer （我是這么理解的，如有誤歡迎指正）。

具體怎么取minimizer，我沒有仔細(xì)研究，大家可以看下面這篇文章：

當(dāng)然minimizer的長(zhǎng)度就自然會(huì)影響到所需要的minimizer的個(gè)數(shù)以及所需要的存儲(chǔ)空間。Kraken默認(rèn)的是15-bp，但是也可以修改。有些人由于設(shè)備的限制可能會(huì)建立一個(gè)MiniKraken，那么在MiniKraken中作者就采用了13-bp，以保證最后建立的數(shù)據(jù)庫(kù)大小在4GB以內(nèi)。

通過(guò)設(shè)置minimizer可以發(fā)現(xiàn)，在鑒定的時(shí)候?qū)τ谝粋€(gè)k-mer的搜索范圍會(huì)大大縮小，從而減少運(yùn)行時(shí)間。

Step4：現(xiàn)在已經(jīng)木有啦
Step5: Sequence ID to taxon map
將SeqID和TxaonID進(jìn)行匹配。

Step6：Set LCA Value
那么我們重點(diǎn)來(lái)講一下Step6?？梢钥匆幌挛恼轮嘘P(guān)于設(shè)定LCA Value的描述。

其實(shí)在這一步就是要將k-mer構(gòu)建成Taxonomy Tree。并給每一個(gè)節(jié)點(diǎn)分配一個(gè)LCA Value。

這個(gè)LCA Value用taxonomic ID number確定的，然后database.kdb中k-mer: count會(huì)變?yōu)?k-mer:lca_taxon_id。

具體地就是在之前的步驟完成后，會(huì)遍歷一遍你下載的數(shù)據(jù)庫(kù)。然后對(duì)于某一物種的序列的所有k-mer，其LCA Value都會(huì)設(shè)置成這個(gè)序列所對(duì)應(yīng)的物種taxonomic ID number。

那么，這個(gè)時(shí)候就會(huì)有一個(gè)問題。如果物種A和物種B共有某一段序列，那么這個(gè)序列的LCA值要怎么確定呢？

這時(shí)候，LCA Value就會(huì)依據(jù)兩個(gè)物種的taxonomic ID number進(jìn)行計(jì)算，而這也就是如何形成Taxonomy Tree的過(guò)程。

舉個(gè)粗暴的例子：
比如菌種a1和菌種a2都屬于菌屬A，并且共有一部分序列，那么假設(shè)序列S是a1和a2共有的。
一開始先遍歷到了a1，所以序列S的k-mer的LCA Value是a1的taxonomic ID number，然后這個(gè)時(shí)候輪到了a2，又出現(xiàn)了這段序列S，這時(shí)候這些k-mer的LCA Value就會(huì)根據(jù)已有的LCA Value和a2的taxonomic ID number進(jìn)行計(jì)算，就會(huì)變成它們對(duì)應(yīng)的屬A的taxonomic ID number。而這個(gè)節(jié)點(diǎn)就會(huì)從最底層往上走一個(gè)分類水平，變成了屬。
那么如果a1和a2不是共同屬于一個(gè)屬A，而是共同屬于科A，那么對(duì)應(yīng)這個(gè)共有序列S的所有k-mer的LCA Value就會(huì)變成科A的taxonomic ID number，這個(gè)節(jié)點(diǎn)也就更接近樹的根節(jié)點(diǎn)了。

這就是建庫(kù)的全部過(guò)程。

鑒定

鑒定的命令十分簡(jiǎn)單。

kraken --db $DBNAME seqs.fa

當(dāng)然，為了加速鑒定速度，kraken也提供了preload指令將數(shù)據(jù)庫(kù)和索引預(yù)先加載到內(nèi)存里去，具體使用的方法是mmap（Box4）

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

Box4：什么是mmap？
mmap是一種內(nèi)存映射文件的方法，即將一個(gè)文件或者其它對(duì)象映射到進(jìn)程的地址空間，實(shí)現(xiàn)文件磁盤地址和進(jìn)程虛擬地址空間中一段虛擬地址的一一對(duì)映關(guān)系。實(shí)現(xiàn)這樣的映射關(guān)系后，進(jìn)程就可以采用指針的方式讀寫操作這一段內(nèi)存，而系統(tǒng)會(huì)自動(dòng)回寫頁(yè)面到對(duì)應(yīng)的文件磁盤上，即完成了對(duì)文件的操作而不必再調(diào)用read,write等系統(tǒng)調(diào)用函數(shù)。相反，內(nèi)核空間對(duì)這段區(qū)域的修改也直接反映用戶空間，從而可以實(shí)現(xiàn)不同進(jìn)程間的文件共享。
常規(guī)文件操作需要從磁盤到頁(yè)緩存再到用戶主存的兩次數(shù)據(jù)拷貝。而mmap操控文件，只需要從磁盤到用戶主存的一次數(shù)據(jù)拷貝過(guò)程。
（Source: https://www.cnblogs.com/huxiao-tee/p/4660352.html）

那么最后我們來(lái)回到最開始的5個(gè)問題：
? 為什么Kraken的分析速度那么快? 使用mmap映射+索引 (C++ &Perl)
? 為什么Kraken的數(shù)據(jù)庫(kù)有幾百G那么大? 原始database大+k-mer大小
? 為什么Kraken建庫(kù)的速度非常慢? 數(shù)據(jù)庫(kù)qsort排序+創(chuàng)建index+遍歷設(shè)置LCA value
? 為什么Kraken數(shù)據(jù)庫(kù)的載入速度非常慢? 數(shù)據(jù)庫(kù)大
? Kraken的數(shù)據(jù)庫(kù)能否拆分使得其能夠分布式運(yùn)行? 原代碼不能，但是可以DIY?

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