3. 比对与变异检测：BWA/GATK标准流程，VCF文件的生成与基本过滤

测序数据下机之后，你手里握着一堆 FASTQ 文件。下一步干嘛？当然是找到这些 reads 在基因组上的位置，然后揪出变异位点。这套流程我跑了不下两百次，从最初的生信小白到现在闭着眼都能写脚本，踩过的坑能绕实验室一圈。

今天咱们就把 BWA 比对和 GATK 变异检测的底裤扒干净。说白了，就是三步：把 reads 怼到参考基因组上 → 找差异 → 生成 VCF 文件。

3.1 BWA 比对：把 reads 放回基因组地图

BWA 是目前最主流的短序列比对工具。我当年刚入行时用过 bowtie2，后来全面转向了 BWA-mem，为什么？因为它对 gap 的处理更聪明，特别是对于长 reads 和存在结构变异的区域。

先建索引。这一步不能省，而且 索引文件必须和参考基因组版本严格对应。我吃过一次亏：用 hg19 的索引匹配了 hg38 的参考序列，结果一堆 reads 比对不上，浪费了三天。

# 下载参考基因组（以 GRCh38 为例）
wget ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 建立 BWA 索引
bwa index -a bwtsw Homo_sapiens.GRCh38.dna.primary_assembly.fa

比对命令很简单，但参数别乱改。我个人习惯用默认参数，除非你明确知道改动的后果。

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

参数	作用	我的建议
-t	线程数	根据服务器核数设定，通常 8-16
-R	Read Group 信息	必须包含 ID、SM、PL，不然 GATK 不认
-M	将 split reads 标记为 supplementary	做 RNA-seq 时用，DNA 一般不推荐

3.2 SAM → BAM：压缩与排序

SAM 文件是文本格式，占地方还慢。所有下游工具都要求 BAM 格式。先转格式，再排序。

# SAM → BAM
samtools view -Sb sample1.sam > sample1.bam

# 排序（按坐标）
samtools sort sample1.bam -o sample1.sorted.bam

# 建索引
samtools index sample1.sorted.bam

你想想看，这里有个很多人忽略的点：排序后的 BAM 必须建立索引 (.bai)，不然 IGV 看不了，GATK 也跑不动。我曾经帮同事排查问题，他的 GATK 总是报 "No index found"，其实就是这步漏了。

3.3 GATK 流程：打磨 BAM 文件

GATK 的变异检测流程，我把核心步骤列出来。这四步是标准套餐，少了哪一步都有可能引入假阳性。

1. MarkDuplicates：标记 PCR 重复。说白了就是同一模板扩出来的多条 reads，只保留一条。
2. FixMateInformation：修复双端 reads 的配对信息，确保两端正确关联。
3. BaseRecalibrator + ApplyBQSR：碱基质量校正。GATK 根据已知 SNP 位点重新校准 base quality score，这一步能大幅提升变异检测准确性。
4. HaplotypeCaller：核心变异检测，输出 GVCF 格式。

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

# 碱基质量校正（已知 SNPs 位点文件）
gatk BaseRecalibrator \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -I sample1.dedup.bam \
    --known-sites dbsnp_151.vcf.gz \
    -O sample1.recal.table

gatk ApplyBQSR \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -I sample1.dedup.bam \
    --bqsr-recal-file sample1.recal.table \
    -O sample1.recal.bam

3.4 HaplotypeCaller：变异检测核心

HaplotypeCaller 是 GATK 最强模块。它不像老方法那样直接比对 reads，而是先在目标区域构建局部单倍型，再判断可能的变异。这样做的好处很明显：对复杂区域（比如 indels 密集区）的检测准确率远高于旧算法。

# 单样本 GVCF 模式
gatk HaplotypeCaller \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -I sample1.recal.bam \
    -O sample1.g.vcf.gz \
    -ERC GVCF

# 如果要做多个样本的 joint genotyping，这个 GVCF 是必须的

参数方面，-ERC GVCF 是推荐选项。它把所有位点的信息都输出（包括非变异位点），这是后续联合分析的基础。如果你只输出 VCF，那只有变异位点的信息，很多统计模型就没法算了。

3.5 VCF 文件生成与基本过滤

单个样本的 GVCF 做完，要汇聚成群体 VCF。这才是大家期待的最终产物。

# 合并 GVCF
gatk GenomicsDBImport \
    -V sample1.g.vcf.gz \
    -V sample2.g.vcf.gz \
    -V sample3.g.vcf.gz \
    --genomicsdb-workspace-path my_database \
    -L intervals.list

# 联合基因分型
gatk GenotypeGVCFs \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V gendb://my_database \
    -O cohort.vcf.gz

VCF 文件生成之后，别急着交差。原始 VCF 里充斥着各种假阳性，必须做过滤。GATK 推荐用 VariantFiltration 做硬过滤，或者用 VQSR 做软过滤。

# 硬过滤示例
gatk VariantFiltration \
    -R Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    -V cohort.vcf.gz \
    -O cohort.filtered.vcf.gz \
    --filter-expression "QD < 2.0" --filter-name "LowQD" \
    --filter-expression "FS > 60.0" --filter-name "HighFS" \
    --filter-expression "MQ < 40.0" --filter-name "LowMQ" \
    --filter-expression "MQRankSum < -12.5" --filter-name "LowMQRankSum" \
    --filter-expression "ReadPosRankSum < -8.0" --filter-name "LowReadPosRankSum"

过滤标准	说明	建议阈值
QD	变异质量值 / 深度，低 QD 通常表示支持证据不足	< 2.0
FS	链偏好 Fisher 检验，高 FS 表示 reads 偏向一条链	> 60.0 (SNP), > 200.0 (INDEL)
MQ	比对质量，低 MQ 表示 reads 比对到重复区域	< 40.0