MUMMER 两个基因组间比较

2020-07-25 本文已影响0人生信师姐

MUMmer出现的历史背景
测序技术刚开始发展的时候，大家得到的序列都是单个基因的长度，所以一般都是逐个基因的比较，用的都是BLAST或FASTA通过逐个基因联配的方式搜索数据库。但是1999年后，越来越多的物种全基因组出现，就需要研究同一物种不同品系进化过程的基因组变化，比如说基因倒置现象。传统的BLAST/FASTA就用不了，就需要用到新的工具。

软件介绍
mummer的名字来源于Maximal Unique Matcher ，最大唯一性比对。mummer这个程序主要是找到参考序列和query序列之间准确匹配的区域。query最大可以有32个。mummer是不容错配的，适合用来画共线性图，但是我们通常的比对都是必须容许一定的错配和gap的，mummer比对完了之后可以使用mummerplot这个程序绘制出共线性图。

Mummer的一个最大特点就是比对速度非常快，对资源的消耗比较少，它是使用一种后缀树的算法使得它的速度比BLASTZ(k-mers)快得多，但是灵敏度低，也就是检测不到比较弱的匹配，但是作者说这都是可以通过修改参数进行改善.

Mummer官网介绍该软件是一个多才多艺的软件包，因为它可以完成生物数据分析中很多的功能。Mummer其实是一个软件包，里面包含了很多工具，这些工具搭配起来使用，可以完成非常多的工作。例如基因组比对，共线性分析，同源序列搜索，重复序列查找，SNP和Indel检测等。

MUMmer能用来研究什么呢？
比如说:

细菌的不同菌株基因组中倒置现象;
人和老鼠的基因组在进化上的重排现象;
还有比较同一物种的不同组装结果等;
等等;

全局比对软件
适用于两个全基因组（基因草图or完整基因组）的快速比对。

那么其比对就分为以下几种情况：
1）完整基因组→完成基因组；
2）草图→完成基因组；
3）草图→草图。

而常见的就是第二种，Mapping a draft sequence to a finished sequence。常用于在草图的结尾工作中将其比对到完整参考基因组序列上，帮助确定每条contig的位置和方向。同时，通过检查这些保守区域，可以提高草图注释结果的准确性。

一、软件安装

MUMmer是开源软件，因此可以通过下载源码编译的方式进行安装，同时biconda上已经有编译好的二进制版本方便用conda进行安装。目前，比较推荐使用源码编译的方式进行安装。

如果在bioconda频道上搜索mummer, 会发现一个pymummer，不要以为这是mummer的源代码用python改写，它仅仅做到了通过调用系统安装的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
#再添加值环境变量即可

二、软件原理

MUMmer的核心基于 Maximal exact matching 算法开发的mummer。其他工具(nucmer,promer,run-mummer1.run-mummer3)都是基于mummer的开发的流程。这些流程的分析策略分为三步：

用mummer在两个输入中找给定长度的极大唯一匹配( Maximal exact matching )
然后将这些匹配区域聚类成较大不完全联配区域, 作为锚定点(anchor)
最后它从每个匹配外部扩展联配, 形成有gap的联配。

1. Maximal exact matching

MUMmer核心是基于后缀树(suffix tree)数据结构的最大匹配路径。根据这个算法开发出来的repeat-match和exact-tandems可以从单个序列中检测重复，mummer则是用于联配两条或两条以上的序列。由于MUMmer的其他工具基本都是基于mummer开发的，于是理解mummer就变得非常重要。

suffix tree:表示一个字符串的所有子字符串的数据结构，比如说abc的所有子字符串就是a,ab,ac,bc,abc.

Maximal Unique Match: 指的是匹配仅在两个比较序列中各出现一次

mummer: 基于后缀树(suffix tree)数据结构，能够在两条序列中有效定位极大唯一匹配(maximal unique matches)，因此它比较适用于产生一组准确匹配(exact matches)以点图形式展示，或者用来锚定从而产生逐对联配(pair-wise alignments)

大部分情况下都不会直接用到mummer，所以只要知道MUMmer历经几次升级，使得mummer可以能够只找在reference和query都唯一的匹配(第一版功能)，也可以找需要在reference唯一的匹配(第二版新增)，甚至不在乎是否唯一的匹配(第三版新增),参数分别为-mum,-mumreference,maxmatch。

repeat-match和exact-tandems比较少用，毕竟参数也不多，似乎有其他更好的工具能用来寻找序列中的重复区。

2. Clustering:聚类

MUMmer的聚类算法能够比较智能地把几个独立地匹配按照顺序聚成一块。分为两种模式gaps和mgaps。这两者差别在于是否允许重排,分别用于run-mummer1,run-mummer3.

image

基于gap和mgaps的输出，第四版还提供了annotate和combineMUMs两个工具增加联配信息。

3. 联配构建工具

基于上述两个工具，作者编写了4个工作流程，方便实际使用。

nucmer: 根据命名我们可以看出，(NUCleotide MUMmer) ，是在核酸水平进行比对的工具。由Perl写的流程，用于联配很相近(closely related)核酸序列。首先找到两条序列之间准确匹配区域，然后进行延伸，在使用mgaps进行cluter程序，最终保留那些满足设定阈值的比对结果。找出全局比对的同源序列。它比较适合定位和展示高度保守的DNA序列。注意，为了提高nucmer的精确性，最好把输入序列先做遮盖(mask)避免不感兴趣的序列的联配，或者修改单一性限制降低重复导致的联配数。
promer：也是Perl写的流程，它以翻译后的氨基酸序列进行联配，工作原理同nucmer.如果两条序列之间亲缘关系比较远的话，我们可以采用promer比对，(PROtein MUMmer)，和nucmer类似，只不过是将两条DNA在氨基酸水平进行比对，promer会将两条序列分别翻译成六种氨基酸模式，正反个三次。然后在氨基酸水平进行比对，这样具有更高的敏感型。注意是DNA在氨基酸水平进行比对，该软件并不支持直接的氨基酸进行比对。使用起来与nucmer类似。但是速度会慢很多，可以用于绘制共线性图，但是不适合来找SNP。
run-mummer1,run-mummer3： 两者是基于cshell写的流程，用于两个序列的常规联配，和promer,nucmer类似，只不过能够自动识别序列类型。它们擅长联配相似度高的DNA序列，找到它们的不同，也就是适合找SNP或者纠错。前者用于1v1无重排，后者1v多有重排

三、nucmer的使用

nucmer的使用也比较简单，后面接选项参数，然后是参考序列reference，query序列即可。都是fasta格式，注意这里面query只能是一个，并不像mummer软件最大支持32个。

nucmer 只适用于DNA比对，有多种比对模式，可以选择比对在参考序列和query上都是唯一的，也就是1对1区域的比对，还可以限定参考序列区域wei1，也就是1对多的比对，还有就是可以对多对的比对。一般选择1对1唯一比对。最终会输出为delta格式文件，这种文件格式非常奇怪，非常的简洁，主要是一堆数字，最开始是程序名，然后是参考序列和query的文件，用于后续程序处理读取文件。对于delta格式的结果，mummer提供了大量工具来进行解析和处理。

nucmer [options] <reference> <query file>

参数我将其分为五个部分，匹配算法，聚类，外延、其他和新增

匹配:

--mum, --mumreference(默认), --maxmatch
--minmatch/-l: 单个匹配最小长度
--forwoard/-f, --reverse/-r: 只匹配正链或只匹配负链。

聚类：

--mincluster/-c: 用于聚类的匹配最低长度，默认65
--maxgap/-g: 两个相邻匹配间的最大gap长度，默认90
--diagdiff/-D: 一个聚类中两个邻接匹配，最大对角差分，默认5
--diagfactor/-d: 也是和对角差分相关参数，只不过和gap长度有关，默认0.12

外延:

--breaklen/-b: 在对联配两端拓展式，在终止后继续延伸的程度，默认200
--[no]extend：是否外延，默认是
--[no]optimize：是否优化，默认是。即在联配分数较低时不会立刻终止，而是回顾整条联配，看能否苟延残喘一会。

其他:

--depend: 显示依赖信息后退出
--coords: 调用show-coords输出坐标信息
--prefix/-p: 输出文件的前缀
--[no]delta: 是否输出delta文件，默认是

新增：

# 在第四版新增的参数
--threads/-t: 多核心
---delta=PATH: 指定位置，而不是当前
--sam-short=PATH：保存为SAM短格式
--sam-long=PATH： 保存为SAM长格式
--save=PREFIX：保存suffix array
--load=PREIFX：加载suffix array

运行后得到一个delta格式的文件，它的作用是记录每个联配的坐标，每个联配中的插入和缺失的距离。下面逐行进行解释

/home/username/reference.fasta /home/username/query.fasta # 两个比较文件的位置
PROMER # 程序运行类型： NUCMER或PROMER
>tagA1 tagB1 3000000 2000000 # 一组联配(可以有多个小匹配)，ref的fastaID，qry的fastaID，ref序列长度，qry序列长度
1667803 1667078 1641506 1640769 14 7 2 # 第一小组 ref起始，ref结束，qry起始，qry结束，错误数(不相同碱基+indel碱基数)，相似错误(非正匹配得分) 终止密码子(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

1. 用法举例

1.1 两个完整度高的基因组

（1）比较常见的用法是把一条连续的序列和另一条连续的序列比。比如说两个细菌的菌株,直接用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: 计算在两个序列中唯一的最大匹配数
# -b: 计算正向和反向匹配数
# -c: 报告反向互补序列相对于原始请求序列的位置

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

run-mummer1 ref.fasta qry.fasta ref_qry

run-mummer1 ref.fasta qry.fasta ref_qry -r    # 仅报告负链匹配序列

（3）高度相似序列，存在重排现象

run-mummer3 ref.fasta qry.fasta ref_qry

（4）以上的run-mummer*比较关注序列的不同之处，那么对于相似度没有那么高的两个序列，就需要用到nucmer。nucmer关注序列的相似之处，所以它允许重排，倒置和重复现象。nucmer允许多对多的比较方式，当然比较常用的是多对一的比较。

nucmer --maxgap=500 --mincluster=100 --prefix=ref_qry ref.fasta qry.fasta
## --maxgap: 两个match间最大gap为500
##--minclster: 至少要有100个match才能算做一簇

注意：第四版中run-mummer1, run-mummer3已经被废弃了，就是尽管保留了，但是没有对它做任何升级的意思。

（5）如果是有点差异的两个序列，可以用翻译的氨基酸序列进行比较

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

1.2 两个基因草图

上面都是两条序列间的比较，但是研究植物的人更容易遇到的是两个物种的基因组都只有scafold级别，甚至是contig级别。那么就可以使用nucmer或promer构建序列间的可能联配。

# 首先过滤低于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和其他数据整理工具，基本上没什么参数，而输出很齐全，非常的人性化。在不知如何开始的时候，可以用这个。

# 已有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. 数据整理

nucmer 最终会输出为delta格式文件。对于delta格式的结果，mummer提供了大量工具来进行解析和处理。

delta格式的结果需要用delta-filter,show-coords和show-tilling进行进一步整理才能用于后续的分析。后续操作基于上面的基因草图和完整基因组比较结果。

首先就可以使用delta-filter工具，看名字就知道这个工具可以对delta格式进行处理。最初的比对结果保留了最多的信息，需要用delta-filter进行一波过滤，除去不太合适的部分。过滤选项有

-i: 最小的相似度 [0,100], 默认0
-l: 最小的匹配长度默认0.
-u: 最小的联配唯一度 [0,100], 默认0
-o: 最大重叠度，针对-r和-q设置。 [0,100], 默认100
-g: 1对1全局匹配，不允许重排
-1: 1对1联配，允许重排，是-r和-q的交集
-m: 多对对联配，允许重排，是-r和-q的合集。
-q: 仅保留每个query在reference上的最佳位置,允许多条query在reference上重叠
-r: 仅保留每个reference在query上的最佳位置,允许多条reference在query上重叠

以上顺序是-i -l -u -q -r -g -m -1.光看参数估计不太明白，来一波图解。referece的一个片段可以联配到query的多个片段上，同样的query的一个片段也可以联配到reference的多个片段上，那么如何取舍呢？

image

通过-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

如何才能验证上面参数运行的结果是符合要求的呢？毕竟数据分析第一原则“不要轻易相信分析结果，需要多次验证才能使用”。

可以先用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

不难发现这个位置锚定的非常不错，至少暂时看起来没有重叠之处

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会尝试根据contig和reference匹配信息构建出tiling path(不好翻译呀。。)，不怎么用得到。

Mummer软件包也提供了更好用的工具show-diff，输入delta格式比对结果，就可以找出两条序列之间的不同，包括gap，重拍，染个题变化，倍增等等信息，非常方便。

show-diff显示大的染色体变化倍增重排或者直接使用dnadiff软件一步生成，结果非常详细，dnadiff可以直接加-d接delta格式的结果，或者更方便直接接两条序列即可，非常方便好用。

四、其他使用示例

下面介绍的是基因组草图比对上完整基因组的情况。

1. 比对

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

nucmer --prefix=ref_qry ref.fasta qry.fasta   #比对
delta-filter -q ref_qry.delta > ref_qry.filter   #过滤，除去不太适合的部分，但结果不适合读
show-coords -rcl ref_qry.delta > ref_qry.coords   #将结果转换为以人类可读的格式显示匹配的坐标
show-aligns ref_qry.delta refname qryname > ref_qry.aligns   #查看某一个匹配的序列的比对情况
show-tiling ref_qry.delta > ref_qry.tiling   #将contig回帖到ref上

show-aligns ：用于显示比对，可以单独数据每个序列的比对情况。
show-coords：用于显示比对坐标
show-tiling ：定位contig的位置，拼接完成图中可以用到。

2. SNP检测
将几个MUMMER组件联合起来还能用于snp的检测。

nucmer --prefix=ref_qry ref.fasta qry.fasta
show-snps -Clr ref_qry.delta > ref_qry.snps 
    # -C指输出唯一匹配的snp 
    # -l 输出结果中包括序列的长度 
    # -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等的结束）的距离
[DIST] ：从这个SNP到最近序列末端的距离
[R] ：覆盖此参考位置的重复alignments数
[Q] ：覆盖此query 位置的重复alignments数
[LEN R] ：参考序列的长度
[LEN Q] ： query序列的长度
[FRM] ：序列（NUCmer）或读码框（PROmer）的方向
[TAGS] ：分别是ref和query的FastA ID。所有位置都与DNA输入序列的正义链有关，而BUFF距离则与排序序列有关。

indel标准：DBUFF为1且P2不变的连续列，总长<50。
BUFF是离最近的mismatch的距离。p2如果是indel只显示起始的碱基位置。

找到python脚本将snp输出文件转成vcf格式

如果比对过程中没有选择唯一比对，这个时候也可以选用delta-filter工具，选择-1选项，过滤掉重复的比对。使用show-snps工具，delta格式文件中记录了序列之间的单碱基错配信息，那么我们只需要使用show-snps工具处理一下即可得到snp结果，而且该工具可以进行多种输出格式的调节，非常方便。

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。如果序列本身碱基非常准确的话，那么得到的SNP就可以继续做后续的分析了。

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

官网：
https://github.com/mummer4/mummer
http://mummer.sourceforge.net/
http://mummer.sourceforge.net/manual/

参考：
http://mummer.sourceforge.net/
http://mummer.sourceforge.net/manual/
https://www.jianshu.com/p/2e184e5c15b7
https://www.jianshu.com/p/33e87e19eb44

MUMMER 两个基因组间比较

一、软件安装

二、软件原理

1. Maximal exact matching

2. Clustering:聚类

3. 联配构建工具

三、nucmer的使用

1. 用法举例

1.1 两个完整度高的基因组

1.2 两个基因草图

1.3 一个基因草图对一个完整基因组

2. 数据整理

四、其他使用示例

猜你喜欢

热点阅读