在使用Purge_dups去冗余時(shí)，用到了Minimaps2，把學(xué)習(xí)的東西整理一下。本文軟件介紹的很多內(nèi)容來自 “生信算法”公眾號(hào)文章，用來作為自己的學(xué)習(xí)記錄。

一.軟件介紹

Minimaps2是李恒大神在2018年發(fā)表在bioinformatics上的一款針對(duì)三代數(shù)據(jù)開發(fā)的比對(duì)工具。網(wǎng)友說這款工具最大的優(yōu)勢(shì)是比對(duì)速度快，我昨天下午跑上了我的數(shù)據(jù)，沒有設(shè)置線程數(shù)，默認(rèn)應(yīng)該是3，結(jié)果今天來了看就生成了一點(diǎn)的結(jié)果，重新設(shè)置到20線程數(shù)，看一下。第一次對(duì)三代數(shù)據(jù)進(jìn)行比對(duì)，所謂的快沒有體驗(yàn)到。

二．minimap2主要思想

minimap2的主要思想是：首先將基因組序列的minimizer存儲(chǔ)在哈希表中(minimizer指一段序列內(nèi)最小哈希值的種子)；然后對(duì)于每一條待比對(duì)序列，找到待比對(duì)序列所有的minimizer，通過哈希表找出其在基因組中的位置，并利用chaining算法尋找待比對(duì)區(qū)域；最后將非種子區(qū)域用動(dòng)態(tài)規(guī)劃算法進(jìn)行比對(duì)，得到比對(duì)結(jié)果。minimap2方法只對(duì)最小哈希值的種子進(jìn)行存儲(chǔ)，可有效降低時(shí)間復(fù)雜度。其比對(duì)實(shí)現(xiàn)主要經(jīng)過以下幾個(gè)步驟：

1.?搜索minimizer

minimizer指的是一段序列內(nèi)最小哈希值的種子，也就是哈希值最小的k-mer。k-mer是長(zhǎng)度為k的序列子片段。DNA序列由A、C、G、T四個(gè)字符組成，按照計(jì)算機(jī)編碼可以看成一個(gè)四進(jìn)制數(shù)。那一個(gè)k-mer就可以看做k位的四進(jìn)制數(shù)。比如GCT的哈希值就是2×4的2次方+1×4的1次方+3×4的0次方=39，所以GCT的哈希值就是39。那么可以算出每一個(gè)k-mer的哈希值，取w窗口內(nèi)最小哈希值的k-mer，就是作者定義的minimizer。

minimap2首先計(jì)算基因組序列的minimizer，存儲(chǔ)到哈希表中。然后計(jì)算待比對(duì)序列的minimizer，通過哈希表就可以查找與基因組中一樣的minimizer在基因組中的位置。這樣每一個(gè)minimizer包含三個(gè)信息：（1）在基因組中的位置；（2）在待比對(duì)序列中的位置；（3）minimizer長(zhǎng)度。

2.?chaining算法

通過chaining就找到一組minimizer后，一個(gè)minimizer就是一個(gè)種子，也是待比對(duì)序列和基因組匹配的區(qū)域。下一步只需將序列的非種子區(qū)域進(jìn)行比對(duì)，與種子區(qū)域連接起來，就是最后的序列比對(duì)結(jié)果。類似于BLAST思想。非種子區(qū)域一般比較短，當(dāng)然是相對(duì)整條待比對(duì)序列來說的。這樣就可以運(yùn)用傳統(tǒng)的NW算法或者SW算法進(jìn)行比對(duì)。

3.minimap2結(jié)果比較

對(duì)于三代PacBio序列（模擬序列），minimap2與其他5個(gè)比對(duì)方法進(jìn)行比較：blasr，bwa，graphmap、minialign和ngmlr。Minima2p在比對(duì)的序列條數(shù)上優(yōu)勢(shì)較大，明顯高于其他5個(gè)三代序列比對(duì)方法。且minimap2也可以用來比對(duì)二代數(shù)據(jù)。

三．下載安裝minimape2

Github地址：

https://github.com/lh3/minimap2/

下載最新的試試

tar -zxvf minimap2-2.24.tar.gz

make

寫進(jìn)環(huán)境變量

運(yùn)行成功

四．用于Purge_dups過程的比對(duì)工作

我安裝這個(gè)軟件是為了配合Purge_dups的使用進(jìn)行比對(duì)過程，

Purge_dups運(yùn)行第一步：根據(jù)覆蓋度計(jì)算分界點(diǎn)(cutoff)

# gzip可以替換成pigz, 進(jìn)行多線程壓縮

minimap2-2.24/minimap2 -x map-hifi prefix.p_ctg.fa hifi_cell_reads.fastq.gz | gzip? > pb_aln.paf.gz

生成的結(jié)果文件為paf格式，第一次接觸到，PAF文件描述了比對(duì)序列的位置、長(zhǎng)度信息，PAF: a Pairwise mApping Format，參考官方鏈接ref2上介紹：

五．其它常用功能、結(jié)果的解讀、FLAG含義

參考 https://www.jianshu.com/p/d1868194b65e

六.參考：

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

ref2:https://github.com/lh3/miniasm/blob/master/PAF.md

本文使用 文章同步助手 同步