3、质量分数与错误率关系：如何解释质量分数

咱们做生信的人，每天跟FASTQ文件打交道，但你真的看懂那一串“质量分数”了吗？

我之前带过一个实习生，拿到FASTQ后只看碱基序列，完全忽略质量行。结果下游比对出来的SNP一半是假的……从那以后，我每次上课都要强调：质量分数就是你在碱基上下的“赌注”，看不懂它，质控就是盲人摸象。

3.1 Phred质量分数：一个数学公式背后的直觉

质量分数的核心公式长这样：

Q = -10 × log10(P)

其中P是该碱基出错的概率。举个例子：

• 如果测序仪认为这个碱基有1%的几率是错的（P=0.01），代入公式：Q = -10 × log10(0.01) = 20。这就是Q20。
• 如果P=0.001，Q=30，也就是大家常说的Q30。

说白了，Q值越大，错误率越低，你对这个碱基就越放心。

3.2 质量分数与错误率的对应表

我习惯把这个对应关系记在脑子里，查起来方便。你也不用死记，看几遍自然就熟了。

Phred质量分数	错误概率 (P)	碱基准确率	常见场景
Q10	0.1 (10%)	90%	早期Sanger测序边缘区域
Q20	0.01 (1%)	99%	很多分析的最低门槛
Q30	0.001 (0.1%)	99.9%	Illumina主流平台典型值
Q40	0.0001 (0.01%)	99.99%	高深度靶向测序要求

你看，从Q20到Q30，错误率降低了10倍。这就是为什么大家总说“Q30是金标准”——多一个数量级的差异，对低频突变检测影响很大。

3.3 在FASTQ里质量分数长什么样？

FASTQ文件里质量分数不是直接写“Q30”这样的数字，而是用ASCII字符编码来表示。最常见的编码方式是Phred+33（也叫Sanger编码）。

@SEQ_ID
AGCTTAGCTAGCTAGCTAGCT
+
IIIIIHHHHGGGGFFFFEEEE

上面第四行就代表质量分数。ASCII字符的十进制值减去33，就是对应的Q值。比如字符'I'的ASCII码是73，73-33=40，所以它代表Q40。字符'!'的ASCII码是33，33-33=0，代表Q0。

3.4 质量分数分布的解读——别被平均值骗了

我记得有一次做全外显子组分析，用FastQC看到平均质量分数有Q35，心想这数据不错。结果一跑GATK BaseRecalibrator，发现尾巴上那些碱基质量虚高得离谱……

为什么？因为质量分数是位置相关的。读长的前几个碱基往往质量不稳定，中间区域最高，末端急速下降。只看一个全局平均值，等于把高点和低点通通抹平了。

所以我在项目中坚持两条原则：

• 要看每个位置的箱线图，特别是读长末尾20bp的质量分布。
• 重点关注低质量碱基的比例，而不是平均Q值。比如Q<20的碱基占比超过10%的话，我会考虑做修剪或者直接弃用该段。

3.5 在实践中应用：一个Python小脚本

你可能需要快速从FASTQ里提取质量信息。下面这段代码是我自己经常用的——读一条记录，计算平均质量分数和低质量碱基比例。

def parse_fastq_quality(fastq_path):
    """
    从FASTQ文件提取质量分数统计
    返回：每条read的平均质量、低质量碱基数
    """
    with open(fastq_path, 'r') as f:
        while True:
            header = f.readline().strip()
            if not header:
                break
            seq = f.readline().strip()
            plus = f.readline().strip()
            qual = f.readline().strip()
            
            # 解析质量分数
            scores = [ord(c) - 33 for c in qual]
            avg_q = sum(scores) / len(scores)
            low_q_count = sum(1 for q in scores if q < 20)
            
            print(f"{header[:20]:20s} avg_Q={avg_q:.1f} low_Q20={low_q_count}")
            
# 使用示例
parse_fastq_quality("sample.fastq")

3.6 一个典型的“避坑”场景

我遇到过一个项目，RNA-seq数据FastQC报告里平均质量Q32，看起来不错。但我试着画了每个位置的箱线图，发现最后30bp的中位数掉到了Q20以下。调用fastp默认参数修剪后，保留了大部分读长，但后期定量发现基因体覆盖度在3'端严重偏低。

后来分析发现：由于3'端质量下降太快，默认参数只剪掉了质量极低的尾部，但那些中等质量区域依然保留，却给比对带来了大量错配。正确做法是按质量分数的分布设定动态阈值，而不是一刀切。

所以，我现在的建议是：

1. 先用FastQC看每个位置的质量分布箱线图，找到质量开始断崖式下降的位置。
2. 再结合read长度分布，决定是裁切尾部还是整条read丢弃。
3. 最后，在报告里展示Q20/Q30的碱基比例变化曲线，而不是只给一个平均Q值。