Kraken軟件已經(jīng)出了非常久，所以網(wǎng)上有不少相關(guān)資料侵佃，因此關(guān)于Kraken的解讀基本參考網(wǎng)上的資料+文獻(xiàn)本身麻昼。資料來源：http://jackywu.site/technology/kraken-code-analysis/
（如有侵權(quán)刪）
從作者說起
這個(gè)實(shí)驗(yàn)室十分鐘情于Genome Biology，Kraken和Kraken2都發(fā)表該雜志上馋辈。Kraken發(fā)表于2014年抚芦，作者和通訊正是上圖所圈選出來的2位。
搜了一下該實(shí)驗(yàn)室之前發(fā)表的論文迈螟，可以看到生信領(lǐng)域常用的軟件Bowtie2叉抡、TopHat、TopHat2答毫、FLASH和HISAT等都出自該實(shí)驗(yàn)室褥民。
當(dāng)然了，正如之前所說烙常，Kraken這個(gè)軟件也是在不斷更新之中。
作者Derrick E Wood在2019年發(fā)布了Kraken的新版本Kraken2。（嗯蚕脏，看了作者的照片侦副，果然是學(xué)習(xí)改變?nèi)松?/p>
那么其實(shí)在2014-2019年之間，該實(shí)驗(yàn)室也有其他人對(duì)Kraken進(jìn)行過迭代更新驼鞭，在2019年發(fā)表了KrakenUniq秦驯，沒錯(cuò)依然是在Genome Biology。
KrakenUniq與Kraken比較主要是對(duì)使用了獨(dú)特的k-mer counts（也就是對(duì)k-mer進(jìn)行了優(yōu)化）挣棕，所以在速度和準(zhǔn)確性上得到了一定的提升译隘。不過，這次我們并不會(huì)深入講解KrakenUniq洛心，有興趣的小伙伴可以自己去讀相關(guān)的paper固耘。
這次講解主要想要回答5個(gè)問題，也就是參考的資料里面提出來的5個(gè)問題：
? 為什么Kraken的分析速度那么快?
? 為什么Kraken的數(shù)據(jù)庫有幾百G那么大?
? 為什么Kraken建庫的速度非常慢?
? 為什么Kraken數(shù)據(jù)庫的載入速度非常慢?
? Kraken的數(shù)據(jù)庫能否拆分使得其能夠分布式運(yùn)行?
The Kraken sequence classification algorithm
首先來看一下Kraken的基本算法词身。簡(jiǎn)單來講厅目，使用Kraken軟件有2步：準(zhǔn)備（建庫）+鑒定。
準(zhǔn)備(建庫)
? 建立k-mer（Box1）對(duì)應(yīng)的taxon數(shù)據(jù)庫
? 將數(shù)據(jù)庫和索引文件映射到內(nèi)存
實(shí)際上建庫的工作只需要在第一次運(yùn)行該軟件時(shí)進(jìn)行即可法严，再次使用的時(shí)候损敷，因?yàn)橐呀?jīng)做好了準(zhǔn)備工作，所以只需要直接對(duì)序列進(jìn)行鑒定即可深啤。
鑒定
? 將待鑒定序列切成k-mer
? 將k-mer比對(duì)到數(shù)據(jù)庫上獲得其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è)樹或者有向無環(huán)圖中同時(shí)擁有x和y作為后代的最深的節(jié)點(diǎn)创倔。
例子:
在下圖Example 1中，x與y的最近公共祖先被標(biāo)記為深綠色焚碌，其他公共祖 先被標(biāo)記為淺綠色畦攘。
Example 1
計(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)的距離摆出。
建庫過程
在正式進(jìn)行建庫之前自然是要下載你所需要的微生物序列朗徊。然后再使用kraken-build命令進(jìn)行建庫。
如果你已經(jīng)進(jìn)行過kraken-build命令偎漫，完成了建庫爷恳，那么再次輸入該命令的時(shí)候，系統(tǒng)就會(huì)提示你步驟已經(jīng)完成象踊。
建庫成功后温亲，我們會(huì)生成下述幾個(gè)文件：
?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īng)不需要進(jìn)行了栈虚。
根據(jù)不同步驟所花費(fèi)的時(shí)間可以發(fā)現(xiàn)，建庫耗時(shí)主要集中在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.
所以可以看到標(biāo)準(zhǔn)的庫下載的是細(xì)菌泌射、古菌以及病毒的RefSeq數(shù)據(jù)粘姜。
但是實(shí)際上我們知道，就人體微生物組數(shù)據(jù)而言熔酷，其實(shí)真菌也是很重要的組成部分孤紧，因此我們可以自主添加真菌數(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就是利用Jellyfish軟件切割k-mer拒秘，生成database.jdb文件, 文件內(nèi)容是 “k-mer: count”号显。
Jellyfish是CBCB（Center for Bioinformatics and Computational Biology）的Guillaume Mar?ais 和 Carl Kingsford 研發(fā)的一款計(jì)數(shù) DNA 的 k-mers 的軟件臭猜。該軟件運(yùn)用 Hash 表來存儲(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ù)庫進(jìn)行了一個(gè)優(yōu)化钾菊，縮減大小。
Step3：Sort set: database.kdb + databse.idx
這就是我們剛才說到的特別慢的一步偎肃。這一步是干啥的呢煞烫？具體可以分為兩步：
Step3.1: 對(duì)database.jdb進(jìn)行排序
Step3.2: 生成索引文件的
由于第一步排序使用了快速排序算法（Box3），因此就特別慢累颂。
Box3：什么是快速排序算法滞详？
快速排序算法是在起泡排序的基礎(chǔ)上進(jìn)行改進(jìn)的一種算法，其實(shí)現(xiàn)的基本思想是：通過一次排序?qū)⒄麄€(gè)無序表分成相互獨(dú)立的兩部分紊馏，其中一部分中的數(shù)據(jù)都比另一部分中包含的數(shù)據(jù)的值小料饥，然后繼續(xù)沿用此方法分別對(duì)兩部分進(jìn)行同樣的操作，直到每一個(gè)小部分不可再分朱监，所得到的整個(gè)序列就成為了有序序列岸啡。
例如，對(duì)無序表{49赫编，38巡蘸，65，97擂送，76悦荒，13，27嘹吨，49}進(jìn)行快速排序搬味，大致過程為：
首先從表中選取一個(gè)記錄的關(guān)鍵字作為分割點(diǎn)（稱為“樞軸”或者支點(diǎn)，一般選擇第一個(gè)關(guān)鍵字）躺苦，例如選取 49身腻；
將表格中大于 49 的放置于 49 的右側(cè)猛蔽，小于 49 的放置于 49 的左側(cè)鸵贬，假設(shè)完成后的無序表為：{27，38绳姨，13愈诚，49她按，65牛隅，97，76酌泰，49}媒佣；
以 49 為支點(diǎn)，將整個(gè)無序表分割成了兩個(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}萍启；
通過以上幾步的排序总珠，最后由子表{13，27勘纯，38}局服、{49}、{49}驳遵、{65}淫奔、{76，97}構(gòu)成有序表：{13堤结，27唆迁，38鸭丛，49，49唐责，65鳞溉，76，97}鼠哥；
那么第二步就是生成索引文件熟菲。具體的我們可以看下面這張圖。
這里Kraken用到了一個(gè)叫做“minimizer”的玩意朴恳。所謂的minimizer就是科盛，由于k-mer的切割方式所以導(dǎo)致臨近的k-mer實(shí)際是很相似的，所以就可以用一個(gè)minimizer來表征一組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ù)庫大小在4GB以內(nèi)乳怎。
通過設(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)來講一下Step6恕出⊙叮可以看一下文章中關(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ù)庫。然后對(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的過程抑进。
舉個(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)了橱野。
這就是建庫的全部過程。
鑒定
鑒定的命令十分簡(jiǎn)單善玫。
kraken --db $DBNAME seqs.fa
當(dāng)然水援，為了加速鑒定速度，kraken也提供了preload指令將數(shù)據(jù)庫和索引預(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)回寫頁面到對(duì)應(yīng)的文件磁盤上，即完成了對(duì)文件的操作而不必再調(diào)用read,write等系統(tǒng)調(diào)用函數(shù)涝开。相反循帐，內(nèi)核空間對(duì)這段區(qū)域的修改也直接反映用戶空間，從而可以實(shí)現(xiàn)不同進(jìn)程間的文件共享舀武。
常規(guī)文件操作需要從磁盤到頁緩存再到用戶主存的兩次數(shù)據(jù)拷貝拄养。而mmap操控文件，只需要從磁盤到用戶主存的一次數(shù)據(jù)拷貝過程。
（Source: https://www.cnblogs.com/huxiao-tee/p/4660352.html）
那么最后我們來回到最開始的5個(gè)問題：
? 為什么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ù)庫能否拆分使得其能夠分布式運(yùn)行? 原代碼不能瘪匿，但是可以DIY?
好了跛梗，今天就講到這里啦～
參考鏈接：http://www.reibang.com/p/eccbc2c616eb
