4. 平均质量分数（Mean Quality Score）计算与整体评估

拿到一个FASTQ文件，你第一眼看什么？

我个人习惯，先扫一眼平均质量分数。这就好比体检报告上的体温——不复杂，但能快速判断这人是不是在发烧。平均质量分数，说白了就是把一条read上每个碱基的Phred质量分数加在一起取平均，得到一个总体的质量标杆。

4.1 平均质量分数怎么算出来的？

先说数学。假设一条read有L个碱基，每个碱基的质量分数是Qi（i从1到L），那么：

Mean Q = (Q1 + Q2 + ... + QL) / L

嗯，就是一个算术平均。但这里有个细节——Q值本身就是对数转换过的，所以平均Q值并不是直接对应平均错误率。你想想看，Q30对应错误率0.1%，Q20是1%，如果一条read前半段是Q30后半段是Q20，平均大概在Q25左右，但这条read的真实错误率并不是0.3%那么简单。

我在项目中遇到过有人把这个概念搞混，直接用平均Q值去推平均错误率，结果质控报告和实际比对结果对不上。后来我给他画了个公式才说明白：

平均错误率 = (1/L) * Σ(10^(-Qi/10))

说白了，均值在Q空间算和在错误率空间算，结果是两码事。我们日常说的平均质量分数，指的就是 算术平均Q值，不是平均错误率再倒推回来的Q值。

4.2 阈值设定——不同场景，不同标准

平均质量分数该设多少？这问题没有标准答案。但我可以给你几个参考线，都是这些年项目里摸出来的：

应用场景 推荐阈值 我个人判断
WGS（全基因组） Q30 以上 低于Q30，变异检测假阳性开始爬升
RNA-seq Q28 以上 比对容忍度稍高，但低质量区会影响表达定量
ChIP-seq Q30 以上 peak calling对质量敏感，低于Q30容易出噪声
单细胞测序 Q25 以上 数据量大，有些平台质量波动大，Q25是个折中
靶向测序（Amplicon） Q28 以上 靶区域覆盖深，但对变异低频检测还是有要求

这些阈值不是拍脑袋定的。我记得有一年做某个肿瘤panel的项目，一开始设了Q25作为平均分阈值，结果下游突变检测的假阳性率高得吓人。后来调成Q30，过滤掉了一批边角料数据，准确率直接提升了两个百分点。

关键原则：阈值设定要跟着你的下游分析走。做低频突变检测，阈值往高调；做物种鉴定，阈值可以适当放低。没有万能阈值，只有适合你实验目的的阈值。

4.3 实现计算——简单脚本书上

你真想算平均质量分数，最常见的做法是用FastQC的"Per sequence quality scores"模块，或者用我自己经常写的Python小脚本：

import gzip

def calc_mean_quality(fastq_file):
    with gzip.open(fastq_file, 'rt') as f:
        total_q = 0
        total_bases = 0
        for i, line in enumerate(f):
            if i % 4 == 3:  # 质量行
                for char in line.strip():
                    q = ord(char) - 33  # Phred+33 编码
                    total_q += q
                    total_bases += 1
    return total_q / total_bases if total_bases else 0

mean_q = calc_mean_quality('sample_R1.fastq.gz')
print(f"平均质量分数: {mean_q:.2f}")

这个脚本速度一般，但在小样本上足够用了。我有个习惯——处理项目前先用这个脚本筛一遍所有样本，看看有没有整体质量偏低到离谱的，直接标记出来，不耽误后面跑完整的FastQC。

4.4 整体评估——不能只看一个数

平均质量分数能反映整条read的"平均健康水平"，但它也有盲区。

举个例子。一条read前半段全是Q40，后半段全是Q10，平均下来是Q25。另一条read整条都是Q25。两条read的平均分一样，但第一条有半截数据几乎不能用，第二条整体还说得过去。平均分没有告诉你这个分布信息。

我建议把平均质量分数和以下几个指标结合起来看：

  • 碱基位置质量分布图——看看质量是不是在末端急剧下滑
  • GC含量偏差——平均分高但GC异常，可能是污染
  • 序列重复率——高重复率加上高平均分，可能是PCR bias
  • N碱基比例——N比例高但平均分还正常，说明测序仪在某些碱基上罢工了

避坑指南：我曾经跟一个合作者吵过一架。他坚持说所有样本平均Q值都在Q33以上，数据质量没问题。结果我拉了一张碱基位置质量图，发现每条read的后20个碱基质量跌到了Q15以下。平均分骗了我们，因为前120个碱基拉高了整体数值。从那以后，我再也不迷信平均分了。

4.5 不同条件、不同测序平台的影响

不同测序平台，质量分布差别挺大。我做过的项目里：

  • Illumina NovaSeq 6000：平均分通常在Q35以上，SBS化学稳定
  • Illumina MiSeq：平均分波动大，有时Q30，有时Q37，跟跑样批次关系很大
  • BGI/MGI平台：平均分略低一点，Q32-Q35，但是一致性不错
  • PacBio HiFi：质量模式不一样，average Q一般在Q20-Q30之间，但这是环形一致测序的正常范围
  • ONT（Nanopore）：平均Q值比较低，Q10-Q15常见，但人家数据结构不同，不能用同一套标准

你想想看，如果把ONT数据的平均Q调到Q30，那基本一条数据都过不了。所以阈值一定要跟测序平台和解码方法对齐。

4.6 质量过滤与mean Q的联动

实际项目中，我们经常用平均质量分数来做过滤条件。比如在fastp里：

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o clean_R1.fastq.gz -O clean_R2.fastq.gz \
      --average_qual 30

这个参数的意思是：扔掉平均质量分数低于Q30的read。加上它之后，数据量一般会减少5%-15%，取决于原始质量。

我遇到过最极端的情况——一个RNA-seq样本，原始数据2.5G read，设了平均Q30过滤后只剩1.2G read。乍一看吓一跳，但一看原始质量报告，原来这个样本的测序质量从中间就开始崩了。去掉那批低质量read后，下游比对率反而从72%升到了89%。所以有时候去掉"有害数据"比保留所有数据更重要。

提示：定阈值时，建议先用5%的数据做个试探性分析。选一个样本跑不同阈值（Q20、Q25、Q30、Q35），看看每个阈值下剩余数据量和下游比对率的关系，找到拐点。这个拐点，往往就是最合适的阈值。

4.7 知识结构总览

下面这张图把平均质量分数在整个质控体系中的位置梳理清楚了。我习惯在做任何大规模项目前，先把这个流程走一遍：

平均质量分数计算与整体评估流程 原始FASTQ文件 Phred质量分数 → 累加 → 取平均 Mean Q = Σ(Qi) / L 设定判据与阈值 场景→推荐阈值 WGS: Q30 / RNA: Q28 结合其他指标 位置质量图 / GC / 重复率 整体质量评估 通过 / 复检 / 重测 平均Q值不能单独看，必须与位置质量分布、GC偏差、重复率联合分析，才能准确评估数据质量

这个流程图的逻辑很简单：从原始数据出发，先计算出平均质量分数，再根据场景设定阈值，同时拉上其他质控指标一起参照，最后做出整体评估。我在实际项目中就是这样一步步走的，很少出大问题。

4.8 常见的理解误区

最后列几个我经常看到的坑：

  • 误区一：平均Q值越高，数据越好。 不一定。高平均分可能是某些区域疯狂好，另一些区域很差，平均掩盖了问题。要结合位置分布看。
  • 误区二：不同平台之间可以直接比平均Q。 别比。Illumina的Q30和PacBio的Q30虽然定义相同，但测序原理不同，错误模式也不一样。拿同一把尺子量不同东西，没有意义。
  • 误区三：过滤后平均Q提升了，说明过滤有效。 嗯，表面上看是这样。但你可能只是把难读的序列扔掉了，而这些序列里可能正好包含你要找的变异。我经历过一次，过滤后平均Q从Q28涨到Q32，结果一个低频突变被一起过滤掉了。
  • 误区四：一个样本平均Q低，整批数据都要重测。 其实可以先把低质量样本单独拿出来分析，看看是文库问题、上机问题还是测序仪问题。我有个项目四个样本平均Q都Q25以下，但其中一个居然是样本本身降解，重新提取后数据就好了。

总结一下我的做法：平均质量分数是第一道安检门。它快速、直观，能帮你筛掉最差的那批数据。但它不是终点。真正靠谱的质控，是把这个指标和其他维度放在一起看，像医生看体检报告一样——单看一项指标看不出什么，综合几项才能下诊断。

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