📚 课程目录第2章测序原理与FASTQ格式

2. FASTQ数据格式详解：质量值编码、序列标识、常见问题

FASTQ这个格式，说白了就是二代测序数据的“通用语言”。你从测序仪拿到的原始数据，十有八九都是这个格式。我最早接触它的时候，觉得不就是四行文本嘛，有啥好研究的？后来在实际项目中吃过亏，才发现里面的门道挺多。

咱们这一章就把FASTQ格式拆开揉碎了讲。你理解了它，后面所有的质控、比对、定量操作才有了根基。

2.1 FASTQ的四行结构

一个标准的FASTQ记录，由且仅由四行组成。别小看这四行，每一行都有自己的规矩。

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTCTCTCCCTCTGCCCCTCNCCGCTCTGC
+
IIIIIIIIIIIIIIIIIIIIIIIIIII8II

咱们一行一行看：

第一行（标识行）：以@开头，后面跟着测序仪的标识符和坐标信息。不同测序平台，这行的格式差异很大。
第二行（序列行）：就是实际的碱基序列，由A、T、C、G和N组成。N代表该位点测序仪无法确定是哪个碱基。
第三行（分隔行）：一般只有一个+符号。有时候后面会重复第一行的标识符，但大多数时候就是个光杆+。
第四行（质量行）：每个字符对应第二行每个碱基的测序质量。字符的数量必须和第二行的碱基数完全一致，多一个少一个都不行。

我的检查习惯：拿到一个新数据集，我第一件事不是跑质控，而是先数一下FASTQ里每个记录是不是都是四行。我遇到过测序仪软件bug，偶尔会输出五行或者三行的坏记录。这种文件直接拿去比对，会直接崩掉。

2.2 质量值编码：Phred分数背后的秘密

质量行那些乱码一样的字符，其实每个都对应一个数字，叫Phred质量分数。公式很简单：

Q = -10 × log₁₀(P)

这里P是碱基测错的概率。Q值越高，说明这个碱基越可靠。

Q值	错误概率	碱基准确率
10	1/10 (10%)	90%
20	1/100 (1%)	99%
30	1/1000 (0.1%)	99.9%
40	1/10000 (0.01%)	99.99%

但是问题来了——质量分数是一个数值，怎么用单个ASCII字符来表示？这就引出了编码方案。

Phred+33 与 Phred+64

最早Illumina用了Phred+64方案，但后来整个行业基本都统一到了Phred+33，也叫Sanger编码。咱们现在看到的FASTQ，99%都是Phred+33。

计算方式：

Phred+33：ASCII字符 = Q值 + 33。Q值范围0-41，对应ASCII码33-74（字符是 ! 到 J）。
Phred+64：ASCII字符 = Q值 + 64。Q值范围0-41，对应ASCII码64-105（字符是 @ 到 i）。

你想想看，如果Q值是0，Phred+33编码出来是!（ASCII 33），而Phred+64编码出来是@（ASCII 64）。这两个字符看着就像正常字符，很容易搞混。

我曾经踩过的坑：有一次我从公共数据库下载了一个老数据集，直接拿FastQC去看。结果质量分布图显示前半段质量奇高，后半段断崖式下跌。我折腾了半天，才发现是编码方案搞错了——数据是Phred+64，而工具默认用了Phred+33。识别的办法很简单：看看质量值的字符范围，如果出现大量低于@的符号（比如 !、"、#），那就是Phred+33；如果最低都是@以上，那很可能是Phred+64。

下面这个简单的Python脚本，可以帮你检查FASTQ的编码类型：

def guess_encoding(qual_string):
    ord_vals = [ord(c) for c in qual_string[:100]]
    min_val = min(ord_vals)
    max_val = max(ord_vals)
    
    if min_val >= 64:
        return "Phred+64"
    elif min_val >= 33 and max_val <= 74:
        return "Phred+33 (Sanger)"
    else:
        return "无法判断，可能编码异常"

2.3 序列标识符：读懂测序仪的语言

第一行的@后面那串东西，不是乱写的。不同测序平台的标识符格式差异很大，但传递的信息都类似。

Illumina平台（最常见）

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG

我习惯把它拆成两个部分：

前半段（冒号分隔）：EAS139:136:FC706VJ:2:2104:15343:197393
- 仪器ID（EAS139）
- 运行编号（136）
- FlowCell ID（FC706VJ）
- Lane号（2）
- Tile号（2104）
- X坐标（15343）
- Y坐标（197393）
后半段（空格后）：1:Y:18:ATCACG
- Pair成员（1表示read1，2表示read2）
- 过滤标记（Y表示通过）
- Index编号（18）
- Index实际序列（ATCACG）

PacBio / ONT 平台

长读长测序的标识符就随意多了，通常只包含read ID和运行信息。

@m54011_190513_220730/5/ccs
@runid_flowcell_sample_read1234

为什么我在乎这个？ 在做多文库合并分析时，如果两个文库的Index序列不一样，你直接从标识符就能看出来。我还见过有人把双端测序的R1和R2混在一起当单端数据跑，一看标识符里的"1"和"2"标记，立刻就能发现问题。

2.4 常见问题与避坑指南

FASTQ看着简单，但实际项目中总会出现各种幺蛾子。我把这些年遇到的典型问题整理了一下。

问题类型	表现	原因	怎么办？
编码方案误判	质量值分布异常	用了错误的Phred编码	用FastQC检查；或用脚本扫描字符范围
质量行与序列长度不匹配	比对工具报错	数据截断或拼接错误	用awk检查每行的长度一致性
行数不是4的倍数	工具无法解析	文件不完整或写入中断	用`wc -l`检查；去除最后不完整记录
标识符包含空格/特殊字符	下游工具分割失败	测序软件或手动修改	用sed替换空格为下划线
碱基序列包含小写字母	比对忽略或报错	某些工具输出时未统一大小写	用`tr 'acgt' 'ACGT'`处理

快速质量检查命令

我一般拿到FASTQ文件，先跑下面这几个命令看看情况：

# 检查记录数是否是4的倍数
cat sample.fastq | wc -l
# 结果除以4，应该得到整数

# 检查每条记录的序列长度是否等于质量行长度
awk 'NR%4==2 {seq_len=length($0)} NR%4==0 {if (length($0) != seq_len) print "长度不匹配, 记录号:", NR/4}' sample.fastq

# 查看前几条记录的标识符
head -n 40 sample.fastq | awk 'NR%4==1'

# 查看质量值的字符范围（判断编码类型）
head -n 400 sample.fastq | awk 'NR%4==0 {for(i=1;i<=length($0);i++){c=substr($0,i,1);if(ord[c]=="") ord[c]=sprintf("%d",ord[c]?ord[c]:0);ord[c]++}} END {for(c in ord) printf("%s:%d ",c,ord[c]); print ""}'

注意：第三个命令里的awk版本可能不支持关联数组的正确初始化。如果报错，直接用Python或Perl写个更稳健的检查脚本更省心。

关于FASTQ格式，我个人有个原则：在开始任何正式分析之前，先花5分钟做一次基础检查。你省掉的这5分钟，可能会在后面变成5小时的debug时间。理解清楚这四行数据，你后面做质控、比对、定量的时候，心里就有底了。