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 值,实际用到时可以直接查:
| 字符 | ASCII | Q (Phred+33) | 准确率 |
|---|---|---|---|
| ! | 33 | 0 | — |
| " | 34 | 1 | — |
| # | 35 | 2 | — |
| …… | …… | …… | …… |
| 5 | 53 | 20 | 99% |
| ? | 63 | 30 | 99.9% |
| I | 73 | 40 | 99.99% |
| J | 74 | 41 | — |
注意:Illumina 1.8+ 的最大 Q 值是 41,所以字符 J 就是上限了。老版本 Illumina 能到 Q62,但那种数据我建议直接重新跑流程——质量标度已经乱套了。
2.4 碱基准确率与 Q 值的关系
把公式反过来,得到错误概率 P = 10^(-Q/10)。我用一个简单表格总结常用阈值:
| Q 值 | 错误概率 | 碱基准确率 | 常见用途 |
|---|---|---|---|
| Q10 | 0.1 (10%) | 90% | 基本不可信,一般过滤掉 |
| Q20 | 0.01 (1%) | 99% | 许多比对工具的最小要求 |
| Q30 | 0.001 (0.1%) | 99.9% | 变异检测推荐阈值 |
| Q40 | 0.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 整个知识体系的逻辑
这张图概括了“原始信号 → 碱基 → 质量分数 → ASCII 存储 → 准确率”的完整链路。你在实际处理 FASTQ 时,就是从右向左解包:先读字符,再算 Q 值,最后映射到准确率。理清这个流程,后面的过滤策略就有了坚实的根基。