tags： Trimmomatic NGS fastq

NGS 原始數(shù)據(jù)過濾對后續(xù)分析至關(guān)重要，去除一些無用的序列也可以提高后續(xù)分析的準確率和效率。Trimmomatic 是一個功能強大的數(shù)據(jù)過濾軟件。

Trimmomatic 介紹

Trimmomatic 發(fā)表的文章至今已被引用了 2810 次，是一個廣受歡迎的 Illumina 平臺數(shù)據(jù)過濾工具。其他平臺的數(shù)據(jù)例如 Iron torrent ，PGM 測序數(shù)據(jù)可以用 fastx_toolkit 、NGSQC toolkit 來過濾。

Trimmomatic 支持多線程，處理數(shù)據(jù)速度快，主要用來去除 Illumina 平臺的 Fastq 序列中的接頭，并根據(jù)堿基質(zhì)量值對 Fastq 進行修剪。軟件有兩種過濾模式，分別對應 SE 和 PE 測序數(shù)據(jù)，同時支持 gzip 和 bzip2 壓縮文件。

另外也支持 phred-33 和 phred-64 格式互相轉(zhuǎn)化，現(xiàn)在之所以會出現(xiàn) phred-33 和 phred-64 格式的困惑，都是 Illumina 公司的鍋（damn you, Illumina!），不過現(xiàn)在絕大部分 Illumina 平臺的產(chǎn)出數(shù)據(jù)也都轉(zhuǎn)為使用 phred-33 格式了。

Trimmomatic 過濾步驟

Trimmomatic 過濾數(shù)據(jù)的步驟與命令行中過濾參數(shù)的順序有關(guān)，通常的過濾步驟如下：

1. ILLUMINACLIP: 過濾 reads 中的 Illumina 測序接頭和引物序列，并決定是否去除反向互補的 R1/R2 中的 R2。
2. SLIDINGWINDOW: 從 reads 的 5' 端開始，進行滑窗質(zhì)量過濾，切掉堿基質(zhì)量平均值低于閾值的滑窗。
3. MAXINFO: 一個自動調(diào)整的過濾選項，在保證 reads 長度的情況下盡量降低測序錯誤率，最大化 reads 的使用價值。
4. LEADING: 從 reads 的開頭切除質(zhì)量值低于閾值的堿基。
5. TRAILING: 從 reads 的末尾開始切除質(zhì)量值低于閾值的堿基。
6. CROP: 從 reads 的末尾切掉部分堿基使得 reads 達到指定長度。
7. HEADCROP: 從 reads 的開頭切掉指定數(shù)量的堿基。
8. MINLEN: 如果經(jīng)過剪切后 reads 的長度低于閾值則丟棄這條 reads。
9. AVGQUAL: 如果 reads 的平均堿基質(zhì)量值低于閾值則丟棄這條 reads。
10. TOPHRED33: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-33。
11. TOPHRED64: 將 reads 的堿基質(zhì)量值體系轉(zhuǎn)為 phred-64。

Trimmomatic 簡單用法

由于 Trimmomatic 過濾數(shù)據(jù)的步驟與命令行中過濾參數(shù)的順序有關(guān)，因此，如果需要去接頭，建議第一步就去接頭，否則接頭序列被其他的過濾參數(shù)剪切掉部分之后就更難匹配更難去除干凈了。

單末端測序模式

在 SE 模式下，只有一個輸入文件和一個過濾之后的輸出文件：

java -jar <path to trimmomatic jar> SE [-threads <threads>] [-phred33 | -phred64] [-trimlog
<logFile>] <input> <output> <step 1> <step 2> ...

-trimlog 參數(shù)指定了過濾日志文件名，日志中包含以下四列內(nèi)容:

• read ID
• 過濾之后剩余序列長度
• 過濾之后的序列起始堿基位置（序列開頭處被切掉了多少個堿基）
• 過濾之后的序列末端堿基位置
• 序列末端處被剪切掉的堿基數(shù)

由于生成的 trimlog 文件中包含了每一條 reads 的處理記錄，因此文件體積巨大（GB 級別），如果后面不會用到 trim 日志，建議不要使用這個參數(shù)。

雙末端測序模式

在 PE 模式下，有兩個輸入文件，正向測序序列和反向測序序列，但是過濾之后輸出文件有四個，過濾之后雙端序列都保留的就是 paired，反之如果其中一端序列過濾之后被丟棄了另一端序列保留下來了就是 unpaired 。

PE 模式下的輸入輸出文件

java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog
<logFile>] >] [-basein <inputBase> | <input 1> <input 2>] [-baseout <outputBase> |
<paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> <step 2> ...

其中 -phred33 和 -phred64 參數(shù)指定 fastq 的質(zhì)量值編碼格式，如果不設(shè)置這個參數(shù)，軟件會自動判斷輸入文件是哪種格式（v0.32 之后的版本都支持），雖然軟件默認的參數(shù)是 phred64，如果不確定序列是哪種質(zhì)量編碼格式，可以不設(shè)置這個參數(shù)。

輸入輸出文件

PE 模式的兩個輸入文件：sample_R1.fastq sample_R2.fastq以及四個輸出文件：sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq

通常 PE 測序的兩個文件，R1 和 R2 的文件名是類似的，因此可以使用 -basein 參數(shù)指定其中 R1 文件名即可，軟件會推測出 R2 的文件名，但是這個功能實測并不好用，因為軟件只能自動識別推測三種種格式的 -basein:

• Sample_Name_R1_001.fq.gz -> Sample_Name_R2_001.fq.gz
• Sample_Name.f.fastq -> Sample_Name.r.fastq
• Sample_Name.1.sequence.txt -> Sample_Name.2.sequence.txt

建議不用 -basein 參數(shù)，直接指定兩個文件名（R1 和 R2）作為輸入。

輸出文件有四個，當然也可以像上文一樣指定四個文件名，但是參數(shù)太長有點麻煩，有個省心的方法，使用 -baseout 參數(shù)指定輸出文件的 basename，軟件會自動為四個輸出文件命名。例如 -baseout mySampleFiltered.fq.gz ，文件名中添加 .gz 后綴，軟件會自動將輸出結(jié)果進行 gzip 壓縮。輸出的四個文件分別會自動命名為：

• mySampleFiltered_1P.fq.gz - for paired forward reads
• mySampleFiltered_1U.fq.gz - for unpaired forward reads
• mySampleFiltered_2P.fq.gz - for paired reverse reads
• mySampleFiltered_2U.fq.gz - for unpaired reverse reads

此外，如果直接指定輸入輸出文件名，文件名后添加 .gz 后綴就是告訴軟件輸入文件是 .gz 壓縮文件，輸出文件需要用 gzip 壓縮。

過濾原理

每一步的過濾如果需要多個參數(shù)，通常用冒號:將各個參數(shù)隔開，當然參數(shù)的先后順序是有要求的。

ILUMINACLIP

從名字可以看出，這一步是為了去除 illumina 接頭的，這個軟件其實就是專為 illumina 平臺數(shù)據(jù)而設(shè)計的。

為了更好理解測序 reads 中為什么會有引物和接頭序列，我畫了一個文庫加上接頭之后的結(jié)構(gòu)示意圖，也把引物結(jié)合部位大概標了出來：

接頭和測序引物結(jié)構(gòu)示意圖

這個文庫結(jié)構(gòu)示意圖理解之后就容易理解測序過程了。

去除接頭以及引物序列看似簡單，但需要權(quán)衡靈敏度（保證接頭和引物去除干凈）和特異性（保證不是接頭和引物的序列不被誤切除），由于測序中可能存在的隨機錯誤讓去接頭這樣一個簡單的操作變的復雜。

雖然理論上接頭序列和引物序列可能出現(xiàn)在 reads 中的任何位置，但實際上序列中出現(xiàn)接頭和引物大部分情況下都是由于文庫插入片段比測序讀長短導致的，這種情況在 reads 的開頭部分是有一段可用序列的，末端包含了接頭的全長或部分序列，如果末端只有接頭的一部分序列，那么去除這殘缺的接頭序列也不是容易的事。

然而，在 PE 測序模式下如果文庫的插入片段比測序讀長短，那么 read1 和 read2 中非接頭序列的那部分會完全反向互補，Trimmomatic 有一個 ‘palindrome’ 模式會利用這個特點進行接頭序列的去除。

下圖中 A、B、C、D 四種情況就是 Trimmomatic 去除接頭和引物的四種模式：

• 紅色條形：被切除的序列
• 綠色條形：保留下來的有效讀長
• 深藍色條形：接頭序列
• 淺藍色條形：引物序列

Trimmomatic 四種模式原理圖解

A 模式：測序 reads 從起始位置開始就包含了完整的接頭序列，那么根據(jù) Illumina 測序原理，這整條 reads 都不可能包含有用序列了，整條 reads 被丟棄。

B 模式：這種相對常見，由于文庫插入片段比測序讀長短，會在 reads 末端包含部分接頭序列，若是這部分接頭序列足夠長是可以識別并去除的，但如果接頭序列太短，比接頭匹配參數(shù)設(shè)置的最短長度還短，那么就無法去除。但是，如果是 PE 測序，可以按照 D 模式去除 reads 末端的很短的接頭序列。

C 模式：PE 測序可能出現(xiàn)這種情況，正向測序和反向測序有部分完全反向互補，但是空載的文庫，兩個接頭直接互連，這樣的 reads 不包含任何有用序列，正反向測序 reads 都被丟棄。

D 模式：是 Trimmomatic 利用 PE 測序進行短接頭序列去除的典范，如果文庫插入片段比測序讀長短，利用正反向測序 reads 中一段堿基可以完全反向互補的特點，將兩個接頭序列與 reads 進行比對，同時兩條 reads 之間也互相比對，可以將 3' 末端哪怕只有 1bp 的接頭序列都可以被準確去除，相對 B 模式去除接頭污染更徹底。

Trimmomatic 使用了一種類似序列比對軟件（例如 Isaac aligner，一個超快速的 alignment 軟件）的兩步策略來搜索潛在的接頭序列。首先，使用接頭序列中的一段種子序列（seed 長度不超過 16bp）與測序 reads 進行比對，如果種子序列在測序 reads 中有足夠好的比對結(jié)果（具體由 seedMismatch 參數(shù)決定），就啟動第二步的接頭全長與 reads 比對。第一步的 seed 搜索速度很快，可以過濾掉沒有接頭污染的 reads ，這種兩步搜索的方法使得接頭序列的查找效率很高。

在第二步的接頭序列和測序 reads 全長比對統(tǒng)計比對分值時，罰分策略考慮了測序堿基的質(zhì)量值Q，每一個比對上的堿基加分 0.6，每一個錯配的堿基減分 Q/10，考慮堿基質(zhì)量值可以降低低質(zhì)量堿基（高測序錯誤率）錯配對整個比對得分的影響。在這個規(guī)則下，一段 12bp 的接頭序列完全比對到 reads 上得分為 7.2, 25bp 的接頭序列完全比對到 reads 上得分為 15。因此在 ILLUMINACLIP 參數(shù)中 simple clip threshold 的值建議為 7-15 之間（即上圖中 A/B 比對模式比對得分閾值）。

對于 palindromic 模式的比對（上圖中 D 模式），可以比對上的序列長度會更長，為了保證識別接頭序列的準確率，比對得分的閾值也更高，例如 reads的 R1 和 R2 中有 50bp 序列可以反向互補匹配，得分為 30。這種模式下，Trimmomatic 可以識別并去除 reads 中非常短的接頭序列。

ILLUMINACLIP 參數(shù)說明：按照規(guī)定順序，ILLUMINACLIP 的參數(shù)列表如下（各個參數(shù)之間以冒號分開），PE 測序需要注意最后一個參數(shù)。對于 SE 測序最后兩個參數(shù)可以不設(shè)置。

ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>

ILLUMINACLIP:TruSeq3-SE:2:30:10 #接頭和引物序列在 TruSeq3-SE 中，第一步 seed 搜索允許2個堿基錯配，palindrome 比對分值閾值 30，simple clip 比對分值閾值 10，palindrome 模式允許切除的最短接頭序列為 8bp（默認值），palindrome 模式去除與 R1 完全反向互補的 R2（默認去除）

fastaWithAdaptersEtc：指定包含接頭和引物序列（所有被視為污染的序列）的 fasta 文件路徑，Trimmomatic 自帶了一個包含 Illumina 平臺接頭和引物序列的 fasta 文件，可以直接用這個。
seedMismatches：指定第一步 seed 搜索時允許的錯配堿基個數(shù)，例如 2。
palindrome clip threshold：指定針對 PE 的 palindrome clip 模式下，需要 R1 和 R2 之間至少多少比對分值（上圖中 D 模式），才會進行接頭切除，例如 30。
simple clip threshold：指定切除接頭序列的最低比對分值（上圖 A/B 模式），通常 7-15 之間。
minAdapterLength：只對 PE 測序的 palindrome clip 模式有效，指定 palindrome 模式下可以切除的接頭序列最短長度，由于歷史的原因，默認值是 8，但實際上 palindrome 模式可以切除短至 1bp 的接頭污染，所以可以設(shè)置為 1 。
keepBothReads：只對 PE 測序的 palindrome clip 模式有效，這個參數(shù)很重要，在上圖中 D 模式下， R1 和 R2 在去除了接頭序列之后剩余的部分是完全反向互補的，默認參數(shù) false，意味著整條去除與 R1 完全反向互補的 R2，當做重復去除掉，但在有些情況下，例如需要用到 paired reads 的 bowtie2 流程，就要將這個參數(shù)改為 true，否則會損失一部分 paired reads。

看一個 PE150 數(shù)據(jù)的測試，就知道 keepBothReads 參數(shù)的重要性了：

$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51

ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 2500 Both Surviving: 1633 (65.32%) Forward Only Surviving: 828 (33.12%) Reverse Only Surviving: 12 (0.48%) Dropped: 27 (1.08%)
TrimmomaticPE: Completed successfully
# 使用 ILLUMINACLIP 默認的第六個參數(shù) false，只有 65.32% paired reads 保留下來

$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51

ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 2500 Both Surviving: 2439 (97.56%) Forward Only Surviving: 22 (0.88%) Reverse Only Surviving: 16 (0.64%) Dropped: 23 (0.92%)
TrimmomaticPE: Completed successfully
# 將 ILLUMINACLIP 第六個參數(shù)改為 true，其余所有參數(shù)均相同，結(jié)果有 97.56% paired reads 保留下來

SLIDINGWINDOW

滑窗剪切，統(tǒng)計滑窗口中所有堿基的平均質(zhì)量值，如果低于設(shè)定的閾值，則切掉窗口。
SLIDINGWINDOW 參數(shù)如下：

SLIDINGWINDOW:<windowSize>:<requiredQuality>

SLIDINGWINDOW:4:15  #設(shè)置 4bp 窗口，堿基平均質(zhì)量值閾值 15

widowSize:設(shè)置窗口大小
requiredQuality：設(shè)置窗口堿基平均質(zhì)量閾值

MAXINFO

包含一個可以自動調(diào)整的過濾條件，在保留盡可能長的序列和保持序列中堿基測序錯誤率盡可能低之間進行平衡，以達到 trim 之后保留序列的價值最大化。

對于不同的應用場景，一條 reads 序列的價值由以下三個因素決定：

1. 允許的最短 read 長度：只有在測序 reads 足夠長的情況下，才能將 reads map 到參考序列的唯一位置上。過短的序列往往能 map 到參考基因組上的多個位置（特異性差），這種非特異性的 map 結(jié)果可以給出的有價值信息很少（對后續(xù)分析沒什么用）。保證 reads 可以 map 到參考序列的唯一位置上需要的 reads 最短長度也與參考序列本身的長度和序列復雜度有關(guān)，但是通常最短 reads 要求不低于 40bp。
2. 需求的額外 reads 長度：在保證 reads 可以 map 到參考序列上的唯一位置的情況下， reads 序列越長，可能對后續(xù)分析越有利。當然，這也需要看應用場景，對于 RNA-seq 來說只需要統(tǒng)計某個位置 map 上的 reads 數(shù)目，這種情況下 reads 長度只要滿足唯一比對要求即可。而對于 de-novo 序列組裝或者 Resequencing 檢測變異，reads 序列長度越長越好，reads 越長可以為后續(xù)分析提供越多的有效證據(jù)。
3. 對測序錯誤的敏感性：即測序錯誤率對后續(xù)分析工作是否造成嚴重影響，有的應用中對 reads 測序錯誤零容忍，就需要非常嚴格的過濾條件，而有的應用對測序錯誤不敏感，就可以使用寬松的過濾條件。

MAXINFO 有兩個參數(shù)，第一個 target read length 控制上面的因素一，即允許的最短 read 長度。第二個參數(shù) strictness 是控制因素二和因素三之間的平衡，即滿足最短 read 長度的情況下，是保留盡可能多的堿基，還是保證盡可能低的測序錯誤率。

MAXINFO 過濾從 reads 3' 端開始進行剪切，在考慮上述三個因素的情況下統(tǒng)計所有可能的 trim 方式的到的 clean reads 的 INFO 分值（即所謂的 reads 價值），這三個因素分別以不同的方式影響最終的 reads INFO 分值：

1. 最短 read 長度：最終保留的 reads 長度與 INFO 分值之間是 logistic 函數(shù)關(guān)系，即 INFO 分值與 reads 長度之間是 S 形增長曲線（類似大腸桿菌生長曲線），大概相當于，當保留下來的 reads 長度小于最短長度時，INFO 分值隨著 reads 長度增加呈指數(shù)級增長，當保留下來的 reads 長度達到最短 read 長度要求之后，INFO 分值的增長會變緩。
2. 超出最短 read 長度的部分：INFO 分值會隨著額外的 reads 長度增加了線性增長，線性系數(shù)與 MAXINFO 的第二個參數(shù)有關(guān)，為 1 - strictness。
3. 測序錯誤敏感性：最終保留下來的所有堿基的質(zhì)量值都被用來計算 reads 包含測序錯誤的概率，出現(xiàn)測序錯誤的概率提升會降低 INFO 分值，錯誤概率對 INFO 分值的影響權(quán)重與 strictness 相關(guān)。

針對一條 read 的任何可能的剪切方式都計算出 INFO 分值，最終的 reads 長度和切除的堿基由 INFO 最大值決定。實際上這三個影響因子作用的方式不同：

1. 當 reads 長度比最短 read 長度還短時 INFO 分值被序列長度主導，因為太短的 reads 根本無法提供足夠的有用信息。
2. 當 reads 長度滿足最短 read 長度要求時，reads 的長度因素就不再是影響 INFO 分值的主導因素了，而且一旦 reads 長度足夠長時，由于 logistic 函數(shù)的關(guān)系，reads 長度就不怎么貢獻 INFO 分值了。
3. 超出最短 read 長度要求的堿基對 INFO 貢獻有限，而另一方面，由于堿基的增加導致 reads 中出現(xiàn)測序錯誤的概率增大，這會導致對 INFO 的罰分。控制這超出最短 read 長度的序列對 INFO 分值的影響到底是罰分還是得分，就看 strictness 參數(shù)了。
4. 對于大部分足夠長的 reads，保留超出最低要求的堿基序列是利還是弊，會依照堿基質(zhì)量值分布和 strictness 參數(shù)來決定如何 trim。

參數(shù)說明：

MAXINFO:<targetLength>:<strictness>

targetLength:使得 reads 可以 map 到參考序列上唯一位置的最短長度（likely）。
strictness：一個介于 0 - 1 之間的小數(shù)，決定如何平衡 最大化 reads 長度 或者 最小化 reads 出現(xiàn)錯誤的概率，當參數(shù)設(shè)置小于 0.2 時傾向于最大化 reads 長度，當參數(shù)設(shè)置大于 0.8 時傾向于最小化 reads 中出現(xiàn)測序錯誤的概率。

LEADING

從 reads 的起始端開始切除質(zhì)量值低于設(shè)定的閾值的堿基，直到有一個堿基其質(zhì)量值達到閾值。

LEADING:<quality>

quality:設(shè)定堿基質(zhì)量值閾值，低于這個閾值將被切除。

TRAILING

從 reads 的末端開始切除質(zhì)量值低于設(shè)定閾值的堿基，直到有一個堿基質(zhì)量值達到閾值。Illumina 平臺有些低質(zhì)量的堿基質(zhì)量值被標記為 2 ，所以設(shè)置為 3 可以過濾掉這部分低質(zhì)量堿基。官方推薦使用 Sliding Window 或 MaxInfo 來代替 LEADING 和 TAILING。

TRAILING:<quality>

quality:設(shè)定堿基質(zhì)量值閾值，低于這個閾值將被切除。

CROP

不管堿基質(zhì)量，從 reads 的起始開始保留設(shè)定長度的堿基，其余全部切除。一刀切，把所有 reads 切成相同的長度。

CROP:<length>

length：reads 從末端除之后保留下來的序列長度

HEADCROP

不管堿基質(zhì)量，從 reads 的起始開始直接切除部分堿基。

HEADCROP:<length>

length：從 reads 的起始開始切除的堿基數(shù)

MINLEN

設(shè)定一個最短 read 長度，當 reads 經(jīng)過前面的過濾之后，如果保留下來的長度低于這個閾值，整條 reads 被丟棄。被丟棄的 reads 數(shù)會被統(tǒng)計在 Trimmomatic 日志的 dropped reads 中。

MINLEN:<length>

length：可被保留的最短 read 長度

TOPHRED33

此選項可以將過濾之后的 Fastq 文件中質(zhì)量值這一行轉(zhuǎn)為 phred-33 格式。

TOPHRED64

此選項可以將過濾之后的 Fastq 文件中質(zhì)量值這一行轉(zhuǎn)為 phred-64 格式。

接頭序列和引物序列

Trimmomatic 也可以自己制作包含接頭和引物序列的 fasta 文件，格式可以參考軟件自帶的 adapters 文件夾中的格式。

adapters 文件夾中包含 illumina 測序 TruSeq2，TruSeq3 針對 SE 和 PE 的通用接頭和引物序列。