测序数据质量控制：FastQC与Trimmomatic实战
你拿到一个原始FASTQ文件，第一件事是什么？
我个人习惯是：先跑FastQC看看数据到底什么成色。好比体检，在抽血化验前总得先拍个片子。
你想想看，如果原始数据里接头没切干净，碱基质量在末端断崖式下跌，那下游的比对和变异检测会被拉下水。我见过太多案例——变异调不出来，最后排查了一圈，发现是3端质量没清理利索。
这一章，我们就来拆解FastQC的质量报告，再实战Trimmomatic的常用参数组合。
先跑一遍FastQC
拿到原始数据后，我的做法是用FastQC扫一圈。命令很简单，但效果很关键：
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc_reports/ -t 4-t 4 开4线程，对上百兆的FASTQ来说够用了。
结果会生成两个文件：一个html报告页，一个zip压缩包。我习惯直接看html，方便快速翻页。
FastQC报告，到底在看什么？
报告里有十多项模块。嗯，这里注意，不是所有图标变红就代表数据报废，要结合测序平台和文库类型判断。我按重要性排个序：
| 模块名称 | 看什么 | 我关注的点 |
|---|---|---|
| Per base sequence quality | 每个碱基位置的Q值箱线图 | 前10bp和末尾10bp波动是否剧烈 |
| Per sequence quality scores | 平均Q值的分布 | 平均值低于Q20的比例是否>20% |
| Adapter Content | 接头序列富集情况 | 是否有任何接头出现在>5%的reads中 |
| Overrepresented sequences | 过度代表序列 | 可能是指标或接头或重复区域 |
| Per base N content | N碱基比率 | 如果在某区域突然升高，说明测序仪没读出来 |
| Per base GC content | GC含量分布 | 偏离理论分布可能提示污染 |
我的判断逻辑：
- 如果只有末尾几个框变红（Q<20），属于正常衰减，Trimmomatic切掉即可。
- 如果整个reads质量都低（平均值
- Adapter Content项只要有黄/红色报警，必须去除接头。
FastQC报告解读实例
我记得有一次，一个WGS项目的FastQC报告里，Per base sequence quality图上绿色条一直在Q25以下摇摆。我后来发现，是文库制备时扩增循环数偏少。像这种系统性偏低，单纯靠Trimmomatic切是救不回来的。
再比如，Adapter Content项，如果R1的3端有Illumina Universal Adapter富集，R2的5端也可能有。这时你就得用TrimmomaticPE把两端的接头都清洗掉。
我曾经踩过的坑：
有一次我对着FastQC报告里的Overrepresented sequences列表发愁，里面显示一个序列占了73%。我以为是接头残留花了半天调Trimmomatic参数，后来发现这居然是目标区域的扩增产物。所以在跑Trimmomatic之前，先看看这个序列的注释，确认是不是真正的外源指标。
Trimmomatic参数实战
Trimmomatic是我处理NGS数据时最常用的工具。我个人习惯用TrimmomaticPE（双端），它比SE模式多了一些配对信息校验。
下面是我常用的一个参数组合，以及每个参数的含义：
java -jar trimmomatic.jar PE \
sample_R1.fastq.gz sample_R2.fastq.gz \
sample_R1_clean.fastq.gz sample_R1_unpaired.fastq.gz \
sample_R2_clean.fastq.gz sample_R2_unpaired.fastq.gz \
ILLUMINACLIP:adapters.fa:2:30:10 \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:36我们来逐个拆解：
| 参数 | 作用 | 我常用的值 |
|---|---|---|
| ILLUMINACLIP | 去除接头与索引序列 | 2:30:10（容错2/30分匹配/10分palindrome模式） |
| LEADING | 从5端切掉Q值低于阈值的碱基 | 3（低于Q3就切） |
| TRAILING | 从3端切掉Q值低于阈值的碱基 | 3 |
| SLIDINGWINDOW | 4碱基窗口内平均Q低于阈值时切除 | 4:15 |
| MINLEN | 过滤掉切完后长度不足的reads | 36（双端150bp我常设36） |
我的建议：
LEADING和TRAILING的阈值别设太高，Q3就已经很宽松。我一般用3，把N碱基切掉就行，因为有些文库在两端本来就是引物段，质量偏低是正常的。如果切太狠，后续的比对率反而会下降。
SLIDINGWINDOW的窗口大小，我个人喜欢用4，不要用3，3窗口太敏感。平均Q值阈值15是个平衡点——切掉Q15以下的区域，但保留Q15以上的信息，不至于把mapped read切没了。
TrimmomaticPE的适配器文件
适配器文件你得根据文库来选。Illumina的通用接头长这样：
>PrefixPE/1
AGATCGGAAGAGCACACGTCTGAACTCCAGTCA
>PrefixPE/2
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT但是，如果你用了Nextera XT文库，接头序列会不一样。我习惯直接去fastqc生成的zip包里提取适配器序列：解压后找到一个叫 fastqc_data.txt 的文件，里面Overrepresented sequences那部分会列出具体序列。
避坑关键：
ILLUMINACLIP的第三个参数（palindrome模式匹配分数阈值），我一般从10开始调。如果接头残留还多，就降到8或6。别一次降太多，否则会把真实reads当成接头切掉。我遇到过一次这样的情况：调成6后，比对率骤降5%，后来发现是Trimmomatic把一部分短插入片段的重复区域当成了palindrome接头切掉了。
去接头与低质量碱基的实战脚本
下面是我给一个项目写的处理脚本，跑完后我看一眼clean数据的FastQC报告，确认效果：
#!/bin/bash
# 双端数据质量控制
FASTQ_DIR=./raw
CLEAN_DIR=./clean
ADAPTER=./adapters/truseq3.fa
for sample in sampleA sampleB sampleC; do
java -jar trimmomatic-0.39.jar PE \
${FASTQ_DIR}/${sample}_R1.fastq.gz \
${FASTQ_DIR}/${sample}_R2.fastq.gz \
${CLEAN_DIR}/${sample}_R1_clean.fastq.gz \
${CLEAN_DIR}/${sample}_R1_unpaired.fastq.gz \
${CLEAN_DIR}/${sample}_R2_clean.fastq.gz \
${CLEAN_DIR}/${sample}_R2_unpaired.fastq.gz \
ILLUMINACLIP:${ADAPTER}:2:30:10 \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:15 MINLEN:36
# 跑clean数据的FastQC检查一下
fastqc ${CLEAN_DIR}/${sample}_R1_clean.fastq.gz -o ${CLEAN_DIR}/qc/
fastqc ${CLEAN_DIR}/${sample}_R2_clean.fastq.gz -o ${CLEAN_DIR}/qc/
done跑完之后，我会去clean目录下的FastQC报告里看一眼Per base sequence quality图和Adapter Content项。如果Adapter Content变绿了，说明接头去干净了。如果Per base sequence quality图前10bp仍然偏低，可能是LEADING设得太高了，我会调回5或直接去掉LEADING。
质量控制的决策流程
下面这张流程图总结了整个质控分析的决策路径。从原始数据出发，先过FastQC，根据报告结果选择不同的Trimmomatic参数组合，最后再验证Clean数据质量：
流程图里那个从Clean FASTQ往回指向FastQC的虚线箭头，是我后来加的。因为每次调完参数后我都习惯重新跑一次FastQC确认接头是否完全去除。你想想看，如果只切了一轮发现还有残留，那第二轮就得降ILLUMINACLIP的匹配阈值，或者加上一个额外的Adapter序列。
写在质量控制的最后
质量控制这件事，说白了就是打磨数据的信噪比。FastQC告诉你信号在哪、噪声在哪，Trimmomatic帮你去掉你不想要的。但质控不是越狠越好——切太大，valid reads不够，变异检测的灵敏度就崩了；切太少，残留的高质量区域被接头污染，假阳性就上来了。
我个人习惯的基线组合就是 ILLUMINACLIP:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36。然后根据FastQC报告的反馈微调。