BWA与GATK全基因组变异检测工作流 - 第1节
为什么要学BWA和GATK？
说实话，我经常被问到这个问题：“生信入门，该先从哪个工具下手？”
我的答案始终是——BWA + GATK。
为什么？
因为它俩组成了目前最主流、最成熟的短读长测序数据分析管线。
你手上的Illumina数据，无论是WGS还是WES，想找人基因组里的SNP、Indel，基本绕不开这个组合。
一句话总结：BWA负责“比对”，GATK负责“找变异”。
整条管线就是：原始FASTQ → BWA比对 → SAM/BAM处理 → GATK变异检测。
我自己最早接触这套流程，是被一个全外显子组项目逼的。
一百多个样本，要批量跑。
那时候还是GATK 3.x版，光跑一个HaplotypeCaller就花了两天两夜。
但跑完之后，看到最后的VCF文件，心里那个踏实——值得。
所以这门课，我会把从原始数据到最终VCF的每一步，手把手带你走一遍。
不会只扔几个命令行就跑，咱们会聊原理、参数、常见坑，还有我自己的实战经验。
背景知识：DNA、测序、参考基因组
先快速过一遍必备背景。如果你对这些已经滚瓜烂熟，可以直接跳到下一节。
DNA序列与测序片段
人基因组大约30亿个碱基对。
测序仪读出来的，是短片段，叫reads。
以Illumina为例，常见的read长度：PE150（双端，每端150个碱基）。
我们拿到的原始数据，一般是FASTQ格式。
里面不光有碱基序列，还有质量分数（Phred score）。
这个分数在后面对比对和变异检测的影响挺大。
一个我踩过的坑：
有次组里同学直接跑BWA，没做QC，结果比对率只有60%。
一看原始FASTQ，低质量碱基太多，接头也没去干净。
所以课前准备：先把FASTQ的质量跑明白。
参考基因组
比对必须有个“地图”——这就是参考基因组。
目前人基因组最新的是GRCh38（也叫hg38）。
你从NCBI、UCSC或Ensembl都能下到FASTA文件。
注意：GATK要求参考基因组必须建立索引。
除了BWA用的bwa index，还要有.dict和.fai文件。
这些索引文件不用自己手动写命令之外的额外操作，后面课程会讲。
| 文件名 | 索引工具 | 用途 |
|---|---|---|
| reference.fasta | - | 参考基因序列本身 |
| reference.fasta.fai | samtools faidx | 快速随机访问任意区域 |
| reference.dict | GATK自带（CreateSequenceDictionary） | 记录染色体名称、长度等元信息 |
| bwa索引文件（.amb, .ann, .bwt, .pac, .sa） | bwa index | BWA自己的索引格式 |
BWA：比对引擎
BWA全称是Burrows-Wheeler Aligner。
它的核心算法基于BWT（Burrows-Wheeler Transform）和FM-index。
说白了，就是用压缩技术加速比对，让一个30亿碱基的基因组能在内存里被快速搜索。
BWA提供了三种比对算法：
- BWA-MEM：最常用，适合70bp-1Mbp长度的reads，默认就是它。
- BWA-ALN：老算法，适合短reads（≤100bp），但现在已经很少人用了。
- BWA-SW：支持长片段和嵌合比对，速度慢。
我个人的习惯：不管reads多长，一律用BWA-MEM。
它的准确性、速度、内存消耗都平衡得很好。
要是你有单端数据，加上-1和-2选项就好。
避坑指南：
BWA-MEM默认输出的是SAM格式（文本）。
如果你不重定向或加管道，文件会直接冲到终端。
我以前犯过这种傻——直接跑了20分钟，屏幕刷了满屏的SAM行，最后啥也没存下来。
所以记得：输出要重定向到文件，或是用管道传给samtools转BAM。
一个典型的BWA比对命令长这样：
# 建立索引（只需一次）
bwa index GRCh38.fasta
# 双端比对：生成 SAM，实时转 BAM，排序
bwa mem -t 8 GRCh38.fasta \
sample_R1.fastq.gz sample_R2.fastq.gz | \
samtools sort -@ 4 -o sample.sorted.bam -这里的关键点：
-t 8：线程数，根据你服务器CPU核数设。samtools sort：比对完直接排序，方便GATK后续处理。-o sample.sorted.bam：指定输出文件。
GATK：变异检测的黄金标准
GATK全称是Genome Analysis Toolkit。
由Broad Institute开发，是目前学术和临床上使用最广的变异检测工具包。
但老实说，它也是一个参数多、流程长、容易出错的工具。
GATK 4把核心流程整合到了几个主要工具里：
- MarkDuplicates：标记PCR重复。
- BaseRecalibrator + ApplyBQSR：碱基质量分数重校正（BQSR）。
- HaplotypeCaller：核心的变异检测工具。
- VariantFiltration / VQSR：硬过滤或训练模型过滤。
我在2018年处理某个千人基因组WGS项目时，没做BQSR这步。
结果HaplotypeCaller出来的变异数量偏多，假阳性率明显偏高。
后来加上BQSR，precision提高了至少5%。
所以——千万别省这一步。
Best Practices Pipeline图示
下面这张图，是整个课程的核心流程：
我整个课程会怎么安排？
下面这张表，是1-8章的路线图：
| 章节 | 核心内容 | 涉及工具 |
|---|---|---|
| 第1章 | 课程介绍与背景 | BWA, GATK |
| 第2章 | 环境搭建：conda、容器、参考基因组 | conda, docker, bwa, GATK |
| 第3章 | BWA索引构建与比对实战 | bwa, samtools |
| 第4章 | 比对后处理：排序、标记重复 | samtools, GATK MarkDuplicates |
| 第5章 | BQSR原理与实操 | GATK BaseRecalibrator, ApplyBQSR |
| 第6章 | HaplotypeCaller单样本检测 | GATK HaplotypeCaller |
| 第7章 | 多样本联合分析 | GATK GenomicsDBImport |
| 第8章 | VCF过滤与质控 | GATK VariantFiltration, bcftools |
你看，整个流程是递进式的。
每一步的输出就是下一步的输入。
所以我的建议是：按顺序来，别跳。
避免后面遇到问题时，搞不清是哪个环节出了错。
小贴士：
如果你准备用自己的数据跑，建议先拿一个小区域（比如chr22）跑通。
我就是这么干的，省了不少时间和CPU资源。
常见疑问与避坑指南
我总结了一些刚开始学这套流程时，容易碰到的问题。
这里统一回答：
- 为什么我的bwa比对率低？
很可能原因：接头污染、低质量reads、选择的参考基因组版本与样本不匹配。 - GATK跑起来太慢，怎么办？
试试只跑目标区域（如exome.bed），或者用内存更大的服务器。
我记得第一次跑30x全基因组时，把实验室的机器活活跑了三天。 - 需不需要做BQSR？
要，强烈建议。我用过很多次，对SNP和Indel的假阳性抑制很明显。
避坑：版本兼容性
我之前遇到过：GATK 4.1 对 Java 版本有要求。
如果你用的是老服务器，先检查java -version。
GATK 4.1需要Java 8+，GATK 4.2需要Java 11+。
不然启动就会报 Java Exception，很烦人。
下一步做什么？
好了，背景和路线图都讲清楚了。
你可以先消化一下上面这些概念。
如果现在让你手动安装这一整套工具，可能还有点懵。