2. 高通量测序原理与 FASTQ 格式

讲 BWA 和 GATK 之前，我觉得还是得先聊清楚一件事：你手上的 FASTQ 文件到底是怎么来的？里面那些符号又代表什么？

我个人习惯是，不管用什么工具，先搞清楚数据源头。你想想看，如果连原始数据长什么样都不清楚，后面做再多的分析，心里也没底。

2.1 测序原理：从 DNA 到 FASTQ

高通量测序的核心，说白了就是“边合成边测序”。我们以 Illumina 平台为例，它是目前最主流的平台。

大致流程是这样的：

  1. 文库构建：把基因组 DNA 打断成几百 bp 的片段，两端接上接头（adapter）。
  2. 簇生成：DNA 片段固定在 flow cell 上，通过桥式 PCR 扩增成一个个簇（cluster）。
  3. 测序：加入带荧光标记的 dNTP，每合成一个碱基就释放一次荧光信号并拍照。
  4. 碱基识别：根据荧光信号的强度和位置，软件（如 RTA）判断每个位置的碱基是什么，同时给出这个判断的可靠程度——也就是质量值。

这一步的结果就是 FASTQ 文件。嗯，这里要注意，你从测序公司拿到的数据，通常就是 FASTQ，而不是原始的图像文件。

关键点：FASTQ 文件里记录的不仅是碱基序列，还有每个碱基的质量分数。质量分数决定了你后续能否信任这个碱基。

2.2 FASTQ 格式：四行一条记录

FASTQ 格式其实不复杂，每一条 read 用四行表示：

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

你看到的这四行，每一行都有它的意义：

行号 标识符 内容
第1行 @ read 的描述信息。包含仪器ID、flow cell编号、坐标、是否过滤等
第2行 （无） 碱基序列，A/T/C/G/N
第3行 + 分隔符。有时 + 后面会重复第1行的描述
第4行 （无） 质量字符串。长度必须等于第2行的序列长度

质量值是用 ASCII 码表示的。我把这个转换关系整理了一下：

ASCII 码值 = Phred 质量值 + 33（现代编码）
Phred 质量值（Q）= -10 * log10(P_error)

举例：
Q=30  -> 错误概率 0.001 -> ASCII 码 63 -> 字符 '?'
Q=20  -> 错误概率 0.01  -> ASCII 码 53 -> 字符 '5'
Q=10  -> 错误概率 0.1   -> ASCII 码 43 -> 字符 '+'

质量值越高，代表你越有信心这个碱基是读对的。我在项目中遇到过，有些旧数据用的是 Phred+64 编码，如果搞混了编码方式，质控报告会看着一团糟。所以拿到 FASTQ 文件的第一件事，就是确认它的质量值编码。

2.3 质量值编码的历史与分辨

为什么会存在不同编码？主要是因为测序平台的发展历史。早期的 Solexa 平台用过一个不同的公式，后来被 Illumina 收购后，又改成了标准的 Sanger 编码。

我建议你用一个简单的办法来区分：

  • Sanger / Illumina 1.8+：ASCII 码范围 33-126，字符最低是 '!'。这是目前最常见的标准。
  • Illumina 1.3-1.7：ASCII 码范围 64-126，字符最低是 '@'。老数据可能会遇到。

你可以用 Linux 的 head 和 od 命令快速看一下质量字符串第一个字符的 ASCII 码：

head -4 your.fastq | tail -1 | od -c

# 如果第一列是 '!' 或 '#'，属于 Sanger 编码
# 如果第一列是 '@' 或 'A'，可能是 Illumina 1.3+ 编码

曾经有一位刚入行的同事，直接拿旧平台的数据跑了整个 GATK 流程，结果变异检测的假阳性率高得吓人。我帮他一查，就是质量编码搞错了。所以，别小看这一步。

2.4 FASTQ 文件的基本质控

拿到 FASTQ 后，不能直接扔给 BWA 就开始跑。你需要先做一次质量评估。最常用的工具就是 FastQC。

运行 FastQC 的命令非常简单：

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o ./qc_reports/

它会生成一个 HTML 报告，里面包含 12 个模块。我个人最看重的是以下几个：

  • Per base sequence quality：看看每个位置的质量分布。如果 3' 端质量掉得厉害，可能需要考虑修剪（trim）。
  • Per sequence GC content：如果 GC 分布出现双峰，可能存在污染或文库扩增偏差。
  • Adapter Content：检查接头是否残留。如果接头含量高，跑 BWA 之前必须先做去接头。
  • Overrepresented sequences：看看有没有异常高的序列出现。

我记得有一次做 ChIP-seq 的数据，FastQC 报告显示 adapter 含量超过了 15%。后来我加了 Trimmomatic 做预处理，结果比对率从 60% 直接升到了 90% 以上。这就是预处理的价值。

2.5 测序数据产生的完整流程

下面这张图，概括了从原始样本到 FASTQ 文件的完整路径。你看一遍，应该能对整个流程有个清晰的印象。

测序数据产生流程 DNA 片段文库 桥式 PCR 簇生成 边合成边测序 荧光信号捕获 碱基识别 (RTA) FASTQ 文件 质量值分配 输出产物 序列 + 质量字符串

这张图把数据产生的关键路径串起来了。从文库构建到最终生成 FASTQ，每一个环节都可能引入误差。你后面用 BWA 做比对时，能不能把 read 正确比对到参考基因组上，很大程度就取决于这些步骤的质量。

2.6 实战小贴士：如何快速检查一个 FASTQ

在你开始正式分析之前，我建议先做三个快速检查：

  • 检查行数：FASTQ 的记录数应该是 4 的倍数。用 wc -l 看一下。
  • 检查序列长度：看看 read 的长度是否一致。如果不一致，可能是测序故障或文库问题。
  • 检查质量字符串长度：第2行和第4行长度必须完全一样。如果不一样，文件可能已损坏。

你可以用下面这个命令快速验证：

zcat sample.fastq.gz | awk 'NR%4==2 {len=length; print len}' | sort | uniq -c

# 输出会告诉你不同长度的 read 各有多少条
# 正常情况下应该只有一种长度
← 上一章 📖 返回目录 下一章 →