📚 课程目录第1章 FASTQ格式解析

第一章：FASTQ格式解析——第一关

大家好，欢迎来到这门课。既然我们要把FASTQ转成H5AD，那第一步总得先搞明白FASTQ里面到底装了什么。这玩意儿看着像天书，其实拆开来看，套路非常固定。

我个人习惯，拿到一个生信文件，先看它的原始文本结构。别急着上工具，先手动打开看一眼（当然，文件别太大）。你一看就知道，FASTQ文件其实就是四行一组的重复堆叠。嗯，就这么简单。

从FASTQ说起

FASTQ格式，说白了就是测序仪吐出来的原始数据。每一段测序片段（我们叫它read），在文件里占4行。这四行是铁打的规矩，一行都不能少。我曾经被一个同事问过，为什么他的脚本总是报错，我让他数数行数，结果发现中间夹了一行空行——测序仪可不会给你空行。

这四行的分工是这样的：

第一行：以@开头，后面跟着序列标识符和可选的描述信息。这条read的名字，就写在这里。
第二行：真正的碱基序列。A、T、C、G、N，干干净净。
第三行：以+开头，后面可以跟上和第一行一样的标识符，也可以什么都不写。就一个+符号。
第四行：质量分数。每一个碱基的测序质量，用ASCII字符编码，长度必须和第二行一样长。

关键点：四行必须严格对应。第二行和第四行的长度必须完全一致，一个是碱基，一个是质量。少一个字符都不行。

四行定乾坤：FASTQ的固定结构

我直接给你看一个真实的例子。你现在打开终端，用head -4看一个FASTQ文件，大概就会看到下面这样：

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG
CTCAAGGTTGTTGCAAGACGGGAGGTAGGTGTCTGGGTAACTCGAGGTAC
+
>>ABBBCBB?@?@?@?@?@?@?@?@?@??????????????????????

你看，是不是很清楚？第一行那个@，你想想看，它就像是个哨兵。什么？看到这个@符号就表示新的一条序列来了。第二行是碱基。第三行就是一个孤零零的+。第四行那些奇奇怪怪的字符，就是质量分数。

这里有一个我经常遇到的坑。有些新手写解析脚本，会用@作为行的开头来判断是否是新的read。逻辑是对的。但要注意，序列行（第二行）也可能以@开头。虽然概率极低，但理论上存在。我建议你用行号模4的方法来精确判断，不要单纯依赖开头的字符。

质量分数：你到底有多靠谱？

质量分数是FASTQ的精髓。它告诉你每一个碱基的测序可信度。怎么算的呢？用的是Phred质量分数。Q = -10 * log10(P)，这里的P是碱基测错的概率。

举个例子：

Q30 → 错误概率0.001，也就是千分之一。99.9%的正确率。
Q20 → 错误概率0.01，百分之一。99%的正确率。

你看，分数越高，越可靠。但问题来了，这个Q值在文件里是怎么存的？直接写数字很浪费空间嘛。所以大家用了ASCII字符来编码。

编码方式主要有两种：

编码体系	公式	ASCII范围	对应Q值范围
Phred+33 (Sanger)	ASCII = Q + 33	33 - 126	0 - 93
Phred+64 (Illumina 1.3+)	ASCII = Q + 64	64 - 126	0 - 62

避坑指南： 我曾经在一个跨平台合作的项目里，被这个编码坑惨了。上游给的FASTQ用的Phred+64，下游的比对软件默认期待Phred+33。结果所有碱基质量分数都高得不正常，导致后续的变异检测几乎全错。所以，拿到FASTQ的第一步，赶紧确认它的质量分数编码。通常看ASCII字符的最小值：如果最小值低于64，是Phred+33；如果都大于等于64，大概率是Phred+64。

怎么区分呢？关键看表头（Header）里的标签。很多软件会在@行后面用Quality encoding之类的字段告诉你。如果没写，就自己算最小ASCII值。我写了一个简单的小脚本，几行Python就能搞定。

import sys

def guess_encoding(file_path):
    min_char = '~'
    with open(file_path, 'r') as f:
        for i, line in enumerate(f):
            # 只检查质量行，即行号模4等于3（0-indexed是3）
            if i % 4 == 3:
                for ch in line.strip():
                    if ch < min_char:
                        min_char = ch
    q_min = ord(min_char)
    if q_min >= 64:
        return "Phred+64 (Illumina 1.3+)"
    else:
        return "Phred+33 (Sanger)"

print(guess_encoding("your.fastq"))

这个脚本很简单，核心就是找到最低的那个ASCII值。跑一遍就知道是什么编码了。

序列标识符：藏在名字里的元数据

@符号后面的那一串，就是序列标识符。它可不是随便写的。每个测序平台都有自己的命名规则，里面藏了大量的元数据。

以Illumina的标识符为例，它一般长这样：

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

拆开来看，用空格分隔成两部分：

第一部分（仪器标识）：
- EAS139：测序仪器名称
- 136：run编号
- FC706VJ：flowcell ID
- 2：lane编号
- 2104：tile编号
- 15343：x坐标
- 197393：y坐标
第二部分（描述信息）：
- 1：是否pair-end，1表示读段1，2表示读段2
- Y：是否通过过滤（Y是pass）
- 18：控制编号
- ATCACG：Index序列或barcode

实用技巧： 我在做多组学整合时，会利用标识符里的坐标信息来判断reads是否来自同一个tile。这可以帮助我做空间转录组的数据质控。你想想看，这个标识符就像是一张身份证，能追溯到你测序时的物理位置。

不过要注意，不同平台的标识符格式差异很大。比如Illumina和PacBio的写法完全不同。我建议你遇到一个新测序平台的数据，先去查一下它的官方命名规范，不要凭经验瞎猜。

一张图看懂FASTQ

光说理论不够直观。我画了张图，把这四行结构、质量编码、标识符的层级关系串起来。你看完就全明白了。

总结一下流程图： 顶部那四个彩色方块就是FASTQ的铁四角。每四行构成一个read。左侧的标识符藏着测序仪信息，右侧的质量分数每条字母对应一个碱基。底部是我让你重点关注的两种质量编码。后续所有的质控、比对、定量，都建立在这个基础之上。

如何用程序处理FASTQ

好，概念讲完了，咱们来看点实际的。你怎么用代码去解析它？我个人习惯用Python的Biopython库，封装得很好。但你如果要做高性能处理，C++或者Rust更合适。

先看一个最简单的解析例子，用Python遍历一个FASTQ文件：

from Bio import SeqIO

for record in SeqIO.parse("sample.fastq", "fastq"):
    print(f"ID: {record.id}")
    print(f"Sequence: {record.seq}")
    print(f"Quality: {record.letter_annotations['phred_quality']}")
    print("---")
    break  # 只看第一条

这段代码会帮你把每条read的ID、碱基序列和Phred质量分数列表取出来。还是那句话，先搞清楚编码。Biopython默认是按Phred+33解析的。如果你的文件是Phred+64，记得要转换。

性能建议： 我处理过几百GB的FASTQ文件，用Python的SeqIO.parse虽然方便，但内存消耗大。我建议你用SeqIO.index来创建索引，或者直接按行读取文本再手动解析。对，用最原始的方式，速度反而更快。

避坑小清单

聊了这么多，我觉得有必要把我在实际项目中踩过的坑整理一下。你遇到这些问题时，能少走很多弯路。

行数不对： 记住，FASTQ文件的行数必须是4的倍数。如果不是，说明文件损坏或不完整。我曾经遇到一个下载了一半的文件，导致下游脚本崩溃了整整一天。
空行陷阱： 大多数工具不处理空行。如果你用脚本逐行读取，碰到空行就直接跳过，否则你的行号索引会乱掉。
@符号误判： 只用line[0] == '@'来判断新纪录是不保险的。更可靠的是用计数器，每4行一组。
质量分数编码混用： 这是一个大坑。我建议无论在哪个环节，都把文件统一转换成Phred+33（Sanger格式）。几乎所有现代工具都支持这个格式。

好了，FASTQ格式这块，该说的都说了。从四行结构到质量编码，再到标识符解读，你应该已经能独立解析一个FASTQ文件了。聊得差不多了，咱们先把这第一关过扎实。