Quality Control II: Trimming （二）：BBDuk

參考：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給予充分的自由