参考基因组与变异检测基础

开始做变异检测之前，有几件事你得先搞清楚。其中之一就是参考基因组——说白了，它就是一条“标准答案”。咱们拿它跟测序数据比，才能知道你在哪个位置上有差异。

参考基因组到底是什么？

参考基因组是一个物种的代表性DNA序列组合。人类参考基因组从2001年发布第一版到现在，已经迭代了好多版本。我个人习惯用GRCh38（也叫hg38），但有些老项目还在用GRCh37（hg19）。这两种版本在染色体坐标、补丁序列上有差异，千万别混用。我在项目中遇到过两次因为hg19和hg38混用造成上下游分析对不上的情况，非常痛苦。

你想想看，如果你用hg19的参考基因组装好了数据，后来想跟别人hg38的结果比对，坐标偏移、字符串差异会让你怀疑人生。所以进入任何一个新项目，第一件事就是确认参考基因组版本。

变异检测的基本流程

从原始FASTQ到最终VCF，中间要过好几道关卡。我一般把它拆成五个步骤：

1. 比对 – 把测序读段怼到参考基因组上，常用BWA-MEM或BWA-MEM2。
2. 排序与去重 – 用samtools或Picard按坐标排序，然后标记或去除PCR重复。
3. 碱基质量重校准（BQSR） – 靠已知SNP对碱基质量值做机器学习修正，GATK4的BaseRecalibrator干这个。
4. 变异检测 – 根据比对结果调用变异，GATK HaplotypeCaller或bcftools mpileup。
5. 过滤与注释 – 根据质量、深度、链偏倚等过滤假阳性，再用VEP/ANNOVAR做功能注释。

嗯，这里要注意：去重和BQSR的顺序不要搞反。我曾经因为先做BQSR后去重，导致重复读段里堆积的假质量问题被BQSR“学到”，结果变异检出的特异性一塌糊涂。正确做法是：比对 → 排序 → 去重 → BQSR → 变异检测。

变异类型你必须要知道

变异不是只有一种。我见过不少新手把SNP和Indel混在一起直接算准确率，结果模型一团糟。常见的分类如下：

咱们这章主要关注SNP和Indel，因为后续做ROC曲线时，这两种变异的假阳性/真阳性比较好定义。结构变异和CNV需要单独评估方法，暂时不做展开。

常用工具与命令

我手边常备这几个工具链。下面给一个最基本的命令行示意（最简版，实际参数会根据测序深度调整）：

# 比对
bwa mem -t 8 reference.fasta sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools sort -@4 -o sample.sorted.bam

# 标记重复
gatk MarkDuplicates \
  -I sample.sorted.bam \
  -O sample.marked.bam \
  -M sample.metrics

# 碱基质量重校准（先建一个已知变异集，如dbSNP）
gatk BaseRecalibrator \
  -R reference.fasta \
  -I sample.marked.bam \
  --known-sites dbsnp_138.vcf \
  -O recal.table

gatk ApplyBQSR \
  -R reference.fasta \
  -I sample.marked.bam \
  --bqsr-recal-file recal.table \
  -O sample.recal.bam

# 变异检测
gatk HaplotypeCaller \
  -R reference.fasta \
  -I sample.recal.bam \
  -O sample.g.vcf.gz \
  -ERC GVCF

注意：GATK4的HaplotypeCaller输出的是gVCF（基因组范围的VCF），下一步还需要联合基因分型。如果你只想快速看单样本变异，可以去掉 -ERC GVCF 参数。

本章的知识体系图

下面这张SVG图概括了变异检测从参考基因组到VCF的完整流程，以及关键节点。你可以把它当作学习路线图。