📚 课程目录第3章数据准备

数据准备：从原始测序到比对（BWA-MEM）

到了这一步，咱们终于要跟真正的测序数据打交道了。说实话，很多新手容易在这一步踩坑——手里拿着 fastq 文件就直接往后跑，结果到 HaplotypeCaller 那一步才发现数据质量不对，又得回头重来。

我自己的习惯是，数据准备阶段宁可多花半小时把每一步都检查清楚，也绝不留着隐患往下传。这条路我走过太多次了。今天就带你从头捋一遍，从原始测序文件到生成合格的 BAM，到底该怎么操作。

1. 原始测序数据格式与质量评估

测序仪下机后，最常见的格式就是 FASTQ。这里面记录着每条读段的碱基序列和对应的质量分数。你打开一看，大概是这个样子：

@EAS139:136:FC706VJ:2:2104:15343:197393
CTCTCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

第一行是读段 ID，包含了机器和坐标信息；第二行是碱基；第三行是个加号分隔；第四行是质量值。质量值用的是 Phred 分数编码，ASCII 字符越高代表置信度越高。

我个人有个铁规矩：拿到 FASTQ 的第一件事，跑 FastQC。这事花不了十分钟，但能省下未来好几天排查问题的功夫。

FastQC 输出必看的几个指标：

每个碱基位置的 Phred 分数箱线图——尾端质量掉到 Q20 以下要小心
GC 含量分布——理论上应该接近正态，有过多的尖峰可能提示污染
超过代表值的 N 含量——N 太多说明测序质量有问题
读段长度分布——如果是变长测序，要了解分布范围

举个例子，有一次我拿到一批肿瘤样本的测序数据，FastQC 一跑发现 GC 含量出现了双峰。我第一反应是样品之间可能存在交叉污染，后来一查果然是索引混了。要不是提前做了 FastQC，那整批分析结果全得废掉。

运行 FastQC 的命令很简单：

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

输出是一份 HTML 报告，绿色是通过、黄色是警告、红色是失败。嗯，这里要注意：警告和失败不意味着数据一定不能用，但你必须清楚它在警告什么，并在后续分析中加以考虑。

2. 参考基因组的准备

比对之前得先把参考基因组准备好。说白了，你得让 BWA 认识你手里的参考序列长什么样。

主流做法是从 NCBI 或 Ensembl 下载参考基因组 FASTA 文件。我个人习惯用 UCSC 的 hg38 或 GRCh38，记得要附带 .fai 索引文件和 .dict 字典文件。

# 下载参考基因组（以 GRCh38 为例）
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
gunzip hg38.fa.gz

# 构建 BWA 索引（这一步会比较耗时）
bwa index -a bwtsw hg38.fa

# 同时需要构建 samtools 索引
samtools faidx hg38.fa

# 以及 GATK 需要的 .dict 文件
gatk CreateSequenceDictionary -R hg38.fa -O hg38.dict

小技巧：构建 bwa index 时指定 -a bwtsw 是针对大基因组（超过 2G）的模式。人类基因组超过 3G，请务必使用这个选项。我曾经看到有人用默认模式跑了一整天还没跑完……那就是选错了算法。

索引构建完成后，你会得到 .amb、.ann、.bwt、.pac、.sa 这五个文件。缺少任何一个，BWA-MEM 都会报错退出。

3. BWA-MEM 比对流程

数据质量确认了，参考基因组也准备好了，那就进入真正的比对环节。BWA-MEM 是目前单样本比对最常用的工具，支持双端测序，算法较新，对长读段和嵌合读段的处理都挺不错。

基本命令格式是这样的：

# 双端测序
bwa mem -t 8 -R '@RG\tID:samplename\tSM:sample01\tPL:ILLUMINA' \
  hg38.fa \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  > sample.sam

说说几个参数：

-t：线程数，我一般设为 CPU 核数减 2，留点资源给系统
-R：Read Group 信息，这是 极其关键 的部分，GATK 的 HaplotypeCaller 依赖 RG 来区分样本和文库
ID、SM、PL 是必填的三项：ID 是测序文库的唯一标识，SM 是样本名称，PL 是测序平台（ILLUMINA、IONTORRENT 等）

重要提醒：不要省略 -R 参数？要是丢了 Read Group 信息，GATK 后期处理会报“缺失 RG”的错误，你还得回头重跑一次比对。我早期做项目时就吃过这个亏，后来养成了习惯，每次写命令前先把 -R 字符串写好再粘贴。

为了便于后续处理，我个人建议一步到位输出压缩格式：

# 直接输出 BAM，省略中间 SAM 文件
bwa mem -t 8 -R '@RG\tID:sample01\tSM:sample01\tPL:ILLUMINA' \
  hg38.fa \
  sample_R1.fastq.gz sample_R2.fastq.gz | \
  samtools view -bS -o sample.raw.bam -

这么做的好处很明显：节省磁盘空间。一个典型全基因组测序的 SAM 文件轻松上几十个 GB，但 BAM 压缩后能缩小到几分之一。

4. 比对后处理：排序、去重

BWA-MEM 输出的 BAM 是 按读段 ID 排序 的，也就是按照测序仪输出的顺序来的。这对于下游分析几乎没用，因为 HaplotypeCaller 需要的是 按基因组坐标排序 的数据。

排序用 samtools：

samtools sort -@ 8 -m 4G -o sample.sorted.bam sample.raw.bam

-m 4G 限制每个线程使用 4G 内存，总内存消耗大约是 4G × 8 = 32G。如果你的服务器内存紧张，可以调低线程数或减小 -m 值。

排序完成后，下一步是 标记重复（Mark Duplicates）。测序文库构建过程中，PCR 扩增会造成部分片段被多次测到，这些重复读段如果不处理，会在变异检测时引入偏向性。

GATK 4 自带了 MarkDuplicates 工具，也可以使用 Picard 或 Sambamba：

# 使用 GATK MarkDuplicates
gatk MarkDuplicates \
  -I sample.sorted.bam \
  -O sample.dedup.bam \
  -M sample.metrics.txt \
  --REMOVE_DUPLICATES false

注意：这里我设置的是 --REMOVE_DUPLICATES false，也就是只标记不删除。为什么？因为很多下游工具（包括 HaplotypeCaller）需要通过标记来判断哪些读段是重复的，但又不希望数据总量的减少影响统计模型。尤其是低频突变检测，保留标记信息比直接删除更稳妥。

标记完成后，建立索引：

samtools index sample.dedup.bam

索引文件 .bai 或 .csi 是必须的，GATK 要求输入 BAM 必须带索引。

5. 数据质控：比对率与覆盖度

到了这一步，你以为数据就可以交付了？别急，先做一轮必要的质控。我通常跑两个命令：

# 比对统计
samtools flagstat sample.dedup.bam

# 覆盖度统计
samtools depth sample.dedup.bam | awk '{sum+=$3} END {print "平均覆盖度:", sum/NR}'

samtools flagstat 的输出会告诉你总读段数、比对上的读段数、成对比对数等。我一般要求：

指标	合格线	说明
总体比对率	> 90%	人类全基因组数据通常可达 95% 以上
成对比对率	> 85%	双端都能比对才算有效配对
重复率	< 20%	超过 30% 要考虑是否有过度扩增
平均覆盖度	30×（WGS）/ 100×（靶向）	随实验设计不同有较大差异

个人经验：如果比对率低于 85%，先别急着调参数。检查一下参考基因组版本是否跟你的数据匹配。我曾遇到过用 hg19 的参考来比对 GRCh38 的测序数据，比对率只有 60%，非常惨。换个正确的参考立马就回到 95% 以上。

这一步做完了，你的 BAM 文件才真正达到了交付 GATK 的标准。我自己在流程管理上习惯用 Makefile 或 Snakemake 把以上所有步骤串起来，这样每次重跑只需要一句命令，也不会因为手误遗漏某个步骤。

最后再啰嗦一句：比对阶段如果发现了异常数据（比方说质量太差导致大量未比对读段），不要硬着头皮往下走。我见过有人用低质量数据跑 HaplotypeCaller，出来的变异列表里一半都是假阳性，后期人工验证全打回去，反而更浪费时间。不如在源头上把数据质量把控好。

至此，你的数据从原始的 FASTQ 经过 FastQC 检查、参考基因组索引构建、BWA-MEM 比对、排序去重，再到最终质控，已经准备好进入变异检测阶段了。