Trimmomatic安装和使用

2022-06-14

Trimmomatic: A flexible read trimming tool for Illumina NGS data

Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics, btu170.

一、Trimmomatic简介

官网：http://www.usadellab.org/cms/?page=trimmomatic

官方文档：TrimmomaticManual_V0.32.pdf

Github网址：https://github.com/usadellab/Trimmomatic

Trimmomatic是一个快速的多线程命令行工具，可以用来处理Illumina (FASTQ)数据以及删除adapter。程序有两种主要模式:Paired-End模式和Single-End模式。Paired-End模式将保持read-pairs的对应关系，并利用read-pairs所包含的附加信息，更好地找到文库制备过程中引入的adapter或PCR引物片段。支持使用gzip或bzip2压缩的文件，并通过使用.gz或.bz2文件扩展名来识别。

Trimmomatic是基于Java开发的，因此需要提前安装Java，才能使用Trimmomatic。目前的Trimmomatic主要包含下面步骤

• ILLUMINACLIP: 从read中切除adapter和Illumina的特殊序列。
• SLIDINGWINDOW: 从read的5'端开始滑动窗口，当窗口内的average quality小于阈值时，切除后边的碱基.。
• MAXINFO: 自适应的trim策略，可以平衡read的长度和错误率保证read的长度最长。
• LEADING: 切除read的开头质量低于阈值的碱基。
• TRAILING: 切除read的末尾质量低于阈值的碱基。
• CROP: 从read的末尾切除特定长度。
• HEADCROP: 从read的开头切除特定长度。
• MINLEN: 舍弃小于最小长度的reads 。
• AVGQUAL: 舍弃平均质量 小于阈值的reads。
• TOPHRED33: 将质量分数转换为Phred-33
• TOPHRED64: 将质量分数转换为Phred-64

二、Trimmomatic安装(trimmomatic-0.39)

1、下载安装

下载编译好的二进制文件

#下载
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
#解压
unzip Trimmomatic-0.39.zip

解压后的文件包括以下几个：

Trimmomatic-0.39
├── adapters
│ ├── NexteraPE-PE.fa
│ ├── TruSeq2-PE.fa
│ ├── TruSeq2-SE.fa
│ ├── TruSeq3-PE-2.fa
│ ├── TruSeq3-PE.fa
│ └── TruSeq3-SE.fa
├── LICENSE
└── trimmomatic-0.39.jar

1 directory, 8 files

• adapters路径下的文件是ILLUMINA常规的接头文件，有单端和双端。
• trimmomatic-0.39.jar是编译之后的文件，需用java执行。

$ java -jar Trimmomatic-0.39/trimmomatic-0.39.jar -h
Usage: 
       PE [-version] [-threads ] [-phred33|-phred64] [-trimlog ] [-summary ] [-quiet] [-validatePairs] [-basein  |  ] [-baseout  |    ] ...
   or: 
       SE [-version] [-threads ] [-phred33|-phred64] [-trimlog ] [-summary ] [-quiet]   ...
   or: 
       -version

2、conda安装

conda install trimmomatic

三、Trimmomatic使用

1、SE模式

对于单端数据，指定一个输入文件和一个输出文件。所需的处理步骤(裁剪、裁剪、适配器裁剪等)作为附加参数指定在输入/输出文件之后。

java -jar  SE [-threads ] [-phred33 | -phred64] [-trimlog
]  
  ...

指定trimlog文件将创建一个记录所有reads 的日志，指出以下细节:

• read的名字，即QNAM
• 保留下的序列长度
• 从起始位置，第一个保留碱基的位置，aka。
• 原始read中最后一个碱基的位置
• 从末尾切除掉的碱基个数

可以根据需要指定多个步骤，方法是在部分处理步骤中描述的在最后使用附加参数。

2、PE模式

对于Pair-End数据，需要两个输入文件，会输出4个文件。其中两个文件对应于"paired"数据，即read1和read2都保留。另外两个对应于"unpaired"数据，在处理的过程中会过滤掉其中一端的reads。

java -jar  PE [-threads ] >] [-basein  |  ] [-baseout  |     ...

• -phred33|-phred64：指定碱基质量编码规则。在0.32版本以后如果缺省，则会根据数据自动选择。默认值是-phred64。
• -threads：表示要使用的线程数，用于提高多核计算机的性能。

## 示例
java -jar trimmomatic-0.30.jar PE \
	# 输入文件
	Sample_R1.fq.gz Sample_R2.fq.gz \
	# 输出文件
	Sample_Trim_paired_R1.fq.gz Sample_Trim_unpaired_R1.fq.gz \
	Sample_Trim_paired_R2.fq.gz Sample_Trim_unpaired_R2.fq.gz \
	# 去接头参数
	ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
	# 切除TruSeq3-PE中提供的Illumina适配器。最初Trimmomatic将寻找种子匹配(16个碱基)，最多允许2个不匹配。如果在配对端读取的情况下达到30分(约50个碱基，50*0.6)，或在单端读取的情况下达到10分(约17个碱基, 17*0.6)，这些"seed"序列将被修剪。
	LEADING:3 \
	# 删除前低质量碱基(低于质量3)
	TRAILING:3 \
	# 删除后低质量碱基(低于质量3)
	SLIDINGWINDOW:4:15 \
	# 扫描4个碱基宽滑动窗口，当每个碱基的平均质量下降到15以下时进行剪切
	MINLEN:36
	# 上述步骤完成后，删除小于36个碱基的reads
	

四、详细的处理细节

1、ILUMINACLIP

这一步用于查找和切除Illumina adapters。理论上虽然adapters和其他实验产生的序列可能位于在reads的任何位置，到目前为止，adapters污染最常见的原因是测序的DNA片段短于读取长度，即segment < read_length。在这种情况下，reads的开头包含有效的数据，但是reads的末尾会读取到到adapters。这将导致在reads的3'端出现部分或全部adapters序列。虽然识别完整的adapters序列比较容易，但识别短的部分adapters序列本质上是困难的。

有趣的是，在Pair-End的数据中，由于reads完全从两端开始，前向和反向reads的非adapters部分将是反向互补。对于Illumina的Pair-End的数据，Trimmomatic包含了第二种adapters识别策略，专门用于adapters测通的情况，并利用Pair-End数据中的证据。这种策略被称为回文(palindrome)模式。

下面的图说明了两种不同的策略。

• A：reads中包含读数据中的整个技术产生的序列，因此标准比对方法足以准确识别并切除。
• B：只有部分的技术序列包含在reads的3'端，因此只能使用短的对齐。这种情况会平衡假阳性和假阴性来进行取舍，对于一些长度低于阈值的的adapters序列片段可能会残留下来。
• C：回文模式检测长的reads。在本例中，根本没有任何有用的片段，并且两条reads都从adapters的序列开始。
• D：使用回文方法来，在reads的3'端有一个非目标序列。然而，由于利用了在配对端模式下可用的额外数据，作为校准的一部分进行比较的区域要长得多。不仅一次比较了两个adapters序列，而且还检查每条reads序列。因此，这种对齐比B中的短对齐要可靠得多，可以去除仅剩有1个碱基的adapters序列。

Trimmomatic使用两步来识别adapters和reads。首先，在reads的每个可能位置测试每个adapters的短段(最大16bp)。这种匹配成为"seed",由seedMismatch参数调控。如果"seed"接近于完全匹配，则adapters和reads之间将计算一个full alignment score。这种两步策略带来了相当大的效率增益，因为"seed"计算较快，同时减少full alignment score的计算。

full alignment score计算方法如下。每匹配一个碱基，评分增加0.6，而每错配一次，比对评分减少Q/10。通过考虑碱基的的质量，读取错误导致的miss-match的影响较小。一个完美匹配的12个碱基序列的score在7个以上，而需要25个碱基才能得到15。因此，我们建议在7 - 15之间的值作为简单对齐模式的阈值。对于回文匹配，可以采用更长的对齐方式，如上所述。因此，在30的范围内这个阈值可以更高。即使这个阈值非常高(需要匹配近50个碱基)，Trimmomatic仍然能够识别非常、非常短的adapters片段。

使用参数

ILLUMINACLIP: :::

ILLUMINACLIP:::::

ILLUMINACLIP::::::

• fastaWithAdaptersEtc：指定包含所有adapters、PCR序列等fasta文件的路径。
• seedMismatches：指定"seed"匹配允许的最大miss-match。
• palindromeClipThreshold：指定PE回文读取对齐时，两个“adapter ligated”读取之间的full alignment score阈值。
• simpleClipThreshold:指定任何adapter序列之间的full alignment score阈值。
• minAdapterLength:除了对齐评分外，回文模式还可以验证已检测到adapter的最小长度。如果未指定，则由于历史原因，默认为8个碱基。然而，由于回文模式有一个非常低的假阳性率，这可以安全地减少，甚至减少到1，以允许更短的adapter片段被删除。
• keepBothReads:在通过回文模式检测到读透后，删除适配器序列，反向读取包含与正向读取相同的序列信息，尽管是反向补码。因此，默认行为是完全删除反向读取。通过为该参数指定true，反向读取也将被保留，这可能会很有用，例如，如果下游工具不能处理成对和非配对读取的组合。一般设置为true。

2、SLIDINGWINDOW

执行滑动窗口修剪，一旦窗口内的平均质量低于阈值就进行切割。通过考虑多个碱基，避免因为一个低质量的碱基导致删除后面的reads中高质量的数据。

使用参数

SLIDINGWINDOW::

• windowSize：滑动窗口的大小
• requiredQuality：窗口内的平均质量。

3、MAXINFO

执行自适应的高质量调整，平衡保持较长reads的好处与保留错误的基础的成本。对于许多应用程序，一次读取的值是三个因素之间的平衡

使用参数

MAXINFO::

• targetLength: This specifies the read length which is likely to allow the location of the read within the target sequence to be determined.
• strictness: 这个值应该设置在0到1之间，它指定尽可能多地保留读长度与删除不正确的碱基之间的平衡。该参数的低值(<0.2)有利于读取较长的数据，而高值(>0.8)有利于读取的正确性。

4、LEADING

从起始位置移除低质量的碱基。当碱基的值低于这个阈值，这个碱基就会被移除，然后研究下一个碱基。

使用参数

LEADING:

• quality:最低质量。

5、TRAILING

从末端移除低质量的碱基。只要一个碱基的值低于这个阈值，这个碱基就会被移除，并且从下一个碱基开始处理(由于从3 '端开始，下一个碱基其实是在刚刚移除的碱基之前)。这种方法可以用来去除特殊的illumina低质量段区域(质量分数为2)，但我们推荐滑动窗口或MaxInfo代替。

使用参数

TRAILING:

• quality:最低质量。

6、CROP

无论读取的质量如何，都将从读取的末尾移除碱基，以便在执行此步骤后读取最大限度地达到指定的长度。在CROP之后执行的步骤可能会进一步缩短reads。

使用参数

CROP:

• length:从reads开始开始所要保留的碱基数。

7、HEADCROP

无论质量如何，从读取开始移除指定数量的碱基。

使用参数

HEADCROP:

• length:从reads起始所要移除的碱基数。

8、MINLEN

此模块删除低于指定最小长度的reads。它通常在所有其他处理步骤之后。这一步删除的读数将被计数，并包括在trimmomatic摘要中。

使用参数

MINLEN:

• length:指定要保留的读取的最小长度。