📚 课程目录第4章数据预处理

数据预处理：Trimmomatic/Cutadapt去接头、过滤低质量序列

数据预处理这步，说白了就是给原始测序数据“洗澡”。

我记得刚开始做RNA-seq时，总觉得这步可有可无。直到有一次，我拿到的数据接头残留高达30%，差异基因列表里全是假阳性。那之后，我再也不敢跳过预处理了。

我们手头的FASTQ文件，除了真正的RNA片段，还混着接头序列、低质量碱基，甚至还有测序仪器的“噪音”。这些东西不清理干净，后面的比对、定量全是白搭。

核心任务就两个

去接头：把测序时遗留的接头序列切掉。
过滤低质量：把质量太差的读段和碱基扔掉。

目前主流的工具就俩：Trimmomatic 和 Cutadapt。我个人习惯是双管齐下——Trimmomatic干粗活，Cutadapt补细活。

我的一点体会： 如果样本质量不错（Q30 > 85%），Trimmomatic的默认参数就能搞定。要是遇到问题样本，比如FFPE组织来源的数据，我会先用FastQC看看接头污染和低质量碱基的位置，再单独调Cutadapt的精准裁切。

Trimmomatic：批量化处理的首选

Trimmomatic用Java写的，支持双端数据（PE）和单端数据（SE）。它的核心优势是滑窗质量扫描，说白了就是滑动窗口算平均质量，低于阈值就咔嚓掉。

常用参数解读

参数	含义	我常用的值
ILLUMINACLIP	去接头，需要提供接头序列文件	`adapters.fa:2:30:10`
LEADING	去掉开头低质量碱基	`3`
TRAILING	去掉结尾低质量碱基	`3`
SLIDINGWINDOW	滑窗质量扫描（窗口大小:平均质量阈值）	`4:15`
MINLEN	保留的最短读段长度	`36`

一个典型的双端处理命令长这样：

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

避坑指南： 我曾经在ILLUMINACLIP参数上栽过跟头。那个 2:30:10 三个数字分别是：最大错配数（2）、比对分值阈值（30）、最短匹配长度（10）。如果接头序列给的不对，或者阈值设太严，接头根本去不干净。我会在第一次跑完后用FastQC检查一下，如果20-30bp处还有峰值，那就是接头残留，得调参数。

Cutadapt：精准活检手术刀

Cutadapt的强项是灵活。它支持正则表达式，可以处理各种奇葩接头。比如RACE-seq、siRNA这种非标准文库，Trimmomatic搞不定的时候，Cutadapt就能上。

它的基本用法很直白：

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
         -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
         -o sample_R1_clean.fastq.gz \
         -p sample_R2_clean.fastq.gz \
         sample_R1.fastq.gz sample_R2.fastq.gz

这里 -a 指定R1的接头，-A 指定R2的接头。实际项目中，我不会每次都手敲接头序列。我习惯把要用的接头写在FASTA文件里，用 -a file:adapters.fa 读入。这样不同文库类型换文件就行，不容易出错。

Cutadapt还有一个很实用的功能：测序质量过滤。用 -q 参数就能直接干：

cutadapt -q 20 -m 36 -o clean.fastq.gz raw.fastq.gz

这个命令会直接把质量低于20的碱基从两端切掉，丢弃长度不足36bp的读段。一气呵成。

选哪个？我的一点心得

场景	推荐工具	原因
标准Illumina PE数据，量大	Trimmomatic	速度快，批量处理方便，滑窗机制靠谱
非标准文库、特殊接头、单端数据	Cutadapt	灵活，支持正则，参数精细可控
需要极低接头残留（如单细胞RNA-seq）	先用Cutadapt再用Trimmomatic	两道关卡，基本把接头和低质量一网打尽
数据质量极差（FFPE、降解样本）	Trimmomatic + 降低滑窗阈值	滑窗机制能保住更多可用数据，但必须配合读段长度分布检查

统计一下： 根据我处理过的100多个RNA-seq项目经验，Trimmomatic和Cutadapt在标准情况下处理结果差异不到5%。真正的差异在于参数调优。我会把两个工具都装在环境里，遇到奇怪的数据就各跑一次对比。

过滤后的评估

处理完了不是直接去比对。我会跑一轮FastQC看看效果，重点关注：

Per base sequence quality：箱线图是否都在绿色区域？
Adapter Content：各位置接头比例是否降到1%以下？
Sequence Length Distribution：长度是否集中在预期区域？有没有意外短于36bp的残留？
Overrepresented sequences：是否还有高丰度未知序列？

如果这三项都过了关，说明预处理是成功的。

参数选择的经验之谈

参数怎么选？我自己的策略是：

接头文件：首先找我做实验时用的接头序列。找不到就用Trimmomatic自带的TruSeq3-PE.fa。如果数据明显有接头但默认文件不对，用FastQC提取过表达序列猜接头序列。
ILLUMINACLIP的三个数字：seedMismatches:palindromeClipThreshold:simpleClipThreshold。我个人从不改种子错配，保持2。后两个值标准库用30:10，长片段文库（>300bp）我会加严到30:15。
滑窗参数：4:15是公认的“甜点”。质量非常好的数据用4:20更干净；FFPE数据我会降到4:10，适当多保长度。
MINLEN：双端数据设36 (8bp 条形码 + 20bp 足够比对的长度)。小RNA文库视情况设18–22。

我的习惯： 每次处理完，我会保留一份过滤日志。Trimmomatic的输出日志里记录了配对保留率、丢弃率。这些数字能帮我在项目总结时快速评估样本整体质量。

好了，预处理这块基本就是这些。记住：数据预处理是RNA-seq分析的第一道坎。这步没做好，后面跑再牛的算法都是白费劲。希望这些经验能帮你少走弯路。