第二章 GATK最佳实践流程概览：从原始FASTQ到VQSR的完整步骤

做变异检测，很多人一上来就想着跑VQSR，但前面的路没走扎实，后面调参数基本是白费功夫。我刚开始接触GATK时也犯过这错：直接把FASTQ丢进HaplotypeCaller，出来的VCF过滤完后一看，假阳性能占一半。说白了，整个流程就像搭积木，每一块都得对齐。

这一章我带你从头过一遍，从FASTQ到VQSR的完整步骤。别嫌啰嗦，每个环节都藏着坑。我当年踩过的，你也能避一避。

2.1 数据预处理：从FASTQ到分析就绪的BAM

原始FASTQ不能直接拿来call变异，得先做三件事：比对、标记重复、碱基质量重校准（BQSR）。我个人习惯在比对后先看一眼统计报告，比如用samtools stats扫一把，确认比对率正常再往下走。

2.1.1 比对 & 排序

主流工具是BWA-MEM。我建议加上-M标记，确保read group信息完整。输出后直接用SortSam或samtools sort排好坐标顺序。

bwa mem -M -t 8 -R '@RG\tID:Sample1\tSM:Sample1\tPL:ILLUMINA' \
  hg38.fasta Sample1_R1.fastq.gz Sample1_R2.fastq.gz \
  | samtools sort -@4 -o sorted.bam -

2.1.2 标记重复 & BQSR

重复标记用MarkDuplicates。BQSR不是强制步骤，但如果你要跑VQSR，我个人建议做——它能降低系统误差，让VQSR的模型更干净。

gatk MarkDuplicates \
  -I sorted.bam \
  -O marked_dups.bam \
  -M metrics.txt

gatk BaseRecalibrator \
  -R hg38.fasta \
  -I marked_dups.bam \
  --known-sites dbsnp_151.hg38.vcf.gz \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  -O recal_data.table

gatk ApplyBQSR \
  -R hg38.fasta \
  -I marked_dups.bam \
  --bqsr-recal-file recal_data.table \
  -O ready.bam

2.2 变异检测与联合基因分型

这一步用HaplotypeCaller对每个样本生成gVCF，然后合并、联合基因分型。我见过有人图快，不生成gVCF直接单样本call，但后续要合并家族或人群数据就很麻烦。

2.2.1 HaplotypeCaller

gatk HaplotypeCaller \
  -R hg38.fasta \
  -I ready.bam \
  -L interval_list.bed \
  --emit-ref-confidence GVCF \
  -O Sample1.g.vcf.gz

2.2.2 合并与基因分型

合并时推荐CombineGVCFs（样本少）或GenomicsDBImport（样本多）。

# 样本少的情况
gatk CombineGVCFs \
  -R hg38.fasta \
  -V Sample1.g.vcf.gz -V Sample2.g.vcf.gz \
  -O cohort.g.vcf.gz

# 再联合基因分型
gatk GenotypeGVCFs \
  -R hg38.fasta \
  -V cohort.g.vcf.gz \
  -O raw_variants.vcf.gz

2.3 硬过滤：VQSR前的“筛子”

VQSR默认需要训练集，如果训练集质量太差，模型会学歪。我建议先用GATK推荐的硬阈值过滤一轮，去掉极低质量的变异。这一步不是最终过滤，而是给VQSR提供干净的种子。

gatk VariantFiltration \
  -R hg38.fasta \
  -V raw_variants.vcf.gz \
  --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
  --filter-name "hard_filter" \
  -O hard_filtered.vcf.gz

2.4 VQSR的输入准备

VQSR需要三个输入：待过滤VCF、已知真集（如HapMap、Omni、1000G）、已知假阳性集（如dbSNP）。这些资源文件下载好后，用VariantAnnotator给VCF加上必要的注释，比如VQSLOD。我一般还会检查已知集和待过滤集的染色体版本是否一致，否则VQSR会直接崩溃。

# 添加必要注释（如果缺的话）
gatk VariantAnnotator \
  -R hg38.fasta \
  -V hard_filtered.vcf.gz \
  --dbsnp dbsnp_151.hg38.vcf.gz \
  -O annotated.vcf.gz

2.5 运行VQSR & 应用过滤阈值

这是整套流程的最后一段。先训练模型，再用训练好的模型对每个变异打分，最后按阈值过滤。我习惯把阈值设在90%的灵敏度点（--sensitivity 90.0），但具体要看你的业务容忍度。

gatk VariantRecalibrator \
  -R hg38.fasta \
  -V annotated.vcf.gz \
  --resource hapmap_3.3.hg38.vcf.gz,known=false,training=true,truth=true,prior=15.0 \
  --resource omni_3.1.hg38.vcf.gz,known=false,training=true,truth=false,prior=12.0 \
  --resource 1000G_phase1.snps.high_confidence.hg38.vcf.gz,known=false,training=true,truth=false,prior=10.0 \
  --resource dbsnp_151.hg38.vcf.gz,known=true,training=false,truth=false,prior=7.0 \
  --max-gaussians 4 \
  -O recalibrate_SNP.recal \
  --tranches-file recalibrate_SNP.tranches

gatk ApplyVQSR \
  -R hg38.fasta \
  -V annotated.vcf.gz \
  --recal-file recalibrate_SNP.recal \
  --tranches-file recalibrate_SNP.tranches \
  --truth-sensitivity-filter-level 90.0 \
  -O vqsr_filtered.vcf.gz

2.6 完整流程概览图

下面这张图把上面提到的步骤串起来。你可以把它贴到工位旁边，每次卡壳时看一眼就知道在哪步断联了。

2.7 步骤汇总表