FASTQ 文件格式详解
FASTQ 这东西,说白了就是给每条 DNA 序列配了一份「质量成绩单」。我在生信行业摸爬滚打这些年,见过不少新手被这个格式坑得够呛——序列截断了、碱基读串了、质量值对不上号…… 这章咱们一步到位,把 FASTQ 的结构、质量分数、常见陷阱都捋清楚。
四行结构 — 一次只读一条
一条标准 FASTQ 记录固定由四行组成,缺一不可。你想想看,如果某个工具把这四行拆散了,后面的分析全都得错位。
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTTAACGGGTTGCACCACGCGGACGCATTGTGAGTA
+
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
- @行(标识行) — 以
@开头,后面跟着测序仪生成的唯一标识符。一般包含仪器名、flowcell 编号、lane、tile、x/y 坐标等。我个人习惯用head -1看一眼 @ 行,就能知道这套数据从哪台机器出来的。 - 序列行 — 就是测出来的碱基序列,大写 A/C/G/T/N。长度必须和质量行一致,否则后面一堆报错等着你。
- + 行 — 通常只是一个
+,但有些格式会把 @ 行的标识再抄一遍(可选)。我一般不管它,但解析时最好别假设内容为空。 - 质量行 — 每个碱基对应的质量分数,用 ASCII 字符表示。字符个数必须跟序列行一样,每个字符对应一个 Phred 分数。
awk 'NR % 4 == 2 {len=length($0)} NR % 4 == 0 {if(length($0)!=len) print "mismatch at line", NR}' sample.fastq
跑一下这个命令,一秒就能发现格式问题。
质量分数 — 把「信心值」变成 ASCII
质量分数是 FASTQ 的灵魂。每条碱基的测序质量用一个数字表示——Phred 质量分数 Q。公式很简单:
Q = -10 × log₁₀(P)
这里的 P 是碱基识别错误的概率。Q 越高,错误概率越低。
| Q 值 | 错误概率 | 碱基准确率 | 实际含义(比如我调过的例子) |
|---|---|---|---|
| 10 | 1/10 = 10% | 90% | 基本不可靠,我会直接过滤掉 |
| 20 | 1/100 = 1% | 99% | 常规过滤下限,很多流程设 Q20 |
| 30 | 1/1000 = 0.1% | 99.9% | 变异的黄金标准,我定过滤器时常用 Q30 |
| 40 | 1/10000 = 0.01% | 99.99% | 非常高,但在 read 末端很难达到 |
但问题来了——计算机怎么存储这个数字?总不能把一个整数直接塞进文本文件吧。于是大家约定:把 Q 值加上一个偏移量,转换成 ASCII 字符。不同世代的测序平台用的偏移不一样:
- Sanger / Illumina 1.8+(最常见):ASCII = Q + 33。质量行内容是 ASCII 33~126 的字符。
- Illumina 1.3~1.7(旧格式):ASCII = Q + 64。质量行内容是 ASCII 64~126。
- Solexa(已淘汰):另有一套映射,遇到概率极低。
B(ASCII 66)。一查才发现,某个老工具用 Illumina 1.3+ 格式写的,而我用 Sanger 编码去解析,导致所有质量值偏大 31。那次比对结果惨不忍睹——错误 SNP 筛出一大堆。所以碰到陌生文件,第一件事用 FastQC 或者 seqtk 检测质量编码。
几类常见「信息」藏在哪里
除了序列和质量,FASTQ 标识行里塞了不少元数据。拿上面那个例子说:
@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
- 仪器 ID:
EAS139 - run ID:
136 - flowcell ID:
FC706VJ - lane 编号:
2 - tile 编号:
2104 - x 坐标:
15343,y 坐标:197393 - pair member:
1(双端测序时的第几端) - 是否通过过滤:
Y - control bits:
18 - index(barcode):
ATCACG
这些信息在多文库拆分、去除 PCR 重复、调试测序问题时非常有用。我每次做变异检测分析前,都会用 awk 统计一下 lane 分布,看看有没有某个 lane 质量特别差——如果有,直接丢弃那个 lane 的数据。
为什么质量分数要用 ASCII 字符?
你想想看,如果直接存数字,比如 37,得占两个字节。但 ASCII 字符 F(对应 70)只占一个字节。在二代测序动辄几百 G 的数据量下,每个碱基省一个字节,整个文件就能小一半。嗯,这里要注意:虽然文件小了,解析却要多一步偏移计算。我用 Python 解析时通常会写一个简单的转换:
def qual_to_phred(ascii_char, offset=33):
return ord(ascii_char) - offset
qual_str = "JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ"
phreds = [qual_to_phred(c, offset=33) for c in qual_str]
print(phreds[:5]) # [42, 42, 42, 42, 42]
一句话总结:FASTQ 本质上是一个四行循环的容器,序列与质量一一对应,质量值藏在 ASCII 背后。掌握它的结构,是走通生信分析 pipeline 的第一步。
这张图把 FASTQ 的四行结构和质量分数映射串在了一起。你看,从左边四行读出来,中间把 ASCII 转成 Phred 值,再对应到错误概率,右边是我给你的检测建议。格式问题不可怕,就怕你不知道自己掉了哪个坑。