Kraken软件可以通过序列对样本进行物种注释，kraken2在该软件基础之上做了一些更新，其中包括注释的加速、支持氨基酸序列的注释等其他特性。
软件安装、基本使用以及专有数据库的下载，可以直接查看官网说明，很详细。链接https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown。
在有部分研究结果的基础上，我们可以把自己积累的数据加到已有的数据库中，用来搭建更适合自己数据分析的数据库。在使用自己序列建库之前，需要先至少安装一个已有的kraken2数据库。
下载物种注释信息
kraken2-build --download-taxonomy --db kraken2_dbtest
下载的文件放在kraken2_dbtest目录
下载物种数据库
kraken2-build --download-library bacteria --db kraken2_dbtest
向数据库中添加自己的序列
kraken2-build --add-to-library test.genomic.fna --db kraken2_dbtest
添加自己的序列时，需要指定taxid等信息，如果没有这些信息，会报以下错误
$ 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)
此时有两种方法可以解决
- 修改序列名
这种方法最简单，直接处理我们提供序列的id，添加对应的taxid信息，示例如下
>seq1|kraken:taxid|999999991 seq1
···
>seq2|kraken:taxid|999999992 seq2
···
>seq3|kraken:taxid|999999993 seq3
···
>seq4|kraken:taxid|999999994 seq4
···
添加的taxid需要注意，不能与原有数据库中的taxid有重复。
- 修改已有的taxonomy文件
这种方法不需要修改自己的序列，可以直接将提供序列的信息加入到taxonomy和library两个目录中的注释文件内。
taxonomy目录下存在多个注释文件，如下载的RefSeq数据库中包含以下文件
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
各文件中的内容及各列含义，可以参照目录下的readme.txt文件。根据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
修改完成后，重新构建数据库索引即可。
$ 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)
数据库构建完成后，建库目录会生成hash.k2d、opts.k2d和taxo.k2d三个文件，其他文件可以删除。
添加自己的序列时，建议直接修改序列名，会方便很多。第二种方法对格式有要求，格式不对建库会失败。另外，自己序列在指定taxid时，不能太大，否则在运行时会超限，返回结果与建库序列的taxid不一致。
物种注释
数据库构建完毕后，可以使用如下命令进行注释。
# 输入序列为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文件会记录每条序列的注释结果，示例如下
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，默认只输出taxid，此处物种的名字+taxid是因为使用了--use-names参数；第四列为序列长度，如果是PE测序的fastq序列，该列为两个数，表示两条配对reads各自的长度，如150|150；其他各列为各种统计结果
- report文件
k2_report.txt文件会记录注释的统计结果，示例如下
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
第一列是序列占比；第二列是比对到数据库目标序列的序列条数(covered to taxid)；第三列为有注释结果的序列条数(assigned to taxid)；第四列是注释结类型，包括 (U)nclassified、(R)oot、(D)omain、(K)ingdom、(P)hylum、(C)lass、(O)rder、(F)amily、(G)enus和(S)pecies；第五列为序列的taxid；最后一列为物种名(该信息如果数据库里有可以获取，没有的话留空)。
参考文章
[1] https://ccb.jhu.edu/software/kraken2/index.shtml?t=manual
网友评论