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

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

從作者說起

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

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

Steven Lloyd Salzberg

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

Publications

當然了，正如之前所說，Kraken這個軟件也是在不斷更新之中。

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

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

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

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

The Kraken sequence classification algorithm

首先來看一下Kraken的基本算法。簡單來講，使用Kraken軟件有2步：準備（建庫）+鑒定。

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

實際上建庫的工作只需要在第一次運行該軟件時進行即可，再次使用的時候，因為已經(jīng)做好了準備工作，所以只需要直接對序列進行鑒定即可。

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

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

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

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

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

建庫過程

具體地，從上圖中我們可以看到建庫過程一共有6步，其中Step4在目前的版本中已經(jīng)不需要進行了。
根據(jù)不同步驟所花費的時間可以發(fā)現(xiàn)，建庫耗時主要集中在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數(shù)據(jù)。
但是實際上我們知道，就人體微生物組數(shù)據(jù)而言，其實真菌也是很重要的組成部分，因此我們可以自主添加真菌數(shù)據(jù)庫。

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ā)的一款計數(shù) DNA 的 k-mers 的軟件。該軟件運用 Hash 表來存儲數(shù)據(jù)，同時能多線程運行，速度快，內(nèi)存消耗小。該軟件只能運行在64位的Linux系統(tǒng)下。其文章于2011年發(fā)表在雜志 Bioinformatics上。

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

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

Box3：什么是快速排序算法？
快速排序算法是在起泡排序的基礎(chǔ)上進行改進的一種算法，其實現(xiàn)的基本思想是：通過一次排序?qū)⒄麄€無序表分成相互獨立的兩部分，其中一部分中的數(shù)據(jù)都比另一部分中包含的數(shù)據(jù)的值小，然后繼續(xù)沿用此方法分別對兩部分進行同樣的操作，直到每一個小部分不可再分，所得到的整個序列就成為了有序序列。
例如，對無序表{49，38，65，97，76，13，27，49}進行快速排序，大致過程為：
首先從表中選取一個記錄的關(guān)鍵字作為分割點（稱為“樞軸”或者支點，一般選擇第一個關(guān)鍵字），例如選取 49；
將表格中大于 49 的放置于 49 的右側(cè)，小于 49 的放置于 49 的左側(cè)，假設(shè)完成后的無序表為：{27，38，13，49，65，97，76，49}；
以 49 為支點，將整個無序表分割成了兩個部分，分別為{27，38，13}和{65，97，76，49}，繼續(xù)采用此種方法分別對兩個子表進行排序；
前部分子表以 27 為支點，排序后的子表為{13，27，38}，此部分已經(jīng)有序；后部分子表以 65 為支點，排序后的子表為{49，65，97，76}；
此時前半部分子表中的數(shù)據(jù)已完成排序；后部分子表繼續(xù)以 65為支點，將其分割為{49}和{97，76}，前者不需排序，后者排序后的結(jié)果為{76，97}；
通過以上幾步的排序，最后由子表{13，27，38}、{49}、{49}、{65}、{76，97}構(gòu)成有序表：{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的個數(shù)以及所需要的存儲空間。Kraken默認的是15-bp，但是也可以修改。有些人由于設(shè)備的限制可能會建立一個MiniKraken，那么在MiniKraken中作者就采用了13-bp，以保證最后建立的數(shù)據(jù)庫大小在4GB以內(nèi)。

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

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

Step6：Set LCA Value
那么我們重點來講一下Step6。可以看一下文章中關(guān)于設(shè)定LCA Value的描述。

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

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

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

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

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

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

這就是建庫的全部過程。

鑒定

鑒定的命令十分簡單。

kraken --db $DBNAME seqs.fa

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

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

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

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

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