Kraken軟件可以通過序列對樣本進(jìn)行物種注釋艰垂，kraken2在該軟件基礎(chǔ)之上做了一些更新，其中包括注釋的加速猜憎、支持氨基酸序列的注釋等其他特性娩怎。
軟件安裝胰柑、基本使用以及專有數(shù)據(jù)庫的下載，可以直接查看官網(wǎng)說明柬讨，很詳細(xì)崩瓤。鏈接https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown踩官。
在有部分研究結(jié)果的基礎(chǔ)上，我們可以把自己積累的數(shù)據(jù)加到已有的數(shù)據(jù)庫中蔗牡，用來搭建更適合自己數(shù)據(jù)分析的數(shù)據(jù)庫颖系。在使用自己序列建庫之前辩越，需要先至少安裝一個(gè)已有的kraken2數(shù)據(jù)庫。
下載物種注釋信息
kraken2-build --download-taxonomy --db kraken2_dbtest
下載的文件放在kraken2_dbtest目錄
下載物種數(shù)據(jù)庫
kraken2-build --download-library bacteria --db kraken2_dbtest
向數(shù)據(jù)庫中添加自己的序列
kraken2-build --add-to-library test.genomic.fna --db kraken2_dbtest
添加自己的序列時(shí)黔攒，需要指定taxid等信息趁啸，如果沒有這些信息，會(huì)報(bào)以下錯(cuò)誤
$ kraken2-build --add-to-library test.genomic.fna --db kraken2_dbtest
Use of uninitialized value $taxid in pattern match (m//) at scan_fasta_file.pl line 43, <> line 1.
Use of uninitialized value $taxid in concatenation (.) or string at scan_fasta_file.pl line 47, <> line 1.
Use of uninitialized value $taxid in pattern match (m//) at scan_fasta_file.pl line 43, <> line 5.
Use of uninitialized value $taxid in concatenation (.) or string at scan_fasta_file.pl line 47, <> line 5.
Masking low-complexity regions of new file... done.
Added "test.genomic.fna" to library (test2)
此時(shí)有兩種方法可以解決
- 修改序列名
這種方法最簡單莲绰，直接處理我們提供序列的id姑丑，添加對應(yīng)的taxid信息蛤签，示例如下
>seq1|kraken:taxid|999999991 seq1
···
>seq2|kraken:taxid|999999992 seq2
···
>seq3|kraken:taxid|999999993 seq3
···
>seq4|kraken:taxid|999999994 seq4
···
添加的taxid需要注意，不能與原有數(shù)據(jù)庫中的taxid有重復(fù)震肮。
- 修改已有的taxonomy文件
這種方法不需要修改自己的序列，可以直接將提供序列的信息加入到taxonomy和library兩個(gè)目錄中的注釋文件內(nèi)留拾。
taxonomy目錄下存在多個(gè)注釋文件，如下載的RefSeq數(shù)據(jù)庫中包含以下文件
citations.dmp
delnodes.dmp
division.dmp
gc.prt
gencode.dmp
images.dmp
merged.dmp
names.dmp
nodes.dmp
nucl_gb.accession2taxid
nucl_wgs.accession2taxid
prelim_map.txt
readme.txt
各文件中的內(nèi)容及各列含義痴柔，可以參照目錄下的readme.txt文件沦偎。根據(jù)readme.txt文件的說明咳蔚，將自己的序列信息添加到文件中豪嚎，需要修改的文件名包括names.dmp谈火、nodes.dmp、accession2taxid和nucl_gb.accession2taxid文件糯耍。
# names.dmp
3059595 | Entrophosporales | | scientific name |
91234567891 | seq1 | | scientific name |
91234567892 | seq2 | | scientific name |
91234567893 | seq3 | | scientific name |
91234567894 | seq4 | | scientific name |
# nodes.dmp
3059595 | 214506 | order | | 4 | 1 code compliant |
91234567891 | 1 | no rank | | 8 | 0 |
91234567892 | 1 | no rank | | 8 | 0 |
91234567893 | 1 | no rank | | 8 | 0 |
91234567894 | 1 | no rank | | 8 | 0 |
# nucl_gb.accession2taxid
Z99999 Z99999.1 77601 3399750
NC_099991 NC_099991.1 91234567891 912345678911
NC_099992 NC_099992.1 91234567892 912345678921
NC_099993 NC_099993.1 91234567893 912345678931
NC_099994 NC_099994.1 91234567894 912345678941
library目錄主要修改map文件，文件名prelim_map.txt温技。
# prelim_map.txt
ACCNUM NC_067270.1 NC_067270
ACCNUM NC_099991.1 NC_099991
ACCNUM NC_099992.1 NC_099992
ACCNUM NC_099993.1 NC_099993
ACCNUM NC_099994.1 NC_099994
修改完成后革为，重新構(gòu)建數(shù)據(jù)庫索引即可舵鳞。
$ kraken2-build --add-to-library test.genomic.fna --db kraken2_dbtest
Masking low-complexity regions of new file... done.
Added "test.genomic.fna" to library (kraken2_dbtest)
數(shù)據(jù)庫構(gòu)建完成后，建庫目錄會(huì)生成hash.k2d系任、opts.k2d和taxo.k2d三個(gè)文件恳蹲，其他文件可以刪除俩滥。
添加自己的序列時(shí)嘉蕾，建議直接修改序列名霜旧，會(huì)方便很多儡率。第二種方法對格式有要求，格式不對建庫會(huì)失敗以清。另外，自己序列在指定taxid時(shí)掷倔，不能太大眉孩，否則在運(yùn)行時(shí)會(huì)超限勒葱，返回結(jié)果與建庫序列的taxid不一致浪汪。
物種注釋
數(shù)據(jù)庫構(gòu)建完畢后凛虽，可以使用如下命令進(jìn)行注釋。
# 輸入序列為fasta
kraken2 --db kraken2_dbtest test.fna --report k2_report.txt --use-names >kraken2_anno.output 2>kraken2_anno.error
# 輸入序列為fastq (PE測序)
kraken2 --db kraken2_dbtest test.fq --paired --report k2_report.txt --use-names >kraken2_anno.output 2>kraken2_anno.error
- output文件
kraken2_anno.output文件會(huì)記錄每條序列的注釋結(jié)果凯旋，示例如下
C NZ_BANN01000056.1 Streptococcus parauberis KCTC 11980BP (taxid 1260132) 353 1301:2 91061:7 1301:35 1260132:36 1301:47 1260132:40 91061:3 1260132:34 91061:5 1239:13 91061:2 1260132:39 1301:9 91061:5 1783272:1 91061:2 186826:5 91061:17 1260132:10 0:7
C NZ_QDLQ01000070.1 Salmonella enterica (taxid 28901) 442 28901:408
C NZ_MOXO01000049.1 Ensifer adhaerens (taxid 106592) 36028 0:2824 133448:2 0:22 28901:1 0:5837 200452:12 0:89 106592:10 0:121 286:5 0:1475 1236:5 656:2 0:9 656:1 0:96 286:3 0:321 106592:5 0:15690 1299326:2 0:2779 595497:5 0:2025 106592:5 0:1513 1224:2 106592:5 0:181 2:3 0:67 1619950:3 0:11 169292:1 0:582 106592:3 0:1 106592:2 0:310 106592:1 0:18 1716172:4 0:78 106592:5 0:1858
U NZ_MOXO01000052.1 unclassified (taxid 0) 34220 0:2408 106592:5 0:31773
U NZ_MOXO01000075.1 unclassified (taxid 0) 29561 0:1499 106592:2 0:28026
第一列表明該序列是否被注釋，C表示注釋成功至非，U表示注釋失斈剖稹睡蟋；第二列為待注釋序列的id踏幻，來源于輸入文件戳杀；第三列是物種的taxid，默認(rèn)只輸出taxid信卡，此處物種的名字+taxid是因?yàn)槭褂昧?code>--use-names參數(shù)隔缀；第四列為序列長度傍菇，如果是PE測序的fastq序列，該列為兩個(gè)數(shù)丢习，表示兩條配對reads各自的長度牵触，如150|150咐低；其他各列為各種統(tǒng)計(jì)結(jié)果
- report文件
k2_report.txt文件會(huì)記錄注釋的統(tǒng)計(jì)結(jié)果揽思，示例如下
95.07 58411 58411 U 0 unclassified
4.93 3028 1408 R 1
0.15 91 91 R1 26517 Salmonella enterica
0.14 89 89 R1 12689
0.13 80 80 R1 11033
第一列是序列占比；第二列是比對到數(shù)據(jù)庫目標(biāo)序列的序列條數(shù)(covered to taxid)见擦；第三列為有注釋結(jié)果的序列條數(shù)(assigned to taxid)羹令；第四列是注釋結(jié)類型，包括 (U)nclassified损痰、(R)oot福侈、(D)omain卢未、(K)ingdom、(P)hylum尝丐、(C)lass显拜、(O)rder爹袁、(F)amily矮固、(G)enus和(S)pecies失息；第五列為序列的taxid档址；最后一列為物種名(該信息如果數(shù)據(jù)庫里有可以獲取，沒有的話留空)守伸。
參考文章
[1] https://ccb.jhu.edu/software/kraken2/index.shtml?t=manual
