Kraken2建库及使用

2024-04-18 本文已影响0人 xiaoji_hb

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

[2] https://github.com/DerrickWood/kraken2