3. 数据预处理：Trimmomatic/Cutadapt去除接头和低质量碱基

拿到测序仪下机的FASTQ文件，你可别急着直接往比对软件里塞。原始数据里藏着接头（Adapter）、低质量碱基、还有那些让人头疼的N碱基。这一章咱们就专门处理这个——把数据洗干净。

我个人习惯先把数据质量用FastQC扫一眼，心里有个底。记得有一次，一个合作项目的样本RNA-seq数据Q30只有60%，我一看就知道是接头残留太多。切掉接头后Q30直接跳到92%——你想想看，这一步能省后面多少麻烦。

3.1 原始FASTQ里到底有什么“脏东西”？

说白了，测序读长往往比插入片段长，测到后面就变成了接头序列。还有，碱基质量从开头到末尾会逐渐下降，末尾那几十个碱基基本可以扔了。

• 接头残留：文库构建时引入的固定序列，比如Illumina的TruSeq、Nextera等。
• 低质量碱基：Q值低于20（错误率1%）或低于15（错误率3%）的碱基。
• N碱基：测序仪无法判断的碱基，长度多了会影响比对。
• 短片段：长度太短的reads，比对时容易多匹配。

3.2 工具选择：Trimmomatic vs Cutadapt

这两个工具在生信圈里几乎人手一份。Trimmomatic是Java写的，专为Illumina数据优化；Cutadapt是Python写的，支持更多平台（包括PacBio/ Nanopore的接头）。我一般这样选：

如果是标准Illumina双端数据，我强烈推荐Trimmomatic。但如果你的文库用了很偏门的接头，或者在做单细胞数据，Cutadapt的自动检测能力会更省心。

3.3 Trimmomatic 实战

先看看我最常用的一条命令：

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

参数说明一下：

• PE：双端模式。单端用SE。
• ILLUMINACLIP：接头文件、允许的最大错配数（2）、palindrome模式阈值（30）、简单模式阈值（10）。这个palindrome模式专门处理双端reads重叠的情况，效率很高。
• LEADING：从5'端开始切掉质量低于3的碱基。
• TRAILING：从3'端开始切。
• SLIDINGWINDOW:4:15：窗口大小为4，如果平均质量低于15，就从窗口位置切掉剩余部分。这个参数最常用，窗口大小和阈值可以按数据调整。
• MINLEN：保留的最短reads长度，短于36的直接丢弃。

3.4 Cutadapt 实战

如果你更习惯Python生态，Cutadapt也很棒。下面是一个双端去除接头的例子：

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
         -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
         -o sample_R1_trimmed.fastq.gz \
         -p sample_R2_trimmed.fastq.gz \
         --minimum-length 36 \
         --quality-cutoff 15,10 \
         sample_R1.fastq.gz sample_R2.fastq.gz

参数解读：

• -a / -A：分别指定R1和R2的接头序列。如果不知道序列，可以用--adapter=auto（需要安装adapterremoval插件）。
• --quality-cutoff 15,10：从3'端往5'端切，直到碱基质量≥15；如果切除后长度低于10则丢弃。这里的“15,10”是cutoff, length_before_trimming的写法，别记反了。
• --minimum-length 36：同Trimmomatic的MINLEN。

3.5 质量评估与双端同步难题

切完数据别急着走，一定要用FastQC再检查一遍。我通常看这几个指标：

• Per base sequence quality：箱线图是否整体在绿色区域？尾部有没有明显下滑？
• Adapter Content：是否降到0%附近？
• Sequence Length Distribution：分布是否集中（滑动窗口裁剪后长度会有差异）？

双端数据的难点在于同步。Trimmomatic会自动处理：如果R1被丢弃，R2也丢弃，并且输出unpaired文件。Cutadapt也有--pair-filter=any等参数。但我见过一位新手把两个文件分别处理，最后R1和R2的行数对不上，跑比对时直接报错。所以，**永远不要在处理双端数据时把两个FASTQ独立出来单独跑**。

3.6 在Nextflow中集成预处理

虽然整门课后面会详细讲Nextflow，但这里先提一句。我会把Trimmomatic或Cutadapt封装成一个process：

process TRIM {
    input:
    tuple val(sample_id), path(reads)

    output:
    tuple val(sample_id), path("${sample_id}_R1_trimmed.fastq.gz"), path("${sample_id}_R2_trimmed.fastq.gz")
    path("${sample_id}_R1_unpaired.fastq.gz"), path("${sample_id}_R2_unpaired.fastq.gz")

    script:
    """
    java -jar trimmomatic-0.39.jar PE -threads ${task.cpus} \
        ${reads[0]} ${reads[1]} \
        ${sample_id}_R1_trimmed.fastq.gz ${sample_id}_R1_unpaired.fastq.gz \
        ${sample_id}_R2_trimmed.fastq.gz ${sample_id}_R2_unpaired.fastq.gz \
        ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
        LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
    """
}

这样一来，每个样本的原始数据都会经过统一的预处理管线，而且自动并行。当然，参数可以放到配置文件中，方便后续调整。

下面这张图总结了整个预处理流程的核心步骤，你可以对照着看：

从图里可以清楚看到：数据先去掉接头，再用滑动窗口剪掉低质量尾巴，然后过滤整体质量，最后处理双端配对关系，输出干净的FASTQ。每一步的参数都会直接影响后续的比对和变异检测质量。