MUMmer出現(xiàn)的歷史背景
測序技術(shù)剛開始發(fā)展的時候，大家得到的序列都是單個基因的長度，所以一般都是逐個基因的比較，用的都是BLAST或FASTA通過逐個基因聯(lián)配的方式搜索數(shù)據(jù)庫。但是1999年后，越來越多的物種全基因組出現(xiàn)，就需要研究同一物種不同品系進化過程的基因組變化，比如說基因倒置現(xiàn)象。傳統(tǒng)的BLAST/FASTA就用不了，就需要用到新的工具。

軟件介紹
mummer的名字來源于 ，最大唯一性比對。mummer這個程序主要是找到參考序列和query序列之間準(zhǔn)確匹配的區(qū)域。query最大可以有32個。mummer是不容錯配的，適合用來畫共線性圖，但是我們通常的比對都是必須容許一定的錯配和gap的，mummer比對完了之后可以使用mummerplot這個程序繪制出共線性圖。

Mummer的一個最大特點就是比對速度非?？?，對資源的消耗比較少，它是使用一種的算法使得它的速度比BLASTZ(k-mers)快得多，但是靈敏度低，也就是檢測不到比較弱的匹配，但是作者說這都是可以通過修改參數(shù)進行改善.

Mummer官網(wǎng)介紹該軟件是一個多才多藝的軟件包，因為它可以完成生物數(shù)據(jù)分析中很多的功能。Mummer其實是一個軟件包，里面包含了很多工具，這些工具搭配起來使用，可以完成非常多的工作。例如基因組比對，共線性分析，同源序列搜索，重復(fù)序列查找，SNP和Indel檢測等。

MUMmer能用來研究什么呢？
比如說:

• 細菌的不同菌株基因組中倒置現(xiàn)象;
• 人和老鼠的基因組在進化上的重排現(xiàn)象;
• 還有比較同一物種的不同組裝結(jié)果等;
• 等等;

全局比對軟件
適用于兩個全基因組（基因草圖or完整基因組）的快速比對。

那么其比對就分為以下幾種情況：
1）完整基因組→完成基因組；
2）草圖→完成基因組；
3）草圖→草圖。

而常見的就是第二種，Mapping a draft sequence to a finished sequence。常用于在草圖的結(jié)尾工作中將其比對到完整參考基因組序列上，幫助確定每條contig的位置和方向。同時，通過檢查這些保守區(qū)域，可以提高草圖注釋結(jié)果的準(zhǔn)確性。

一、軟件安裝

MUMmer是開源軟件，因此可以通過下載源碼編譯的方式進行安裝，同時biconda上已經(jīng)有編譯好的二進制版本方便用conda進行安裝。目前，比較推薦使用源碼編譯的方式進行安裝。

如果在bioconda頻道上搜索mummer, 會發(fā)現(xiàn)一個pymummer，不要以為這是mummer的源代碼用python改寫，它僅僅做到了通過調(diào)用系統(tǒng)安裝的MUMmer的工具的方式運行而已，并且功能目前實在是太弱了。

# MUMmer3.23
wget https://gigenet.dl.sourceforge.net/project/mummer/mummer/3.23/MUMmer3.23.tar.gz
tar -xf MUMmer3.23.tar.gz
cd  MUMmer3.23
make install


# MUMmer4.00-beta2
wget https://github.com/mummer4/mummer/releases/download/v4.0.0beta2/mummer-4.0.0beta2.tar.gz
tar xf mummer-4.0.0beta2.tar.gz
cd mummer-4.0.0beta2
./configure --prefix=$HOME/biosoft/mummer-4.0.0beta2 && make && make install
#再添加值環(huán)境變量即可

二、軟件原理

MUMmer的核心基于 Maximal exact matching 算法開發(fā)的mummer。其他工具(nucmer,promer,run-mummer1.run-mummer3)都是基于mummer的開發(fā)的流程。這些流程的分析策略分為三步：

• 用mummer在兩個輸入中找給定長度的極大唯一匹配( Maximal exact matching )
• 然后將這些匹配區(qū)域聚類成較大不完全聯(lián)配區(qū)域, 作為錨定點(anchor)
• 最后它從每個匹配外部擴展聯(lián)配, 形成有g(shù)ap的聯(lián)配。

1. Maximal exact matching

MUMmer核心是基于后綴樹(suffix tree)數(shù)據(jù)結(jié)構(gòu)的最大匹配路徑。 根據(jù)這個算法開發(fā)出來的repeat-match和exact-tandems可以從單個序列中檢測重復(fù)，mummer則是用于聯(lián)配兩條或兩條以上的序列。由于MUMmer的其他工具基本都是基于mummer開發(fā)的，于是理解mummer就變得非常重要。

suffix tree:表示一個字符串的所有子字符串的數(shù)據(jù)結(jié)構(gòu)，比如說abc的所有子字符串就是a,ab,ac,bc,abc.

Maximal Unique Match: 指的是匹配僅在兩個比較序列中各出現(xiàn)一次

mummer: 基于后綴樹(suffix tree)數(shù)據(jù)結(jié)構(gòu)，能夠在兩條序列中有效定位極大唯一匹配(maximal unique matches)，因此它比較適用于產(chǎn)生一組準(zhǔn)確匹配(exact matches)以點圖形式展示，或者用來錨定從而產(chǎn)生逐對聯(lián)配(pair-wise alignments)

大部分情況下都不會直接用到mummer，所以只要知道MUMmer歷經(jīng)幾次升級，使得mummer可以能夠只找在reference和query都唯一的匹配(第一版功能)，也可以找需要在reference唯一的匹配(第二版新增)，甚至不在乎是否唯一的匹配(第三版新增),參數(shù)分別為-mum,-mumreference,maxmatch。

repeat-match和exact-tandems比較少用，畢竟參數(shù)也不多，似乎有其他更好的工具能用來尋找序列中的重復(fù)區(qū)。

2. Clustering:聚類

MUMmer的聚類算法能夠比較智能地把幾個獨立地匹配按照順序聚成一塊。分為兩種模式gaps和mgaps。這兩者差別在于是否允許重排,分別用于run-mummer1,run-mummer3.

image

image

基于gap和mgaps的輸出，第四版還提供了annotate和combineMUMs兩個工具增加聯(lián)配信息。

3. 聯(lián)配構(gòu)建工具

基于上述兩個工具，作者編寫了4個工作流程，方便實際使用。

• nucmer: 根據(jù)命名我們可以看出，(NUCleotide MUMmer) ，是在核酸水平進行比對的工具。由Perl寫的流程，用于聯(lián)配很相近(closely related)核酸序列。首先找到兩條序列之間準(zhǔn)確匹配區(qū)域，然后進行延伸，在使用mgaps進行cluter程序，最終保留那些滿足設(shè)定閾值的比對結(jié)果。找出全局比對的同源序列。它比較適合定位和展示高度保守的DNA序列。注意，為了提高nucmer的精確性，最好把輸入序列先做遮蓋(mask)避免不感興趣的序列的聯(lián)配，或者修改單一性限制降低重復(fù)導(dǎo)致的聯(lián)配數(shù)。

• promer：也是Perl寫的流程，它以翻譯后的氨基酸序列進行聯(lián)配，工作原理同nucmer.如果兩條序列之間親緣關(guān)系比較遠的話，我們可以采用promer比對，(PROtein MUMmer)，和nucmer類似，只不過是將兩條DNA在氨基酸水平進行比對，promer會將兩條序列分別翻譯成六種氨基酸模式，正反個三次。然后在氨基酸水平進行比對，這樣具有更高的敏感型。注意是DNA在氨基酸水平進行比對，該軟件并不支持直接的氨基酸進行比對。使用起來與nucmer類似。但是速度會慢很多，可以用于繪制共線性圖，但是不適合來找SNP。

• run-mummer1,run-mummer3： 兩者是基于cshell寫的流程，用于兩個序列的常規(guī)聯(lián)配，和promer,nucmer類似，只不過能夠自動識別序列類型。它們擅長聯(lián)配相似度高的DNA序列，找到它們的不同，也就是適合找SNP或者糾錯。前者用于1v1無重排，后者1v多有重排

三、nucmer的使用

nucmer的使用也比較簡單，后面接選項參數(shù)，然后是參考序列reference，query序列即可。都是fasta格式，注意這里面query只能是一個，并不像mummer軟件最大支持32個。

nucmer 只適用于DNA比對，有多種比對模式，可以選擇比對在參考序列和query上都是唯一的，也就是1對1區(qū)域的比對，還可以限定參考序列區(qū)域wei1，也就是1對多的比對，還有就是可以對多對的比對。一般選擇1對1唯一比對。最終會輸出為delta格式文件，這種文件格式非常奇怪，非常的簡潔，主要是一堆數(shù)字，最開始是程序名，然后是參考序列和query的文件，用于后續(xù)程序處理讀取文件。對于delta格式的結(jié)果，mummer提供了大量工具來進行解析和處理。

nucmer [options] <reference> <query file>

參數(shù)我將其分為五個部分，匹配算法，聚類，外延、其他和新增

匹配:

--mum, --mumreference(默認), --maxmatch
--minmatch/-l: 單個匹配最小長度
--forwoard/-f, --reverse/-r: 只匹配正鏈或只匹配負鏈。

聚類：

--mincluster/-c: 用于聚類的匹配最低長度，默認65
--maxgap/-g: 兩個相鄰匹配間的最大gap長度，默認90
--diagdiff/-D: 一個聚類中兩個鄰接匹配，最大對角差分，默認5
--diagfactor/-d: 也是和對角差分相關(guān)參數(shù)，只不過和gap長度有關(guān)，默認0.12

外延:

--breaklen/-b: 在對聯(lián)配兩端拓展時，在終止后繼續(xù)延伸的程度，默認200
--[no]extend：是否外延，默認是
--[no]optimize：是否優(yōu)化，默認是。即在聯(lián)配分數(shù)較低時不會立刻終止，而是回顧整條聯(lián)配，看能否茍延殘喘一會。

其他:

--depend: 顯示依賴信息后退出
--coords: 調(diào)用show-coords輸出坐標(biāo)信息
--prefix/-p: 輸出文件的前綴
--[no]delta: 是否輸出delta文件，默認是

新增：

# 在第四版新增的參數(shù)
--threads/-t: 多核心
---delta=PATH: 指定位置，而不是當(dāng)前
--sam-short=PATH：保存為SAM短格式
--sam-long=PATH： 保存為SAM長格式
--save=PREFIX：保存suffix array
--load=PREIFX：加載suffix array

運行后得到一個delta格式的文件，它的作用是記錄每個聯(lián)配的坐標(biāo)，每個聯(lián)配中的插入和缺失的距離。下面逐行進行解釋

/home/username/reference.fasta /home/username/query.fasta # 兩個比較文件的位置
PROMER # 程序運行類型： NUCMER或PROMER
>tagA1 tagB1 3000000 2000000 # 一組聯(lián)配(可以有多個小匹配)，ref的fastaID，qry的fastaID，ref序列長度，qry序列長度
1667803 1667078 1641506 1640769 14 7 2 # 第一小組 ref起始，ref結(jié)束，qry起始，qry結(jié)束，錯誤數(shù)(不相同堿基+indel堿基數(shù))，相似錯誤(非正匹配得分) 終止密碼子(NUCMER為0)。 如果結(jié)束大于起始，表示在負鏈。
-145 # qry的145有插入
-3   # qry的145+3=148有插入
-1   # qry的145+3+1=149有插入
40   # qry的145+3+1+40=149有缺失
0 # 表示當(dāng)前匹配結(jié)束
1667804 1667079 1641507 1640770 10 5 3 # 第二小組
-146
-1
-1
-34
0

1. 用法舉例

1.1 兩個完整度高的基因組

（1）比較常見的用法是把一條連續(xù)的序列和另一條連續(xù)的序列比。比如說兩個細菌的菌株,直接用mummer

wget http://mummer.sourceforge.net/examples/data/H_pylori26695_Eslice.fasta
wget http://mummer.sourceforge.net/examples/data/H_pyloriJ99_Eslice.fasta
mummer -mum -b -c H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta > 26695_J99.mums
# -mum: 計算在兩個序列中唯一的最大匹配數(shù)
# -b: 計算正向和反向匹配數(shù)
# -c: 報告反向互補序列相對于原始請求序列的位置

（2）高度相似序列，不含重排

run-mummer1 ref.fasta qry.fasta ref_qry

run-mummer1 ref.fasta qry.fasta ref_qry -r    # 僅報告負鏈匹配序列

（3）高度相似序列，存在重排現(xiàn)象

run-mummer3 ref.fasta qry.fasta ref_qry

（4）以上的run-mummer*比較關(guān)注序列的不同之處，那么對于相似度沒有那么高的兩個序列，就需要用到nucmer。nucmer關(guān)注序列的相似之處，所以它允許重排，倒置和重復(fù)現(xiàn)象。nucmer允許多對多的比較方式，當(dāng)然比較常用的是多對一的比較。

nucmer --maxgap=500 --mincluster=100 --prefix=ref_qry ref.fasta qry.fasta
## --maxgap: 兩個match間最大gap為500
##--minclster: 至少要有100個match才能算做一簇

注意： 第四版中run-mummer1, run-mummer3已經(jīng)被廢棄了，就是盡管保留了，但是沒有對它做任何升級的意思。

（5）如果是有點差異的兩個序列，可以用翻譯的氨基酸序列進行比較

promer --mum --maxgap=500 --mincluster=100 --prefix=promer …/Data/O157.fna …/Data/K12.fna

1.2 兩個基因草圖

上面都是兩條序列間的比較，但是研究植物的人更容易遇到的是兩個物種的基因組都只有scafold級別，甚至是contig級別。那么就可以使用nucmer或promer構(gòu)建序列間的可能聯(lián)配。

# 首先過濾低于1kb的序列
bioawk -c fastx '{if (length($seq) > 1000) print ">"$name "\n"$seq}' ~/reference/genome/rice_contigs/HP1 > HP103_1kb.fa
bioawk -c fastx '{if (length($seq) > 1000) print ">"$name "\n"$seq}' ~/reference/genome/rice_contigs/HP119.fa > HP119_1kb.fa
# 處理，時間會比較久，因為分別有20109，17877條contig
nucmer --prefix HP103_HP119 HP103_1kb.fa HP119_1kb.fa &

1.3 一個基因草圖對一個完整基因組

這里可以比較一下水稻日本晴基因組和其他地方品種

nucmer  --prefix IRGSP1_DHX2 ~/reference/genome/IRGSP1.0/IRGSP-1.0_genome.fasta ~/reference/genome/rice_contigs/DHX2.fa

在第四版中新增了一個dnadiff，進一步封裝nucmer和其他數(shù)據(jù)整理工具，基本上沒什么參數(shù)，而輸出很齊全，非常的人性化。在不知如何開始的時候，可以用這個。

# 已有delta文件
dnadiff -d IRGSP1_DHX2.delta


# 未有delta文件
dnadiff IRGSP1_DHX2 ~/reference/genome/IRGSP1.0/IRGSP-1.0_genome.fasta ~/reference/genome/rice_contigs/DHX2.fa

2. 數(shù)據(jù)整理

nucmer 最終會輸出為delta格式文件。對于delta格式的結(jié)果，mummer提供了大量工具來進行解析和處理。

delta格式的結(jié)果需要用delta-filter,show-coords和show-tilling進行進一步整理才能用于后續(xù)的分析。后續(xù)操作基于上面的基因草圖和完整基因組比較結(jié)果。

首先就可以使用delta-filter工具，看名字就知道這個工具可以對delta格式進行處理。最初的比對結(jié)果保留了最多的信息，需要用delta-filter進行一波過濾，除去不太合適的部分。過濾選項有

• -i: 最小的相似度 [0,100], 默認0
• -l: 最小的匹配長度 默認0.
• -u: 最小的聯(lián)配唯一度 [0,100], 默認0
• -o: 最大重疊度，針對-r和-q設(shè)置。 [0,100], 默認100
• -g: 1對1全局匹配，不允許重排
• -1: 1對1聯(lián)配，允許重排，是-r和-q的交集
• -m: 多對對聯(lián)配，允許重排，是-r和-q的合集。
• -q: 僅保留每個query在reference上的最佳位置,允許多條query在reference上重疊
• -r: 僅保留每個reference在query上的最佳位置,允許多條reference在query上重疊

以上順序是-i -l -u -q -r -g -m -1.光看參數(shù)估計不太明白，來一波圖解。referece的一個片段可以聯(lián)配到query的多個片段上，同樣的query的一個片段也可以聯(lián)配到reference的多個片段上，那么如何取舍呢？

image

通過-i,-l可以先過濾一些比較短，并且相似度比較低的匹配情況。進一步，計算長度和相似度的乘積(加權(quán)最長增加子集)，對于-q而言就是保留左2，對于-r則是保留右3. 這就是傳說中的三角關(guān)系，這種關(guān)系可以用-m保留或者用-q消滅。

比如說我想看contig和reference兩者唯一匹配，并且長度在1000，相似度大于90.

delta-filter -i 89 -l 1000 -1 IRGSP1_DHX2.delta > IRGSP1_DHX2_i89_l1000_1.delta.filter

如何才能驗證上面參數(shù)運行的結(jié)果是符合要求的呢？畢竟數(shù)據(jù)分析第一原則“不要輕易相信分析結(jié)果，需要多次驗證才能使用”。

可以先用show-coord以人類可讀的格式顯示匹配的坐標(biāo)。

show-coords -r IRGSP1_DHX2_i89_l1000_1.delta.filter > IRGSP1_DHX2_i89_l1000_1.coord      # -r：以refID排序，相對的，還有-q，以queryID排序
less IRGSP1_DHX2_i89_l1000_1.coord

不難發(fā)現(xiàn)這個位置錨定的非常不錯，至少暫時看起來沒有重疊之處

image

用show-aligns看某一個匹配的序列比對情況。

show-aligns IRGSP1_DHX2_i89_l1000_1.delta.filter chr01 DHX2_00006753 | less

image

針對reference有很長的組裝序列的情況，還可以用show-tilling將contig回貼到reference上，如果裝了gnuplot還能用mummerplot可視化點圖.show-tiling會嘗試根據(jù)contig和reference匹配信息構(gòu)建出tiling path(不好翻譯呀。。)，不怎么用得到。

Mummer軟件包也提供了更好用的工具show-diff，輸入delta格式比對結(jié)果，就可以找出兩條序列之間的不同，包括gap，重拍，染個題變化，倍增等等信息，非常方便。

show-diff顯示大的染色體變化 倍增 重排或者直接使用dnadiff軟件一步生成，結(jié)果非常詳細，dnadiff可以直接加-d接delta格式的結(jié)果，或者更方便直接接兩條序列即可，非常方便好用。

四、其他使用示例

下面介紹的是基因組草圖比對上完整基因組的情況。

1. 比對

• ref.fasta :完成圖序列；
• qry.fasta：接近完成圖的contig序列。

nucmer --prefix=ref_qry ref.fasta qry.fasta   #比對
delta-filter -q ref_qry.delta > ref_qry.filter   #過濾，除去不太適合的部分，但結(jié)果不適合讀
show-coords -rcl ref_qry.delta > ref_qry.coords   #將結(jié)果轉(zhuǎn)換為以人類可讀的格式顯示匹配的坐標(biāo)
show-aligns ref_qry.delta refname qryname > ref_qry.aligns   #查看某一個匹配的序列的比對情況
show-tiling ref_qry.delta > ref_qry.tiling   #將contig回帖到ref上

show-aligns ：用于顯示比對，可以單獨數(shù)據(jù)每個序列的比對情況。
show-coords：用于顯示比對坐標(biāo)
show-tiling ：定位contig的位置，拼接完成圖中可以用到。

2. SNP檢測
將幾個MUMMER組件聯(lián)合起來還能用于snp的檢測。

nucmer --prefix=ref_qry ref.fasta qry.fasta
show-snps -Clr ref_qry.delta > ref_qry.snps 
    # -C指輸出唯一匹配的snp 
    # -l 輸出結(jié)果中包括序列的長度 
    # -r 按照 ref的ID和snp位置信息進行排序

$head ref_qry.snps 

    [P1]  [SUB]  [P2]      |   [BUFF]   [DIST]  |  [LEN R]  [LEN Q]  | [FRM]  [TAGS]
========================================================================================
   11820   G A   11820     |       30    11820  |  7367045  7342681  |  1  1  utg0_segment0 utg0_segment0
   11850   A G   11850     |        8    11850  |  7367045  7342681  |  1  1  utg0_segment0 utg0_segment0
   11858   T C   11858     |        8    11858  |  7367045  7342681  |  1  1  utg0_segment0 utg0_segment0
   11921   C T   11921     |       40    11921  |  7367045  7342681  |  1  1  utg0_segment0 utg0_segment0
   11961   T C   11961     |       40    11961  |  7367045  7342681  |  1  1  utg0_segment0 utg0_segment0

[P1] ：SNP在參考序列中的位置。
For indels, this position refers to the 1-based position of the first character before the indel, e.g. for an indel at the very beginning of a sequence this would report 0. For indels on the reverse strand, this position refers to the forward-strand position of the first character before indel on the reverse-strand, e.g. for an indel at the very end of a reverse complemented sequence this would report 1.
對于indels，此位置是指indel前一個堿基的位置，例如，對于位于序列開始處的indel ，它將報告0。對于反向鏈上的indel，此位置是指反向鏈上indel之前第一個字符的前鏈位置，例如，對于反向補碼序列末尾的索引，該位置將報告1。
[SUB] ：reference 和query中此位置的堿基或gap ;
[P2]：SNP在query序列中的位置 ;
[BUFF]： 在相同的比對中（alignment），從該SNP到最近的不匹配（alignment、indel、SNP等的結(jié)束）的距離
[DIST] ：從這個SNP到最近序列末端的距離
[R] ：覆蓋此參考位置的重復(fù)alignments數(shù)
[Q] ：覆蓋此query 位置的重復(fù)alignments數(shù)
[LEN R] ： 參考序列的長度
[LEN Q] ： query序列的長度
[FRM] ：序列（NUCmer）或讀碼框（PROmer）的方向
[TAGS] ：分別是ref和query的FastA ID。所有位置都與DNA輸入序列的正義鏈有關(guān)，而BUFF距離則與排序序列有關(guān)。

indel標(biāo)準(zhǔn)：DBUFF為1且P2不變的連續(xù)列，總長<50。
BUFF是離最近的mismatch的距離。p2如果是indel只顯示起始的堿基位置。

找到python腳本將snp輸出文件轉(zhuǎn)成vcf格式

如果比對過程中沒有選擇唯一比對，這個時候也可以選用delta-filter工具，選擇-1選項，過濾掉重復(fù)的比對。使用show-snps工具，delta格式文件中記錄了序列之間的單堿基錯配信息，那么我們只需要使用show-snps工具處理一下即可得到snp結(jié)果，而且該工具可以進行多種輸出格式的調(diào)節(jié)，非常方便。

nucmer --mum --maxgap=500 --mincluster=100 --prefix=nucmer …/…/…/Data/O157.fna …/…/…/Data/K12.fna
delta-filter -1 -q -r nucmer.delta > nucmer.filter
show-snps -C -H -I -T -r -l nucmer.filter >nucmer.snp

那么這個就是我們最終得到的參考序列與query序列之間的SNP。如果序列本身堿基非常準(zhǔn)確的話，那么得到的SNP就可以繼續(xù)做后續(xù)的分析了。

注意事項：
1、Mummer的安裝需要make、sed、g++，awk等工具；
2、promer只是用于核酸在氨基酸水平比對，并不能直接用于氨基酸比對。

官網(wǎng)：
https://github.com/mummer4/mummer
http://mummer.sourceforge.net/
http://mummer.sourceforge.net/manual/

參考：
http://mummer.sourceforge.net/
http://mummer.sourceforge.net/manual/
http://www.itdecent.cn/p/2e184e5c15b7
http://www.itdecent.cn/p/33e87e19eb44