三代序列比对算法minimap2，就是快，快，快！

本文来自“生信算法”公众号。

        由于二代测序序列长度短（~300 bp），错误率低（<0.02%），其测序错误主要为替换错误，可通过哈希查找或数据压缩等方法直接查询每条序列在基因组中的准确比对区域，得到比对结果。而三代测序序列一方面错误率较高（~15%），且错误类型主要为插入或删除错误，难以通过文本查找直接找到序列在基因组中的比对区域，需要采用针对性的查询策略找到比对区域。另一方面，在比对阶段，由于三代序列读段长（平均长度超过10 kbp），直接采用全局比对（Needleman-Wunsch，NW）或局部比对（Smith-Waterman，SW）比对算法需构建L×L大小的得分矩阵（其中L表示序列长度），计算空间复杂度较高。基于以上两方面分析，目前绝大多数针对二代测序的比对方法不适合处理三代测序数据。因此， 需要开发针对三代测序数据的序列比对算法。

        minimap2算法是Heng Li大牛2018年在Bioinformatics期刊上发表的一篇三代序列比对算法文章。目前谷歌引用次数已经达到431次（截止时间2019.11.3日），可谓是三代序列处理一大利器。当时正好我也在研究三代序列比对，发现minimap2方法确实好用，最大的优势就是速度快，是真的快。而且比对结果比较不错，要说缺点，那主要的缺点就是耗内存，是牺牲内存换速度，毕竟二者不可兼得嘛。

minimap2主要思想

        minimap2的主要思想是：首先将基因组序列的minimizer存储在哈希表中(minimizer指一段序列内最小哈希值的种子)；然后对于每一条待比对序列， 找到待比对序列所有的minimizer，通过哈希表找出其在基因组中的位置， 并利用chaining算法寻找待比对区域；最后将非种子区域用动态规划算法进行比对，得到比对结果。minimap2方法只对最小哈希值的种子进行存储，可有效降低时间复杂度。 

1. 搜索minimizer

        minimizer指的是一段序列内最小哈希值的种子，也就是哈希值最小的k-mer。k-mer是长度为k的序列子片段。DNA序列由A、C、G、T四个字符组成，按照计算机编码可以看成一个四进制数。那一个k-mer就可以看做k位的四进制数。比如GCT的哈希值就是2×4的2次方+1×4的1次方+3×4的0次方=39，所以GCT的哈希值就是39。那么可以算出每一个k-mer的哈希值，取w窗口内最小哈希值的k-mer，就是作者定义的minimizer。

        minimap2首先计算基因组序列的minimizer，存储到哈希表中。然后计算待比对序列的minimizer，通过哈希表就可以查找与基因组中一样的minimizer在基因组中的位置。这样每一个minimizer包含三个信息：（1）在基因组中的位置；（2）在待比对序列中的位置；（3）minimizer长度。

2. chaining算法

        通过chaining就找到一组minimizer后，一个minimizer就是一个种子，也是待比对序列和基因组匹配的区域。下一步只需将序列的非种子区域进行比对，与种子区域连接起来，就是最后的序列比对结果。类似于BLAST思想。非种子区域一般比较短，当然是相对整条待比对序列来说的。这样就可以运用传统的NW算法或者SW算法进行比对。

minimap2结果比较

        对于三代PacBio序列（模拟序列），minimap2与其他5个比对方法进行比较：blasr，bwa，graphmap、minialign和ngmlr。从图（a）可以看出，minimap在比对的序列条数上优势较大，明显高于其他5个三代序列比对方法。图（b）是对于二代测序序列数据的比对结果，由于minimap2主要针对三代序列，图（b）主要说明minimap2也可以比对二代测序数据。

        下表是minimap2对三代Nanopore测序数据的比对结果（人类基因组测序数据），从图中可以看出minimap2速度简直太快了，15.9分钟，至少速度提升了50%。

minimap2安装

        minimap2安装非常方便，直接在GitHub网站上进行下载后即可安装，（https://github.com/lh3/minimap2）。用C语言编写。下载后直接make就可以完成安装。

        安装及运行步骤如下：

git clone https://github.com/lh3/minimap2
cd minimap2 && make
# long sequences against a reference genome
./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam
# create an index first and then map
./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
# use presets (no test data)
./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio CCS genomic reads
./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
# man page for detailed command line options
man ./minimap2.1

 总结

        序列比对序列数据分析的重要研究内容。序列比对是为得到序列之间的相似性或同源性， 将两条序列或多条序列按照一定规则进行排列的操作。序列比对是序列分析最为基本和重要的步骤。通过序列比对， 可以确定序列间的相似性，方便序列聚类。二代测序平台输出的序列片段一般较短（~200 bp）， 而最新的三代测序技术可以产生长度超过10 kbp的序列数据， 目前二代测序序列比对算法不能快速有效处理三代测序序列， 需要研究者开发快速高效的三代序列比对算法。 而minimap2就是针对三代序列设计的一款非常好用的比对软件。希望在你以后的测序数据分析中可以用到。

参考文献：

1. Li H. Minimap2: pairwise alignment for nucleotide sequences[J]. Bioinformatics, 2018, 34(18): 3094-3100.