02 比对与变异检测流程

从 FASTQ 到 VCF 的标准流水线。这一章给出完整的 bash 命令序列，适合在服务器上跑。

BWA-MEM 比对

# 建索引（一次性）
bwa index reference.fa

# 比对（每个样本）
bwa mem -t 8 -R "@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA" \
  reference.fa sample1_R1.fq.gz sample1_R2.fq.gz \
  | samtools sort -@ 4 -o sample1.sorted.bam

samtools index sample1.sorted.bam

-R 参数加 read group 信息，GATK 后续步骤必须有。

去重复 + BQSR

# 标记 PCR 重复
gatk MarkDuplicates \
  -I sample1.sorted.bam \
  -O sample1.dedup.bam \
  -M sample1.dedup_metrics.txt

# Base Quality Score Recalibration
gatk BaseRecalibrator \
  -R reference.fa \
  -I sample1.dedup.bam \
  --known-sites dbsnp.vcf.gz \
  --known-sites mills_indels.vcf.gz \
  -O sample1.recal_table

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

HaplotypeCaller（胚系变异）

# 单样本模式（输出 gVCF）
gatk HaplotypeCaller \
  -R reference.fa \
  -I sample1.recal.bam \
  -O sample1.g.vcf.gz \
  -ERC GVCF

# 多样本联合分型
gatk CombineGVCFs -R reference.fa \
  -V sample1.g.vcf.gz -V sample2.g.vcf.gz \
  -O cohort.g.vcf.gz

gatk GenotypeGVCFs -R reference.fa \
  -V cohort.g.vcf.gz \
  -O cohort.vcf.gz

Mutect2（体细胞变异）

gatk Mutect2 \
  -R reference.fa \
  -I tumor.recal.bam \
  -I normal.recal.bam \
  -normal normal_sample \
  --germline-resource gnomad.vcf.gz \
  -O somatic_unfiltered.vcf.gz

gatk FilterMutectCalls \
  -R reference.fa \
  -V somatic_unfiltered.vcf.gz \
  -O somatic_filtered.vcf.gz

流水线工具

手动跑上面这些命令容易出错。推荐用 nf-core/sarek：

nextflow run nf-core/sarek \
  --input samplesheet.csv \
  --genome GRCh38 \
  --tools mutect2,haplotypecaller \
  --outdir results/

为什么 GATK Best Practices 这么麻烦

新手第一次跑这套常常想问"为什么要这么多步？" 每一步都是为了解决一个具体的假阳性来源：

步骤	解决什么	跳过会怎样
MarkDuplicates	PCR 扩增产生的重复 reads	duplicate 让 VAF 估计偏高，假阳性 ↑
BQSR	测序仪 quality score 系统偏差	假阳性 ↑10-30%，漏报 ↑
HaplotypeCaller local assembly	重复区域 / indel 处比对错误	indel 假阴性高
gVCF + GenotypeGVCFs	多样本联合分型校准信息	"我有变异但只在一个样本里" 漏报
VQSR / hard filter	工具固有假阳性	报告里 90% 是噪声

每一步都不是装饰。理解每一步在防什么，比死记参数更有价值。

常见坑

坑 1：忘了 read group 信息

GATK 后续步骤强制要求 BAM 里有 @RG 标签。bwa mem -R "@RG\tID:..." 一开始就要加。否则 MarkDuplicates 直接报错。后期补 RG 用 gatk AddOrReplaceReadGroups。

坑 2：BWA 用了 -M 标志

老教程用 -M 让多比对 reads 标记为 secondary，但新版 GATK + Picard 期望 supplementary。用 -M 在新版 pipeline 里会让 MarkDuplicates 行为异常。不加 -M 是默认正确。

坑 3：GVCF 文件多样本时手动 cat

GVCF 是结构化文件，不能直接 cat。必须用 gatk CombineGVCFs 或 GenomicsDBImport，否则下游报错。

坑 4：Mutect2 不指定 germline-resource

--germline-resource gnomad.vcf.gz 让 Mutect2 把已知胚系变异从体细胞结果里去掉。不加这个参数 tumor-normal 配对的结果还能凑合，tumor-only 模式直接错。

坑 5：BAM 文件没 sort 就 dedup

MarkDuplicates 要求 BAM 已按坐标排序。bwa mem 输出的 SAM 是 read-name 顺序，要先 samtools sort -@ 4。直接 dedup 会报错或给出错误结果。

下一步

接着深入：

• 03 VCF 注释与可视化 — VCF 拿到手后的第一步
• 05 变异过滤与质量评估 — 决定结果可信度的关键

横向延伸：

• 04 maftools 肿瘤突变分析 — Mutect2 输出的 MAF 用 maftools 出图
• nf-core/sarek 教程 — 一键跑完整 pipeline

参考资源

• BWA 手册
• GATK Best Practices
• nf-core/sarek
• DeepVariant