比对:把 reads 放回基因组
测序仪输出的是 FASTQ:一条条没有位置信息的序列和质量值。比对把每条 read 放回参考基因组上它最可能来自的位置,得到 BAM 文件。之后还有两道清理:标记 PCR 重复,校正碱基质量。变异检测读取的就是这份清理过的 BAM。
比对:给每条 read 找到位置
GATK 最佳实践用 BWA-MEM 把 reads 比对到参考基因组。新项目多用 GRCh38,较早的项目用 GRCh37,TCGA-LAML 的 MAF 就是 GRCh37 坐标。双端 reads 成对比对:两条方向相对、间距合理,就是一对正常配对的 reads。结果写成 SAM/BAM(或更紧凑的 CRAM),按坐标排序并建立索引,之后才能快速取出任一区域的 reads。
每条 read 都有一个比对质量 MAPQ,表示比对位置可信的程度。基因组里的重复序列、假基因和同源基因家族会让短 read 有好几个同样合适的位置,BWA 对这类 read 给出 MAPQ 0。变异检测会忽略 MAPQ 很低的 reads,所以这些区域在 WES 里实际上是“看不见”的。
标记重复与碱基质量校正
标记重复。 MarkDuplicates 找出两端比对起点和方向都相同的 read 对,认为它们来自同一个原始分子,只保留质量最好的一对,其余打上重复标记(并不删除)。光学重复也在这一步识别。重复率就是在这里算出来的。带 UMI 的文库可以区分“恰好同起点的不同分子”,标准 WES 大多不用 UMI。
碱基质量校正。 测序仪报告的碱基质量有系统偏差。BQSR 在已知变异位点之外统计实际的错配,按读段组、测序循环和相邻碱基建模,再重写每个碱基的质量值。Mutect2 依据碱基质量给每条 read 的证据加权,所以校正后的质量值更可信。
读较早的文献时还会看到 indel 局部重比对这一步。GATK4 已不再需要它,因为 Mutect2 会在每个活跃区域内重新做局部组装。
BAM 里与 QC 有关的信息
| BAM 中的信息 | 由哪一步产生 | 对应的 QC 指标 |
|---|---|---|
| 比对位置与 CIGAR | BWA-MEM 比对 | 比对率、soft-clip 比例 |
| MAPQ | 比对位置是否唯一 | 低 MAPQ reads 比例 |
| TLEN(模板长度) | 双端 reads 的间距 | 插入片段长度分布 |
| 重复标记(flag 1024) | MarkDuplicates | 重复率 |
| 碱基质量 | 测序仪给出,BQSR 校正 | ≥Q30 碱基比例 |
| 读段组 RG | 样本、文库、lane 信息 | 区分批次,供去重和 BQSR 使用 |
soft-clip 比例偏高,常见原因是插入片段太短、读进了接头,也可能是结构变异的断点。
常见误区
- 参考基因组版本混用。 GRCh37 与 GRCh38 的坐标不同,染色体名也可能有、无 “chr” 前缀。BAM、BED、PoN 和注释必须使用同一版本。
- 肿瘤和正常用不同流程处理。 两者应使用同一参考基因组、同一比对软件和同样的参数,否则流程差异本身就会制造假的“体细胞”信号。
- 以为低 MAPQ 区域没有突变。 那里只是看不见。PMS2 这类有高度同源假基因的基因需要专门的方法。