第4章 数据预处理：Trimmomatic进行adapter去除和低质量过滤

拿到FASTQ数据之后，很少有人会直接拿去做比对。为什么？你想想看，测序仪出来的数据，多多少少带着接头残留、低质量碱基，还有那些说不清道不明的测序错误。直接拿这些“毛坯房”去建分析流程，结果基本就是一团浆糊。

我个人习惯，拿到数据的第一件事就是跑一遍Trimmomatic。这工具我用了快十年，它虽然界面不花哨，但干活踏实。说白了，它就像一把瑞士军刀——切接头、裁低质量、筛短片段，一把搞定。

4.1 为什么要做数据预处理？

我见过不少新手，觉得预处理就是“多此一举”。直到有一次，我帮一个团队排查差异基因分析结果异常的问题，查了两天，最终发现是接头残留导致大量read比对不上参考基因组。从那以后，我再也不敢跳过这步了。

数据预处理的核心目标其实就三个：

• 去掉接头序列——测序片段比文库长度短时，测序仪会把接头的序列也读出来
• 过滤低质量碱基——3'端质量往往越往后越差，需要切除
• 去除过短的read——太短的序列比对时容易产生假阳性

4.2 Trimmomatic简介

Trimmomatic由Usadellab实验室开发，支持SE（单端）和PE（双端）模式。它基于Java，所以跑之前要先确认JVM装好了。

它的运行模式很有意思——你可以把各种处理步骤按照顺序拼成一个“管道”。每一步处理完，结果传给下一步。我刚开始用的时候总觉得参数太多记不住，后来用的久了，发现常用的也就那么几个。

# 基本检查：验证Java环境
java -version

# 下载Trimmomatic（以v0.39为例）
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip

注意Trimmomatic需要两个文件：可执行JAR包和接头序列文件。接头文件在adapters/目录下，根据你的建库试剂盒挑选对应的。我个人习惯把常用接头文件统一整理到一个目录，方便随时调用。

4.3 核心参数解读

Trimmomatic的参数很多，但真正核心的就这5个。我列一张表，方便你对照：

参数	功能	常用值	避坑点
ILLUMINACLIP	去除接头和引物	2:30:10	seed mismatch过高会漏接头
SLIDINGWINDOW	滑动窗口质量过滤	4:20	窗口大小越小，过滤越激进
LEADING	切除5'端低质量碱基	3	一般不需要改
TRAILING	切除3'端低质量碱基	3	和SLIDINGWINDOW有重叠
MINLEN	保留的最小read长度	36或50	设太高会丢掉太多数据

4.4 运行Trimmomatic：双端示例

PE模式是分析中最常见的。嗯，这里要注意——PE模式的输出文件有4个（两对paired + 两个unpaired），很多人第一次跑的时候会懵。

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  sample_R1_paired.fastq.gz sample_R1_unpaired.fastq.gz \
  sample_R2_paired.fastq.gz sample_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:20 \
  MINLEN:36

为什么会输出unpaired文件？因为双端测序中，一条read被过滤掉了，它的配对标对另一条read就变成“单身”了。这些单身read虽然不能做双端比对，但如果单独用，仍然包含可用信息。我曾经在分析低质量样本时靠这些unpaired read挽回了7%的有效数据量。

4.5 单端模式

SE模式简单得多，只有两个输出文件（一个clean，一个trash）。参数逻辑和PE一致，但不需要考虑配对关系。

java -jar trimmomatic-0.39.jar SE \
  sample.fastq.gz \
  sample_clean.fastq.gz \
  ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:20 \
  MINLEN:36

4.6 参数调优经验

说实话，没有放之四海皆准的参数组合。我一般会根据数据质量先跑一个小样本试一遍，看日志输出的SURVIVING比例来决定怎么调。

举个真实案例。有个RNA-seq数据，我按默认参数跑完，发现SURVIVING只有58%。检查日志看到大量read被MINLEN:36过滤。我怀疑是接头去除不干净，导致read被切得太短。后来我把ILLUMINACLIP的第二个参数从30改成25（降低比对的严谨度），SURVIVING直接提升到79%。

4.7 检查过滤结果

跑完之后，别急着往下走。Trimmomatic会在终端输出摘要信息，记得截图或重定向到文件。我习惯这样收集日志：

java -jar trimmomatic-0.39.jar PE ... 2>&1 | tee trimmomatic_sample.log

典型的日志摘要长这样：

Input Read Pairs: 1000000
Both Surviving: 850000 (85.00%)
Forward Only Surviving: 30000 (3.00%)
Reverse Only Surviving: 25000 (2.50%)
Dropped: 95000 (9.50%)

如果Both Surviving低于70%，我个人建议先排查问题，而不是直接进入下游分析。

4.8 常见陷阱

4.9 批量处理脚本

跑单样本测试没问题后，就该批量跑了。我习惯写一个简单的Shell循环：

for r1 in *R1.fastq.gz; do
  sample=${r1%_R1.fastq.gz}
  r2=${sample}_R2.fastq.gz
  java -jar trimmomatic-0.39.jar PE \
    $r1 $r2 \
    ${sample}_R1_paired.fastq.gz ${sample}_R1_unpaired.fastq.gz \
    ${sample}_R2_paired.fastq.gz ${sample}_R2_unpaired.fastq.gz \
    ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:20 MINLEN:36
done

这段代码会处理当前目录下所有*R1.fastq.gz文件，并对应对找到对应的R2文件。我个人觉得这个结构最简单，改起来也方便。

4.10 本章流程总览

下面这张图概括了Trimmomatic预处理的完整数据流向，从原始FASTQ到过滤后的输出，每个环节我都标了主要的参数作用点。你可以把它贴在工位旁做参考：

4.11 核心参数速查表

最后，我把ILLUMINACLIP的三个数字含义拆开讲清楚，因为这是最容易搞混的部分：

参数位置	名称	解释	推荐值
2:30:10 中的 2	seed mismatches	比对接头seed时允许的错配数	2（太大会漏，太小会误切）
2:30:10 中的 30	palindrome clip threshold	PE模式下双端read重叠匹配得分阈值	30（质量好的数据可降到20）
2:30:10 中的 10	simple clip threshold	单端模式下接头比对的最低得分	10（和palindrome一起调）