参考：BBDuk Guide - Archive

在我们了解了如何使用trimmomatic之后，我们开始进一步了解另外一种trim工具BBDuk

首先小编要声明：如果想要完全掌握一个工具是需要较长时间的钻研和学习的，这里呢只是提供BBDuk处理数据的基本逻辑和基本的一些操作，还有很多高级的操作可以在BBDuk脚本文件查看，或者等到实际运用的时候再进行调用~

对比Trimmomatic

小编认为BBDuk的灵活度还是较高的，和trimmomatic相比的话，BBDuk其实是采取了很多不一样的思维去处理数据，这可以体现在接头去除和质量控制多个环节。

咱们先按照BBDuk的官方文档的介绍顺序对主要的一些参量的使用进行介绍，然后再用实际代码串联运用~

关键参量介绍

Kmer：特定长度K的连续核苷酸片段

Kmer的应用可以说是贯穿了BBDuk始终的，其实Kmer的概念是很好理解的，解释长度为kbp的DNA片段（在此碱基测序语境下，如果在蛋白质会有不一样），那么通常我们也就是从read或者是reference序列里面取子集来时生成Kmer的。

但是如果我们真的想要深入了解Kmer的含义，我们必须得先了解hdist（Hamming Distance）和qhdist（Query Hamming Distance）的区别。

Hamming Distance的定义如下：两个等长字符串之间对应位置上不同字符的数量（豆包）

所以在这个语境下Hamming Distance表示了两个碱基序列的差异程度，其的基本单位可以理解为碱基个数，BBDuk实际比对的原则是（在hdist模式下）根据我们设定的hdist大小来生成reference所有的可能变异序列Kmers，然后再利用滑动窗口来逐个取read片段与我们的Kmer库进行比对，如果精确比对（必须完美匹配）我们认为匹配成功。

我们来细致分析上面的机理：首先比如我们设定Kmer的长度是5（实际是会更长的，一般是15到到31），那么对于referenceATTCGAGT这样一个8碱基的序列，我们可以首先分割出的原始Kmer数就有4个ATTCG,TTCGA,TCGAG,CGAGT，然后我们假设此时的hdist为1，也就是我们允许read片段和reference片段有1个碱基的差异。我们的操作是（对于ATTCG这个Kmer而言），一个一个碱基来看，如果是第一个碱基变异有三种，第二个碱基变异也有3中，那么一共有5个碱基，所以此片段在hdist=1的条件下就有15中可能，在加上其自己就一共有16个，再乘上原始片段数4，所以reference序列在Kmer长度为5，hdist=1的条件下有64个Kmer。那么如果我们的read序列是ATGCGTGAT9个碱基，那么从5'(A那一端）开始就可以依次截取5个Kmer，我们截取第一个Kmer:ATGCG，然后把此序列依次和这64个Kmer对比，如果有序列一致的，说明就匹配成功。对比完之后我们接着对下一个片段对比（TGCGT）

这里有必要说一下为什么我们不选择直接拿reference原始的Kmer去和read作对比，因为这一步骤耗时比较长，但是生成和精确判断这一步骤耗时会较短（具体时间复杂度是多少小编也没有研究）

然后我们把条件调回更加符合普遍的情况：如果read有Mbp，reference有Nbp，设定Kmer长为Kbp，hdist简记为h，那么我们来计算一下内存占用和时间花费（如果生成一个Kmer记为T1，精确比对Kmer记为T2）

内存占用其实是和总共生成了多少Kmer是成正比的（这里的比例系数我们就不写了）我们就直接用Kmer的数量来表示内存占用

一共有（N-K+1）个原始Kmer，每个Kmer挑出h个位点进行变异，一共可以有 $\binom{K}{h}$ 种组合方法，每个位点可以有3种变异方式，那么一个Kmer就有 $3^{h}$ 个变异，所以reference共产生了（N-K+1）* $\binom{K}{h}$ * $3^{h}$ +（N-K+1）个Kmer，然后我们计算时间就是用前面的【（N-K+1）* $\binom{K}{h}$ * $3^{h}$ +（N-K+1）】*T1+（M-K+1）*【（N-K+1）* $\binom{K}{h}$ * $3^{h}$ +（N-K+1）】*T2

其实这个数量级的增长是很快的，读者可以自己计算一下就可以知道

上面是hdist的内存占用和时间花费，如果我们需要去减少内存占用的话我们可以尝试另外一种比对方法：qhdist

这种方法是这么操作的：我们在read当中截取Kmer然后根据qhdist（其实也就是Hamming Distance，只是在这个模式里这么写而已）来生成read的变异Kmer，然后将read的变异Kmer与原始的referenceKmer进行比对，同样我们举一个例子（我们就还是借用hdist的例子吧）

referenceATTCGAGT，readATGCGTGAT，Kmer=5，qhdist=1，首先我们生成reference原始Kmer（也就是只截取不变异）：ATTCG,TTCGA,TCGAG,CGAGT，然后我们先截取read的第一个KmerATGCG，我们对read的这个Kmer进行变异，一共是由15中变异方式，那么我们用一共16个（含不变的本身）Kmer去和reference的原始Kmer对比，如果精确匹配那么说明就匹配成功。这一个read的Kmer结束之后接着对下一个Kmer进行截取变异。

那现在我们来分析为什么这个可以节省内存，但却要付出时间的代价。

同样我们设定：如果read有Mbp，reference有Nbp，设定Kmer长为Kbp，qhdist简记为h，那么我们来计算一下内存占用和时间花费（如果生成一个Kmer记为T1，精确比对Kmer记为T2）

那么我们每一次比较一个read的Kmer需要生成多少变异？（暂且不计算reference的Kmer，因为是线性增长，所以当做小量）： $\binom{K}{h}$ * $3^{h}$ +1，那么我们每次其实仅仅用到这么多内存，当我们再次比较的时候我们就把这些read的Kmer全部清除再次生成即可，所以对比之前的内存（reference所有的Kmer一直储存到程序结束）和现在的储存（当前的read的一个Kmer的变异，唯一一直储存的仅仅是reference的原始Kmer），储存的优势高下立判

然后我们再来分析时间的花费为什么第二种更多。我们来分析生成的Kmer总共要耗时：

【 $\binom{K}{h}$ * $3^{h}$ +1】*（M-K+1）*T1，比较耗时：【 $\binom{K}{h}$ * $3^{h}$ +1】*（M-K+1）*（N-K+1）*T2

所以看到其实在比较耗时上两者是没有差别的，但是在分析耗时上两者一个是与M-K+1成正比一个是与N-K+1乘正比，对于庞大的DNA数据而言，M通常是远远大于N的，所以在内存上qhdist与hdist开销之比接近M-K+1/(N-K+1)

上面用到了一点点数学推导，但是总的结论就是如果我们追求速度当然可以使用hdist，但是如果我们想要减小内存消耗我们使用qhdist

除此之外，小编也相信，通过对hdist和qhdist工作机制的了解，我们肯定对Kmer究竟是何物也有一定的认识了，我们接着介绍其他的参数

题外话：如果选用Kmer长度大于31如何操作

我们一般设定Kmer的长度是15-31bp之间，如果Kmer长度设定超过31比如40，那么判断的方式会有一点变化，可以认为reference生成的Kmer仍然是31bp的，但是read截取的Kmer是40bp的，如果要求匹配，实际上需要使得40bp的Kmer在reference的Kmer当中有10个匹配的31bp的Kmer，和31bp的其实很像，但是其实要求更高了因为这是直接用40bp作为一个单位，但凡其中有一个31bp没有匹配上那么就失败了，感觉乍一看似乎hdist失去效用了其实不然，我们可以简单推导一下，假设40bp的readKmer和40bp的referenceKmer，如果hdist=2，这40bp的referenceKmer说明只可以有2个位点变异，但是如果我们差分为10个31bp的Kmer，似乎可以有20个位点变异，但是你仔细想一下如果在40bp的非起始和末尾端，每个中间碱基至少被31bp重复两次，也就是如果31bp中除了第一个31bp的5'和最后一个31bp的3'变异不会影响别的31bp，夹在其中的有任何一个碱基变异立马要求与之有相同位点的Kmer必须在这个点变异，所以其实10个31bp总的变异次数也不会超过2次，最根本原因就是40<2*31(官方文档仅仅注解了40bp和31bp操作方式很接近，为了教程的完整性我在这里说下我对这个操作的理解，不对的话请评论纠正~）

Mink：序列端降低Kmer长度匹配接头

我们在进行质量处理当中有很大一个任务就是处理接头残存，但是一般来说接头残存一般可能是8bp等比较小的片段，我们一般设定Kmer的长度是15-31bp之间，所以我们是无法直接完美匹配这种仅有部分残留的序列的，为此BBduk引入了另外一种参数mink，mink表示最小的Kmer长度

比如我们设定Kmer长度是25bp，mink=10bp，那么在判断端碱基的时候，如果没有匹配成功，Kmer长度就会从25递减至10bp，如果10bp也没有成功，那么说明匹配失败。首先我们需要去确定是从哪端开始读取，ktrim=r即设定是从右端也就是3'端开始读取（反之ktrim=l），这里的递减的含义是：比如最开始进行比对的是末端25bp的Kmer，如果比对不成功，那么就会截取末端24bp的Kmer，同时reference的Kmer可以在原有25bp中截取24bp（但是截取是仅仅在两端，简单例子就是ATGCG，现在Mink=3，那么新的reference就是ATG，GCG，中间的不会截取），但是这里不会继续占用内存，产生mink比较完即删。总的来说设置mink的作用就是为了避免接头碱基残留过少所以所以降低kmer长度以提高灵敏度，但是请注意mink不可以调整得过低，不然很有可能是假阳性。

hdist2，qhdist2：与Mink搭配使用调节准确度

咱们来做这样一个数学题：如果Kmer长度是25，那么随机匹配成功的概率是多少？答案就是 $1/4^{25}$ ，如果允许hdist=3，那么随机匹配的概率是多少呢？（64* $\binom{25}{3}$ / $4^{25}$ ）但是由于分子很大，所以随机匹配的概率仍然很低，这里随机匹配的概率可以认为是假阳性的概率。

那么如果按照我们刚才的操作，如果我们使用了mink，那么在原本的maxk（即为Kmer原始长度）没有匹配成功的时候，会递减降低长度，如果我们降低到12bp的时候仍然保持hdist=3会发生什么？这里就按照我们刚才的计算方式可以得出 $\binom{12}{3}*64/4^{12}$ ，其实此时的假阳性概率相较于之前是6,730,923,356,124倍，那么为了保证我们仍然有较好的准确度，我们可以再设置一个hdist记为hdist2（同理qhdist2），那么在全长的时候用hdist，在开始递减之后用hdist2.根据我们刚才的分析也可以知道我们使用的hdist2应当是小于hdist的。

当然了这其中的设置也有其粗糙的一面，比如在全长25bp的时候的准确度其实是不如24bp的准确度的，原因就在于hdist是离散分布的，当降低hdist就会对假阳性概率出现很大变动，所以可能会出现极值。但是也不用慌张，就如我们刚才分析，hdist2更多适配的是更小的Kmer，一方面如果提高准确度的情况下仍然能够检测出来当然最好，另一方面即使出现一段较高水平的准确度（甚至高出我们的要求），但是会因为Kmer长度的减小逐渐回归预期的（可以自己计算一下）

（这里再强调一遍Kmer适度选择的重要性：Kmer选择的越大，那么specificity准确度就越高，如果Kmer选择的越小那么sensity灵敏度越高，至于究竟选择如何的Kmer长度需要根据实际的运用需求和一定的试验，但是值得注意的是不管Kmer如何选择，一定要保证Kmer的长度是短于reference的长度的（如果Kmer长度选择大于31其实前面也提及了，可以满足更大的准确度要求））

Kmer Skipping：准确度与运行速度、内存占用的tradeoff

我们已经基本介绍完BBDuk核心的处理数据的思维，也了解了不同算法之间的优劣，但是有一个小小问题，如果我们的dataset真的太大怎么办呢，比如人类基因组有3GB左右的数据量，我们如何快速处理这么多数据，如果我们需要去追求一定的处理速度和尽可能减少内存占用的话，我们可以使用另外一个参数Kmer Skipping

Kmer Skipping有三种模式：rskip（reference skip），qskip（query skip），speed（前两种结合）

首先是rskip，比如我们设定rskip=4，那么就说明我们仅仅会产生1/4个reference序列，即每四个Kmer记录一个Kmer，那么这个方式主要可以减少内存的占用，所以一般和hdist模式搭配使用（当然了生成的Kmer少了，运行速度也自然更快）

其次是qskip，比如我们设定qskip=4，那么说明我们对于每四个read的Kmer序列仅仅会比对一个，那么这个运行方式对于内存的占用一定程度上较少（但是影响不大，因为本来基数不大），更多影响的是生成这些Kmer的时间总和，所以一般是使用在qhist模式当中

最后是speed，speed可以取0-15之间的任意整数，表示会同时跳过speed/16部分的Kmer，比如我们设定speed=10，那么如果在hdist模式当中，我们就在read序列16个Kmer我们仅仅会选择6个出来比对，在reference也仅仅会产生原本数量的6/16个Kmer，如果在qhdist模式当中同理。那么speed参数可以说是同时减小了内存占用（主要体现在reference上）和加快了运行速度（体现在Kmer生成数减少和比对次数减少）

但是有得就有失，追求速度和内存必然会付出一定的准确度代价，有些Kmer没有发生比对，所以有一定的概率产生假阴性，我们选取的skip值越大，那么产生假阴性的概率也就会越大

Numbers of Kmer Hits

我们知道增大准确度的很重要的方式就是增大Kmer长度和降低Hamming Distance，但是除此之外还有别的方式可以提高准确度，就是对一个read被匹配的次数等进行限制。这里介绍三种提高准确度的方式：首先是 “mkh” (minkmerhits)，最少匹配次数：这要求一个read当中至少有mkh个Kmer可以正确匹配，如果mkh=2，这里需要注意的是，不是read的一个Kmer片段需要再reference中找到两个Kmer匹配，而是一个read至少有2个Kmer可以正确匹配。比如说read长10bp，Kmer长5bp，那么read一共可以产生6个Kmer，mkh=2要求这六个Kmer中至少有两个可以匹配

其次是“mkf” (minkmerfraction)，最少匹配比率，同样是是上面这个例子，如果我们设定mkf=0.5，就要求有一半的Kmer可以匹配，那么也就是6个中有3个可以可以匹配

最后是“mcf” (mincoveredfraction)，最少覆盖比率，这个相对上面两个有些不同，如果我们上面分别是一头一尾Kmer实现了匹配，那么就说明read当中的所有碱基实际上都被成功匹配的Kmer覆盖了，所以此时的mcf=1，那么如果是第一和第二个Kmer实现匹配呢？那么mcf=60%，但是mkf=33.3%

我们接着来总结一下这三个参数有何使用上的不同：首先mkh是绝对数量，所以要求mkh是要根据read的长度（当然还有准确度的要求）进行调配的，那么如果mkh越大自然也就说明在当前的read准确度的要求也越大。其次是mkf和mcf的对比，虽然这两者全部都是比率，但是也有很大的不同。读者不妨试想一下，如果read和reference是同一序列（均为80bp），那么我们现在要求hdist=0（其实仍然是可以匹配的），Kmer长度为30，如果read第40个碱基处出现了读取错误，但是reference不变，我们来计算一下此时的mkf=21/51（约等于0.41），此时的mcf=0.99，差距很明显，但是在我们得出结论之前，我们再来看如果第一个碱基出错，mkf=50/51（约等于0.98），此时mcf=0.99，那么这两个例子首先告诉我们一点mkf是和碱基的位置高度相关的，在这一点上mcf会更加稳定

但是如果我们允许Hamming Distance=1，对于mcf而言，碱基出错在端碱基和中间碱基有什么不同（这一次我们不再使用read和reference只差一个碱基，但是长度仍然借用），我们采用半定量的方式进行考虑，首先对于端碱基而言，有且只有一个Kmer是覆盖到它的，但是对于中间碱基而言是有30个Kmer覆盖的，那么也就是说只有在第一个Kmer当中出现了其他任何一个碱基也同时出错，那么第一个碱基就没有被覆盖，但是对于中间的碱基而言，它有30次机会，只要这30个kmer当中存在一个Kmer使得它是唯一出错的碱基，那么就可以被覆盖。所以可以发现其实mcf在端碱基处是更加敏感的，这也正好可以和mkf实现互补。

Maskmiddle 和 rcomp

我们刚才提到如何提高准确度，但有些时候我们想要提高灵敏度，这个时候我们需要降低一些准确度，除了我们刚才提到的可以调节的参量，我们还可以调节这两个参量：Maskmiddle和rcomp

首先需要声明的事，这两个量类似于开关，只有开和关两种操作，而不是可以赋值的变量

Maskmiddle是是否忽略中间碱基是否配对，那么很自然地如果是偶数如Kmer长度是10，那么就忽略第五第六个碱基（Maskmiddle=2），如果是奇数如11，那么就忽略第六个碱基（Maskmiddle=1）。（那么实际上可以理解为：将中间的碱基替换为N，那么任何一个碱基都可完成配对）。但是读者不妨想想为什么我们把中间的碱基降低准确度来换取灵敏度，为什么不替换其他位置的碱基比如端碱基？

其实在上面提到过，处于中间位置的碱基对于Kmer是否可以匹配影响是大于端碱基的，所以我们牺牲一个单位的（或者两个）碱基是否可以正常配对来换来最大的灵敏度自然是最好的。如果我们把端位碱基替换N，那么能够影响到的也即是一个Kmer而已，“利益”没有最大化。

除此之外，一般端位碱基的保守性较好，生物学意义也更大，从这一点出发也应该替换中间碱基。

rcomp是是否匹配反向序列，这一点其实在Illumina测序平台体现的不明显：因为Illumina测序结果一般是不会出现read其中插入实际序列的方向序列，但是在其他的平台比如nanopole等可能就会出现这种特殊情形，当我们把rcomp模式打开，那么不仅会生成reference的正向序列，同时也有反向互补序列，每一个read 的Kmer都会在这两部分序列里面扫描，有任一匹配成功即成功

实际代码操作

在我们基本介绍完BBDuk的基本思路和控制准确度和灵敏度的基本思想之后，我们来进行代码学习

首先安装bbmap（BBDuk是其中的一个工具）

conda install bbmap -y

首先输入有两种方式，一种是单端的一种是双端的（和Trimmomatic类似）

bbduk.sh in=input.fastq.gz out=output.fastq.gz #单端
bbduk.sh in1=read1.fq.gz in2=read2.fq.gz out1=clean1.fq.gz out2=clean2.fq.gz #双端

（输入输出文件可压缩，但是注意拓展名）

输出文件有3种方式（针对于有筛选条件的）：对于单端测序第一种是out，表示通过了筛选的条件，第二种是outm表示没有通过筛选条件的，那么对于双端测序而言，out表示表示两个都通过了筛选，outm在默认情况下表示至少有其一没有通过筛选（即使通过筛选的一部分也放在里面）（也可以调节outm的参数rieb=f，这样就可以使得只有在两个全部没有通过筛选的才会放在outm里面（但是这个不会影响out）），outs一般是用于双端情况下的，表示有一个通过但是另外一个没有通过的时候，那个通过了的序列会存在outs（同样rieb参数也不会影响outs）

下面是豆包生成的教学代码，小编会给予进一步解释

#!/bin/bash# 假设BBDuk已安装，且所在路径已加入环境变量
# 若未加入环境变量，需使用完整路径，如：/path/to/bbmap/bbduk.sh# 示例1：去除测序数据中的接头序列
# 输入：双端测序文件 (read1.fq.gz, read2.fq.gz)
# 输出：去接头后的文件 (clean1.fq.gz, clean2.fq.gz)
bbduk.sh \in1=read1.fq.gz \in2=read2.fq.gz \out1=clean1.fq.gz \out2=clean2.fq.gz \ref=adapters.fa \  # 接头序列数据库（BBMap自带）ktrim=r \          # 从右端（3'端）修剪接头k=23 \             # k-mer长度为23mink=11 \          # 最小k-mer长度为11hdist=1 \          # 允许1个错配tbo \              # 对于双端数据，同时修剪成对的readstpe \              # 对于双端数据，确保成对reads长度一致# 示例2：过滤低质量序列和短序列
# 输入：去接头后的文件
# 输出：过滤后的高质量文件 (filtered1.fq.gz, filtered2.fq.gz)
bbduk.sh \in1=clean1.fq.gz \in2=clean2.fq.gz \out1=filtered1.fq.gz \out2=filtered2.fq.gz \qtrim=rl \         # 从两端（5'和3'）修剪低质量碱基trimq=20 \         # 质量阈值为20（Phred评分）minlen=50 \        # 保留长度≥50bp的序列maq=20 \           # 平均质量≥20的序列才保留# 示例3：去除序列中的污染（如宿主基因组序列）
# 输入：高质量序列文件
# 输出：去污染后的文件 (decontam1.fq.gz, decontam2.fq.gz)
bbduk.sh \in1=filtered1.fq.gz \in2=filtered2.fq.gz \out1=decontam1.fq.gz \out2=decontam2.fq.gz \outm1=contam1.fq.gz \  # 输出被移除的污染序列（可选）outm2=contam2.fq.gz \ref=host_genome.fa \   # 宿主基因组参考序列k=31 \                 # 使用31-mer进行比对hdist=1 \              # 允许1个错配maskmiddle=5 \         # 掩蔽k-mer中间5个碱基（减少对中间区域的依赖）filter# 示例4：单端序列处理（仅处理单个文件）#此示例不再讲解
bbduk.sh \in=single_end_reads.fq.gz \out=clean_single_end.fq.gz \ref=adapters.fa \ktrim=r \k=23 \qtrim=r \trimq=15 \minlen=30

1. 双端测序去接头

# 示例1：去除测序数据中的接头序列
# 输入：双端测序文件 (read1.fq.gz, read2.fq.gz)
# 输出：去接头后的文件 (clean1.fq.gz, clean2.fq.gz)
bbduk.sh \in1=read1.fq.gz \in2=read2.fq.gz \out1=clean1.fq.gz \out2=clean2.fq.gz \ref=adapters.fa \  # 接头序列数据库（BBMap自带）ktrim=r \          # 从右端（3'端）修剪接头k=23 \             # k-mer长度为23mink=11 \          # 最小k-mer长度为11hdist=1 \          # 允许1个错配tbo \              # 对于双端数据，同时修剪成对的readstpe \              # 对于双端数据，确保成对reads长度一致

ktrim之前提到过是规定从哪一边开始检查接头，ktrim=r是从3'，ktrim=l是从5'，ktrim=rl是同时

ref=adapter.fa这个是自带文件，如果我们想要自定义接头序列替换即可（注意地址写清楚）

k=23是指正常（full-length相对于mink）kmer的长度是23bp

mink，hdist不再解释

tbo是对于双端数据，由于有效DNA片段部分互补，所以在某一段检测到接头，自动在另一pair剪除相应部分接头

tpe是将双端数据的长度剪切一直，保证互补配对

注意这里只是对接头进行剪除，没有涉及过滤某个序列，所以简单用out就好（认为全部满足筛选条件）

（这里其实开启去接头模式的关键参数是ktrim）

2. 双端测序过滤

bbduk.sh \in1=clean1.fq.gz \in2=clean2.fq.gz \out1=filtered1.fq.gz \out2=filtered2.fq.gz \qtrim=rl \         # 从两端（5'和3'）修剪低质量碱基trimq=20 \         # 质量阈值为20（Phred评分）minlen=50 \        # 保留长度≥50bp的序列maq=20 \           # 平均质量≥20的序列才保留filter

trimq=20是剪切低质量碱基的参数，将质量低于20的碱基定义为低质量碱基，在qtrim的方向控制下剪切低质量碱基（直到遇到高质量碱基）

minlen是表示对碱基进行修剪之后，仅仅保留长度大于50的序列

maq=20表示对碱基进行修剪之后，仅仅保留平均质量大于20的序列

那么这里就有筛选了，所以out就表示通过筛选的序列，但是后面的分析没有用到没有通过的序列，所以没有用outm和outs来储存这些序列

# 示例3：去除序列中的污染（如宿主基因组序列）
# 输入：高质量序列文件
# 输出：去污染后的文件 (decontam1.fq.gz, decontam2.fq.gz)
bbduk.sh \in1=filtered1.fq.gz \in2=filtered2.fq.gz \out1=decontam1.fq.gz \out2=decontam2.fq.gz \outm1=contam1.fq.gz \  # 输出被移除的污染序列（可选）outm2=contam2.fq.gz \ref=host_genome.fa \   # 宿主基因组参考序列k=31 \                 # 使用31-mer进行比对hdist=1 \              # 允许1个错配maskmiddle=5 \         # 掩蔽k-mer中间5个碱基（减少对中间区域的依赖）filter #开启过滤模式

这里开启了kmer-filtering模式，如果有kmer匹配那么就过滤该read，ref记为reference序列

这里是我们经常使用的代码，需要再强调一下，我们看到这几个代码其实是非常相像的，所以每个代码执行的功能是需要一定的关键参数控制的，比如开启去接头就是参数ktrim控制，开启长度和质量过滤就是minlen和maq控制，开启kmer-filtering就是filter参数控制（如果没有filter那么就不会进行过滤）

一般我们可以将BBDuk操作分成多个部分来执行，就如上面所示，但是请注意一点，操作的先后会严重影响到结果，一般来说我们先进行接头处理（因为接头越完整越好），然后进行一系列过滤操作（如kmer-filtering），最后才会在minlen和map进行操作，当然了只是建议，BBDuk给予充分的自由