4. 长读长比对工具选型与实战
长读长测序（PacBio HiFi、ONT）产生的读段长度可达 10–100 kb，错误模式以插入/缺失为主，且存在局部高错误区。把这类读段比对到参考基因组，不能沿用短读长时代的 BWA-MEM 思路，而应选择面向长读长的种子链式比对器。本节聚焦 minimap2 与 Winnowmap 两类主流工具，讲清算法原理、参数调优、结果评估与可视化。
一、算法原理：从种子到链再到比对
minimap2 的核心是“种子—链—比对”三段式。它先用 minimizer 采样参考与读段，得到短精确匹配种子；再对种子按参考坐标与读段坐标做共线性链化（chaining），把落在同一对角带上的种子串成锚点链；最后对高分链做碱基级比对（base-level alignment），输出 CIGAR。对长读长，minimap2 默认使用 -x map-ont 或 -x map-hifi 预设，前者允许更高错误率与更长间隙，后者因 HiFi 准确率高而收紧参数。
Winnowmap 的关键改进在种子采样：它用“最小化加权采样”（weighted minimizer）对参考中高频出现的 minimizer 降权，从而抑制重复序列产生的假种子。对高度重复的基因组（如着丝粒、rDNA 簇），Winnowmap 的链化更干净，比对到重复区的读段不会因种子爆炸而错配到旁系同源位点。代价是它需要预先构建参考的重复权重文件（meryl 计数 + winnowmap 生成），流程比 minimap2 多一步。
二、参数调优：按数据类型与目标选择
选型的第一步是确定预设。ONT 原始读段用 -x map-ont，PacBio HiFi 用 -x map-hifi，若做剪接感知比对（RNA）用 -x splice。预设之外，最常调的是 -k（种子长度）、-w（minimizer 窗口）、-A（匹配得分）、-B（错配罚分）、-O（间隙开放）、-E（间隙延伸）。对高错误率 ONT，可适当降低 -A 或增大 -O，让比对器更容忍长间隙。
重复区比对是长读长最易出问题的地方。若发现读段在着丝粒或 rDNA 区出现大量次优比对（MAPQ 低、SA 标签多），可改用 Winnowmap，或对 minimap2 增加 --secondary=no 减少次优输出，并用 -N 控制保留的次优比对数量。对结构变异检测，建议保留次优比对（-N 5 左右），因为断裂点附近的读段常需要次优链来定位。
# minimap2 比对 ONT 读段到参考
minimap2 -t 16 -ax map-ont -N 5 --secondary=no \
ref.fa reads.fq.gz | samtools sort -@ 8 -o aln.bam
# Winnowmap 比对（需先建重复权重）
meryl count k=15 ref.fa output ref.meryl
winnowmap -W ref.meryl -ax map-ont ref.fa reads.fq.gz \
| samtools sort -@ 8 -o aln.winnowmap.bam三、比对结果评估：看什么指标
比对完成后，先看整体比对率与 MAPQ 分布。用 samtools flagstat 看 mapped 比例，用 samtools stats 看插入片段长度、错配率、覆盖度。长读长比对率通常应在 95% 以上（HiFi 可达 99%），若明显偏低，需检查参考版本、读段质量或预设是否匹配。
更细的评估看 CIGAR 与 MAPQ。高 MAPQ（≥30）且 CIGAR 以 M 为主，说明比对可信；若大量读段 MAPQ=0 且 CIGAR 含长 I/D，往往落在重复区或存在结构变异。可用 samtools view 抽取 MAPQ=0 的读段，统计其参考位置是否集中在重复注释区。对结构变异，重点看 SA 标签与软剪切（S）长度，软剪切超过 1 kb 常提示存在大片段插入或参考缺失。
| 指标 | 观察命令 | 健康范围（示例） | 异常提示 |
|---|---|---|---|
| 比对率 | samtools flagstat | HiFi ≥99%，ONT ≥95% | 参考版本错、读段污染 |
| MAPQ 分布 | samtools view | cut -f5 | 主峰 ≥30 | 重复区、多拷贝基因 |
| 软剪切长度 | samtools view | CIGAR 解析 | 多数 <500 bp | 大插入、参考缺失 |
| SA 标签数 | samtools view -h | grep SA | 少量 | 结构变异、重复 |
四、可视化：把比对结果画出来
评估之后需要可视化。常用工具包括 IGV（看单读段比对与覆盖度）、samtools depth（覆盖度曲线）、mosdepth（快速覆盖度统计）以及 IGV.js 做网页嵌入。对结构变异，可用 sniffles 或 cuteSV 从 BAM 中调用 SV，再用 IGV 查看断裂点附近的读段堆叠。
一个实用的可视化流程是：先用 mosdepth 生成覆盖度 BED，再用 bedtools genomecov 转成 bedGraph，最后用 IGV 加载 BAM 与 bedGraph。若发现某区域覆盖度骤降且读段软剪切集中，可判定为潜在缺失；若覆盖度骤升且读段出现大量 SA 标签，可判定为潜在重复或插入。
ref.meryl，否则降权会失效，导致重复区比对错误。最后，比对不是终点，而是结构变异检测、甲基化调用、组装纠错等下游分析的起点。评估时不要只看比对率，还要结合 MAPQ、CIGAR、SA 标签与覆盖度曲线综合判断。只有把比对结果放到具体生物学问题中看，才能判断参数是否调对、工具是否选对。