Kraken的原理淺析

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

特別說明
Kraken軟件已經出了非常久底循，所以網上有不少相關資料，因此關于Kraken的解讀基本參考網上的資料+文獻本身槐瑞。資料來源：http://jackywu.site/technology/kraken-code-analysis/
（如有侵權刪）

從作者說起

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

可能有人已經知道這位通訊作者了，沒錯是個大佬了悼沿。

Steven Lloyd Salzberg

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

Publications

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

作者Derrick E Wood在2019年發(fā)布了Kraken的新版本Kraken2非驮。（嗯交汤，看了作者的照片，果然是學習改變人生）

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

KrakenUniq與Kraken比較主要是對使用了獨特的k-mer counts（也就是對k-mer進行了優(yōu)化）纵顾，所以在速度和準確性上得到了一定的提升伍茄。不過，這次我們并不會深入講解KrakenUniq施逾，有興趣的小伙伴可以自己去讀相關的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）對應的taxon數(shù)據(jù)庫
? 將數(shù)據(jù)庫和索引文件映射到內存

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

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

舉個例子，如上圖掉伏，如果輸入的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é)點對應的物種。

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)的距離给僵。

建庫過程

在正式進行建庫之前自然是要下載你所需要的微生物序列。然后再使用kraken-build命令進行建庫详拙。
如果你已經進行過kraken-build命令想际，完成了建庫，那么再次輸入該命令的時候溪厘，系統(tǒng)就會提示你步驟已經完成胡本。
建庫成功后，我們會生成下述幾個文件：
?database.kdb: Contains the k-mer to taxon mappings ?database.idx: Contains minimizer offset locations in database.kdb
?taxonomy/nodes.dmp: Taxonomy tree structure + ranks
?taxonomy/names.dmp: Taxonomy names

具體地畸悬，從上圖中我們可以看到建庫過程一共有6步侧甫，其中Step4在目前的版本中已經不需要進行了。
根據(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文件, 文件內容是 “k-mer: count”。

Jllyfish是CBCB（Center for Bioinformatics and Computational Biology）的Guillaume Mar?ais 和 Carl Kingsford 研發(fā)的一款計數(shù) DNA 的 k-mers 的軟件泰偿。該軟件運用 Hash 表來存儲數(shù)據(jù)熄守，同時能多線程運行，速度快，內存消耗小裕照。該軟件只能運行在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：什么是快速排序算法餐屎？
快速排序算法是在起泡排序的基礎上進行改進的一種算法，其實現(xiàn)的基本思想是：通過一次排序將整個無序表分成相互獨立的兩部分玩祟，其中一部分中的數(shù)據(jù)都比另一部分中包含的數(shù)據(jù)的值小腹缩，然后繼續(xù)沿用此方法分別對兩部分進行同樣的操作，直到每一個小部分不可再分空扎，所得到的整個序列就成為了有序序列藏鹊。
例如，對無序表{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}诸衔，繼續(xù)采用此種方法分別對兩個子表進行排序；
前部分子表以 27 為支點颇玷，排序后的子表為{13笨农，27，38}帖渠，此部分已經有序谒亦；后部分子表以 65 為支點，排序后的子表為{49空郊，65份招，97，76}狞甚；
此時前半部分子表中的數(shù)據(jù)已完成排序锁摔；后部分子表繼續(xù)以 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的個數(shù)以及所需要的存儲空間。Kraken默認的是15-bp也殖，但是也可以修改土思。有些人由于設備的限制可能會建立一個MiniKraken，那么在MiniKraken中作者就采用了13-bp忆嗜，以保證最后建立的數(shù)據(jù)庫大小在4GB以內己儒。

通過設置minimizer可以發(fā)現(xiàn)，在鑒定的時候對于一個k-mer的搜索范圍會大大縮小捆毫，從而減少運行時間闪湾。

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

Step6：Set LCA Value
那么我們重點來講一下Step6绩卤⊥狙可以看一下文章中關于設定LCA Value的描述。

其實在這一步就是要將k-mer構建成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都會設置成這個序列所對應的物種taxonomic ID number宏胯。

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

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

舉個粗暴的例子：
比如菌種a1和菌種a2都屬于菌屬A，并且共有一部分序列鹰祸，那么假設序列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進行計算街图，就會變成它們對應的屬A的taxonomic ID number。而這個節(jié)點就會從最底層往上走一個分類水平懒构，變成了屬餐济。
那么如果a1和a2不是共同屬于一個屬A，而是共同屬于科A胆剧，那么對應這個共有序列S的所有k-mer的LCA Value就會變成科A的taxonomic ID number絮姆，這個節(jié)點也就更接近樹的根節(jié)點了。

這就是建庫的全部過程秩霍。

鑒定

鑒定的命令十分簡單篙悯。

kraken --db $DBNAME seqs.fa

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

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

Box4：什么是mmap？
mmap是一種內存映射文件的方法颠悬，即將一個文件或者其它對象映射到進程的地址空間矮燎，實現(xiàn)文件磁盤地址和進程虛擬地址空間中一段虛擬地址的一一對映關系。實現(xiàn)這樣的映射關系后赔癌，進程就可以采用指針的方式讀寫操作這一段內存诞外，而系統(tǒng)會自動回寫頁面到對應的文件磁盤上，即完成了對文件的操作而不必再調用read,write等系統(tǒng)調用函數(shù)灾票。相反浅乔，內核空間對這段區(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+遍歷設置LCA value
? 為什么Kraken數(shù)據(jù)庫的載入速度非常慢? 數(shù)據(jù)庫大
? Kraken的數(shù)據(jù)庫能否拆分使得其能夠分布式運行? 原代碼不能贤壁，但是可以DIY?

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