Trimmomatic

Trimmomatic

Trimmomatic 支持多线程，处理数据速度快，主要用来去除 Illumina 平台的 Fastq 序列中的接头，并根据碱基质量值对 Fastq 进行修剪。软件有两种过滤模式，分别对应 SE 和 PE 测序数据，同时支持 gzip 和 bzip2 压缩文件。另外也支持 phred-33 和 phred-64 格式互相转化。

#安装：

不建议直接试用conda安装，调用命令时会提醒找不到 TruSeq3-SE

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
java -jar trimmomatic-0.39.jar -h

Usage:

由于 Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关，因此，如果需要去接头，建议第一步就去接头，否则接头序列被其他的过滤参数剪切掉部分之后就更难匹配更难去除干净了。

单末端测序模式

在 SE 模式下，只有一个输入文件和一个过滤之后的输出文件：

java -jar <path to trimmomatic jar> SE [-threads <threads>] [-phred33 | -phred64] [-trimlog
<logFile>] <input> <output> <step 1> <step 2> ...

-trimlog 参数指定了过滤日志文件名，日志中包含以下四列内容:

• read ID
• 过滤之后剩余序列长度
• 过滤之后的序列起始碱基位置（序列开头处被切掉了多少个碱基）
• 过滤之后的序列末端碱基位置
• 序列末端处被剪切掉的碱基数

由于生成的 trimlog 文件中包含了每一条 reads 的处理记录，因此文件体积巨大（GB 级别），如果后面不会用到 trim 日志，建议不要使用这个参数。

双末端测序模式

在 PE 模式下，有两个输入文件，正向测序序列和反向测序序列，但是过滤之后输出文件有四个，过滤之后双端序列都保留的就是 paired，反之如果其中一端序列过滤之后被丢弃了另一端序列保留下来了就是 unpaired 。

java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog
<logFile>] >] [-basein <inputBase> | <input 1> <input 2>] [-baseout <outputBase> |
<paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> <step 2> ...

其中 -phred33 和 -phred64 参数指定 fastq 的质量值编码格式，如果不设置这个参数，软件会自动判断输入文件是哪种格式（v0.32 之后的版本都支持），虽然软件默认的参数是 phred64，如果不确定序列是哪种质量编码格式，可以不设置这个参数。

Example:

cat ../list.txt |while read id ;do (arr=($id);java -jar /public/home/djs/software/Trimmomatic-0.39/trimmomatic-0.39.jar PE -phred33 ../${arr[0]} ../${arr[1]} -baseout ${arr[0]}.trim.fastq.gz ILLUMINACLIP:/public/home/djs/software/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:150 );done

#Trimmomatic 过滤步骤:
Trimmomatic 过滤数据的步骤与命令行中过滤参数的顺序有关，通常的过滤步骤如下：

1. ILLUMINACLIP: 过滤 reads 中的 Illumina 测序接头和引物序列，并决定是否去除反向互补的 R1/R2 中的 R2。
2. SLIDINGWINDOW: 从 reads 的 5’ 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗。
3. MAXINFO: 一个自动调整的过滤选项，在保证 reads 长度的情况下尽量降低测序错误率，最大化 reads 的使用价值。
4. LEADING: 从 reads 的开头切除质量值低于阈值的碱基。
5. TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基。
6. CROP: 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度。
7. HEADCROP: 从 reads 的开头切掉指定数量的碱基。
8. MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads。
9. AVGQUAL: 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads。
10. TOPHRED33: 将 reads 的碱基质量值体系转为 phred-33。
11. TOPHRED64: 将 reads 的碱基质量值体系转为 phred-64。

输入输出文件

PE 模式的两个输入文件：sample_R1.fastq sample_R2.fastq以及四个输出文件：sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq sample_paired_R1.clean.fastq sample_unpaired_R1.clean.fastq

通常 PE 测序的两个文件，R1 和 R2 的文件名是类似的，因此可以使用 -basein 参数指定其中 R1 文件名即可，软件会推测出 R2 的文件名，但是这个功能实测并不好用，因为软件只能自动识别推测三种种格式的 -basein:

• Sample_Name_R1_001.fq.gz -> Sample_Name_R2_001.fq.gz
• Sample_Name.f.fastq -> Sample_Name.r.fastq
• Sample_Name.1.sequence.txt -> Sample_Name.2.sequence.txt
  建议不用 -basein 参数，直接指定两个文件名（R1 和 R2）作为输入。

输出文件有四个，当然也可以像上文一样指定四个文件名，但是参数太长有点麻烦，有个省心的方法，使用 -baseout 参数指定输出文件的 basename，软件会自动为四个输出文件命名。例如 -baseout mySampleFiltered.fq.gz ，文件名中添加 .gz 后缀，软件会自动将输出结果进行 gzip 压缩。输出的四个文件分别会自动命名为：

• mySampleFiltered_1P.fq.gz - for paired forward reads
• mySampleFiltered_1U.fq.gz - for unpaired forward reads
• mySampleFiltered_2P.fq.gz - for paired reverse reads
• mySampleFiltered_2U.fq.gz - for unpaired reverse reads
  此外，如果直接指定输入输出文件名，文件名后添加 .gz 后缀就是告诉软件输入文件是 .gz 压缩文件，输出文件需要用 gzip 压缩。

ILUMINACLIP

每一步的过滤如果需要多个参数，通常用冒号:将各个参数隔开，当然参数的先后顺序是有要求的。

ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>

ILLUMINACLIP:TruSeq3-SE:2:30:10 #接头和引物序列在 TruSeq3-SE 中，第一步 seed 搜索允许2个碱基错配，palindrome 比对分值阈值 30，simple clip 比对分值阈值 10，palindrome 模式允许切除的最短接头序列为 8bp（默认值），palindrome 模式去除与 R1 完全反向互补的 R2（默认去除）

fastaWithAdaptersEtc：指定包含接头和引物序列（所有被视为污染的序列）的 fasta 文件路径，Trimmomatic 自带了一个包含 Illumina 平台接头和引物序列的 fasta 文件，可以直接用这个。
seedMismatches：指定第一步 seed 搜索时允许的错配碱基个数，例如 2。
palindrome clip threshold：指定针对 PE 的 palindrome clip 模式下，需要 R1 和 R2 之间至少多少比对分值（上图中 D 模式），才会进行接头切除，例如 30。
simple clip threshold：指定切除接头序列的最低比对分值（上图 A/B 模式），通常 7-15 之间。
minAdapterLength：只对 PE 测序的 palindrome clip 模式有效，指定 palindrome 模式下可以切除的接头序列最短长度，由于历史的原因，默认值是 8，但实际上 palindrome 模式可以切除短至 1bp 的接头污染，所以可以设置为 1 。
keepBothReads：只对 PE 测序的 palindrome clip 模式有效，这个参数很重要，在上图中 D 模式下， R1 和 R2 在去除了接头序列之后剩余的部分是完全反向互补的，默认参数 false，意味着整条去除与 R1 完全反向互补的 R2，当做重复去除掉，但在有些情况下，例如需要用到 paired reads 的 bowtie2 流程，就要将这个参数改为 true，否则会损失一部分 paired reads。

看一个 PE150 数据的测试，就知道 keepBothReads 参数的重要性了：

$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51

ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 2500 Both Surviving: 1633 (65.32%) Forward Only Surviving: 828 (33.12%) Reverse Only Surviving: 12 (0.48%) Dropped: 27 (1.08%)
TrimmomaticPE: Completed successfully
# 使用 ILLUMINACLIP 默认的第六个参数 false，只有 65.32% paired reads 保留下来

$ java -jar trimmomatic-0.36.jar PE -phred33 F-2-test_R1.fastq.gz F-2-test_R2.fastq.gz -baseout F-2.fastq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:51

ILLUMINACLIP: Using 1 prefix pairs, 0 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 2500 Both Surviving: 2439 (97.56%) Forward Only Surviving: 22 (0.88%) Reverse Only Surviving: 16 (0.64%) Dropped: 23 (0.92%)
TrimmomaticPE: Completed successfully
# 将 ILLUMINACLIP 第六个参数改为 true，其余所有参数均相同，结果有 97.56% paired reads 保留下来

SLIDINGWINDOW

滑窗剪切，统计滑窗口中所有碱基的平均质量值，如果低于设定的阈值，则切掉窗口。
SLIDINGWINDOW 参数如下：

SLIDINGWINDOW:<windowSize>:<requiredQuality>

SLIDINGWINDOW:4:15  #设置 4bp 窗口，碱基平均质量值阈值 15

LEADING

从 reads 的起始端开始切除质量值低于设定的阈值的碱基，直到有一个碱基其质量值达到阈值。

LEADING:<quality>

quality:设定碱基质量值阈值，低于这个阈值将被切除。

TRAILING

从 reads 的末端开始切除质量值低于设定阈值的碱基，直到有一个碱基质量值达到阈值。Illumina 平台有些低质量的碱基质量值被标记为 2 ，所以设置为 3 可以过滤掉这部分低质量碱基。官方推荐使用 Sliding Window 或 MaxInfo 来代替 LEADING 和 TAILING。

TRAILING:<quality>

quality:设定碱基质量值阈值，低于这个阈值将被切除。

CROP

不管碱基质量，从 reads 的起始开始保留设定长度的碱基，其余全部切除。一刀切，把所有 reads 切成相同的长度。

CROP:<length>

length：reads 从末端除之后保留下来的序列长度

HEADCROP

不管碱基质量，从 reads 的起始开始直接切除部分碱基。

HEADCROP:<length>

length：从 reads 的起始开始切除的碱基数

MINLEN

设定一个最短 read 长度，当 reads 经过前面的过滤之后，如果保留下来的长度低于这个阈值，整条 reads 被丢弃。被丢弃的 reads 数会被统计在 Trimmomatic 日志的 dropped reads 中。

MINLEN:<length>

length：可被保留的最短 read 长度

TOPHRED33

此选项可以将过滤之后的 Fastq 文件中质量值这一行转为 phred-33 格式。

TOPHRED64

此选项可以将过滤之后的 Fastq 文件中质量值这一行转为 phred-64 格式。

接头序列和引物序列

Trimmomatic 也可以自己制作包含接头和引物序列的 fasta 文件，格式可以参考软件自带的 adapters 文件夹中的格式。
adapters 文件夹中包含 illumina 测序 TruSeq2，TruSeq3 针对 SE 和 PE 的通用接头和引物序列。