3. 参考基因组与比对策略

长读长结构变异检测的第一步，不是测序，而是把读段（reads）放到一个可信的坐标参照系上。这个参照系就是参考基因组，而把读段映射到参考上的过程就是比对。参考选错、比对参数不当，后续无论用 Sniffles、cuteSV 还是 SVIM，都会得到大量假阳性或漏检。本章围绕三件事展开：参考基因组的选择与预处理、长读长比对算法的原理、以及比对参数调优与结果评估。

一、参考基因组的选择与预处理

参考基因组的选择要匹配样本的物种、亚种和倍性。以人类为例，GRCh38 与 T2T-CHM13 是当前两个主流选择。GRCh38 生态成熟、注释完善、与大量公共数据集兼容；T2T-CHM13 补齐了着丝粒、近端着丝粒短臂等约 8% 此前缺失的区域，对检测这些区域的结构变异更有利，但部分工具链和注释资源仍在适配。若研究对象是特定族群，还应考虑是否使用带 alternate contig 或族群特异组装的参考，避免把族群特异序列误判为变异。

预处理的核心是建立索引并屏蔽干扰区域。以 minimap2 为例，常用命令如下：

# 为 GRCh38 建立 minimap2 索引（-d 输出 .mmi）
minimap2 -d GRCh38.mmi GRCh38.fa

# 建立 faidx 索引，供后续工具按坐标取序列
samtools faidx GRCh38.fa

# 生成 .dict，供 GATK 等工具使用
gatk CreateSequenceDictionary -R GRCh38.fa

此外，应准备 BED 文件屏蔽高重复、低可比对区域（如 ENCODE blacklist），并在结果评估时排除这些区域，否则着丝粒、端粒附近的比对堆积会制造大量假 SV。参考序列中的 N 碱基比例、contig 命名一致性（chr1 与 1 混用）也需提前统一，否则会导致坐标错位。

二、长读长比对算法原理

长读长比对与短读长比对有本质差异。短读长依赖精确 k-mer 种子，而长读长（PacBio HiFi、ONT）错误率虽已降至 0.1%～5%，但读段长达数千至数十 kb，且含插入、缺失，因此需要能容忍错误的种子与链化策略。

minimap2 采用“minimizer 种子 + 链化 + 碱基级比对”三段式。它先在参考和读段上取 minimizer（局部最小的 k-mer），用哈希表找到匹配种子；再对种子做链化（chaining），把共线性的种子连成锚链；最后对锚链之间的空隙做动态规划比对。minimap2 提供多种预设：-ax map-ont 用于 ONT 读段，-ax map-hifi 用于 HiFi，-ax asm5/asm10/asm20 用于组装比对。其优势是速度快、内存低，适合大规模样本。

NGMLR 则针对结构变异检测做了优化。它使用“种子—链化—局部比对”流程，但在链化阶段采用更宽松的 gap 惩罚，并引入“SV-aware”的比对模型，使跨越 SV 断点的读段更容易被正确比对，而不是被切成软剪切（soft clip）。因此 NGMLR 在检测大片段插入、缺失时召回率往往更高，但速度慢于 minimap2。实际工程中常见策略是：用 minimap2 做快速初筛和覆盖度评估，用 NGMLR 对候选区域或全样本做精细比对。

长读长比对算法流程与 SV 断点处理对比 minimap2（速度优先） ① minimizer 种子提取（k-mer 局部最小） ② 种子链化（chaining，共线性锚链） ③ 锚链间动态规划碱基级比对 SV 断点处：易产生 soft clip 或 分裂比对，需下游工具重组 NGMLR（SV 敏感优先） ① 种子提取 + 宽松 gap 惩罚 ② SV-aware 链化（容忍大 gap） ③ 局部比对 + 断点精修 SV 断点处：读段跨越断点被完整比对， 插入/缺失信号更清晰 工程策略：minimap2 初筛覆盖度 → NGMLR 精细比对候选区域

三、比对参数调优

参数调优的目标是让真实 SV 信号在比对结果中保留，同时抑制重复区域的假信号。以下参数最值得关注：

  • -x 预设：决定 k-mer 长度、gap 惩罚等基础参数。ONT 用 map-ont，HiFi 用 map-hifi，切勿混用。
  • -k / -w：k-mer 长度与 minimizer 窗口。k 越大越特异，但低复杂度区域种子越少；w 越小种子越密、越慢。
  • -a：输出 SAM/BAM，SV 检测必须加。
  • --MD：在 BAM 中保留 MD 标签，便于后续工具还原错配。
  • -t：线程数，影响速度不影响结果。
  • -R：读组标签，多样本合并时用于区分来源。

一个典型的 ONT 比对命令如下：

minimap2 -ax map-ont -t 16 --MD \
  -R '@RG\tID:sample1\tSM:sample1' \
  GRCh38.mmi sample1.fastq.gz \
  | samtools sort -@ 8 -o sample1.sorted.bam
samtools index sample1.sorted.bam

若发现某区域覆盖度异常低或 soft clip 比例异常高，可尝试降低 -k 或增大 -w 的密度，或改用 NGMLR 重比对。对于已知的高重复区域，应结合 BED 屏蔽，而不是一味调参。

四、比对结果评估

比对完成后，不能直接进入 SV calling，必须先评估质量。关键指标包括：

指标含义合理范围（示例）异常处理
比对率成功比对读段占比ONT ≥ 90%，HiFi ≥ 98%过低检查参考版本、污染
平均覆盖度有效比对碱基数 / 基因组长度SV 检测建议 ≥ 15×不足则补测序
MAPQ 分布比对质量值分布MAPQ≥20 占比高低 MAPQ 多提示重复区
soft clip 比例被剪切碱基占比通常 < 10%过高提示 SV 或参考缺失
插入片段长度读段比对跨度与文库一致异常提示嵌合或结构变异

可用 samtools flagstat、samtools stats 和 mosdepth 快速获得上述指标：

samtools flagstat sample1.sorted.bam
samtools stats sample1.sorted.bam | grep ^SN
mosdepth -t 4 -b blacklist.bed sample1 sample1.sorted.bam
评估时建议同时查看覆盖度曲线和 soft clip 分布。若某区域覆盖度骤降且 soft clip 骤升，往往正是候选 SV 所在，应保留该区域进入下游 calling，而不是简单过滤掉。
不要用短读长比对的默认参数直接跑长读长数据。短读长参数 k 小、gap 惩罚严，会把长读段中的真实插入/缺失切成大量 soft clip，导致 SV 漏检。

总结来说，参考基因组决定坐标基准，比对算法决定信号保留程度，参数调优决定信噪比，结果评估决定是否值得继续。四者环环相扣，任何一环失守，后续 SV 检测都会事倍功半。建议在正式 calling 前，先用一个已知 SV 的阳性样本（如 NA12878 或 HG002）跑通全流程，确认比对指标和 SV 召回率符合预期，再进入大规模样本分析。

← 上一章 📖 返回目录 下一章 →