過(guò)濾軟件的比較與選擇：cutadapt/fastp/trimmomatic

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

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

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

CutAdapt，2010

1. 基于Python，作者是個(gè)德國(guó)人，長(zhǎng)得還挺帥氣(???) 不過(guò)都9年過(guò)去了，嗯。。
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長(zhǎng)度大于50才保留
-o -p #生成文件：過(guò)濾后的R1 R2
read1.fastq read2.fastq #輸入文件

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

fastp，2018

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

能進(jìn)行質(zhì)控，生成比f(wàn)astqc美觀、全面的報(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. 最簡(jiǎn)單的命令：
   fastp -i r1.fq -o rr1.fq -i r2.fq -o rr2.fq
5. 這篇介紹寫的不錯(cuò)：知乎
   但他說(shuō)一般下機(jī)數(shù)據(jù)要經(jīng)過(guò)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)！
   所以建議先去接頭，否則接頭被剪更無(wú)法有效去除。

3. PE數(shù)據(jù)常用參數(shù)：
   ILLMINACLIP: 注意以下參數(shù)以：隔開(kāi)
   <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 模式下可以切除的接頭序列最短長(zhǎng)度，默認(rèn)值是 8。但實(shí)際上 palindrome 模式可以切除短至 1bp 的接頭污染，所以可以設(shè)置為 1。
   <keepBothReads> 重要參數(shù)！第一次做的時(shí)候沒(méi)加這個(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的起始端開(kāi)始切除質(zhì)量值低于設(shè)定的閾值的堿基，直到有一個(gè)堿基其質(zhì)量值達(dá)到閾值。一般用LEADING:3???
TRAILING:20 #一般用3，因?yàn)镮llumina 平臺(tái)有些低質(zhì)量的堿基質(zhì)量值被標(biāo)記為2，所以設(shè)置為 3 可以過(guò)濾掉這部分低質(zhì)量堿基
SLIDINGWINDOW:4:20 #滑窗剪切，統(tǒng)計(jì)滑窗口中所有堿基的平均質(zhì)量值，如果低于設(shè)定的閾值，則切掉窗口。此處設(shè)置4bp窗口，閾值20，一般閾值用15。
MINLEN:50 #可被保留的最短 read 長(zhǎng)度

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

幾種軟件比較

以上?？梢詔est一下trimmomatic的true參數(shù)，還有fastp試一下到底強(qiáng)大在哪里。