3. 数据输入与质量控制：从原始FASTQ到干净数据

这一章咱们聊聊NGS分析的第一步——数据入场和清洗。很多人觉得质量控制是体力活，不就是跑个FastQC、截个接头嘛。但我告诉你，这一步做不好，后面GATK跑出来的变异调用全是废的。我在项目中就吃过这个亏，所以今天把话说透。

3.1 FASTQ文件格式——你的数据长什么样？

先搞清楚我们的原材料。FASTQ不是普通的文本文件，它是序列+质量的打包袋。每一份测序read占4行：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTACGT
+
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ

第一行是标识符，包含仪器ID、流道坐标、barcode等信息。第二行是碱基序列。第三行只有一个加号，偶尔跟重复的标识符。第四行是碱基质量分数，每个符号对应一个ASCII码，代表该碱基的测序错误概率。

质量分数用的是Phred体系。早期版本有Phred+33和Phred+64，现在主流Illumina数据都用Phred+33。简单换算：ASCII值减去33就是Q值。Q20意味着错误概率1%，Q30是0.1%。打个比方，Q30的碱基你读1000个，大概只有1个是错的。

我的经验：千万别假设所有FASTQ都是Phred+33。有些旧数据或特殊平台会用Phred+64。我曾遇到一个项目，所有比对率都异常低，折腾两天才发现是质量编码搞错了。跑FastQC时看一眼"Per base sequence quality"图的纵坐标范围——如果是0-40就是Phred+33，如果0-60就要警惕。

3.2 原始数据质量评估：FastQC怎么看？

拿到FASTQ先别急着跑流程。我的习惯是：先扔给FastQC把所有样品都看一眼。命令行很简单：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -t 4 -o fastqc_reports/

-t指定线程数，-o输出目录。FastQC会生成HTML报告，里面塞了11个模块。你不需要每个都深究，但有三块必须盯紧：

模块 红灯阈值 怎么办？
Per base sequence quality 任何位置的中位质量低于Q28 考虑末端裁剪（trimming）
Adapter Content 任何位置adapter比例 >5% 必须去除接头
Overrepresented sequences 任何序列占比 >0.1% 可能是接头、rRNA或扩增产物

看报告有个诀窍：不要只盯着Summary里的红叉叉。我曾经见过一个样品，所有条目都是绿色，但Per base GC content出现了双峰——后来发现是两种不同物种的混样。你会习惯地右滑鼠标，看看每个模块的详细曲线。尤其是Kmer Content那一块，异常峰值往往暗示PCR扩增偏好或污染。

避坑：FastQC对压缩文件（.fastq.gz）不是不能读，但单个文件太大时会内存溢出。我一般先解压到临时目录，或者用zcat管道传过去：

zcat sample_R1.fastq.gz | head -n 400000 | fastqc stdin: -o fastqc_reports/

这招在服务器上尤其好用，解压100G的文件不占用磁盘空间。

3.3 接头去除与低质量过滤：Trimmomatic vs cutadapt

质检完了，就该动刀子了。市面上工具一堆，我常用的就两个：Trimmomatic和cutadapt。它们各有长短，我来对比一下。

对比项 Trimmomatic cutadapt
接头去除方式 滑动窗口比对（ILLUMINACLIP） 精确匹配或模糊匹配（允许错配）
滑动窗口质量裁剪 内置SLIDINGWINDOW，非常方便 需要结合其他工具（如seqtk）
双端模式 优秀，支持PE文件同步裁剪 也支持，但需要额外参数
速度 快（Java实现，多线程） 更快（C实现，单线程但可并行）
输出文件管理 自动生成paired/unpaired四个文件 手动指定输出

我个人习惯：跑illumina双端数据首选Trimmomatic。原因很简单——它把质量裁剪和接头去除打包在一起，一条命令搞定，而且处理的paired/unpaired逻辑非常清晰。看一个典型命令：

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  R1_paired.fastq.gz R1_unpaired.fastq.gz \
  R2_paired.fastq.gz R2_unpaired.fastq.gz \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:True \
  LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

拆解一下参数：

  • ILLUMINACLIP后面跟着接头文件、允许的最大错配数（2）、配对时得分阈值（30）、单端得分阈值（10）、最小接头长度（2）、是否保持反向互补（True）。
  • LEADING:3——去掉read开头质量低于3的碱基。想省事就设5。
  • TRAILING:3——同样处理尾部。
  • SLIDINGWINDOW:4:15——以4个碱基为窗口，若窗口平均质量低于15就截断。这个参数我一般调到4:20，更保守一点。
  • MINLEN:36——小于36bp的read直接丢弃。

为什么我偏爱滑动窗口？你想想看，测序仪越往后读质量越差，但如果整条read的平均质量还行，直接丢弃太可惜了。滑动窗口就是帮你切掉坏尾巴，保留好头。我用这个策略处理过一批深度很高的RNA-seq数据，本来只有60%的mapping率，加了这个参数后直接跳到85%。

如果数据比较特殊，比如有大量未知接头，或者想精细控制错配，我会换cutadapt。它的接头检测更灵活，支持正则和双端内部匹配。举个例子：

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
         -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
         -o R1_trimmed.fastq.gz -p R2_trimmed.fastq.gz \
         -j 4 -q 15 -m 36 \
         sample_R1.fastq.gz sample_R2.fastq.gz

-a和-A分别是R1和R2的接头序列，-j线程数，-q是质量阈值（类似Trimmomatic的滑动窗口？其实不一样，cutadapt的-q是全序列裁剪末端低质量碱基）。注意cutadapt的-q不如Trimmomatic的滑动窗口精细，所以如果你需要精确控制，还是用Trimmomatic。

3.4 过滤决策的黄金法则

说了这么多工具，到底什么时候必须过滤？我总结了几条硬指标：

  • 接头污染 >1%：必须去接头。我见过一个全外显子项目，接头残留达到15%，结果GATK的HaplotypeCaller把接头序列当成了变异——冤不冤？
  • 中位质量 < Q30：强烈建议裁剪末端，至少切除最后10bp。如果整个read质量都低（比如中位Q20），直接丢弃。
  • GC含量异常：检查是否为文库偏差，如果是PCR扩增偏好，可以考虑用--clip-overlap处理双端重叠区域。
  • 重复率过高：在FastQC的Sequence Duplication Levels里如果出现30%以上重复，且不是扩增产物，先怀疑是不是测序深度不够或者PCR过度。这一步可以在后续用MarkDuplicates处理。

3.5 一张图看懂质控流程

下面是我平时走完一个WGS/WES样品时的标准质量流程示意图。它把今天讲的东西串起来了。

NGS数据质量控制流程图 GATK上游数据质量控制流程图 原始FASTQ（R1/R2） FastQC质量评估 通过？ 是 否 Trimmomatic/Cutadapt 干净FASTQ → BWA比对 说明：如果FastQC未通过，则进入Trimmomatic/Cutadapt清洗，然后重新跑FastQC验证。 迭代直到质量达标。clean FASTQ才能安全送进BWA和GATK管线。

3.6 一些你可能会踩的坑

最后分享几个实战教训。

  • 不要把所有样品用一个参数清洗。不同文库制备方法、不同测序批次，接头残留和质量分布差异很大。我习惯给每个样品单独跑FastQC，然后根据具体情况微调Trimmomatic参数。比如有的样品GC含量偏低，可能是文库捕获出了偏差，这时候裁剪窗口就不能设太高，否则雪上加霜。
  • 双端数据要同步处理。如果用了cutadapt，记得加上--pair-filter=any，否则R1和R2的条数会不对齐。我自己因为忘加这个参数，导致下游BWA报错，排查了整整一下午。用Trimmomatic默认就处理好，省心。
  • 接头文件要选对。Trimmomatic的ILLUMINACLIP需要指定包含接头序列的FASTA文件。不同的试剂盒对应的接头不同——TruSeq、Nextera、NEB Ultra II都各有专用文件。我犯过一个低级错误：用TruSeq的接头文件去处理Nextera数据，结果一个接头都没去干净，因为序列完全不一样。
  • 别指望一刀切能解决所有问题。如果某个样品的平均质量低于Q20，再剪也无法挽回。这时候要考虑是测序失败还是文库问题，直接标记为不合格，不要在分析管线里强塞。我见过有人强行过滤后保留了一堆低质量数据，导致变异检测结果里有大量假阳性——这种人我称之为“数据强迫症患者”，得不偿失。

好了，数据清洗这一步说到这。记住：干净的数据是GATK流程的基石，花再长时间在上面都不为过。

← 上一章 📖 返回目录 下一章 →