如何使用MUMmer比对大片段序列测序技术刚开始发展的时候,大家得到的序列都是单个基因的长度,所以一般都是逐个基因的比较,用的都是BLAST或FASTA通过逐个基因联配的方式搜索数据库。但是1999年后,越来越多的物种全基因组出现,比如说在1999年出现了Helicobacter pylori的第二类菌株的基因组序列,就需要研究同一物种不同品系进化过程的基因组变化,比如说基因倒置现象。传统的BLAST/FASTA就用不了,就需要用到新的工具,这就是MUMmer出现的历史背景。 那么MUMmer能用来研究什么呢?比如说细菌的不同菌株基因组中倒置现象,人和老鼠的基因组在进化上的重排现象。还有比较同一物种的不同组装结果等。MUMmer的算法基础(suffix tree)使得它的速度比BLASTZ(k-mers)快得多,但是灵敏度低,也就是检测不到比较弱的匹配,但是作者说这都是可以通过修改参数进行改善. 安装MUMmer是开源软件,因此可以通过下载源码编译的方式进行安装,同时biconda上已经有编译好的二进制版本方便用conda进行安装。目前,我个人比较推荐使用源码编译的方式进行安装。目前MUMmer已经更新到第四版,但是还在测试中,所以文章也没有发,求稳还是用3.23.
# MUMmer3.23wget https://gigenet.dl./project/mummer/mummer/3.23/MUMmer3.23.tar.gztar -xf MUMmer3.23.tar.gzcd MUMmer3.23make install# MUMmer4.00-beta2wget https://github.com/mummer4/mummer/releases/download/v4.0.0beta2/mummer-4.0.0beta2.tar.gztar xf mummer-4.0.0beta2.tar.gzcd mummer-4.0.0beta2./configure --prefix=$HOME/biosoft/mummer-4.0.0beta2 && make && make install 为了方便使用记得将软件路径加入PATH。 MUMmer使用方法MUMmer的核心基于 Maximal exact matching 算法开发的
Maximal exact matchingMUMmer核心是基于后缀树(suffix tree)数据结构的最大匹配路径。 根据这个算法开发出来的
mummer: 基于后缀树(suffix tree)数据结构,能够在两条序列中有效定位极大唯一匹配(maximal unique matches),因此它比较适用于产生一组准确匹配(exact matches)以点图形式展示,或者用来锚定从而产生逐对联配(pair-wise alignments) 大部分情况下都不会直接用到 repeat-match和exact-tandems比较少用,毕竟参数也不多,似乎有其他更好的工具能用来寻找序列中的重复区。 Clustering:聚类
联配构建工具基于上述两个工具,作者编写了4个工作流程,方便实际使用。
重点介绍一下 nucmer [options] 参数我将其分为五个部分,匹配算法,聚类,外延、其他和新增 匹配: --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-340 用法举例
比较常见的用法是把一条连续的序列和另一条连续的序列比。比如说两个细菌的菌株,直接用 wget http://mummer./examples/data/H_pylori26695_Eslice.fastawget http://mummer./examples/data/H_pyloriJ99_Eslice.fastamummer -mum -b -c H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta > 26695_J99.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 以上的 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
上面都是两条序列间的比较,但是研究植物的人更容易遇到的是两个物种的基因组都只有scafold级别,甚至是contig级别。那么就可以使用 # 首先过滤低于1kb的序列bioawk -c fastx '{if (length($seq) > 1000) print '>'$name '\n'$seq}' ~/reference/genome/rice_contigs/HP1 > HP103_1kb.fabioawk -c fastx '{if (length($seq) > 1000) print '>'$name '\n'$seq}' ~/reference/genome/rice_contigs/HP119.fa > HP119_1kb.fa# 处理,时间会比较久,因为分别有20109,17877条contignucmer --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
# 已有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 数据整理之前得到的数据还需要用 最初的比对结果保留了最多的信息,需要用
以上顺序是 通过 比如说我想看contig和reference两者唯一匹配,并且长度在1000,相似度大于90. delta-filter -i 89 -l 1000 -1 IRGSP1_DHX2.delta > IRGSP1_DHX2_i89_l1000_1.delta.filter 如何才能验证上面参数运行的结果是符合要求的呢?毕竟数据分析第一原则“不要轻易相信分析结果,需要多次验证才能使用”。 可以先用 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 不难发现这个位置锚定的非常不错,至少暂时看起来没有重叠之处 用 show-aligns IRGSP1_DHX2_i89_l1000_1.delta.filter chr01 DHX2_00006753 | less 针对reference有很长的组装序列的情况,还可以用 |
|