4. 数据准备

拿到测序数据，第一件事不是直接跑比对。
你想想看，如果原始数据里一堆接头污染、低质量碱基，后面分析做得再花哨也没用。

这章咱们就聊聊怎么把FASTQ「洗」干净。
我个人的习惯是：先质检，再清洗，再质检。这样才能确保数据质量可追溯。

为什么要做质控？

测序仪不是完美的机器。
它会引入测序错误、接头残留、低质量碱基。尤其Illumina平台，读长越长末端质量越差。

我遇到过一个大样本项目，直接跑HaplotypeCaller，结果出来大量假阳性SNP。
回头一查——原始数据Q20比例不到85%。白跑了一周时间。

FastQC：质量检查第一关

FastQC是最常用的质检工具。它生成的HTML报告让我能一眼看出数据状况。
我个人习惯先对整个批次做一个快速扫描，找出明显有问题的样本。

运行FastQC

# 单样本
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./fastqc_results/

# 批量处理，所有样本
fastqc *.fastq.gz -o ./fastqc_results/ -t 8

参数说明：
-o 指定输出目录，-t 使用的线程数。

如何读报告？

FastQC报告里我最关注这几个模块：

multiqc ./fastqc_results/ -o ./multiqc_report/

Trimmomatic：数据清洗利器

FastQC查出问题后，就该Trimmomatic上场了。
它能完成三件事：切除接头、滑窗裁剪低质量碱基、过滤短读。

基础命令框架

java -jar trimmomatic-0.39.jar PE \
  -threads 8 \
  -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 \
  LEADING:3 \
  TRAILING:3 \
  SLIDINGWINDOW:4:15 \
  MINLEN:36

我来拆解一下这些参数实际含义：

• PE：paired-end模式。如果是SE，换成SE即可。
• -phred33：质量值编码。现在基本都是Phred+33，旧数据可能是+64。
• ILLUMINACLIP：接头切除。TruSeq3-PE.fa是接头库文件，后面三个数字分别代表：种子错配数、palindrome score、simple score。
• LEADING:3 / TRAILING:3：切除首尾质量低于Q3的碱基。
• SLIDINGWINDOW:4:15：每4个碱基一个窗口，平均质量低于Q15则切除。
• MINLEN:36：切除后读长短于36bp则丢弃。

实际项目中的参数调整

不同项目参数差别挺大。我列一个经验表，供你参考：

质控后的验证

清洗完后，一定要再做一次FastQC。
对比清洗前后的报告，确认三个指标：

1. 接头残留接近0
2. 低质量区域明显减少
3. 读长分布集中在目标长度范围内

我常用的一个快速对比命令：

# 统计清洗前后读长数量变化
echo "原始数据:" && zcat sample_R1.fastq.gz | awk 'NR%4==1' | wc -l
echo "清洗后:" && zcat sample_R1_paired.fastq.gz | awk 'NR%4==1' | wc -l

如果清洗后读长数量减少超过40%，那就得回头检查参数是不是太激进了。

常见错误处理

项目里最常遇到的问题，我简单列一下：

• 内存不足：Trimmomatic是Java程序，用-Xmx指定堆内存，比如java -Xmx16g -jar trimmomatic-0.39.jar
• 接头库不存在：去Trimmomatic安装目录下的adapters/找对的文件
• 两个read文件名顺序反了：PE模式必须R1在前，R2在后，否则会出现mating信息错乱
• 质量编码错误：老数据用Phred+64，Trimmomatic不认识会报错。加上-phred64即可

我记得有一次，一个合作方给的数据明明是最新的Novaseq产出，结果Trimmomatic死活报Phred编码错误。排查了半天，原来是他们压缩时用了pigz多线程，不同线程写入顺序乱了，导致一个FASTQ里混着新旧两种质量值编码。最后只能让他们重新从bcl解压。