数据清洗：从混乱的reads到干净的起点

数据清洗，说白了就是给原始测序数据“洗澡”。

你拿到的FASTQ文件里，不仅有目标DNA片段，还混着接头(adapter)、低质量碱基和过短的序列。这些杂质不清理掉，下游分析全是坑。我个人习惯在每次项目开始时都把数据清洗放在第一步，绝不跳过。

为什么要做数据清洗？

测序仪吐出来的数据，质量参差不齐。你想想看，一个典型的Illumina测序实验会产生几百万条reads，其中大概有5%-15%包含接头污染。如果直接拿去比对或组装，会有什么后果？

• 比对率降低 – 包含接头的read会匹配不到参考基因组上
• 假阳性SNP – 低质量碱基被错误识别为变异位点
• 组装断裂 – 短read或带接头的read打断contig连续性
• 定量偏差 – 在RNA-seq或ChIP-seq中，残留接头会导致错误计数

Trimmomatic：我的主力清洗工具

Trimmomatic是我用得最多的清洗工具。它把去接头、质量修剪和长度过滤整合在一次运行里，效率很高。我记得第一次用Trimmomatic时，被它的参数搞晕了，后来发现只要理解几个核心参数，剩下的都是套路。

安装和基本用法

Trimmomatic是Java程序，下载jar包就能跑。我习惯放在/usr/local/bin/或项目根目录下：

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
cd Trimmomatic-0.39

java -jar trimmomatic-0.39.jar PE \
  -phred33 \
  input_R1.fastq.gz input_R2.fastq.gz \
  output_R1_paired.fastq.gz output_R1_unpaired.fastq.gz \
  output_R2_paired.fastq.gz output_R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

PE表示双端测序，单端用SE。-phred33是质量值编码方式，现在绝大多数数据都用这个。

去接头：ILLUMINACLIP参数

去接头是清洗中最关键的一步。Trimmomatic通过ILLUMINACLIP参数组实现，格式为：

ILLUMINACLIP:接头文件:seedMismatches:palindromeClip:simpleClip:minAdapterLength:keepBothReads

参数	说明	我常用的值
接头文件	包含接头序列的FASTA文件	TruSeq3-PE.fa
seedMismatches	接头匹配允许的错配数	2
palindromeClip	回文模式匹配的错配数	30
simpleClip	简单模式匹配的错配数	10
minAdapterLength	最小接头匹配长度	2
keepBothReads	保留两条read的接头信息	keepBothReads

质量修剪：LEADING、TRAILING和SLIDINGWINDOW

低质量碱基会引入假阳性。Trimmomatic提供了三种修剪策略：

• LEADING:3 – 从read开头开始剪，遇到质量值低于3的碱基就切掉
• TRAILING:3 – 从read末尾往前剪，同样以质量值3为阈值
• SLIDINGWINDOW:4:15 – 用4bp窗口滑动，如果窗口内平均质量低于15，就从这个窗口位置切开

这里我踩过坑。以前我把SLIDINGWINDOW的窗口设成1:20，结果每条read被剪得稀碎，长度都不够用了。后来改成4:15，平衡多了。你想想看，窗口太小会让修剪过于激进，窗口太大会放过局部低质量区域。4bp窗口是个不错的折中。

长度过滤：MINLEN

最后一步，去掉太短的read。Trimmomatic通过MINLEN参数过滤，默认36bp。我个人建议这个值设成你后续分析的最小需求。比如做RNA-seq的差异表达，通常要求read长度大于50bp；做16S扩增子，可能100bp就够了。

# 常见用法
MINLEN:50   # 保留长度≥50bp的read，长度不够的直接丢掉

cutadapt：另一种去接头利器

cutadapt在很多场景下比Trimmomatic更灵活。它专注于去接头和质量修剪，语法也更直观。当你有多种接头类型需要处理时，我会优先选择cutadapt。

基本去接头命令

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
         -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
         -o output_R1.fastq.gz -p output_R2.fastq.gz \
         input_R1.fastq.gz input_R2.fastq.gz

• -a 指定R1的接头序列（3'端）
• -A 指定R2的接头序列（3'端）
• -o 和 -p 分别指定R1和R2的输出文件

质量修剪和长度过滤

cutadapt也可以做质量修剪和长度过滤，但需要搭配额外参数：

cutadapt -q 20 -m 50 \
         -a ADAPTER_R1 -A ADAPTER_R2 \
         -o clean_R1.fastq.gz -p clean_R2.fastq.gz \
         raw_R1.fastq.gz raw_R2.fastq.gz

-q 20表示从3'端修剪碱基，直到遇到质量值大于等于20的碱基。-m 50表示丢弃长度小于50bp的read。

两种工具对比：什么时候用哪个？

场景	推荐工具	理由
标准Illumina双端数据	Trimmomatic	参数集成度高，一次运行搞定所有
多种接头混合	cutadapt	支持多接头列表，匹配更灵活
单端数据	两者均可	cutadapt语法更简洁
需要精确控制接头匹配	cutadapt	支持锚定(anchored)和通配符匹配
批量处理大量样本	Trimmomatic	并行处理效率高，可写脚本一次跑完

数据清洗的完整流程（SVG图）

下面这张图展示了从原始FASTQ到干净数据的完整清洗流程。每个步骤都会输出统计信息，方便你检查清洗效果。

避坑指南：我踩过的那些坑

做数据清洗这么多年，有些问题几乎每批数据都会遇到。这里分享几个实用经验。

接头文件要选对

我曾经在一个项目中用了错误的接头文件。数据是Nextera文库制备，但我不小心用了TruSeq的接头。结果清洗后大部分read都被切掉了接头，实际是目标片段被误认为接头。检查一下：

# 查看数据中是否还有接头污染
# 取清洗后未配对的前100条read检查
zcat clean_R1_unpaired.fastq.gz | head -100 | grep "AGATCGGAAG"

如果还能看到大量接头序列，说明接头文件或参数不对。

质量值编码问题

老数据（比如2014年之前的Illumina数据）可能用-phred64编码。现在大部分数据是-phred33。混用会导致质量值解释错误，修剪结果会完全乱掉。我一般先用FastQC看一眼编码格式，再决定Trimmomatic的-phred参数。

双端数据配对问题

Trimmomatic处理双端数据时，会丢掉配对失败的单端read。这可能导致R1和R2的read数不一致。解决方法是始终使用_paired文件做下游分析，_unpaired文件在必要时再用。

清洗后的质量检查

清洗完成不等于万事大吉。我习惯再做一轮FastQC，和清洗前的报告对比：

• 接头污染：清洗后应该在FastQC的“Adapter Content”模块中几乎看不见接头序列
• 质量分布：清洗后的read质量箱线图应该整体右移，尾部变短
• read长度：长度分布图应该更集中，短尾巴被剪掉了
• read数：清洗前后read数比例正常，一般损失5%-20%属于合理范围

如果损失比例太高（比如超过40%），说明你的参数设得太严，或者原始数据质量确实很差。嗯，这时候需要回头调整参数，或者考虑是否要重新测序。