multiple whole genome alignment

写在前面 ：在上亿年的进化历程中，基因组经历了大大小小的改变。从小的核苷酸突变、插入、缺失到大的基因缺失、重复、基因组重排和水平基因转移。基因组比对可以通过比对序列中的同源区域，找到DNA中的改变。理解每种改变的速率和模式将有助于了解多种多样的生物学过程。

经过google，高分推荐的multiple whole genome alignment软件有Mavue, Mugsy和Multiz三个。这三个软件都可以通过conda安装。

但相比于conda，源码附带的readme是我们了解这个软件的重要参考，因此我更偏爱源码。

Mugsy是Angiuoli SV和Salzberg SL于2010年，开发的一款用于multiple whole genome alignment 的工具。

特点 ：

文件要求

选项

MAF格式

mugsy 脚本可以完成multiple alignment 的所有过程。

其中，核心程序是：

例子：

选择10个Lp genome，进行预实验，记录时间。

话说好久啊，10个基因组花了3个小时。那我有500多个···

Mugsy安装目录下有转换文件：

MAF格式：

FASTA格式：

把比对信息抽出来，抽成一块一块的，不再是原来摞在一起的：

第一个参数是fasta文件，第二个是单一文件的输出目录

然后传到MEGAX里就可以看了。

Mauve分为原始Mauve算法和progressive Mauve算法。

original Mauve 的劣势：

progressiveMauve算法的优势：

文件格式

输入文件可以是FastA，Multi-FastA或 GenBank。如果 一个文件中有多个序列，即Multi-FastA，Mauve会先将它们合并起来 ，再将合并后的序列和其他序列比对。

选项汪纤

重要参数的默认设置

比对完成后，Mauve会产生多个比对相关文件。下面介绍一下每个文件包含的信息及对应的格式。

.alignment文件以 eXtended Multi-FastA格式存储Mauve产生的比对数据。

XMFA格式中存储了多个共线性块的比对情况，不同共线性块以 = 分割。

每个共线性块中，一个基因组有一条对应的fasta格式的序列。其中，定义行给出这条序列位于基因组的位置和方向（正负链）。

这些共线性块共同组成了基因组比对结果。

Mauve略微改造了下xmfa格式，要求基因组序列中的核苷酸能且仅能出现一次。

.island文件中存储了比对中发现的genomic islands，以tab键分割。

Island指的是比对中一部分基因组有，而另一部分基因组没有的区域。

一个迹郑island由一个序列的基因组坐标定义，其中另一个基因组在比对的那部分中包含长度为n或更长的缺口。缺口的长度n可以人为定义。是不是很拗口，看下例子就明白了。

设定n=5，则称基因组0和2在5-10位置上有个岛

对应的.islands文件中，一行记录一个岛。格式：基因组A➡️岛最左侧在A中的位置➡️岛最右侧在A中的位置➡️基因组B➡️岛最左侧在B中的位置➡️岛最右侧在B中的位置。

第一行记录在基因组0中，核苷酸4至11与基因组1中的核苷酸4至5对齐。

岛的长姿陵颂度=｜(rightA - leftA) - (rightB - leftB)｜。

什么是骨架，支撑性的东西。

.backbone文件存储了在所有基因组中都是保守的比对区域。

保守区域定义为长度大于等于x，期间gap不超过y个核苷酸的alignment片段。

来自第一个基因组的片段[22,256-22,371]分别在第二个[20,147-20,299]和第三个基因组[22,255-22,370]中保守。

跟original mauve不同的是，progressive mauve的backbone文件不再要求在所有基因组中都是保守的， align regions conserved among two or more genomes 。backbone文件按照 seq_0_leftend 列排序。

可以从backbone文件中推出island位置

行和列都是基因组文件，按照输入顺序排序。

0代表没有一个相同的核苷酸，1代表每个位置的核苷酸都是相同的。

记录了SNP模式（按照输入顺序排序）以及SNP在每个基因组中的位置。

第一行列出了4个同源基因，每个基因分别来自一个基因组。0:Z03:2818-3750中：

第2行中，2::2801-3733代表该区域没有注释到的基因，因此没有locus_tag，但是是同源的。

第3行中，基因组1上注释到了两个基因。

第4行中，基因组2上没有注释到基因。

相同的测试文件，Mauve用了大约45分钟。这明明比Mugsy快好吧···

有LASTZ+ Multiz和LAST+ Multiz两种方案。鉴于UCSC用lastz，我就先试试第一种方案。

https://github.com/mcfrith/last-genome-alignments/issues/7

因为用的是UCSC的自动化流程-doBlastzChainNet.pl，所以需要配置一些东西

ucsc上的doBlastzChainNet.pl可以自动完成pairwise alignement的所有步骤，包括partition, blastz, cat, chainRun, chainMerge, net, load, download, cleanup, syntenicNet。最后直接得到可用于Multiz的结果文件。

基因组文件需要先压缩成gz格式，然后再用faToTwoBit将基因组文件转为.2bit格式，用twoBitInfo获取各染色体的信息。

创建一个名为DEF的文件，里面为lastz所用参数。

参数意义详见： http://genomewiki.ucsc.edu/index.php/DoBlastzChainNet.pl

为了缩短比对时间，将genome切成了chunk。为了充分利用多核系统，UCSC开发了parasol系统，用于管理任务。

a. 比对前，要先开启parasol

b. 比对中

c. 比对后

最后 BASE/axtChain/target.query.synNet.maf 为结果文件，作为Multiz的输入文件。

合并pairwise的结果：

写在最后：

---