過濾軟件的比較與選擇：cutadapt/fastp/trimmomatic

注：還沒有完全搞明白，先總結(jié)一下特點(diǎn)和使用栓拜，之后再慢慢體會(huì)座泳、總結(jié)經(jīng)驗(yàn)
本次只針對(duì)雙端PE
算法都沒好好讀，因?yàn)榭床欢?=

首先幕与，我們對(duì)數(shù)據(jù)進(jìn)行過濾钳榨，是為了：

去掉接頭
去掉低質(zhì)量reads
去掉污染序列
在盡量去掉上述序列的同時(shí)，保留盡可能多的有用數(shù)據(jù)纽门，減少損失

CutAdapt，2010

1. 基于Python营罢，作者是個(gè)德國人赏陵，長得還挺帥氣(???) 不過都9年過去了，嗯饲漾。蝙搔。
2. 不僅支持illumina，還支持SOLID考传，454等平臺(tái)產(chǎn)出的數(shù)據(jù)
3. 支持輸入.gz
4. 需要自己先檢測(cè)接頭類型（fastqc等）吃型，然后搜索接頭序列是啥，手動(dòng)輸入到參數(shù)里僚楞。但是有一個(gè)參數(shù) -n勤晚，若是兩種接頭枉层，也可以指定然后去除：-n 2
5. 一般命令：

cutadapt -a -A #a是read1的3'接頭，A是read2的3'接頭(5'接頭的反向互補(bǔ)序列)
-e 0.1 -0.5 -m 50 #去除接頭后read長度大于50才保留
-o -p #生成文件：過濾后的R1 R2
read1.fastq read2.fastq #輸入文件

4. 本次分析沒用赐写，所以詳細(xì)參數(shù)可以閱讀--help

fastp鸟蜡，2018

1. 基于c++這種強(qiáng)大的語言所以算法比較高效，中科院深圳所發(fā)的挺邀。還沒用過揉忘，不過身邊做RNA-Seq的倆師兄強(qiáng)烈推薦，有空可以test一下端铛。
2. 主題就是ultra-fast泣矛，all-in-one，而且是只處理FASTQ也就是illuminate下機(jī)數(shù)據(jù)
3. 特點(diǎn)：

能進(jìn)行質(zhì)控禾蚕，生成比fastqc美觀您朽、全面的報(bào)告，但是我看了一遍夕膀，不如fastqc直觀虚倒、fresh-friendly
號(hào)稱去除低質(zhì)量序列的方法類似于trimmomatic但是更快
自動(dòng)識(shí)別序列并去除
支持illuminate short read，也一定程度支持Pacbio/Nanopore long reads产舞，具體支持到什么程度魂奥，需要試驗(yàn)。
參數(shù)眾多易猫，但是挺有條理的耻煤，而且很多都是默認(rèn)不是必需參數(shù)，不會(huì)“新手退散”

4. 最簡單的命令：
   fastp -i r1.fq -o rr1.fq -i r2.fq -o rr2.fq
5. 這篇介紹寫的不錯(cuò)：知乎
   但他說一般下機(jī)數(shù)據(jù)要經(jīng)過fastqc+cutadapt+trimmomatic准颓，有點(diǎn)不太理解哈蝇，要這么麻煩嗎？

Trimmomatic攘已，2014

1. 也是很好用的炮赦，引用量超高，good at去除低質(zhì)量reads样勃，只針對(duì)illuminate數(shù)據(jù)

2. 最重要的特點(diǎn)：對(duì)數(shù)據(jù)的處理步驟與參數(shù)的順序有關(guān)吠勘！
   所以建議先去接頭，否則接頭被剪更無法有效去除峡眶。

3. PE數(shù)據(jù)常用參數(shù)：
   ILLMINACLIP: 注意以下參數(shù)以：隔開
   <fastaWithAdaptersEtc>: 指定包含接頭和引物序列（所有被視為污染的序列）的 fasta 文件
   <seed mismatches>: 第一步seed搜索時(shí)允許的mismatch個(gè)數(shù)剧防，一般2。
   <palindrome clip threshold>: 指定針對(duì) PE的palindrome clip模式下辫樱，需要R1和 R2之間至少多少比對(duì)分值峭拘，才會(huì)進(jìn)行接頭切除，例如30。
   <simple clip threshold>: 指定切除接頭序列的最低比對(duì)分值鸡挠，一般7-15之間辉饱。
   <minAdapterLength>: 只對(duì) PE 測(cè)序的 palindrome clip 模式有效，指定 palindrome 模式下可以切除的接頭序列最短長度宵凌，默認(rèn)值是 8鞋囊。但實(shí)際上 palindrome 模式可以切除短至 1bp 的接頭污染，所以可以設(shè)置為 1瞎惫。
   <keepBothReads> 重要參數(shù)溜腐！第一次做的時(shí)候沒加這個(gè)參數(shù)，結(jié)果20%+的數(shù)據(jù)Unpaired瓜喇，扔掉不現(xiàn)實(shí)挺益，比對(duì)處理太麻煩！正確用法：只對(duì) PE 測(cè)序的 palindrome clip 模式有效乘寒，R1 和 R2 在去除了接頭序列之后剩余的部分是完全反向互補(bǔ)的望众，默認(rèn)參數(shù) false，意味著整條去除與 R1 完全反向互補(bǔ)的 R2伞辛，當(dāng)做重復(fù)去除掉烂翰，但在有些情況下，例如需要用到 paired reads 的 bowtie2 流程蚤氏，就要將這個(gè)參數(shù)改為 true甘耿，否則會(huì)損失一部分 paired reads。

4. 本次所用命令：（也是公司報(bào)告中所用的）

java -jar trimmomatic-0.38.jar PE -threads 2 #雙端模式竿滨，兩個(gè)線程
ILLUMINACLIP: #顧名思義佳恬，去illumina接頭
TruSeq3-PE.fa: #接頭文件，需要指定全路徑
2:30:10 # 默認(rèn)格式為 2:30:10:8:false于游，可改做：2:30:10:8:true
LEADING:20 #從reads的起始端開始切除質(zhì)量值低于設(shè)定的閾值的堿基毁葱，直到有一個(gè)堿基其質(zhì)量值達(dá)到閾值。一般用LEADING:3???
TRAILING:20 #一般用3贰剥，因?yàn)镮llumina 平臺(tái)有些低質(zhì)量的堿基質(zhì)量值被標(biāo)記為2倾剿，所以設(shè)置為 3 可以過濾掉這部分低質(zhì)量堿基
SLIDINGWINDOW:4:20 #滑窗剪切，統(tǒng)計(jì)滑窗口中所有堿基的平均質(zhì)量值蚌成，如果低于設(shè)定的閾值前痘，則切掉窗口。此處設(shè)置4bp窗口笑陈，閾值20，一般閾值用15葵袭。
MINLEN:50 #可被保留的最短 read 長度

trimmomatic PE模式默認(rèn)處理2個(gè)文件涵妥，也就是說，sh腳本中使用本辦法只能一次列舉R1 R2兩個(gè)文件，不能 In_R1 In_R2 IP_R1 IP_R2這樣四個(gè)文件都列出來蓬网，事實(shí)證明會(huì)報(bào)錯(cuò)窒所，trimmomatic有點(diǎn)傻傻的不知道第三個(gè)開始的文件該干嘛。
所以要批量做帆锋，需要寫循環(huán)吵取，或者是認(rèn)真閱讀使用說明的參數(shù)。

5. trimmomatic的更多解讀可以參考這個(gè)锯厢，寫得很詳細(xì)皮官。目前我理解的是以上。

最后附一個(gè)圖：
出自：Chen et al. Source Code for Biology and Medicine 2014, 9:8. Software for pre-processing Illumina nextgeneration sequencing short read sequences.

幾種軟件比較

以上实辑∞嗲猓可以test一下trimmomatic的true參數(shù)，還有fastp試一下到底強(qiáng)大在哪里剪撬。