从kraken2注释结果中提取指定物种序列

对二代测序的数据，可以使用一些软件进行物种注释，从而了解样本环境中的微生态状况。常用的物种注释软件包括MetaPhlAn系列、Kraken系列等等，不好意思我就记得这俩，因为平时就用着俩。
MetaPhlAn是在HUMAnN中附带的，而且使用过程中存在两个问题：一个是检出率不够高，除非是丰度或者匹配程度特别高的物种，否则很难检测出来；另一个是数据库更新不够方便，目前的特征序列数据库应该还是201901的版本，所以有些新物种不会得到注释。而Kraken采用哈希表作为检索方式，虽然特异性尚需讨论，但检出率非常高，而且数据库构建相对比较开放，基本上以ncbi的物种数据库作为依托，应该会比较全面，更新也会很及时。我总结其特点就是“建库一两周、载入三分钟、检索十秒钟”，哈希表简直就是用空间换时间的代表。
使用Kraken系列进行检测的过程中有一个问题，就是如果注释出了一个感兴趣的特殊物种，就需要将这些序列从原始数据当中提取出来进一步分析。虽然有些工具能够进行提取，比如samtools之类的，但考虑到平台的通用性（比如我有时会用Mac来处理一部分数据和代码），所以写一段代码来处理这个问题，并与大家交流一下。

文件准备

文件准备需要参照一下我之前的文章宏基因组单个样本数据处理流程当中的部分内容。
**seqs.hostfree.fastq：**去除宿主序列之后的fastq序列文件，其实就是所有的（其实还是含有部分未能去除宿主序列）微生物序列。
**seqs.hostfree.kra：**对seqs.hostfree.fastq序列文件使用kraken2进行注释比对之后生成的output文件，每行是对应的每一条序列信息，总共为5列：第一列只有U和C两个值，标记该序列是Unclassified（未能匹配）还是Classified（已被匹配）；第二列是该序列的测序信息，与seqs.hostfree.fastq文件中的测序信息完全吻合；第三列是该序列对应的物种kraken-id，是Kraken对物种的编码id，需要在建库文件当中寻找一下；第四列是该序列的长度；第五列是该序列中每个k-mers的LCA映射。而本次操作主要需要第二列和第三列。
**seqs.hostfree.mpa：**这是使用Kraken进行注释之后采用MetaPhlAn格式表述的物种分类，这里用处并不大，只是在处理之前观察一下目标物种。

代码运行

接下来的处理代码采用R语言进行编写。

# 读入mpa格式的物种表格。
tax.seqs.count <- read.csv("seqs.hostfree.mpa", sep = "\t", header = F)
# 这个时候就可以自己看看tax.seqs.count里注释了哪些物种以及有没有自己感兴趣的物种。
# 最好能够定位到“种”，因为其他级别我还没试过。

# 读取序列匹配表格。
reads.table <- read.csv("seqs.hostfree.kra", sep = "\t", header = F)
# 保留已经匹配的序列，就是把第一列为“C”的挑出来。
classified.reads <- reads.table[grep("C", reads.rna.table[ , 1], invert = F), ]
# 可以看看里面都有哪些序列，以及不同物种id的大致数量。
summary(as.factor(classified.reads[ , 3]))

# 现在寻找特定物种的序列，比如咱们随便挑一个krakenid为“2697049”的序列吧。
# 读入fastq文件，基本思路是先把每一行读入，再按照所在行数对4的余数进行区分：余1为测序信息，余2为序列本身，余3为符号“+”，余0为质控数据。
# 记得读取时千万要使用comment.char = ""的标记去除读取注释，因为fastq文件当中的质控行可能以字符“#”开头，读取的时候整行会被略过，影响后续内容的读入。
reads.in <- read.table("seqs.hostfree.fastq", sep = "\n", comment.char = "")
# 把数据的四列分离出来，先分出第一列，随后每生成一列就拼接入数据表中。
reads.data <- reads.in[as.numeric(rownames(reads.out))%%4 == 1, ]
reads.data <- cbind(reads.data, reads.in[as.numeric(rownames(reads.out))%%4 == 2, ])
reads.data <- cbind(reads.data, reads.in[as.numeric(rownames(reads.out))%%4 == 3, ])
reads.data <- cbind(reads.data, reads.in[as.numeric(rownames(reads.out))%%4 == 0, ])

# 从reads数据中选出符合注释要求的序列。
# 这段代码的意思就是提取reads.data（就是序列文件）中第一列（就是序列信息）属于classified.reads中第三列（物种id）为指定物种（就是之前随便挑选的“2697049”）的序列信息（就是第二列序列信息）。反正我刚打完疫苗，已经被这段的描述绕晕了……
reads.selected <- reads.data[reads.data[ , 1] %in% paste0('@', classified.reads[classified.reads$V3 == "2697049", 2]), ]

# 将这些序列写入文档。
write.table(reads.selected, file = "seqs.selected.fastq", sep = '\n', row.names = F, col.names = F, quote = F)

现在已经把fastq格式的序列弄出来，那么就用以下这段Linux命令将fastq转换为fasta格式。

$ awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' seqs.selected.fastq >seqs.selected.fasta

物种选择

这里有个新问题：我怎么知道物种的kraken-id呢？其实答案就在kraken2的建库过程当中。当我们使用kraken2-build命令进行检索数据库构建时，会下载ncbi的物种数据库进行构建，在数据库目录下会生成子目录，比如virus、bacteria、archaea、fungi等。每个子目录下又会生成相应物种的序列总和文件，一般文件名是library.fna。在这些目录下面用文本工具打开该文件，搜索感兴趣的物种名称，就能看到对应的krakenid。
当然，也可以使用一段命令来提取这些id信息，我稍后补充。

结果验证

一般严谨起见，对于提取出来的序列需要进行验证。可以直接使用blast官方网站进行验证，或者下载对应物种的序列由比对软件在本地进行验证。虽然开头下机数据特别多，但是一般商业测序也就20×-30×的深度，所以提取出来的序列数据量一般不会特别大，用官方网站进行验证应该足够。
提取出来的序列除了进行验证，也可以做很多事情。比如通过比对，观察所处的基因位置和整体基因组的覆盖情况，这些就看各自的发挥了。