数据质量指标体系
质量控制这件事，说白了就是建立一套标准。
标准不清晰，后头的分析全是空中楼阁。我最早做RNA-seq项目时，上来就直接跑比对，结果一大堆垃圾数据混进去，后续差异表达算出来全是假阳性。嗯，从那以后我再也不敢跳过质控了。
多组学数据的核心指标，其实就三大类：测序质量、比对率、覆盖度。这三者是递进关系，一个环节出问题，后面的指标必然受影响。
一、测序质量 —— 一切的基础
测序质量得分，行业内通行标准是 Phred 质量分数 (Q-score)。公式很简单：Q = -10 × log₁₀(P)，P 是碱基识别错误的概率。
举个例子：Q30 意味着错误率 1/1000，也就是 99.9% 的准确率。Q20 对应 1% 的错误率。我个人习惯，对于常规全基因组测序，底线设在 Q30 > 80%。低于这个数，我会直接要求重测。
核心阈值参考：
- Q20 (99% 准确率)： 至少 90% 碱基达到
- Q30 (99.9% 准确率)： 至少 80% 碱基达到
- GC 含量偏差： 理论值与实测值偏差 < 5%
- N 占比： 未识别碱基占比 < 1%
一个小技巧： 看 FASTQC 报告时，别只看平均值。我记得有一次项目里 Q30 均值有 88%，但后 50bp 的序列直接掉到了 Q15。这种情况，截断处理比直接扔掉更划算。
我个人建议把 FASTQC 跑出来的结果做成这样一份表：
| 样本 | 总 reads | Q20 占比 | Q30 占比 | GC 偏差 | 结论 |
|---|---|---|---|---|---|
| 正常样本 | 50M | 93.2% | 86.7% | +1.2% | 通过 |
| 问题样本 | 45M | 81.4% | 68.3% | +8.5% | 不通过 |
| 边缘样本 | 48M | 89.1% | 78.9% | -3.4% | 需截断 |
下面是一份典型的 FASTQC 质控报告截图，注意看末尾的警告：
Filename: sample_01.fastq
Total Sequences: 52678342
Filtered Sequences: 0
Sequence length: 150
%GC: 42.5
>>Per base sequence quality pass
>>Per tile sequence quality pass
>>Per sequence quality scores pass
>>Per base N content warn
>>Sequence Duplication Levels fail
>>Overrepresented sequences fail
>>Adapter Content fail嗯，这里要注意：FASTQC 里的 "fail" 并不一定代表这组数据完全不能用。它只是告诉你某个指标偏离了期望分布。比如 Adapter Content 出现 fail，你用 cutadapt 或 Trimmomatic 切一下就行。
二、比对率 —— 映射的质量衡量
测序质量过关了，下一步就是比对率。你想想看，如果有一半的 reads 都怼不到参考基因组上，那这组数据基本上就是废的。
| 数据类型 | 总比对率期待值 | 唯一比对率期待值 | 说明 |
|---|---|---|---|
| WGS 人全基因组 | ≥ 95% | ≥ 80% | 参考基因组完整度高 |
| WES 全外显子 | ≥ 90% | ≥ 75% | 探针捕获效率影响很大 |
| RNA-seq | ≥ 70% | ≥ 55% | 转录组有重复区域 |
| 植物 / 非模式生物 | ≥ 60% | ≥ 45% | 参考基因组质量差 |
我曾经踩过一个坑： 有一个植物 RNA-seq 项目，比对率只有 53%。我当时第一反应是转录组重复高，但后来发现是因为参考基因组里线粒体序列被污染了。比对时，多拷贝 reads 全怼到了线粒体区域。所以，比对率不达标先查参考基因组质量。
比对率异常的情况，我通常按以下步骤排查：
- 检查重复率： 如果重复率超过 40%，说明可能是 PCR 扩增偏差或建库有问题。
- 检查 rRNA 占比： 对于 RNA-seq，rRNA 如果超过 10%，说明 rRNA 去除不彻底。
- 检查比对参数： STAR 或 HISAT2 的参数设置，特别是 --outFilterMultimapNmax 和 --outFilterMismatchNmax。
# STAR 比对参数示例（我常用的）
STAR --genomeDir /ref/GRCh38 \
--readFilesIn sample_01_R1.fastq sample_01_R2.fastq \
--outSAMtype BAM SortedByCoordinate \
--outFilterMultimapNmax 10 \
--outFilterMismatchNmax 3 \
--outFilterMismatchNoverLmax 0.05 \
--alignIntronMin 20 \
--alignIntronMax 1000000这里的几个参数为什么要这么设？
outFilterMultimapNmax 10：允许每个 read 最多比对到 10 个位置。太高容易引入多拷贝基因的噪音，太低又会丢掉真正的同源基因。outFilterMismatchNmax 3：最多允许 3 个错配。这是 150bp 读长的常见设置，错配率控制在 2% 以内。
三、覆盖度 —— 深度与广度
覆盖度这个指标，其实包含两个维度：
深度： 平均每个碱基被多少条 reads 覆盖（一般说“30x”就指平均 30 条 reads）。
广度： 有多少比例的基因组区域被覆盖到了（通常要求 ≥ 99%）。
我举个例子你就会明白。假如一个基因组的覆盖度是 30x，但某个关键区域只有 1x 覆盖，那么在这个区域里，一个碱基的测序错误就能让你得出错误的结论。
| 应用场景 | 推荐覆盖深度 | 最低广度要求 |
|---|---|---|
| 全基因组测序 (WGS) 变异检测 | 30-40x | ≥ 99% 区域 ≥ 10x |
| 全外显子组测序 (WES) | 100-200x | ≥ 95% 靶向区域 ≥ 50x |
| RNA-seq 基因表达定量 | 20-50M reads | ≥ 70% 基因 ≥ 10 counts |
| 单细胞 RNA-seq | 50000 reads/cell | ≥ 80% 细胞 ≥ 2000 genes |
查看覆盖度时，我习惯用 samtools depth 统计全基因组分布，再用一个简单的 Python 脚本做分位数统计。你总不能光看平均值就结束了吧？
# 获取染色体 1 上所有位点的覆盖深度
samtools depth -r chr1 sample_01.bam > coverage_chr1.txt
# 覆盖深度分布统计
awk '{
if($3 >= 1) covered++
total++
sum+=$3
if($3 >= 10) covered_10x++
if($3 >= 20) covered_20x++
if($3 >= 30) covered_30x++
} END {
print "平均深度:", sum/total
print "覆盖率 (≥1x):", (covered/total)*100 "%"
print "覆盖率 (≥10x):", (covered_10x/total)*100 "%"
print "覆盖率 (≥20x):", (covered_20x/total)*100 "%"
print "覆盖率 (≥30x):", (covered_30x/total)*100 "%"
}' coverage_chr1.txt输出示例：
平均深度: 32.4
覆盖率 (≥1x): 99.8%
覆盖率 (≥10x): 97.2%
覆盖率 (≥20x): 91.5%
覆盖率 (≥30x): 80.3%一个实用经验： 我处理 WGS 数据时，要求 ≥ 10x 的覆盖率必须达到 99% 以上。如果达不到这个数，CNV 的检测质量会变得非常差。对于外显子捕获的数据，重点看 捕获效率 (PCT_TARGET_BASES_20X)，这个值低于 90% 说明探针设计或杂交效率有问题。
四、核心指标速查表
最后，我整理了一份各平台常见指标的合格范围。你可以直接把它贴在工位旁当参考。
| 平台类型 | 关键指标 | 合格范围 | 警戒阈值 |
|---|---|---|---|
| Illumina WGS | Q30 比例 | > 85% | < 80% |
| 总比对率 | > 95% | < 90% | |
| 覆盖均一性 (均匀度) | Gini < 0.2 | Gini > 0.3 | |
| PacBio HiFi | 平均精度 | > Q20 | < Q15 |
| 总产量 | > 15Gb / cell | < 10Gb / cell | |
| 次读长 | > 10kb | < 5kb | |
| ONT 超长读长 | 中位 read 长度 | > 20kb | < 10kb |
| 质量得分 | > Q10 | < Q7 |
这些指标之间的关系不是孤立的。测序质量差，直接导致比对率下降；比对率下降，有效覆盖度自然就缩水。所以，你做质控时一定要按顺序检查：先看 FASTQC，再看比对日志，最后算覆盖度。跳一步，就可能让后续分析徒劳无功。
我个人习惯把每个样本的质控结果汇总成一个 JSON 文件，方便自动化判断。这样项目跑完后，直接扫一眼就知道哪些样本需要重测或剔除。
{
"sample_01": {
"Q30_pct": 87.3,
"mapping_rate_pct": 96.5,
"unique_mapping_pct": 82.1,
"mean_depth": 33.2,
"coverage_10x_pct": 99.1,
"gc_deviation_pct": 0.8,
"status": "PASS"
},
"sample_02": {
"Q30_pct": 72.4,
"mapping_rate_pct": 88.7,
"unique_mapping_pct": 74.3,
"mean_depth": 22.5,
"coverage_10x_pct": 94.6,
"gc_deviation_pct": -6.2,
"status": "FAIL"
}
}这样一份 JSON，配合一个简单的条件判断脚本，就能自动化完成大样本量的质控分诊。嗯，效率高了不少。