FASTQ预处理：接头去除、低质量碱基修剪、长度过滤

拿到原始FASTQ数据之后，第一件事不是急着跑比对。你得先看看数据质量，把该去掉的杂质去掉。我个人习惯把这一步叫做"数据洗澡"——洗得干净，后面分析才舒坦。

这一章我们专门聊Trimmomatic和cutadapt这两个工具。说实话，我入行那几年，这两个家伙帮我挡了不少坑。有一次项目赶得紧，我图省事跳过了预处理，结果跑出来的表达矩阵跟预期差了一大截。嗯，从那以后我再也不敢偷懒了。

为什么非做预处理不可？

你想想，测序仪出来的数据，其实混着不少"垃圾"：

• 接头序列：建库时人为加上的，不剪掉就会污染比对结果
• 低质量碱基：测序末尾质量直线下降，留着会引入错误变异
• 过短的读段：短到没法可靠比对，留着也是浪费计算资源

说白了，预处理就是个去伪存真的过程。我见过不少新人，觉得多组学分析高大上，结果卡在最基础的预处理上。我的建议是：把这一步做扎实，后面90%的问题都不会出现。

接头去除：cutadapt vs Trimmomatic

接头去除是预处理的重头戏。我们先来看看两个主流工具怎么选。

特性	cutadapt	Trimmomatic
输入格式	FASTQ, FASTA	FASTQ（压缩/不压缩）
接头检测方式	对齐匹配，支持模糊匹配	精确匹配 + 允许错配
双端数据支持	支持，但需手动配置	原生双端协同处理
运行速度	单线程较快	多线程并行
语法简洁性	简洁，一条命令搞定	参数较多，但功能全

我个人更推荐cutadapt做接头去除，原因很简单——它的模糊匹配算法在处理非标准接头时更灵活。当然，如果你已经知道接头序列是标准的Illumina接头，Trimmomatic也没问题。

cutadapt实战

这是我常用的单端数据命令：

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
         -o clean_sample1.fastq.gz \
         sample1.fastq.gz \
         --minimum-length 36 \
         --trim-n

参数说明：

• -a：指定接头序列（我这里写的是Illumina通用接头，你要根据建库试剂盒调整）
• -o：输出文件名
• --minimum-length：去除接头后，读段短于这个长度的直接丢弃
• --trim-n：修剪读段两端的N碱基

双端数据稍微复杂一点：

cutadapt -a ADAPTER_FWD -A ADAPTER_REV \
         -o R1_trimmed.fastq.gz -p R2_trimmed.fastq.gz \
         R1.fastq.gz R2.fastq.gz \
         --minimum-length 36 \
         --pair-filter=any

--pair-filter=any这个参数很有意思——只要一端读段被丢弃，另一端也被丢弃。我在做单细胞数据时就吃过这个亏。一开始没加这个参数，结果R1和R2数量对不上，后面合并的时候报错。嗯，亏吃多了就记住了。

Trimmomatic的接头去除方式

Trimmomatic用的是另一种逻辑。它不是直接在命令行写接头序列，而是用适配器文件：

java -jar trimmomatic-0.39.jar SE \
     sample1.fastq.gz clean_sample1.fastq.gz \
     ILLUMINACLIP:adapters.fa:2:30:10 \
     LEADING:3 TRAILING:3 \
     SLIDINGWINDOW:4:15 MINLEN:36

这里ILLUMINACLIP:adapters.fa:2:30:10的冒号参数含义是：

• 2：种子匹配允许的最大错配数
• 30：匹配得分阈值（越高越严格）
• 10：比对得分阈值（越高越容忍模糊匹配）

我在项目里遇到过一种情况：接头序列不太标准，Trimmomatic死活剪不干净。后来换成cutadapt的-a加--overlap=5参数，问题就解决了。所以我的建议是——两个工具都装好，遇到问题随时切换。

低质量碱基修剪：滑动窗口法

接头去掉了，接下来处理低质量碱基。这里我强烈推荐滑动窗口法，说白了就是从读段5'端开始，每次取一个固定大小的窗口，计算窗口内碱基的平均质量。如果低于阈值，就把整个窗口及之后的部分全部剪掉。

Trimmomatic的滑动窗口参数：

SLIDINGWINDOW:4:15

意思是用4个碱基的窗口，平均质量低于15时开始裁剪。这个15的阈值怎么定？我个人经验：

• 常规RNA-seq：15-20
• 单细胞：10-15（因为单细胞数据本身质量波动大，剪太狠会丢太多细胞）
• 长读长ONT数据：7-10（质量体系不同）

有一次我处理一个单细胞数据集，用默认的20做阈值，结果过滤掉了60%的 reads。后来降到12，保留率恢复到80%，而且下游聚类结果反而更好。所以说，参数不是固定的，要根据数据质量微调。

cutadapt没有内置滑动窗口，但可以配合质量修剪参数：

cutadapt -q 15,10 \
         -o output.fastq.gz \
         input.fastq.gz

-q 15,10的意思是：5'端质量低于15的剪掉，3'端质量低于10的剪掉。这种两端修剪的方式比滑动窗口简单，但效果也不错。

长度过滤：守住底线

剪完接头和质量后，你会发现有些读段变得很短——短到没法用了。长度过滤就是设置一个底线，太短的就扔掉。

Trimmomatic里用MINLEN:36：

java -jar trimmomatic-0.39.jar PE \
     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:adapters.fa:2:30:10 \
     LEADING:3 TRAILING:3 \
     SLIDINGWINDOW:4:15 \
     MINLEN:36

cutadapt里用--minimum-length 36：

cutadapt -a ADAPTER -o output.fastq.gz \
         input.fastq.gz --minimum-length 36

我建议的长度底线：

应用场景	建议最小长度
标准RNA-seq（75bp/100bp）	30-36 bp
单细胞10x（28+91）	20-25 bp
ChIP-seq（36-50bp）	20-25 bp
长读长ONT/PacBio	100-200 bp

完整的预处理流程：一张图说清楚

下面这张图展示了从原始FASTQ到干净FASTQ的完整预处理流程。我自己画的时候特意把关键参数标出来了，方便你对照。

双端数据处理的特别提醒

双端数据的预处理有一个非常容易踩的坑——成对性。处理过程中，如果R1被丢弃但R2保留，或者反过来，就会导致两条文件的行数不一致。后续做比对时，比对软件直接崩溃。

总结几个关键点：

• 始终保持成对：Trimmomatic的PE模式会自动处理，cutadapt需要加--pair-filter=any
• 检查输出行数：处理完后用wc -l检查R1和R2的行数
• 日志文件别删：Trimmomatic会输出处理摘要，里面有存活reads的比例，留着做质量控制报告