2、Phred质量值体系:从ASCII到碱基准确率

测序数据里最基础也最绕不开的概念,就是 Phred 质量值。你可能在 FASTQ 文件里见过一堆带符号的字符串,比如 IIIIII 或 AAAAAA。嗯,这些符号背后藏着一个精确的数字,告诉你每个碱基的可靠程度。我刚开始接触时也觉得这层转换挺绕,但弄懂了之后,一切就变得很直白。

2.1 质量值是怎么来的?

Phred 质量值最早是为 Sanger 测序开发的,后来 NGS 数据也沿用了这套体系。它的核心公式很简单:

Q = -10 × log10(P)

其中 P 是碱基识别错误的概率。打个比方,如果某个碱基的 Q 值是 20,那么它的错误概率就是 1/100,准确率 99%。Q30 就是错误概率 1/1000,准确率 99.9%。

我在项目中遇到过刚入行的同学,看到 FASTQ 里一串 III 能直接看出质量最高只有 40。其实你只要记住几个关键阈值就行:

  • Q10 → 90% 准确率,一般直接扔掉
  • Q20 → 99% 准确率,可接受
  • Q30 → 99.9% 准确率,很多分析要求这个我建议作为过滤底线

2.2 从 ASCII 到质量值的转换

FASTQ 文件里存的是字符,不是数字。每一个质量字符对应一个 ASCII 码,再减去一个固定的偏移就能得到 Q 值。这就是所谓的 Phred+偏移 编码。

早期 Illumina 用过 Phred+64,不过现在绝大多数平台(Sanger、Illumina 1.8+)都统一用 Phred+33。转换就是一句话:

Q = ord(quality_char) - 33

举个例子:字符 I 的 ASCII 是 73,减去 33 得 40,这就是 Q40。要是你看到 !(ASCII 33),减去 33 就是 0,说明这个碱基完全不可信。

不同编码方案对比

编码名称 偏移量 Q 范围 ASCII 范围 常见平台
Sanger / Illumina 1.8+ 33 0 – 41 ! – J 大多数现代 NGS 数据
Illumina 1.3 – 1.7 64 0 – 62 @ – ~ 旧 Illumina 数据
Illumina 1.0 – 1.2 64 0 – 40 B – h 非常老的 GA 数据

我的个人习惯:拿到 FASTQ 之后第一件事就是看质量行第一个字符的 ASCII 码。如果第一个字符是 !,基本就是 Phred+33;如果是 @,那就很可能是旧 Illumina 格式。这招让我跳过不少坑。

2.3 ASCII 范围速查表

下面这张表帮你快速定位质量字符对应的 Q 值,实际用到时可以直接查:

字符ASCIIQ (Phred+33)准确率
!330—
"341—
#352—
……………………
5532099%
?633099.9%
I734099.99%
J7441—

注意:Illumina 1.8+ 的最大 Q 值是 41,所以字符 J 就是上限了。老版本 Illumina 能到 Q62,但那种数据我建议直接重新跑流程——质量标度已经乱套了。

2.4 碱基准确率与 Q 值的关系

把公式反过来,得到错误概率 P = 10^(-Q/10)。我用一个简单表格总结常用阈值:

Q 值错误概率碱基准确率常见用途
Q100.1 (10%)90%基本不可信,一般过滤掉
Q200.01 (1%)99%许多比对工具的最小要求
Q300.001 (0.1%)99.9%变异检测推荐阈值
Q400.0001 (0.01%)99.99%超高精度,少见

避坑指南:我曾经在分析公共数据时遇到一个 "混合编码" 的 FASTQ,同一条序列的质量行里既有小于 ! 的字符(ASCII 32),又有 @。一查发现是旧版 Illumina 未正确转换合并的结果。这时候千万不能直接减 33 或 64,得先识别编码。我写过一个简单的探测函数,用第一个字符的 ASCII 做判断:

def guess_phred_offset(first_char):
    ascii_val = ord(first_char)
    if ascii_val >= 64:
        return 64   # 旧 Illumina
    elif ascii_val >= 33:
        return 33   # 新格式
    else:
        raise ValueError("Unsupported encoding")

2.5 实战:在 FASTQ 过滤中应用质量值

我在做 RNA-seq 时经常写一个小脚本,把质量行转换成数组,然后统计低于某个阈值的碱基比例。如果一条 read 有超过 20% 的碱基 Q<20,直接扔掉。下面是最核心的转换片段:

import gzip

def parse_fastq_quality(line_bytes):
    # 假设编码为 Phred+33
    quality_str = line_bytes.decode().strip()
    quals = [ord(c) - 33 for c in quality_str]
    return quals

# 使用时,按 4 行一组读取 FASTQ
with gzip.open('sample.fastq.gz', 'rb') as f:
    for name, seq, _, qual in zip(*[f]*4):
        quals = parse_fastq_quality(qual)
        low_q = sum(1 for q in quals if q < 20)
        if low_q / len(quals) > 0.2:
            continue   # 丢弃
        # 保存高质量的 read

关键点:不要直接对字符做大小比较,一定要先转换到整数 Q 值。另外,记得解码时忽略末尾的换行符,不然 ord('\n') 会把最后一位算错。

2.6 整个知识体系的逻辑

Phred 质量值数据流与转换逻辑 原始信号 碱基识别 质量分数 Q ASCII 字符 Phred+33 转换公式 Q = ASCII(字符) - 33 碱基准确率 准确率 = 1 - 10^(-Q/10) Q20 → 99%

这张图概括了“原始信号 → 碱基 → 质量分数 → ASCII 存储 → 准确率”的完整链路。你在实际处理 FASTQ 时,就是从右向左解包:先读字符,再算 Q 值,最后映射到准确率。理清这个流程,后面的过滤策略就有了坚实的根基。