minimap2长读长序列比对——病毒三代测序数据的通用比对引擎

原理

三代测序（Nanopore、PacBio）读长轻松超过10 kb，但这东西的错误率曾经高达~15%。传统短读长比对工具（比如BWA-MEM）碰到这种又长又吵的数据，要么跑不动，要么结果一塌糊涂。minimap2就是为这个场景设计的——它用minimizer种子采样的方式，把比对拆成三步：收集种子（seeding）、共线性链接（chaining）、碱基级比对（alignment）。核心思路很简单：长读长意味着可跳过的重复区更多，没必要像短读长那样逐个碱基去死磕。

跟BWA-MEM比一下：同样是比对10 kb的Nanopore读长，minimap2快30倍以上，灵敏度还更高。因为BWA-MEM的算法假设读长就几百bp，碰上长片段结构变异就撒手不管了，而minimap2的凹形gap罚分函数天然能处理长插入/缺失。在病毒溯源里，minimap2是ARTIC扩增子分析流程的默认比对引擎——从Zika到SARS-CoV-2，凡是用Nanopore做现场病毒基因组测序的，几乎都跑过minimap2。在埃博拉实时基因组监测等溯源案例中，它把测序读长比对到参考基因组，为下游变异检测打下基础。

分析步骤

1. 安装

源码编译或通过Conda安装。源码编译两步完事：git clone下来后make即可，依赖简单，仅需系统装有zlib开发文件（zlib1g-dev）。Conda一行：conda install -c bioconda minimap2。Windows用户建议走WSL。

git clone https://github.com/lh3/minimap2
cd minimap2 && make

2. 构建参考基因组索引

比对前先把参考基因组（比如病毒的GenBank参考序列）建成.mmi索引文件。索引只需构建一次，反复使用。病毒基因组通常很小（几kb到几十kb），索引构建秒级完成。

minimap2 -d ref.mmi reference.fasta

3. 读长比对：按平台选预设

minimap2最省心的地方是-x预设参数。不用自己调k-mer大小、窗口宽度，直接按测序平台和数据类型选预设就行。Nanopore数据用lr:hq（Q20+高质量读长，Dorado basecaller，错误率<1%）或map-ont（传统Guppy数据）。PacBio HiFi用map-hifi；PacBio CLR老数据用map-pb。如果是cDNA或直接RNA测序，加上splice模式做剪接感知比对。病毒扩增子数据走ARTIC官方流程的话，minimap2默认用map-ont预设——别被读长迷惑，虽然ARTIC扩增子只有~400 bp，但因为数据来自Nanopore平台，错误特征跟Illumina完全不同，map-ont的k-mer和gap罚分配置才是正确匹配。sr预设面向~1%错误率的Illumina类短读长数据，别用在Nanopore扩增子上。

# Nanopore Q20+ 高质量读长（推荐，ONT 官方博客已改用 lr:hq）
minimap2 -ax lr:hq ref.mmi ont_reads.fastq.gz > aln.sam

# PacBio HiFi
minimap2 -ax map-hifi ref.mmi hifi_reads.fastq.gz > aln.sam

# Nanopore 传统 Guppy 数据 / ARTIC 扩增子（fieldbioinformatics 管线默认）
minimap2 -ax map-ont ref.mmi ont_reads.fastq.gz > aln.sam

4. SAM to BAM，排序，建索引

比对出来的SAM是人类可读的文本格式，体积大。用samtools转成BAM、按坐标排序、建索引，后续才能快速提取特定区间的比对。

minimap2 -ax lr:hq ref.mmi reads.fastq.gz | samtools sort -O BAM -o sorted.bam
samtools index sorted.bam

5. 覆盖度评估与变异检测

病毒基因组分析的核心输出有两个：覆盖深度分布（看哪些区域扩增子dropout了）和变异列表（VCF）。覆盖度用samtools depth或bedtools genomecov生成。变异检测方面，ARTIC管线从v1.4.5起已将默认工具从medaka切换为Clair3（基于深度神经网络，对长INDEL支持更好，是当前官方唯一推荐的变异检出工具）。medaka曾是ARTIC的默认变异检出工具，用人类数据训练，病毒场景下可能出现训练偏差，现已从ARTIC主线彻底移除。通用替代方案是bcftools mpileup，无训练偏差问题，但精度通常不如Clair3。

# 生成覆盖度统计
samtools depth -a sorted.bam > coverage.txt

# ARTIC 管线（v1.4.5+）内部调用 Clair3 做变异检测
artic minion --scheme-name artic-inrb-mpox --scheme-version v1.0.0 --read-file reads.fastq sample_name

# bcftools 独立替代方案
bcftools mpileup -f ref.fasta sorted.bam | bcftools call -mv -Ov -o variants.vcf

参数选择建议

常见问题

minimap2 和 BWA-MEM 怎么选？
看数据。Nanopore/PacBio长读长用minimap2（快30倍以上且支持长INDEL）。Illumina短读长用BWA-MEM（比对精度略高，soft clipping少）。如果你的数据是ARTIC扩增子的~400 bp短读长但来自Nanopore平台——用minimap2的map-ont预设（跟ARTIC官方fieldbioinformatics管线一致），别用sr预设。sr预设是为Illumina ~1%错误率数据设计的，不匹配Nanopore的错误特征。两边都跑一遍看比对率也行，不费事。

map-ont、map-hifi 和 lr:hq 预设怎么区分？
这几个预设对应不同测序化学和碱基识别版本的错误特征。map-ont适用于Guppy 3.x-5.x的Nanopore数据（错误率~5-10%），也是ARTIC fieldbioinformatics管线的默认预设；lr:hq是v2.27版本引入的，针对Q20+（Dorado basecaller，错误率<1%）数据，参数更激进，ONT官方博客已于2025年推荐lr:hq作为Nanopore数据的首选预设；map-hifi是PacBio HiFi专供（错误率<1%）。用错预设不会报错，但比对率和INDEL精度会打折。

为什么比对病毒基因组时比对率很低？
三个常见原因。一是参考序列选错了——比如用A系参考株去比对B系病毒，差异超过15%时minimap2会丢很多读长，换更近缘的参考株就行。二是数据里有大量宿主reads——临床样本的宿主比例可超99%，先跑Kraken2或hostile去宿主再做比对。三是忘记去引物——扩增子测序的引物序列不比对到参考基因组会被soft clip掉，看起来比对率低但其实是正常的，用iVar或cutadapt先trim引物再比对。

参考文献

1. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34(18):3094-3100. DOI: 10.1093/bioinformatics/bty191
2. Li H. New strategies to improve minimap2 alignment accuracy. Bioinformatics. 2021;37(23):4572-4574. DOI: 10.1093/bioinformatics/btab705
3. Quick J, Grubaugh ND, Pullan ST, et al. Multiplex PCR method for MinION and Illumina sequencing of Zika and other virus genomes directly from clinical samples. Nature Protocols. 2017;12(6):1261-1276. DOI: 10.1038/nprot.2017.066