GATK 最佳实践流程：从原始测序数据到 GVCF 的完整链路

说到从原始数据到GVCF的链路，我个人习惯先把整体流程在脑子里过一遍。这个流程不是拍脑袋想出来的，是Broad Institute经过多年沉淀定下来的标准做法。说白了，就是把一堆乱七八糟的原始测序数据，洗干净、对齐、打磨好，最后输出一份高质量的中间结果。

你先记住一句话：GVCF的质量，直接取决于前面每一步做没做干净。我在项目中见过太多因为跳过某一步导致最终结果无法使用的例子。嗯，下面我们就把这条链路拆开来看。

完整链路一览

从原始FASTQ到GVCF，一共需要走5个主要阶段。我习惯把它们分成数据准备、预处理、变异检测三大块：

阶段	步骤	核心输出	我用到的关键工具
数据准备	原始FASTQ → 质量控制	Clean FASTQ	FastQC, Trimmomatic
比对与排序	Clean FASTQ → BAM	Sorted BAM	bwa mem, samtools sort
标记重复	Sorted BAM → Dedup BAM	Marked BAM	GATK MarkDuplicates
碱基质量校正	Dedup BAM → Recal BAM	Recalibrated BAM	GATK BaseRecalibrator, ApplyBQSR
GVCF生成	Recal BAM → GVCF	GVCF文件	GATK HaplotypeCaller

第一步：从原始测序数据到干净数据

拿到原始FASTQ，第一件事不是直接怼进去跑bwa。你先想想看，测序仪给出来的数据里带了多少杂质？接头污染、低质量碱基、还有那些读得太短的片段，这些东西不处理干净，后续比对全是坑。

我通常这样处理：

• 跑FastQC看一下整体质量分布
• 用Trimmomatic或cutadapt去掉接头
• 滑窗切掉质量低于Q20的尾巴
• 扔掉长度小于36bp的短读

第二步：比对到参考基因组

干净数据拿来，下一步就是比对。这里我用bwa mem，几乎是行业标配了。命令很简单：

bwa mem -t 8 -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \
  reference.fa sample1_R1.fastq.gz sample1_R2.fastq.gz \
  > sample1.sam

注意这里的Read Group一定要写对。我在项目中见过太多因为忘了加SM标签导致联合分析时样本混淆的情况。你想想看，如果每个样本的RG信息不完整，到后面GenotypeGVCFs阶段根本识别不了样本来源。

比对完的SAM文件太大，我会直接走管道转成BAM并排序：

bwa mem -t 8 -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \
  reference.fa sample1_R1.fastq.gz sample1_R2.fastq.gz \
  | samtools sort -@ 8 -o sample1.sorted.bam -

第三步：标记重复序列

为什么一定要做这一步？因为PCR扩增会引入重复读段，如果不标记掉，这些重复会严重高估覆盖深度，导致后续变异检测时假阳性率飙升。

我用GATK的MarkDuplicates来操作：

gatk MarkDuplicates \
  -I sample1.sorted.bam \
  -O sample1.dedup.bam \
  -M sample1.metrics.txt \
  --REMOVE_DUPLICATES false

嗯，这里要注意，我一般只标记不删除。为什么？因为HaplotypeCaller会自动处理这些标记，但如果你真删掉了，某些区域的深度信息就丢了。曾经有个项目，团队为了省空间直接删了重复，结果低覆盖区域的变异检测几乎全废了。

第四步：碱基质量分数重校准（BQSR）

这一步是GATK最佳实践里让我觉得最巧妙的设计。测序仪给的碱基质量分数有系统偏差，BQSR就是拿已知的SNP位点做校准，把这些偏差纠正过来。

分两步走。先建立校正模型：

gatk BaseRecalibrator \
  -R reference.fa \
  -I sample1.dedup.bam \
  --known-sites known_snp.vcf \
  --known-sites known_indel.vcf \
  -O sample1.recal.table

再应用到BAM上：

gatk ApplyBQSR \
  -R reference.fa \
  -I sample1.dedup.bam \
  --bqsr-recal-file sample1.recal.table \
  -O sample1.recal.bam

第五步：生成GVCF

重头戏来了。HaplotypeCaller是GATK里最核心的工具，它不像传统方法那样一个个位点半对齐，而是对整个区域做局部组装，输出每个位点的全部可能性。

命令其实很直接：

gatk HaplotypeCaller \
  -R reference.fa \
  -I sample1.recal.bam \
  -O sample1.g.vcf.gz \
  -ERC GVCF

-ERC GVCF这个参数是开关。不加的话，它只输出有变异的位点；加上之后，它会为每一个位点都输出一个参考置信度模型。你想想看，这意味着什么？意味着后面联合分析时，即使某个样本在某位点没测到变异，我们也能知道它这里到底有多"自信"是纯合参考。

整条链路的流程图

下面的SVG把整个流程串起来了。你可以看到从原始FASTQ到最终GVCF的每一步，以及每个阶段的关键产出文件。我个人喜欢这样用图来指导实际工作，每一步都对应着明确的输入输出。

每条链路上的避坑总结

在实际项目中跑这条链路，有几个点我反复栽过跟头，专门列出来：

• Read Group遗漏：@RG里的SM标签必须和样本名称严格对应，否则联合分析阶段报错
• 重复标记策略：标记不删除，保留reads上的duplicate flag给下游判断
• BQSR已知位点：只能用人种匹配的dbSNP，别混用版本
• HaplotypeCaller线程数：-Xmx给足内存，单个样本至少4G起步，WGS建议8G以上
• 中间文件保存：dedup.bam和recal.bam不要删，后续分析可能还要回溯