MUMmer共线性分析与SNP检测

系统发育相关的基因组之间既存在保守性又存在可变性。有些序列片段的数目以及顺序具有保守性，这种保守性可以使用共线性（synteny）或同线性（colinearity）来进行描述。共线性主要强调两方面，一是序列的同源性，二是序列片段的排列顺序。同时即使很近缘的基因组也可能存在大量的变异和多态性，这种变异可能构成了不同个体与群体性状差异的基础。单核苷酸多态性(single-nucleotide polymorphism，SNP)是指由于单个核苷酸位置上存在转换或颠换等变异所引起的DNA序列多态性，常用来研究近缘物种基因组的进化。

MUMmer（http://mummer.sourceforge.net/）是一个大尺度的长片段DNA和蛋白序列全局比对工具，可以方便的比较两个物种完整基因组之间的差别。

MUMmer3安装方法如下所示（http://mummer.sourceforge.net/manual/）：

tar -zxvf MUMmer3.23.tar.gz
cd MUMmer3.23
make install

MUMmer4的安装（https://github.com/mummer4/mummer/releases）：

tar -zxvf mummer-4.0.0beta2.tar.gz
cd mummer-4.0.0beta2
./configure --prefix=/data/tengwenkai/software/MUMmer4.0
make
make install

MUMmer原理

MUMmer的核心基于Maximalexact matching算法开发的mummer。其他工具（nucmer、promer、run-mummer1、run-mummer3）都是基于mummer的开发的流程。这些流程的分析策略分为三部：

①用mummer在两个输入中找给定长度的极大唯一匹配（Maximal exact matching）

②然后将这些匹配区域聚类成较大不完全联配区域，作为锚定点（anchor）

③最后它从每个匹配外部扩展联配，形成有gap的联配。

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

概念1：suffix tree: 表示一个字符串的所有子字符串的数据结构，比如说abc的所有子字符串就是a、ab、ac、bc、abc。

概念2：Maximal Unique Match指的是匹配仅在两个比较序列中各出现一次。

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

大部分情况下都不会直接用到主程序mummer，所以只要知道MUMmer历经几次升级，使得mummer可以能够只找在reference和query都唯一的匹配(第一版功能)，也可以找在reference唯一但不一定在query唯一的匹配(第二版新增)，甚至不在乎是否唯一的匹配(第三版新增)，参数分别为-mum、-mumreference、-maxmatch。repeat-match和exact-tandems比较少用，毕竟参数也不多，似乎有其他更好的工具能用来寻找序列中的重复区。

MUMmer的聚类算法能够比较智能地把几个独立的匹配按照顺序聚成一组匹配，分为两种模式gaps和mgaps（如下图所示）。这两者差别在于是否允许重排（例如倒置现象），分别用于run-mummer1、run-mummer3。

mummer最适合生成可以显示为dot plot的精确匹配列表，或者用作生成成对比对中的锚点。基于mummer，作者编写了以下4个pipeline，方便实际使用：

nucmer：由Perl写的流程，用于联配很相近(closely related)核酸序列。它比较适合定位和展示高度保守的DNA序列。注意，为了提高nucmer的精确性，最好把输入序列先做遮盖(mask)避免不感兴趣的序列的联配，或者修改单一性限制降低重复导致的联配数。

promer：也是Perl写的流程，工作原理类似nucmer。其在进行任何精确匹配之前，将输入序列被翻译成所有六种读框的氨基酸。这使得promer能够鉴定在DNA水平上可能不保守的保守蛋白质序列的区域，并因此使其具有比nucmer更高的灵敏度。注意，灵敏度的增加将导致大量输出高度相似的序列，因此建议仅当输入序列太分散难以产生合理数量的nucmer输出时才使用promer。

run-mummer1和run-mummer3：两者是基于cshell写的流程，用于两个序列的常规联配，和promer、nucmer类似，只不过能够自动识别序列类型。它们擅长联配相似度高的DNA序列，找到它们的不同，也就是适合找SNP或者纠错。前者用于1v1无重排，后者1v多有重排。

MUMmer使用

由于MUMmer有一个主程序和4个主流程，因此决定每种情况下最佳的MUMmer比对程序十分重要。MUMmer的使用情况可能有以下几种：

①两个完成序列的全局比对，例如两个细菌基因组的比较。独立的mummer程序，与mummerplot结合，可能是可视化两个序列的全局比对所必需的，有助于确定两个序列之间的差异，其使用如下所示：

./mummer [options] <reference-file> <query-files>
-mum：只寻找在reference和query都唯一的匹配
-mumreference：寻找在reference唯一但不一定在query唯一的匹配（默认）
-maxmatch：寻找所有匹配，不在乎是否唯一
-n：只匹配字符a、c、g、t，可以大写或小写，忽略掉被mask的序列
-l：匹配的最短长度，默认为20
-b：同时查找正向链和反向互补链的匹配
-r：只查找反向互补链的匹配
-s：显示匹配的子字符串
-c：汇报与原始链对应的反向互补匹配的query-position
-F：不管输入序列的数目，强制4列的输出结果格式
-L：显示query序列的长度

使用mummer对两个基因组进行分析：

MUMmer4.0/bin/mummer -mum -b -c -n 1171_armatimo.fasta 142_armatimo.fasta > 1171_142.mums

结果如下所示（第一列为查询基因组中的位置，第二列为参考基因组中的位置，第三列为匹配长度）：

Mummerplot使用方法如下所示：

mummerplot [options] <match file>
match file：匹配文件，由mummer、nucmer、promer等程序生成(后缀.out、.cluster、.delta、.tiling)
-f, --filter：只展示.delta比对中best匹配（在一对多模式中）
--fat：只展示使用fattest比对的序列
-p|prefix：设置输出结果的文件前缀，默认为'out'
-rv：x11格式结果背景颜色反转
-r|IdR：指定X轴绘制的序列ID
-q|IdQ：指定Y轴绘制的序列ID
-R|Rfile：通过文件Rfile指定参考序列的绘制顺序
-Q|Qfile：通过文件Qfile指定查询序列的绘制顺序，Rfile/Qfile可以是fasta序列文件，也可以是序列ID的列表
-s|size：输出图片的大小，可选small、medium、large，相当于--small、--medium、--large，默认为small。
-S, --SNP：在比对中标出SNP位点
-t|terminal：输出结果为x11、postscript、png，相当于--x11、--postscript、--png，默认为x11，x11是一种互动展示，postscript为矢量格式
-t|title：设置图片的标题，默认为none
-x|xrange：设置x轴的范围[min:max]
-y|yrange：设置y轴的范围[min:max]

注意，MUMmer3中的mummerplot只适合gnuplot 5.0以下的版本。使用mummerplot对分析结果作图：

MUMmer4.0/bin/mummerplot -t postscript -p 1171_142 1171_142.mums

作图结果如下所示：

不同颜色分别代表不同方向（正向链/反向互补链）上的匹配。

②没有重排的高度相似序列，例如同一个属或种的基因组。当比较两个几乎相同的序列，比对的目的通常是SNP和small InDel的鉴定。原MUMmer 1.0的pipeline仍然是这类分析的一个方便的工具，使用方法如下所示（注意，MUMmer 4.0的scripts/run-mummer1.sh脚本中bindir须修改为包含mumer和gap命令的路径，由于4.0版安装后的bin中没有gap命令，因此可设置为MUMmer3.23的路径；此外MUMmer3.23中的run-mummer1脚本有一点错误，需要在21行tail命令后面添加-n，此处可参考run-mummer1.sh）：

MUMmer3.23/run-mummer1 142_armatimo.fasta 391_armatimo.fasta 142_391

结果如下所示（第一列为查询基因组）：

可以看出，由于是很近缘的物种的基因组，匹配区域很长，但有一个长序列的重排（此重排很可能是由于基因组序列未作复制起点校正所致）。Gaps文件给出了匹配之间的gap长度，如下所示（第五列为连续匹配之间的gap长度）：

如果正向链匹配效果不好，还可以查询反向互补链的匹配与gap：

MUMmer3.23/run-mummer1 142_armatimo.fasta 391_armatimo.fasta 142_391 -r

③有重排的高度相似序列，有时候两个序列是高度相似的，但是会出现大片段的序列重排、颠倒或插入。为了比对这些序列应该使用run-mummer3，它使用一种聚类方法，允许这些大规模突变的类型，但保留了run-mummer1的许多其他功能。为了更准确地寻找SNP，您可以编辑脚本，并将-D选项添加到combineMUMs命令行，从而产生一个仅两个序列之间差异位置的简明文件。用法如下所示：

MUMmer3.23/run-mummer3 142_armatimo.fasta 391_armatimo.fasta 142_391

重排之后gap减少。在脚本里添加-D后的align文件给出了gap处的碱基差异，如下所示：

④较相似序列的比对，run-mummer1和run-mummer3更多地关注两个序列之间的区别，而nucmer关注的是什么是相同的。它对一致性对齐的限制很少，所以重新排列，反转和重复都将被nucmer识别。nucmer的使用方法如下所示（需要预先安装Foundation.pm）：

nucmer [options] <Reference> <Query>
Reference：参考基因组，含有多条序列的FASTA文件名
Query：要匹配的基因组，含有多条序列的FASTA文件名
--mum, --mumreference(默认), --maxmatch：与mumer相同
-b, --breaklen：一个比对尝试延伸的最大距离，默认为200
-c, --mincluster：一个匹配聚类簇的最短长度，默认为65
-D, --diagdiff：一个聚类中两个邻接匹配的最大对角差分，默认5
-d, --diagfactor一个聚类中两个邻接匹配的最大对角差分与gap长度的比值，默认为0.12
--noextend：不执行聚类簇延长步骤，默认关闭
-f, --forward：只使用查询序列的正向链
-g, --maxgap：一个聚类中两个邻接匹配的最大gap长度，默认为90
-l, --minmatch：一个匹配的最短长度，默认为20
-L, --minalign：一个聚类延伸后比对的最短长度，默认为0
-r, --reverse：只使用查询序列的反向互补链
--nosimplify：不简化比对，当使用序列与自身比对来寻找重复时可以选此选项，默认关闭
-p, --prefix：输出结果delta文件的前缀，默认为out
--sam-short：保存SAM短格式到文件路径
--sam-long：保存SAM长格式到文件路径
-t, --threads：程序运行使用的核数

使用nucmer对两个基因组进行比较分析：

MUMmer4.0/bin/nucmer --mum -g 500 -c 100 -p 1171_142 142_armatimo.fasta 1171_armatimo.fasta

运行后得到一个delta格式的文件，它的作用是记录每个联配的坐标，每个联配中的插入和缺失的距离。使用show-coords脚本可以将delta文件转换为易读的匹配坐标：

MUMmer4.0/bin/show-coords -r 1171_142.delta > 1171_142.coords

其中-r表示按照参考序列的ID以及坐标进行分类，结果如下所示：

使用show-aligns可以查看具体的序列比对情况，如下所示：

MUMmer4.0/bin/show-aligns -r 1171_142.delta 142_armatimo1 1171_armatimo1 > 1171_142.aligns

结果如下所示：

使用mummerplot进行可视化，如下所示：

MUMmer4.0/bin/mummerplot -t postscript -p 1171_142 1171_142.delta

作图结果如下所示：

⑤较不相似的序列比对，很多基因的DNA序列差异较大，但蛋白序列是保守的，因此比较蛋白序列能寻找到更多的匹配，promer可以将DNA序列翻译成蛋白序列进行比对，其使用参数与nucmer类似，如下所示：

MUMmer4.0/bin/promer --mum -p 1171_142 142_armatimo.fasta 1171_armatimo.fasta

使用delta-filter过滤掉Identity低于50%的匹配：

MUMmer4.0/bin/delta-filter -q -r -i 50 1171_142.delta > 1171_142.filter

进行可视化作图：

MUMmer4.0/bin/mummerplot -t postscript -p 1171_142 1171_142.filter

作图结果如下所示：

⑥检测SNP，SNP主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态性，因此在检测SNP时需要对基因组进行比对，排除插入缺失、基因重排的影响，寻找匹配聚类簇中的单核苷酸变异位点，如下所示：

MUMmer4.0/bin/nucmer -p 142_391 142_armatimo.fasta 391_armatimo.fasta

重复序列可能会掩盖可能的SNP，因此使用delta-filter去除一对多、多对多中的冗余匹配：

MUMmer4.0/bin/delta-filter -r -q 142_391.delta > 142_391.filter

使用show-snps显示匹配中的SNP，如下所示：

MUMmer4.0/bin/show-snps -Clr 142_391.filter > 142_391.snps

结果如下所示：

END