如何使用Trimmomatic对NGS数据进行质量过滤

引言

随着高通量测序技术（Next-Generation Sequencing, NGS）的快速发展，生物信息学领域面临着海量数据的处理挑战。NGS数据通常包含大量的短读长（short reads），这些读长在测序过程中可能会受到各种因素的影响，导致数据质量下降。为了提高后续分析的准确性，对NGS数据进行质量过滤是必不可少的一步。Trimmomatic是一款广泛使用的工具，专门用于对NGS数据进行质量控制和过滤。本文将详细介绍如何使用Trimmomatic对NGS数据进行质量过滤。

Trimmomatic简介

Trimmomatic是一款由Java编写的开源工具，专门用于处理Illumina平台的NGS数据。它能够对原始测序数据进行多种质量控制操作，包括去除低质量碱基、去除接头序列、裁剪读长等。Trimmomatic支持单端和双端测序数据的处理，并且具有较高的灵活性和效率。

主要功能

• 去除低质量碱基：根据用户设定的质量阈值，去除读长中质量较低的碱基。
• 去除接头序列：通过比对已知的接头序列，去除读长中的接头污染。
• 裁剪读长：根据用户设定的长度阈值，裁剪读长的起始或末尾部分。
• 滑动窗口过滤：通过滑动窗口的方式，动态评估读长的质量，并进行相应的裁剪。
• 保留或丢弃读长：根据过滤结果，保留或丢弃不符合质量要求的读长。

安装Trimmomatic

在使用Trimmomatic之前，首先需要确保系统上已经安装了Java运行环境（JRE）。Trimmomatic可以通过以下步骤进行安装：

1. 下载Trimmomatic：访问Trimmomatic的官方网站（http://www.usadellab.org/cms/?page=trimmomatic）下载最新版本的Trimmomatic。
2. 解压缩：将下载的压缩包解压到指定目录。
3. 设置环境变量：为了方便使用，可以将Trimmomatic的路径添加到系统的环境变量中。

# 假设Trimmomatic解压到/opt/Trimmomatic目录
export TRIMMOMATIC_HOME=/opt/Trimmomatic
export PATH=$PATH:$TRIMMOMATIC_HOME

使用Trimmomatic进行质量过滤

基本命令格式

Trimmomatic的基本命令格式如下：

java -jar trimmomatic.jar [PE/SE] [-threads <threads>] [-phred33/-phred64] [-trimlog <logfile>] <input1> <input2> <output1_paired> <output1_unpaired> <output2_paired> <output2_unpaired> <options>

• PE/SE：指定输入数据的类型，PE表示双端测序数据，SE表示单端测序数据。
• -threads：指定使用的线程数，默认为1。
• -phred33/-phred64：指定质量分数的编码方式，phred33是Illumina 1.8+版本的标准，phred64是早期版本的标准。
• -trimlog：指定日志文件的路径。
• input1/input2：输入文件的路径，对于双端测序数据，input1和input2分别对应两个读长文件。
• output1_paired/output1_unpaired：输出文件的路径，paired表示保留的成对读长，unpaired表示未成对的读长。
• options：具体的过滤选项。

常用过滤选项

Trimmomatic提供了多种过滤选项，以下是一些常用的选项：

• ILLUMINACLIP：去除接头序列。
• SLIDINGWINDOW：滑动窗口过滤。
• LEADING：去除读长起始部分的低质量碱基。
• TRLING：去除读长末尾部分的低质量碱基。
• CROP：裁剪读长到指定长度。
• HEADCROP：去除读长起始部分的指定长度。
• MINLEN：保留长度大于指定值的读长。

示例：双端测序数据的质量过滤

假设我们有一对双端测序数据文件input_R1.fastq和input_R2.fastq，我们希望使用Trimmomatic对其进行质量过滤。具体步骤如下：

1. 去除接头序列：使用ILLUMINACLIP选项去除接头序列。
2. 滑动窗口过滤：使用SLIDINGWINDOW选项进行滑动窗口过滤。
3. 去除起始和末尾的低质量碱基：使用LEADING和TRLING选项去除起始和末尾的低质量碱基。
4. 保留长度大于50的读长：使用MINLEN选项保留长度大于50的读长。

具体的命令如下：

java -jar trimmomatic.jar PE -threads 4 -phred33 \
  input_R1.fastq input_R2.fastq \
  output_R1_paired.fastq output_R1_unpaired.fastq \
  output_R2_paired.fastq output_R2_unpaired.fastq \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  SLIDINGWINDOW:4:20 \
  LEADING:3 \
  TRLING:3 \
  MINLEN:50

参数解释

• ILLUMINACLIP:TruSeq3-PE.fa:2:30:10：使用TruSeq3-PE.fa文件中的接头序列进行过滤，允许最多2个错配，接头序列的最小匹配分数为30，接头序列的最小匹配长度为10。
• SLIDINGWINDOW:4:20：使用滑动窗口过滤，窗口大小为4，平均质量分数阈值为20。
• LEADING:3：去除读长起始部分质量分数低于3的碱基。
• TRLING:3：去除读长末尾部分质量分数低于3的碱基。
• MINLEN:50：保留长度大于50的读长。

示例：单端测序数据的质量过滤

对于单端测序数据，命令格式与双端测序数据类似，只是不需要指定第二个输入文件和输出文件。假设我们有一个单端测序数据文件input.fastq，我们希望使用Trimmomatic对其进行质量过滤。具体步骤如下：

1. 去除接头序列：使用ILLUMINACLIP选项去除接头序列。
2. 滑动窗口过滤：使用SLIDINGWINDOW选项进行滑动窗口过滤。
3. 去除起始和末尾的低质量碱基：使用LEADING和TRLING选项去除起始和末尾的低质量碱基。
4. 保留长度大于50的读长：使用MINLEN选项保留长度大于50的读长。

具体的命令如下：

java -jar trimmomatic.jar SE -threads 4 -phred33 \
  input.fastq \
  output_paired.fastq \
  ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 \
  SLIDINGWINDOW:4:20 \
  LEADING:3 \
  TRLING:3 \
  MINLEN:50

参数解释

• ILLUMINACLIP:TruSeq3-SE.fa:2:30:10：使用TruSeq3-SE.fa文件中的接头序列进行过滤，允许最多2个错配，接头序列的最小匹配分数为30，接头序列的最小匹配长度为10。
• SLIDINGWINDOW:4:20：使用滑动窗口过滤，窗口大小为4，平均质量分数阈值为20。
• LEADING:3：去除读长起始部分质量分数低于3的碱基。
• TRLING:3：去除读长末尾部分质量分数低于3的碱基。
• MINLEN:50：保留长度大于50的读长。

结果分析

Trimmomatic运行完成后，会生成多个输出文件。对于双端测序数据，通常会生成四个文件：

• output_R1_paired.fastq：保留的成对读长文件，对应第一个读长文件。
• output_R1_unpaired.fastq：未成对的读长文件，对应第一个读长文件。
• output_R2_paired.fastq：保留的成对读长文件，对应第二个读长文件。
• output_R2_unpaired.fastq：未成对的读长文件，对应第二个读长文件。

对于单端测序数据，通常会生成一个文件：

• output_paired.fastq：保留的读长文件。

用户可以根据需要选择使用这些文件进行后续分析。通常情况下，成对的读长文件用于后续的比对和组装分析，而未成对的读长文件可以用于其他分析或直接丢弃。

总结

Trimmomatic是一款功能强大且灵活的工具，能够有效地对NGS数据进行质量过滤。通过去除低质量碱基、接头序列和裁剪读长等操作，Trimmomatic能够显著提高后续分析的准确性。本文详细介绍了Trimmomatic的安装、基本命令格式、常用过滤选项以及双端和单端测序数据的质量过滤示例。希望本文能够帮助读者更好地理解和使用Trimmomatic进行NGS数据的质量过滤。

参考文献

• Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 30(15), 2114-2120.
• Trimmomatic官方文档：http://www.usadellab.org/cms/?page=trimmomatic