基于kmer快速搜索的read分类工具
写在前面的
以前在做宏基因组数据分析时经常会有快速鉴定测序数据的reads包含的微生物种类的需求。那时常用到基于kmer快速搜索算法的软件如kraken2。这类基于kmer搜索的算法无需序列比对，因此相较于比对算法的显著优势是速度快。我曾经认为，只有在微生物组领域，可能才会大量用到基于kmer搜索的reads快速分类算法。因为微生物组的参考基因组数据相当大，目前NCBI已测序的包含病毒在内的微生物参考基因组总计已经超过了1T。与微生物组参考序列相比，人类基因组3G的大小似乎很难会用到kmer快速搜索。然而这个孤陋寡闻的认知，在最近被打破。
背景
最近对人类基因组的个性化分析任务中，我们恰好构建了一个由若干长度小于30bp的随机序列组成的文件大小约为30G的参考序列数据。而我们的其中一个分析需求是将公共数据库中大量人源性NGS深度测序的reads快速与这些随机序列比较。如果采用传统的比对方式，在分析速度上会大打折扣。因此我突然意识到基于kmer搜索的算法也可以应用到这个任务中。尽管kraken2软件可以支持自建参考序列的kmer数据库。但是kraken2毕竟是为微生物序列注释而设计，其中在构建kmer数据库时需要输入固定格式的taxonomy分类信息文件。碍于时间限制，我懒得为了迎合kraken2的格式刻意去生成随机序列的taxonomy分类信息文件。于是我开始检索信息，查询是否存在一些现成的工具，可以实现通过kmer搜索的方式基于一类参考序列从另一类序列中划分出能够匹配到这类参考序列和不能匹配的这类参考序列的文件，最后通过查询锁定到bbduk.sh这个脚本。
bbduk.sh：基于kmer搜索的reads修剪及分类工具
bbduk.sh是归属于BBMap的一个工具，其实我之前的推文分享过BBMap工具下的repair.sh脚本，该脚本用于检查双端fastq文件中是否存在不成对的reads并去除这些不成对reads。bbduk则负责用于检查双端reads中是否存在包含接头、引物、污染序列等之类在下游分析中不希望出现的序列，并去除这些序列。其实这个功能本质上是上面我们写到的需求。由于bbduk.sh和repair.sh都属于BBMap，有了上次的repair.sh使用经验，bbduk.sh也遵循类型的命令语法。
bbduk.sh in=sample_R1.fastq in2=sample_R2.fastq ref=kmers_reference.fasta k=24 hdist=0 fixjunk=T outm=sample_matchkmer_R1.fq.gz outm2=sample_matchkmer_R2.fq.gz out=sample_unmatchkmer_R1.fq.gz out2=sample_unmatchkmer_R2.fq.gz -Xmx90g threads=3#-Xmx90g: 设定使用内存为90g，默认为4G#threads：线程数量，注意-Xmx是指程序运行时的总内存 threads参数指定线程后每个线程分配的内存数会被线程数均分，要确保根据读取的参考序列和reads数量保证单独每个线程的内存够用#hdist=0: 设定序列距离，0代表不允许kmer搜索中存在错配#k=24：设定搜索kmer长度为24bp#fixjunk=T：搜索中自动识别并忽视reads或参考序列中的可能存在问题的碱基#outm, outm2：双端模式下match到reference kmer序列数据集的reads集合，因为bbduk.sh设计初衷是用于数据前处理时从fastq文件中去除一些序列，因此默认match的reads属于污染物，是被去除的，unmatch是被保留的，所以outm和outm2的结果需要特别声明才能输出，而我们这里需要的是这部分信息。#程序运行结果Version 39.26Set threads to 3Allocating kmer table: 0.025 seconds.Initial:Memory: max=96636m, total=96636m, free=96552m, used=84mAdded 731499509 kmers; time: 426.895 seconds.Memory: max=96636m, total=96636m, free=8433m, used=88203mInput is being processed as pairedStarted output streams: 0.339 seconds.Processing time: 3.596 seconds.Input: 2065220 reads 105206212 bases.Contaminants: 247926 reads (12.00%) 11834358 bases (11.25%)Total Removed: 247926 reads (12.00%) 11834358 bases (11.25%)Result: 1817294 reads (88.00%) 93371854 bases (88.75%)Time: 431.637 seconds.Reads Processed: 2065k 4.78k reads/secBases Processed: 105m 0.24m bases/sec
等待上面的命令运行完成后，输出结果中则包含了我们划分好的包含了目标序列的reads及不包含目标序列的reads，然后我们可以继续开展下游分析。