Mummer 用法简析

本文详细介绍了MUMmer工具中的nucmer、promer以及run-mummer1和run-mummer3的工作流程,用于核酸序列比对。nucmer适用于高度保守的DNA序列定位,promer则基于氨基酸序列比对。此外,文章还涵盖了MUMmer的参数设置,如--mum、--mincluster等,并提供了比对实例和如何使用dnadiff和delta-filter进行结果过滤和分析。最后,提到了如何绘制共线性图谱。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

动动发财的小手,点个赞吧!

4个工作流程

nucmer

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

promer

​ 以翻译后的氨基酸序列进行联配,工作原理同nucmer.

rum-mummer1 run-mummer3

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

nucmer 参数详解

匹配:

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

实例

比较常见的用法是把一条连续的序列和另一条连续的序列比。比如说两个细菌的菌株,直接用mummer

两个完整度高的基因组

mummer -mum -b -c A.fasta B.fasta > out.mums

# -mum: 计算在两个序列中唯一的最大匹配数
# -b: 计算正向和反向匹配数
# -c: 报告反向互补序列相对于原始请求序列的位置

「高度相似」序列,不含重排

run-mummer1 ref.fasta qry.fasta ref_qry

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

「高度相似」序列,存在重排

run-mummer3 ref.fasta qry.fasta ref_qry

相似度不高

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

nucmer --maxgap=500 --mincluster=100 --prefix=ref_qry ref.fasta qry.fasta

# --maxgap: 两个match间最大gap为500
#--minclster: 至少要有100个match才能算做一簇

有点差异

可以用翻译的氨基酸序列进行比较

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

基因草图(scaffold, contig级别)

使用nucmerpromer构建序列间的可能联配。

# 首先过滤低于1kb的序列

awk -c fastx '{if (length($seq) > 1000) print ">"$name "\n"$seq}' A.fasta > A_1kb.fa

awk -c fastx '{if (length($seq) > 1000) print ">"$name "\n"$seq}' B.fa > HB_1kb.fa

# 处理,时间会比较久,因为分别有20109,17877条contig

nucmer --prefix A_B A_1kb.fa B_1kb.fa &

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

nucmer  --prefix A_B A.fa B.fa

在第四版中新增了一个dnadiff,进一步封装nucmer和其他数据整理工具,基本上没啥参数,而输出很齐全,非常的人性化。在不知如何开始的时候,可以无脑用这个。

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

# 未有delta文件
dnadiff A_B A.fasta B.fa

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上重叠

show-coord 显示匹配坐标

show-coords -r A_B.delta.filter > A_B.coord

# -r:以refID排序,相对的,还有-q,以queryID排序less IRGSP1_DHX2_i89_l1000_1.coord


「实际应用」

nucmer

nucmer --maxmatch -c 100 -b 500 -l 50 ref.fasta qur.fasta

#--maxmatch 使用所有锚匹配,而不管它们的唯一性
#--mincluster/-c: 用于聚类的匹配最低长度,默认65
#--breaklen/-b: 在对联配两端拓展式,在终止后继续延伸的程度,默认200
#--minmatch/-l: 单个匹配最小长度

delta-filter -1 -m -i 90 -l 100 out.delta > out.filtered.delta   

#-1: 1对1联配,允许重排,是-r和-q的交集
#-m: 多对对联配,允许重排,是-r和-q的合集
#-i: 最小的相似度 [0,100], 默认0
#-l: 最小的匹配长度 默认0.

show-coords -THrd out.filtered.delta > out.filtered.coords   

#-T 将输出切换为制表符分隔格式 
#-H 不打印输出头
#-r 按参考 ID 和坐标对输出行进行排序
#-d 在附加中显示对齐方向

sed ':a; N;s/\n/ /; ta' a.txt  ## 利用sed的跳转功能
#绘制共线性图谱
#需要提前安装gunplot ps2pdf程序(sudo apt-get install xxx)
mummerplot --postscript test.delta

#将图像转换为pdf文件
ps2pdf out.ps out.pdf

本文由 mdnice 多平台发布

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值