Trimmomatic

前述

Trimmomatic工具是用于illumina二代测序数据的reads处理，主要对接头(adapter)序列和低质量序列进行过滤。下面是使用该工具处理双端测序(PE)数据时，常用参数的一些说明。

参考文档

Trimmomatic工具的参考文献
Trimmomatic工具官网
Trimmomatic工具使用手册

软件使用

执行命令

## 双端测序数据使用方法
# 使用v0.32版本:
1.  java -jar trimmomatic-0.32.jar PE \
2.		[-threads <threads>] \
3.		[-phred33|-phred64] \
4.		[-trimlog <trimLogFile>] \
5.		[<inputFile1> <inputFile2>] \
6.		[<outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] \
7.		[ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>] \
8.		[SLIDINGWINDOW:<windowSize>:<requiredQuality>] \
9.		[LEADING:<quality>] \
10.		[TRAILING:<quality>] \
11.		[MINLEN:<length>]

# 举例:
1. java -jar trimmomatic-0.32.jar PE \
2.		-threads 16 \
3.		-phred33 \
4.		-trimlog trim.log \
5.		input_forward_R1.fq.gz input_reverse_R2.fq.gz \
6.		output_forward_paired_R1.fq.gz  output_forward_unpaired_R1.fq.gz output_reverse_paired_R2.fq.gz output_reverse_unpaired_R2.fq.gz \
7.		ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true \
8.		SLIDINGWINDOW:5:20 \
9.		LEADING:3 \
10.		TRAILING:3 \
11.		MINLEN:36

参数说明

PE 设置使用trimmomatic处理双端数据，单端数据用(‘SE’)
-threads 16 设置线程数为16
-phred33 设置碱基的质量格式(默认-phred64，自v0.32版本之后可自动识别是phred33还是phred64)
-trimlog trim.log 设置trimmommatic工具处理的日志文件为’trim.log’，每两行为一对reads信息

生成的日志文件’trim.log’包含5列信息：

  1) read名称
  2) 切除后剩余的read长度
  3) 切除后第一个碱基所在位置(0-base)，也就是起始位置切除了几个碱基
  4) 切除后最后一个碱基所在位置
  5) read末尾切除了几个碱基
  # 由于该生成文件较大，如后续步骤不需该文件信息，可考虑不设置

input_forward_R1.fq.gz 输入的forward正向链R1对应的fastq文件
input_reverse_R2.fq.gz 输入的reverse反向链R2对应的fastq文件
output_forward_paired_R1.fq.gz 处理后输出的R1对应reads成对fastq文件
output_forward_unpaired_R1.fq.gz 处理后输出的R1对应reads不成对的fastq文件
output_reverse_paired_R2.fq.gz 处理后输出的R2对应reads成对fastq文件
output_reverse_unpaired_R2.fq.gz 处理后输出的R2对应reads不成对的fastq文件
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true 切除illumina接头参数设置。

说明ILLUMINACLIP的各参数之前，先解释一下要使用的两种切除接头的模式和比对分值计算方法

1. Simple mode

simple模式，它可用于切除任何序列(any technical sequence, 暂且称之为污染序列)。该模式的方法是将污染序列直接与一条read进行比对（局部比对)，将进行下面4个步骤(情形)：

 Simple mode是根据与污染序列的比对情况判断切除序列
 A) 污染序列3'端的一部分与read5'端有overlap，去除整条read；
 B) 污染序列从read的5'端开始能与read完全overlap, 去除整条read；
 C) 污染序列在read中间有overlap，5'端一部分没有overlap的保留，overlap的部分及之后的进行切除；
 D) 污染序列3'端的一部分与read有overlap，切除overlap部分.

2. Palindrome mode
   palindrome模式，它仅用于read测穿(read-through)的情形， 该模式的方法是：在识别接头序列之前，先将接头序列加在成对reads起始（‘adapter ligated’ read，read-with-adapter, 暂且称加adapter后的这对reads为a-R1和a-R2），再将a-R1和a-R2进行比对（全局比对)。 该模式将进行如下4个步骤(情形)：
   Palindrome mode
   

 Palindrome mode: 是根据a-R1, a-R2的overlap判断切除的序列
 A) overlap是在接头序列与reads的起始之间，即overlap只有接头序列，没有有用的信息，这种情形将这两条read都整条去除；
 B) overlap有一部分是DNA插入片段，一部分是整个接头序列，该情形会将这对reads从overlap序列及剩余序列(3'端)切除；
 C) overlap中有接头序列的一小部分，与B)相同方式切除；
 D) overlap中不存在接头序列.

再说明一下比对分值的计算方法：

match：匹配一个碱基加分，分值为$log_10{4}$(约等于0.602)；
 mismatch: 错配一个碱基减分，分值为Q/10，
 	   其中，Q为对应碱基的碱基质量值，由于一般碱基质量值区间为: [0-40]，即错配的罚分值区间[0-4]。
 	   该分值公式说明碱基质量越高的碱基出现错配时罚分越多.
 ## 例如：
 # 12bp的接头序列能够完全比对到read上的分值为：12 * 0.602 ≈ 7；
 # 25bp的接头序列完全比对到read上的分值为：25 * 0.602 ≈ 15；
 # a-R1和a-R2有50bp完全匹配的分值为50 * 0.602 ≈ 30.

下面是ILLUMINACLIP的各参数说明：

 ":<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>"
 - ":TruSeq3-PE.fa"  # <fastaWithAdaptersEtc>是fasta格式的接头序列文件路径
 - ":2"  # <seed mismatches>是将接头序列的一段(不超过16bp)作为seed，与reads进行比对，能够容忍的最大错配(mismatch)数，这里是最多2个错配
 - ":30"  # <palindrome clip threshold>是 a-R1, a-R2的比对分值阈值，达到阈值，才进行切除，这里设置阈值为30（大约比对50bp碱基）
 - ":10"  # <simple clip threshold>是任意(接头)序列与read比对最低分阈值，大于这个阈值，才进行切除，这里设置阈值为10（大约比对17bp碱基）
 - ":8"  # <minAdapterLength>只作用于Palindrome模式，是设置检测到接头序列的最小长度（默认为8，甚至可设置为1）
 - ":true"  # <keepBothReads>只作用于palindrome模式，是设置是否保留反向链，这里是说去除接头序列后，由于正反链包含相同的序列信息(尽管序列是反向互补的)，默认情况(":false")下会去除反向链，设置为":True"则保留反向链

SLIDINGWINDOW:5:20 设置滑动窗口阈值，以为5bp为窗口，这5bp碱基的平均质量值低于20，要进行切除
LEADING:3 设置从reads起始开始，去除质量低于阈值或为’N’的碱基，直到达到阈值不再去除，这里设置阈值为3。
TRAILING:3 设置从reads末尾开始，去除质量低于阈值的碱基或为’N’的碱基直到达到阈值不再去除，这里设置阈值为3，这种方法是去除特定的illumina平台低质量区域(由于illumina会将低质量碱基标记为2)，官方操作文档中建议使用 SLIDINGWINDOW 或 MAXINFO代替[这里未给出MAXINFO参数说明] ）。
MINLEN:36 设置read切除后的最短长度，这里设置长度至少为36bp，长度小于36bp的reads被去除。

除了这些参数外，如果我们想要切除reads开头的碱基，可以使用HEADCROP 参数。

HEADCROP:9表示切除reads开头9个碱基。

如果想切除reads末端的序列，仅保留前一部分，可以使用CROP 参数。
假设reads长度为150，CROP:130代表切除reads后20个碱基，仅保留前130个碱基。