3. 测序数据质量评估与预处理
讲变异检测，第一关就是测序数据本身的质量。你想想看，如果原始数据里全是错误碱基，后续的变异准确性评估本质上就是「垃圾进，垃圾出」。我个人习惯把数据质量评估看作整个分析流程的「体检环节」，这一步省不得。
你在弄变异准确性评估或者是ROC曲线教程之前，首先要搞清楚：你的数据到底能不能用？质量评估是整个变异检测 pipeline 的基石。今天我就把这块内容掰开揉碎了讲清楚。
3.1. 为什么要做质量评估？
测序仪出来的原始数据，看着是一堆碱基（A/T/C/G），实际上里面有不少噪声。这些噪声来自测序本身，比如平台误差、试剂衰减、样本降解等等。我刚开始做生物信息分析时，有一次直接拿没过滤的 FASTQ 跑变异检测，结果假阳性率高得吓人。嗯，从那以后我再也不敢跳过这一步了。
简单来总结，质量评估的目的就是两件事：
- 判断数据整体质量是否达标：碱基准确度够不够？读长分布是否异常？
- 识别需要剪裁或丢弃的测序读段：接头污染、低质量尾部、测序错误过多的片段。
3.2. 核心评估指标
我们一般用 Phred 质量分数来量化每个碱基的准确性。它的定义很简单：
Q = -10 × log₁₀(P)其中 P 是该碱基判定错误的概率。Q30 意味着错误概率千分之一，也就是碱基准确率 99.9%。
| 质量分数 | 错误概率 | 碱基准确率 | 含义 |
|---|---|---|---|
| Q20 | 1% | 99% | 通常可接受 |
| Q30 | 0.1% | 99.9% | 推荐标准 |
| Q40 | 0.01% | 99.99% | 高质量 |
关键经验：在我经手的项目中，一般要求 >85% 的碱基达到 Q30。如果低于这个标准，变异检测的假阳性会明显增加。特别是低频突变检测，质量不够的话根本没法做。
3.3. 流行工具：FastQC
我评估数据最常用的工具是 FastQC。它不是唯一的，但可以说是工业标准。你线上有 Java 环境就能跑，简单直接。
基本命令长这样：
java -jar fastqc_v0.12.1.jar sample_R1.fastq.gz -o ./fastqc_reports它会生成一个 HTML 报告，包含10多个模块。但说句实话，不是每个模块都值得你去细看。我一般会优先关注下面这4张图：
- Per base sequence quality：碱基在测序循环上的质量分布箱线图。如果末端（后 50bp）质量掉到 Q20 以下，我通常会考虑截短。
- Per sequence quality scores：所有读段的平均质量分数分布。峰值如果低于 Q27，说明数据整体偏差。
- Adapter content：接头残留比例。超过 5% 的话，必须去除接头。
- Sequence duplication levels：重复序列比例。高重复率可能提示 PCR 扩增偏差或文库复杂性不足。
避坑小技巧：我曾经遇到过一个问题，FastQC 明明报告「Per base sequence quality」显示了一片红，但其实是该批次测序仪质量编码格式用了 Illumina 1.3 的 Phred+64 体系而非默认的 Sanger Phred+33。先把编码格式搞清楚再下结论，否则你会白白扔掉很多好数据。
3.4. 质量预处理三步走
数据评估完，发现问题怎么办？我主要分三步来处理：
- 去除接头（Adapter Trimming）：用工具识别读段中包含的接头序列并剪掉。注意双端测序的话通常两端都要处理。推荐
cutadapt或Trimmomatic。 - 质量截短（Quality Trimming）：从读段 3' 端开始，逐碱基切除质量低于阈值的尾部。一般我用
cutadapt -q 20。 - 长度过滤（Length Filtering）：截短后，太短的读段（比如 <30 bp）意义不大，直接丢弃。
下面是一个标准处理脚本的示例：
# 使用 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:adapters.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:36参数说明：
| 参数 | 作用 |
|---|---|
ILLUMINACLIP |
比对并切除接头序列 |
LEADING:3 |
切除 5' 端质量低于 3 的碱基 |
TRAILING:3 |
切除 3' 端质量低于 3 的碱基 |
SLIDINGWINDOW:4:20 |
以 4bp 为窗口扫描，如果平均质量 <20，切除整个窗口及后续序列 |
MINLEN:36 |
丢弃处理后长度小于 36bp 的读段 |
注意：生产中不要盲目使用统一参数。不同的测序平台（HiSeq vs NovaSeq）、不同的文库构建方法、甚至不同的样本类型，最优参数都会变化。我曾经在某个大型肿瘤测序项目中，发现 Trimmomatic 的 SLIDINGWINDOW:4:20 参数太保守，导致大量短读段被丢掉，最后不得不把阈值放宽到 4:15 才保住数据量。灵活调整才是王道。
3.5. 质量评估与预处理流程图
下面这张图把上面讲的全流程串起来了。从原始 FASTQ 到高质量读段，每一步都很关键。
3.6. 处理后的检查
预处理完了，我建议再做一次 FastQC，然后前后对比一下。
- 看质量箱线图是否明显改善：特别是读段尾部，应该从红色变成绿色。
- 检查接头残留是否已清除：一般要在 1% 以下。
- 统计数据总量损失：一般预处理会丢掉 10~20% 的数据。如果超过 30%，你就得回去查文库构建或者测序过程是不是出了什么问题。
讲到这里，你基本知道怎么对一个测序数据集做质量评估和预处理了。后面无论你是跑 GATK 还是做 ROC 曲线评估变异准确性，这一步都是前提。磨刀不误砍柴工，数据干净了，后面的分析才会靠谱。