【bioinfo】学习使用一款数据质控软件（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’)

• -thread 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文件

• 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被去除。

• 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: 是根据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的各参数说明：

"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"则保留反向链

关于使用的接头序列文件，在Trimmomatic工具使用手册中也给出说明：

• 如果接头文件中序列的名称是Prefix开头，并且有"/1"和"/2"分别对应正向和反向的接头序列，才会使用Palindrome的模式进行trimming（该模式会同时考虑PE read的overlap去接头序列）;
• 否则使用Simple模式，如果使用Simple模式按默认的分值阈值（ILLUMINACLIP: Adaper:2:30:10:8:true），需要大约17bp match的adapter比对上才去除。
• 另外，Palindrome模式是在R1和R2可以match的情况下，使用Prefix接头去除，如果R1/R2不能match，Prefix的接头会按Simple模式去除。

修改日志

• 2024-03-25: 增加adapter文件序列名说明
• 2024-01-15: 修改图
• 2019-06-24: 新增