写在前面
Trimmomatic工具是用于illumina二代测序数据的reads处理,主要对接头(adapter)序列和低质量序列进行过滤。下面是使用该工具处理双端测序(PE)数据时,常用参数的一些说明。
参考文档
软件使用
执行命令
## 双端测序数据使用方法
# 使用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: 新增