📚 课程目录第4章变异检测流程简介

4. 变异检测流程简介

这几年带团队做变异检测项目，我发现自己经常需要跟新人解释同一件事——变异检测到底是怎么一步步跑下来的。流程看起来就那么几步，但每一步踩过的坑，真是说都说不完。这一章我把整个流程拆开来讲，顺便聊聊我这些年攒下来的一些实战经验。

4.1 整个流程长什么样

先说个大框架。变异检测不是上来就找变异，你得先把原始测序数据洗得干干净净，再比对到参考基因组上，然后才能开始检测。我习惯把流程分成四大块：原始数据预处理 → 序列比对与排序 → 变异检测 → 结果过滤与输出。

下面这张图是我自己画的一个流程概览，你把它贴在工位旁边，走流程时碰到什么问题，瞄一眼就知道卡在哪一步了。

你仔细看这个图，左侧贴了两个提示——那都是我真实碰到过的事。尤其是预处理阶段，你以为只是去掉接头、修一下低质量碱基？没那么简单，参数设不对，后面全是垃圾数据。

4.2 每一步在做什么

4.2.1 原始数据预处理

这一步的目标很单纯：把机器吐出来的原始信号变成干净、可用的序列。具体包括去接头（adapter trimming）、质量修剪（quality trimming）、以及去除PCR重复。

我之前有个项目，测序深度做到60×，结果检测出来的变异数量比预期多了两倍。排查到最后才发现——接头残留太多，导致比对时大量read没比对上去，软件就把这些未比对的序列当成了"新变异"。你说冤不冤？

我的建议：预处理阶段多花10分钟调参数，后面能省出10小时的重跑时间。尤其是去接头这一步，一定要用实际数据的接头序列做验证，不要只看软件默认配置。

4.2.2 序列比对与排序

把清洗后的read映射到参考基因组上，这一步的工具有很多选择。我个人用得最多的是 BWA-MEM 和 Bowtie2。选哪个？主要看你的数据类型。

工具	适用场景	速度	准确率
BWA-MEM	长read（≥70bp）、人类全基因组	中	高
Bowtie2	短read、RNA-seq、细菌基因组	快	中
minimap2	长读长测序（PacBio / ONT）	中	高

比对完成后，必须做排序。很多人觉得排序就是排个序嘛，随便做做就好。嗯，我刚开始也是这么想的。直到有一次我跳过了排序直接跑变异检测，GATK 直接报了个让人摸不着头脑的错误。后来才发现它要求输入的BAM文件必须按坐标排序——不然死活不认。

# 用SAMtools排序
samtools sort -@ 8 -o sample.sorted.bam sample.bam

# 建立索引
samtools index sample.sorted.bam

4.2.3 标记重复与碱基质量校正

这一步特别容易被忽视，但恰恰是决定变异检测准确性的关键环节。PCR扩增会产生大量重复read，如果不标记掉，这些重复会误导变异频率的计算，让你以为某个变异位点的支持度很高——实际上只是PCR偏好造成的假象。

GATK 的 MarkDuplicates 是社区标配。至于 BQSR（碱基质量分数校正），我建议不是必需但别省。我曾经在一个肿瘤-正常配对项目里，没跑BQSR直接做变异检测，结果假阳性率比跑过的高了将近一倍。

# 标记重复（GATK）
gatk MarkDuplicates \
    -I sample.sorted.bam \
    -O sample.markdup.bam \
    -M sample.markdup.metrics.txt

# 碱基质量校正（GATK BQSR）
gatk BaseRecalibrator \
    -R reference.fa \
    -I sample.markdup.bam \
    --known-sites known_variants.vcf \
    -O recalibration.table

gatk ApplyBQSR \
    -R reference.fa \
    -I sample.markdup.bam \
    --bqsr-recal-file recalibration.table \
    -O sample.recalibrated.bam

注意：BQSR需要提供已知变异位点文件（比如dbSNP）。如果没有已知位点，这一步可以跳过，但要清楚这样做的后果——质量分数校准不准，变异检测的准确率会有所下降。

4.3 核心检测策略的对比

现在主流的变异检测策略分两类：基于比对（pileup-based） 和 基于组装（assembly-based）。怎么选？说说我的理解。

策略	代表工具	优势	劣势
基于比对	GATK UnifiedGenotyper、bcftools mpileup	速度快、内存低、适合SNP检测	对结构变异不敏感
基于组装	GATK HaplotypeCaller、Platypus	能检测复杂变异、InDel准确率高	计算资源消耗大
混合策略	FreeBayes、Streika2	兼顾速度和灵敏度	配置复杂

我个人目前的主力是 GATK HaplotypeCaller。它在局部区域做de novo组装，既能检出SNP又能检出InDel，而且对结构变异的检测能力比传统的pileup方式强很多。当然，代价就是慢。曾经有一个全基因组项目，32核服务器跑了整整两天才跑完。

实际项目中的选择经验：

如果是全基因组SNP筛查，用bcftools mpileup + call，快且够用。
如果是肿瘤-正常配对、稀有变异筛查，用GATK HaplotypeCaller或Streika2。
如果是长读长测序，直接上minimap2 + longshot或Clair3。

4.4 流程质量控制：你得盯着这几个数

流程跑完了，你怎么知道结果靠不靠谱？光看VCF文件里的数据量是不够的。我每次跑完流程都会盯着下面这几个指标，一条一条过。

比对率（Mapping Rate）：正常情况下，人类全基因组数据的比对率应该在95%以上。低于90%的，基本可以断定数据质量或比对参数有问题。
覆盖度（Coverage Depth）：检测SNP至少需要10×的覆盖度，检测稀有变异建议20×以上。覆盖度太低，检出率会断崖式下跌。
重复率（Duplicate Rate）：PCR重复率超过15%就要警惕了。我遇到过重复率高达40%的样本，最后只能重新建库测序。
过渡到颠换比（Ti/Tv Ratio）：人类全基因组的Ti/Tv比大约是2.0~2.1。如果远低于这个值，说明检测流程可能引入了大量假阳性。
已知位点检出率：用dbSNP中的已知位点做对照，看看你的流程能检出多少。低于80%的，得回去检查参数了。

一个实用的小脚本：跑完流程后，用下面的命令快速获取关键QC指标。

# 比对率
samtools flagstat sample.markdup.bam | grep "mapped (%"

# 平均覆盖度
samtools depth sample.recalibrated.bam | awk '{sum+=$3} END {print sum/NR}'

# 重复率
grep "PERCENT_DUPLICATION" sample.markdup.metrics.txt | tail -1

# Ti/Tv比
bcftools stats sample.vcf.gz | grep "TSTV"

4.5 输出格式与数据交付

变异检测的最终结果一般是 VCF（Variant Call Format） 格式。你也可能会遇到BCF（二进制的VCF）或者gVCF（基因组范围的VCF，包含所有位点的信息）。

VCF文件的核心字段包括：CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO、FORMAT。其中QUAL和FILTER两个字段直接影响后续的变异筛选。QUAL是变异质量分数，值越高代表可信度越高；FILTER是过滤标记，PASS表示通过过滤。

# 一个典型的VCF记录行
#CHROM  POS   ID      REF   ALT   QUAL   FILTER   INFO
chr21   941   rs123   G     A     523.77 PASS     DP=128;AF=0.445
chr22   1589  .       T     G     89.23  LowQual  DP=32;AF=0.0625

这里有个容易混淆的点：QUAL高分不等于变异一定真实。它只是代表该变异在当前数据下的支持度比较高。但如果是系统误差（比如测序偏好、比对错误），QUAL再高也是假阳性。所以后续的过滤和校准是必不可少的步骤——这也是我们后面几章要重点讲的内容。

总结一下这一章的核心内容：

变异检测流程分四步：预处理 → 比对 → 处理（排序/去重/BQSR） → 检测
每一步都有对应的工具和参数，别跳步骤，也别盲信默认值
检测策略的选择要根据数据类型和项目目标来定，没有"万能方案"
QC指标是判断流程是否跑偏的"仪表盘"——养成跑完就看的好习惯

流程搭好之后，下一步就是看检测结果准不准了。ROC曲线、混淆矩阵、F1-score……这些评估手段才是你判断流程好坏的"裁判"。别急，接下来几章我会带你一个一个过。