测序基础：FASTQ、BAM 与参考基因组

做生信分析，说白了就是跟这三个文件格式打交道。FASTQ、BAM 和参考基因组，它们贯穿了整个流程的始终。我刚开始接触这一行时，觉得不就是几个文件嘛，有什么好学的？结果后面做 HaplotypeCaller 时，发现很多奇奇怪怪的错误，源头都是对这几个基础格式理解不到位。

所以咱们这一章，我把这些年踩过的坑和总结的经验一起讲给你听。

FASTQ：你看到的第一个测序文件

测序仪下机后，你拿到的第一个文件就是 FASTQ。它存的是每条 read 的序列和对应的碱基质量。

核心要点：FASTQ 中每 4 行为一条完整记录。第 1 行是 @ 开头的标识符，第 2 行是序列，第 3 行是 +，第 4 行是质量值（与序列等长）。

一个典型的 FASTQ 记录长这样：

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

质量值那行看起来像乱码，其实每个字符编码了一个数字。这个数字是 Phred 质量分数，公式是 Q = -10 * log10(P)，其中 P 是该碱基测错的概率。举个例子，Q30 表示错误率千分之一，Q20 就是百分之一。

我的习惯：拿到一个 FASTQ，我第一件事就是检查质量分布。用 FastQC 看一眼，如果后半段质量掉得厉害，我会用 fastp 或 Trimmomatic 做一下修剪。别怕丢数据，烂数据进去，烂结果出来。

不过你要注意，不同测序平台给的质量编码体系不一样。Illumina 1.8+ 用的是 Phred+33，Sanger 也是这个。但有些老版本用的是 Phred+64。我以前接过一个项目，没留意这个差异，结果下游分析全乱了。嗯，这里有个口诀：「33 常见，64 少用，转换用 seqtk 或 BioPython。」

编码体系 ASCII 范围 质量范围 常见平台
Phred+33 (Sanger/Illumina 1.8+) 33~126 0~93 主流 Illumina 数据
Phred+64 (Illumina 1.3/1.5) 64~126 0~62 老版本 Illumina

BAM：比对后的核心数据格式

把 FASTQ 比对到参考基因组上后，你就得到了 SAM/BAM 文件。SAM 是文本格式，体积大；BAM 是它的压缩二进制版，生信分析中用的基本是 BAM。

一个 BAM 文件里，每一行代表一条 read 的比对结果。关键字段我列一下：

  • QNAME：read 的名字，对应 FASTQ 中 @ 后面的部分。
  • FLAG：一个整数，用位运算表示比对属性。比如 0x4 表示 unmapped，0x10 表示 reverse strand。
  • RNAME：参考序列的名字，比如 chr1 或 1。
  • POS：比对到参考序列的起始位置（1-based）。
  • MAPQ：比对质量，越高越可靠。
  • CIGAR：描述比对情况的字符串，比如 10M1D20M 表示 10 个匹配、1 个缺失、20 个匹配。
  • SEQ 和 QUAL：read 的序列和质量，跟 FASTQ 一样。

避坑指南：我曾经在处理全外显子数据时，发现 BAM 文件里有很多 read 被标记为 duplicate。后来一查，是 PCR 扩增导致。如果你不做标记（mark duplicate），HaplotypeCaller 会高估这些区域的覆盖度，导致假阳性变异。记住：BAM 的排序和去重是必须步骤。

查看 BAM 内容，最常用的工具是 samtools：

# 查看 BAM 头信息
samtools view -H sample.bam

# 查看 BAM 记录（前 10 条）
samtools view sample.bam | head -10

# 统计比对情况
samtools flagstat sample.bam

参考基因组：比对的标准答案

参考基因组，说白了就是一条「标准」的模板序列。你把测序得到的片段往这个模板上「贴」，就知道它们来自哪里。

人类参考基因组最有名的是 GRCh38（也叫 hg38）和 GRCh37（hg19）。GRCh38 比 GRCh37 多了很多修正，比如修复了某些 gaps，加入了替代序列。

关键：选择哪个版本，取决于你的下游分析工具和已有数据。比如一些公共数据库的注释文件是基于 hg19 的，你用 GRCh38 去跑，就得先做 liftover 转换。

参考基因组文件通常是这样组织的：

  • FASTA 文件：包含所有染色体的序列。
  • FAI 索引：让工具能快速定位到染色体的某一段。
  • dict 文件：Picard 或 GATK 需要的字典文件，记录了染色体名字和长度。

构建索引的命令：

# 创建 FASTA 索引
samtools faidx reference.fasta

# 创建 GATK 格式的字典
gatk CreateSequenceDictionary -R reference.fasta

我个人习惯是，把参考基因组放到一个统一目录，所有项目都链接过去。不要每个项目都复制一份，太占空间，也容易搞混版本。我曾经犯过一个错：在一个项目里混淆了 hg19 和 GRCh38，结果跑出来的变异列表跟数据库对不上，排查了大半天。自那以后，我每次建项目第一步就是确认参考基因组的 MD5 值。

三者如何串联起来

FAQST → BAM → 参考基因组，它们的关系是这样：

  1. 测序仪产出 FASTQ。
  2. 比对工具（如 BWA）把 FASTQ 的 read 与参考基因组匹配，产出 SAM/BAM。
  3. BAM 作为 HaplotypeCaller 的输入，结合参考基因组，进行变异检测。

没有参考基因组，比对无从谈起；没有高质量的 BAM，变异检测就是空中楼阁。每一步的质控都至关重要。

我的经验：拿到新数据，我会先跑一个 mini pipeline 测数据质量：FASTQ 用 FastQC，BAM 用 samtools stats，再画一个覆盖度分布图。整个流程走下来，数据能不能用、哪部分需要清理，心里就八九不离十了。

测序数据分析核心数据流 测序仪 FASTQ 文件 产出 BAM 文件 比对 (BWA/STAR) 参考基因组 (FASTA) 参考序列 HaplotypeCaller 变异检测 质控节点：FastQC | samtools stats | flagstat

回到我们这门课的主题——HaplotypeCaller。说句实话，如果你不把 FASTQ、BAM 和参考基因组吃透，后面的分析就像在沙子上盖房子。所以这一章虽然讲的是「基础」，但我觉得它比后面任何一章都重要。

← 上一章 📖 返回目录 下一章 →