Minimap2 用户手册
概述
- 名称
- minimap2 - DNA序列集合之间的映射和比对
内容
- 简介
- 描述
- 选项
- 索引选项
- 映射选项
- 对齐选项
- 输入/输出选项
- 预设选项
- 其他选项
- 输出格式
- 限制
- 参考
简介
Minimap2是一个快速的序列映射和比对程序,能够找到长噪声读段之间的重叠,或者将长读段或它们的组装映射到参考基因组上,并且可以选择性地进行详细比对(即CIGAR)。目前,它能够高效地处理从几千碱基到约100兆碱基长度的查询序列,错误率约为15%。Minimap2以PAF或SAM格式输出。
选项
索引选项
k INT
: 最小化k-mer长度 [15]w INT
: 最小化窗口大小 [k-mer长度的2/3]。最小化是w个连续k-mer窗口中的最小k-mer。H
: 使用同聚物压缩(HPC)最小化。HPC序列是通过将同聚物运行压缩为单个碱基来构建的。HPC最小化是HPC序列上的最小化。I NUM
: 索引时最多加载NUM个目标碱基到RAM [4G]。--idx-no-seq
: 不在索引中存储目标序列。节省磁盘空间和内存,但这样生成的索引将不适用于-a或-c选项。
映射选项
U INT1[,INT2]
: k-mer出现的下限和上限 [10,1000000]。e INT
: 每隔INT碱基对采样一个高频最小化 [500]。g NUM
: 如果在NUM-bp内没有最小化,则停止链的延伸 [10k]。r NUM1[,NUM2]
: 链化和基础对齐的带宽 [500,20k]。
对齐选项
A INT
: 匹配得分 [2]。B INT
: 不匹配的惩罚 [4]。O INT1[,INT2]
: 间隙开放惩罚 [4,24]。
输入/输出选项
a
: 生成CIGAR并在SAM格式中输出比对。Minimap2默认以PAF输出。o FILE
: 将比对输出到FILE [stdout]。Q
: 在输入文件中忽略碱基质量。
预设选项
x STR
: 预设 []。这个选项同时应用多个选项。它应该在其他选项之前应用,因为后面应用的选项将覆盖-x设置的值。
输出格式
Minimap2默认以成对映射格式(PAF)输出映射位置。PAF是一个制表符分隔的文本格式,每行至少包含12个字段。
限制
Minimap2在通过长低复杂性区域时可能产生次优比对,因为在这些区域种子位置可能是次优的。
Minimap2需要SSE2或NEON指令来编译。可以添加非SSE2/NEON支持,但这会使Minimap2的速度慢几倍。
以下是GitHub上lh3/minimap2页面的中文翻译,按照Markdown格式整理:
# 开始使用
## 获取Minimap2
git clone https://github.com/lh3/minimap2
cd minimap2 && make
长序列与参考基因组比对
./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam
先创建索引再映射
./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
./minimap2 -a MT-human-ont.mmi test/MT-orang.fa > test.sam
使用预设(无测试数据)
./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam # PacBio CLR基因组读取
./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam # Oxford Nanopore基因组读取
./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS基因组读取 (v2.19或更高版本)
./minimap2 -ax lr:hq ref.fa ont-Q20.fq.gz > aln.sam # Nanopore Q20基因组读取 (v2.27或更高版本)
./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam # 短的基因组成对末端读取
./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam # 剪接长的读取(链不明确)
./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam # 噪声的Nanopore直接RNA-seq
./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam # 最终PacBio Iso-seq或传统的cDNA
./minimap2 -ax splice --junc-bed anno.bed12 ref.fa query.fa > aln.sam # 优先考虑注释的连接点
./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf # 种内asm-to-asm比对
./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf # PacBio读取重叠
./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf # Nanopore读取重叠
详细命令行选项
man ./minimap2.1
目录
- 开始使用
- 用户指南
- 安装
- 一般用法
- 使用案例
- 映射长噪声基因组读取
- 映射长mRNA/cDNA读取
- 寻找长读取之间的重叠
- 映射短准确基因组读取
- 全基因组/组装比对
- 高级特性
- 处理>65535个CIGAR操作
- 可选的cs标签
- 使用PAF格式
- 算法概览
- 获取帮助
- 引用minimap2
- 开发者指南
- 限制
用户指南
Minimap2是一个多功能序列比对程序,可将DNA或mRNA序列与大型参考数据库比对。典型用例包括:(1)将PacBio或Oxford Nanopore基因组读取映射到人类基因组;(2)寻找错误率约为15%的长读取之间的重叠;(3)针对参考基因组的PacBio Iso-Seq或Nanopore cDNA或直接RNA读取进行剪接感知比对;(4)比对Illumina单端或成对末端读取;(5)组装到组装比对;(6)两个密切相关物种之间的全基因组比对,差异度低于15%。
对于约10kb的噪声读取序列,minimap2的速度是BLASR、BWA-MEM、NGMLR和GMAP等主流长读取映射器的数十倍。在模拟长读取上更准确,产生的比对具有生物学意义,可为下游分析做好准备。对于大于100bp的Illumina短读取,minimap2的速度是BWA-MEM和Bowtie2的三倍,在模拟数据上同样准确。详细的评估可从minimap2论文或预印本获得。
安装
Minimap2针对x86-64 CPU进行了优化。你可以从发布页面获取预编译的二进制文件:
curl -L https://github.com/lh3/minimap2/releases/download/v2.28/minimap2-2.28_x64-linux.tar.bz2 | tar -jxvf -
./minimap2-2.28_x64-linux/minimap2
如果你想从源代码编译,你需要安装C编译器、GNU make和zlib开发文件。然后在源代码目录中输入make
进行编译。如果出现编译错误,尝试make sse2only=1
禁用SSE4代码,这将使minimap2速度略慢。
Minimap2也适用于支持NEON指令集的ARM CPU。要为32位ARM架构(如ARMv7)编译,使用make arm_neon=1
。要为64位ARM架构(如ARMv8)编译,使用make arm_neon=1 aarch64=1
。
Minimap2可以使用SIMD Everywhere (SIMDe)库将实现移植到不同的SIMD指令集。要使用SIMDe编译,使用make -f Makefile.simde
。要为ARM CPU编译,请使用上面给出的与ARM相关的命令行和Makefile.simde
。
一般用法
不加任何选项时,minimap2以参考数据库和查询序列文件为输入,生成近似映射(即坐标仅是近似值,输出中没有CIGAR),以PAF格式呈现:
minimap2 ref.fa query.fq > approx-mapping.paf
你可以要求minimap2在PAF的cg标签生成CIGAR:
minimap2 -c ref.fa query.fq > alignment.paf
或以SAM格式输出比对:
minimap2 -a ref.fa query.fq > alignment.sam
Minimap2可以无缝处理gzip压缩的FASTA和FASTQ格式输入。你不需要先将FASTA和FASTQ之间转换或解压缩gzip文件。
对于人类参考基因组,minimap2在映射前需要几分钟生成参考的最小化索引。为了减少索引时间,你可以使用选项**-d**保存索引,并在minimap2命令行上用索引文件替换参考序列文件:
minimap2 -d ref.mmi ref.fa # 索引
minimap2 -a ref.mmi reads.fq > alignment.sam # 比对
重要的是,需要注意的是,一旦你构建了索引,如**-k**、-w、-H和**-I**这样的索引参数在映射期间不能更改。如果你为不同类型的数据运行minimap2,你可能需要保留多个用不同参数生成的索引。这使得minimap2与BWA不同,BWA不管查询数据类型如何总是使用相同的索引。
使用案例
Minimap2对所有应用使用相同的基础算法。然而,由于它支持不同类型的数据(例如,短与长读取;DNA与mRNA读取),minimap2需要调整以获得最佳性能和准确性。通常建议使用选项**-x**选择一个预设,该选项同时设置多个参数。默认设置与map-ont
相同。
映射长噪声基因组读取
minimap2 -ax map-pb ref.fa pacbio-reads.fq > aln.sam # PacBio CLR读取
minimap2 -ax map-ont ref.fa ont-reads.fq > aln.sam # Oxford Nanopore读取
minimap2 -ax map-iclr ref.fa iclr-reads.fq > aln.sam # Illumina Complete Long Reads
map-pb
和map-ont
的区别在于map-pb
使用同聚物压缩(HPC)最小化为种子,而map-ont
使用普通最小化为种子。经验评估表明,HPC最小化在对齐PacBio CLR读取时提高了性能和灵敏度,但在对齐Nanopore读取时却有所损害。map-iclr
使用调整后的对齐计分矩阵,考虑到读取中整体错误率较低,其中转换错误比颠换错误少。
映射长mRNA/cDNA读取
minimap2 -ax splice:hq -uf ref.fa iso-seq.fq > aln.sam # PacBio Iso-seq/传统cDNA
minimap2 -ax splice ref.fa nanopore-cdna.fa > aln.sam # Nanopore 2D cDNA-seq
minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam # Nanopore直接RNA-seq
minimap2 -ax splice --splice-flank=no SIRV.fa SIRV-seq.fa # 针对SIRV对照的映射
有不同种类的长读取RNA-seq技术,包括传统的全长cDNA、EST、PacBio Iso-seq、Nanopore 2D cDNA-seq和直接RNA-seq。它们产生的数据质量和属性各不相同。默认情况下,-x splice
假定读取相对于转录本链的方向未知。它尝试两轮比对以推断方向,如果可能的话,将链写入ts
SAM/PAF标签。对于Iso-seq、直接RNA-seq和传统全长cDNA,应用-u f
强制minimap2只考虑正向转录本链是可取的。这加快了比对速度,并略微提高了准确性。对于噪声的Nanopore直接RNA-seq读取,建议使用较小的k-mer大小,以增加对首尾外显子的灵敏度。
Minimap2根据最大得分子段的得分对齐比对,不包括内含子,并在SAM中将最佳比对标记为主比对。当一个剪接基因也有未剪接的假基因时,minimap2不会有意地偏好剪接比对,尽管在实践中它更经常将剪接比对标记为主比对。默认情况下,minimap2输出多达五个次要比对(即在RNA-seq映射的背景下可能是假基因)。这可以通过选项**-N**进行调整。
对于长的RNA-seq读取,minimap2可能会产生由基因融合/结构变异引起的嵌合比对,或者由于内含子长度超过最大内含子长度**-G**(默认为200k)而引起的嵌合比对。目前,不建议应用过大的**-G**,因为这会减慢minimap2的速度,有时会导致错误的比对。
值得注意的是,默认情况下-x splice
更倾向于GT[A/G]…[C/T]AG而不是GT[C/T]…[A/G]AG,然后是其他剪接信号。考虑一个额外的碱基可以提高噪声读取的连接点准确性,但在对齐广泛使用的SIRV对照数据时会降低准确性。这是因为SIRV不遵循进化保守的剪接信号。如果你正在研究SIRV,你可以应用--splice-flank=no
让minimap2只模拟GT…AG,忽略额外的碱基。
自v2.17以来,minimap2可以可选地接受注释基因作为输入,并优先考虑注释的剪接连接点。要使用此功能,你可以:
paftools.js gff2bed anno.gff > anno.bed
minimap2 -ax splice --junc-bed anno.bed ref.fa query.fa > aln.sam
这里,anno.gff
是GTF或GFF3格式的基因注释(gff2bed
自动测试格式)。gff2bed
的输出是12列的BED格式,或BED12格式。使用--junc-bed
选项时,如果对齐的连接点与注释中的连接点匹配,minimap2会增加一个奖励得分(由--junc-bonus
调整)。选项--junc-bed
也接受5列BED,包括链场。在这种情况下,每一行表示一个定向连接点。
寻找长读取之间的重叠
minimap2 -x ava-pb reads.fq reads.fq > ovlp.paf # PacBio CLR读取重叠
minimap2 -x ava-ont reads.fq reads.fq > ovlp.paf # Oxford Nanopore读取重叠
同样,ava-pb
使用HPC最小化,而ava-ont
使用普通最小化。通常不建议在重叠模式中执行基础比对,因为它很慢,可能会产生假阳性重叠。然而,如果性能不是问题,你可以尝试添加-a或-c。
映射短准确基因组读取
minimap2 -ax sr ref.fa reads-se.fq > aln.sam # 单端比对
minimap2 -ax sr ref.fa read1.fq read2.fq > aln.sam # 配对末端比对
minimap2 -ax sr ref.fa reads-interleaved.fq > aln.sam # 配对末端比对
当指定两个读取文件时,minimap2依次从每个文件中读取,并将它们合并为内部的交错流。如果两个读取在输入流中相邻并且具有相同的名称(如果存在/[0-9]后缀,则会进行修剪),则认为它们是配对的。可以混合单端和配对末端读取。
Minimap2对短剪接读取的处理效果不佳。对于短读取,有许多能够处理RNA-seq的映射器。
全基因组/组装比对
minimap2 -ax asm5 ref.fa asm.fa > aln.sam # 组装到组装/参考比对
对于跨物种全基因组比对,需要根据序列差异调整计分系统。
高级特性
处理>65535个CIGAR操作
由于设计缺陷,BAM不适用于具有>65535个操作的CIGAR字符串(SAM和CRAM可以)。然而,对于超长纳米孔读取,minimap2可能会比BAM的能力对齐约1%的读取碱基,其CIGAR很长。如果你将这样的SAM/CRAM转换为BAM,Picard和最近的samtools会抛出错误并中止。旧版本的samtools和其他工具可能会创建损坏的BAM。
为了避免这个问题,你可以在minimap2命令行添加选项**-L**。此选项将长CIGAR移动到CG标签,并将完全剪辑的CIGAR留在SAM CIGAR列。当前不读取CIGAR的工具(例如合并和排序)仍然可以处理这样的BAM记录;读取CIGAR的工具将有效忽略这些记录。已决定未来的工具将无缝识别由选项**-L**生成的长CIGAR记录。
TL;DR:如果你使用超长读取并且使用仅处理BAM文件的工具,请添加选项**-L**。
可选的cs标签
cs SAM/PAF标签对不匹配和INDEL处的碱基进行编码。它匹配正则表达式/(:[0-9]+|\*[a-z][a-z]|[=\+\-][A-Za-z]+)+/
。像CIGAR一样,cs由一系列操作组成。每个前导字符指定操作;后面的序列是参与操作的序列。
cs标签由命令行选项--cs
启用。例如,以下比对:
CGATCGATAAATAGAGTAG---GAATAGCA
|||||| |||||||||| |||| |||
CGATCG---AATAGAGTAGGTCGAATtGCA
表示为:6-ata:10+gtc:4*at:3
,其中:[0-9]+
表示相同区块,-ata
表示删除,+gtc
表示插入,*at
表示参考碱基a
被查询碱基t
替换。它类似于MD SAM标签,但是独立的,更容易解析。
如果使用--cs=long
,cs字符串还将包含比对中的相同序列。上述示例将变为=CGATCG-ata=AATAGAGTAG+gtc=GAAT*at=GCA
。cs的长格式在一条字符串中编码了参考和查询序列。cs标签还编码内含子位置和剪接信号(详见minimap2手册页)。
使用PAF格式
Minimap2还附带了一个(java)script paftools.js
,用于处理PAF格式的比对。它从组装到参考的比对中调用变体,基于比对提升BED文件,转换格式,并提供各种评估的实用工具。详情请参阅misc/README.md。
算法概览
在以下内容中,minimap2命令行选项前面有一个短划线,并且用粗体突出显示。描述可能有助于调整minimap2参数。
- 读取**-I** [= 4G] 参考碱基,提取(-k, -w)-最小化并将其索引在哈希表中。
- 读取**-K** [= 200M] 查询碱基。对于每个查询序列,执行步骤3到7:
- 对于查询上的每个(-k, -w)-最小化,检查对参考索引。如果一个参考最小化不是在顶部**-f** [= 2e-4] 最频繁的,收集它在参考中的出现,称为种子。
- 按参考中的位置对种子进行排序。使用动态规划将它们串联。每个串联表示一个潜在的映射。对于读取重叠,报告所有串联然后转到步骤8。对于参考映射,执行步骤5到7:
- 让_P_是主映射集合,最初是一个空集。对于每个串联,根据它们的串联得分从最好到最差:如果在查询上,该串联与_P_中的串联重叠
–mask-level [= 0.5] 或更高比例的较短串联,将该串联标记为_P_中的串联的次要;否则,将该串联添加到_P_。 - 保留所有主映射。还要保留多达**-N** [= 5] 个顶级次要映射,如果它们的串联得分高于它们对应的主映射的**-p** [= 0.8]。
- 如果请求了比对,过滤出可能导致既长插入又长删除的内部种子。从最左边的种子扩展。在内部种子之间执行全局比对。如果全局比对沿累积得分下降**-z** [= 400],则拆分串联,不考虑长间隔。从最右边的种子扩展。输出串联及其比对。
- 如果输入中还有更多查询序列,转到步骤2,直到没有更多查询为止。
- 如果还有更多参考序列,从一开始就重新打开查询文件,然后转到步骤1;否则停止。
开发者指南
Minimap2不仅是一个命令行工具,还是一个程序库。它提供了C API来构建/加载索引和对比对索引的序列。文件example.c演示了C API的典型用法。头文件minimap.h提供了更详细的API文档。Minimap2的目标是保持此头文件中的API稳定。文件mmpriv.h包含额外的私有API,可能会频繁更改。
此存储库还提供了对C API子集的Python绑定。文件python/README.rst提供了完整的文档;python/minimap2.py显示了一个例子。这个Python扩展,mappy,也可以通过PyPI pip install mappy
或通过BioConda conda install -c bioconda mappy
获得。
限制
- minimap2可能会在长低复杂性区域产生次优比对,其中种子位置可能是次优的。这不应是一个大问题,因为即使最优比对在这些区域也可能是错误的。
- minimap2需要在x86 CPU上的SSE2指令或ARM CPU上的NEON。可以添加非SIMD支持,但这将使minimap2慢几倍。
- minimap2不适用于单个查询或数据库序列约20亿碱基或更长(确切地说是2,147,483,647)。所有序列的总长度可以远远超过这个阈值。
- minimap2经常错过小外显子。