如何使用MUMmer比對大片段序列

測序技術剛開始發(fā)展的時候障陶，大家得到的序列都是單個基因的長度炫狱，所以一般都是逐個基因的比較脸候，用的都是BLAST或FASTA通過逐個基因聯(lián)配的方式搜索數(shù)據庫顽素。但是1999年后粱胜，越來越多的物種全基因組出現(xiàn)飒硅，比如說在1999年出現(xiàn)了Helicobacter pylori的第二類菌株的基因組序列砂缩，就需要研究同一物種不同品系進化過程的基因組變化，比如說基因倒置現(xiàn)象三娩。傳統(tǒng)的BLAST/FASTA就用不了庵芭，就需要用到新的工具，這就是MUMmer出現(xiàn)的歷史背景雀监。

那么MUMmer能用來研究什么呢双吆？比如說細菌的不同菌株基因組中倒置現(xiàn)象眨唬，人和老鼠的基因組在進化上的重排現(xiàn)象。還有比較同一物種的不同組裝結果等好乐。MUMmer的算法基礎(suffix tree)使得它的速度比BLASTZ(k-mers)快得多匾竿，但是靈敏度低，也就是檢測不到比較弱的匹配蔚万，但是作者說這都是可以通過修改參數(shù)進行改善.

安裝

MUMmer是開源軟件岭妖，因此可以通過下載源碼編譯的方式進行安裝，同時biconda上已經有編譯好的二進制版本方便用conda進行安裝反璃。目前区转，我個人比較推薦使用源碼編譯的方式進行安裝。目前MUMmer已經更新到第四版版扩，但是還在測試中，所以文章也沒有發(fā)侄泽，求穩(wěn)還是用3.23.

多說一句礁芦，如果在bioconda頻道上搜索mummer, 會發(fā)現(xiàn)一個pymummer，不要以為這是mummer的源代碼用python改寫悼尾，它僅僅做到了通過調用系統(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

為了方便使用記得將軟件路徑加入PATH闺魏。

MUMmer使用方法

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

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

Maximal exact matching

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

概念1：suffix tree: 表示一個字符串的所有子字符串的數(shù)據結構娩怎，比如說abc的所有子字符串就是a,ab,ac,bc,abc.
概念2：Maximal Unique Match: 指的是匹配僅在兩個比較序列中各出現(xiàn)一次

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

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

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

Clustering:聚類

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

gaps

mgaps

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

聯(lián)配構建工具

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

• nucmer: 由Perl寫的流程蒋院，用于聯(lián)配很相近(closely related)核酸序列亏钩。它比較適合定位和展示高度保守的DNA序列。注意欺旧，為了提高nucmer的精確性姑丑，最好把輸入序列先做遮蓋(mask)避免不感興趣的序列的聯(lián)配篮条，或者修改單一性限制降低重復導致的聯(lián)配數(shù)对碌。
• promer：也是Perl寫的流程，它以翻譯后的氨基酸序列進行聯(lián)配而账，工作原理同nucmer.
• run-mummer1,run-mummer3： 兩者是基于cshell寫的流程称龙，用于兩個序列的常規(guī)聯(lián)配留拾，和promer,nucmer類似，只不過能夠自動識別序列類型鲫尊。它們擅長聯(lián)配相似度高的DNA序列痴柔，找到它們的不同，也就是適合找SNP或者糾錯疫向。前者用于1v1無重排咳蔚，后者1v多有重排

重點介紹一下nucmer的使用。reference和query文件都需要時fasta格式搔驼，每個都可以有多條序列屹篓。

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: 也是和對角差分相關參數(shù)，只不過和gap長度有關焊刹，默認0.12

外延:

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

其他:

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

新增：

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

運行后得到一個delta格式的文件以清，它的作用是記錄每個聯(lián)配的坐標，每個聯(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結束障贸，qry起始，qry結束吟宦，錯誤數(shù)(不相同堿基+indel堿基數(shù))篮洁，相似錯誤(非正匹配得分) 終止密碼子(NUCMER為0)。 如果結束大于起始殃姓，表示在負鏈袁波。
-145 # qry的145有插入
-3   # qry的145+3=148有插入
-1   # qry的145+3+1=149有插入
40   # qry的145+3+1+40=149有缺失
0 # 表示當前匹配結束
1667804 1667079 1641507 1640770 10 5 3 # 第二小組
-146
-1
-1
-34
0

用法舉例

兩個完整度高的基因組

比較常見的用法是把一條連續(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: 報告反向互補序列相對于原始請求序列的位置

或者是高度相似序列蜗侈，不含重排

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

或者是高度相似序列篷牌，存在重排現(xiàn)象

run-mummer3 ref.fasta qry.fasta ref_qry

以上的run-mummer*比較關注序列的不同之處，那么對于相似度沒有那么高的兩個序列踏幻，就需要用到nucmer枷颊。nucmer關注序列的相似之處，所以它允許重排该面，倒置和重復現(xiàn)象夭苗。nucmer允許多對多的比較方式，當然比較常用的是多對一的比較隔缀。

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

注意一點： 第四版中run-mummer1, run-mummer3已經被廢棄了题造，就是盡管保留了，但是沒有對它做任何升級的意思猾瘸。

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

promer --prefix=ref_qry ref.fasta qry.fasta

兩個基因草圖

上面都是兩條序列間的比較丢习，但是研究植物的人更容易遇到的是兩個物種的基因組都只有scafold級別，甚至是contig級別淮悼。那么就可以使用nucmer或promer構建序列間的可能聯(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 &

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

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

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

在第四版中新增了一個dnadiff，進一步封裝nucmer和其他數(shù)據整理工具瞧挤，基本上沒啥參數(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

數(shù)據整理

之前得到的數(shù)據還需要用delta-filter,show-coords和show-tilling進行進一步整理才能用于后續(xù)的分析癌刽。后續(xù)操作基于上面的基因草圖和完成基因組比較結果役首。

最初的比對結果保留了最多的信息，需要用delta-filter進行一波過濾显拜，除去不太合適的部分衡奥。過濾選項有

• -i: 最小的相似度 [0,100], 默認0
• -l: 最小的匹配長度 默認0.
• -u: 最小的聯(lián)配唯一度 [0,100], 默認0
• -o: 最大重疊度，針對-r和-q設置远荠。 [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的多個片段上尼摹，那么如何取舍呢？

多對多

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

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

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

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

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

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)這個位置錨定的非常不錯吏颖，至少暫時看起來沒有重疊之處

coord信息

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

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

alignment

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

順便放一下自己的知識星球，如果你覺得我對你有幫助的話衬吆。

知識星球