multiple whole genome alignment
写在前面：在上亿年的进化历程中，基因组经历了大大小小的改变。从小的核苷酸突变、插入、缺失到大的基因缺失、重复、基因组重排和水平基因转移。基因组比对可以通过比对序列中的同源区域，找到DNA中的改变。理解每种改变的速率和模式将有助于了解多种多样的生物学过程。
经过google，高分推荐的multiple whole genome alignment软件有Mavue, Mugsy和Multiz三个。这三个软件都可以通过conda安装。
但相比于conda，源码附带的readme是我们了解这个软件的重要参考，因此我更偏爱源码。
Mugsy: fast multiple alignment of closely related whole genomes
1. Mugsy
Mugsy是Angiuoli SV和Salzberg SL于2010年，开发的一款用于multiple whole genome alignment 的工具。
特点：
- 可以接受组装基因组草图文件（contig/scaffold，一个文件中有多条fasta序列）
- 不需要参考基因组
- 非常适合closely related genomes
- 可以识别到duplication、rearrangement和大规模的gain/loss。
- 目前测试的是几十个细菌基因组。
1.1 安装
wget https://sourceforge.net/projects/mugsy/files/latest/download
tar xvzf download
#将Mugsy安装路径加到PATH里
vim ~/.bashrc
export PATH="/picb/evolgen/users/gushanshan/software/mugsy/mugsy_1.2.2:$PATH"
source ~/.bashrc
1.2. 运行
1.2.1 输入
文件要求
- 一个或多个FASTA文件
- 每个文件应该包括单个物种的所有序列
- FASTA的header中不能包括：或-
- 模糊字符将被转换成N
选项
-p: 输出文件前缀
--directory：用于储存输出和临时文件的目录，必须是全路径，否则会报错
1.2.2 输出
MAF格式
1.2.3 例子
mugsy脚本可以完成multiple alignment 的所有过程。
其中，核心程序是：
mugsyWGA - whole genome aligner based on Seqan::TCoffee
synchain-mugsy - segmentation program to produce locally collinear
blocks LCBs from a set of anchors
nucmer - 3.20 release bundled for convenience with new utility
delta2maf and modified delta-filter to add support for reporting
duplications
例子：
mugsy --directory /local/scratch --prefix mygenomes genome1.fsa genome2.fsa genome3.fsa
1.3 预实验
选择10个Lp genome，进行预实验，记录时间。
10 genomes
Starting Nucmer: Wed Sep 23 12:43:43 CST 2020
.........
Finished Nucmer Wed Sep 23 12:51:10 CST 2020
Starting MUGSYWGA: Wed Sep 23 12:51:10 CST 2020
Finished MUGSYWGA: Wed Sep 23 15:55:35 CST 2020
Final output (MAF format): path/output/pre_experiment.maf
Finished Wed Sep 23 15:55:37 CST 2020
话说好久啊，10个基因组花了3个小时。那我有500多个···
1.4 maf格式转为fasta
Mugsy安装目录下有转换文件：
./maf2fasta.pl < a.maf > a.fasta
MAF格式：
##maf version=1 scoring=tba.v8
# tba.v8 (((human chimp) baboon) (mouse rat))
a score=23262.0
s hg18.chr7 27578828 38 + 158545518 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s panTro1.chr6 28741140 38 + 161576975 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s baboon 116834 38 + 4622798 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
s mm4.chr6 53215344 38 + 151104725 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
s rn3.chr4 81344243 40 + 187371129 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG
a score=5062.0
s hg18.chr7 27699739 6 + 158545518 TAAAGA
s panTro1.chr6 28862317 6 + 161576975 TAAAGA
s baboon 241163 6 + 4622798 TAAAGA
s mm4.chr6 53303881 6 + 151104725 TAAAGA
s rn3.chr4 81444246 6 + 187371129 taagga
a score=6636.0
s hg18.chr7 27707221 13 + 158545518 gcagctgaaaaca
s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca
s baboon 249182 13 + 4622798 gcagctgaaaaca
s mm4.chr6 53310102 13 + 151104725 ACAGCTGAAAATA
FASTA格式：
>hg18.chr7 27578828 38 + 158545518
AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
>panTro1.chr6 28741140 38 + 161576975
AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
>baboon 116834 38 + 4622798
AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
>mm4.chr6 53215344 38 + 151104725
-AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
>rn3.chr4 81344243 40 + 187371129
-AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG
=
>hg18.chr7 27699739 6 + 158545518
TAAAGA
>panTro1.chr6 28862317 6 + 161576975
TAAAGA
>baboon 241163 6 + 4622798
TAAAGA
>mm4.chr6 53303881 6 + 151104725
TAAAGA
>rn3.chr4 81444246 6 + 187371129
taagga
=
>hg18.chr7 27707221 13 + 158545518
gcagctgaaaaca
>panTro1.chr6 28869787 13 + 161576975
gcagctgaaaaca
>baboon 249182 13 + 4622798
gcagctgaaaaca
>mm4.chr6 53310102 13 + 151104725
ACAGCTGAAAATA
=
把比对信息抽出来，抽成一块一块的，不再是原来摞在一起的：
Args <- commandArgs()
a<-file(Args[6],"r")
output<-Args[7]
lines=readLines(a,n=1)
filename="1"
while(length(lines)!=0){
if(lines !="="){
cat(paste0(lines,"\n"),file=paste0(paste(output,filename,sep = "/"),".fasta"),append = T)
lines=readLines(a,n=1)
}else{
filename=as.character(as.numeric(filename)+1)
lines=readLines(a,n=1)
}
}
close(a)
第一个参数是fasta文件，第二个是单一文件的输出目录
mkdir a #存储输出文件
Rscript ./singleFasta.R ./a.fasta a/
然后传到MEGAX里就可以看了。
2. Mauve
Mauve分为原始Mauve算法和progressive Mauve算法。
original Mauve 的劣势：
- 比较适合closely-related 物种
- 不比对部分基因组共享的大区域
- 识别不到部分基因组共享的重排区域
- 为了准确估计基因组重排，必须人工调整最小LCB权重
- 很难比对经常出现片段重复的基因组
progressiveMauve算法的优势：
- 能够比对更多的基因组
- 能够比对分歧更大的基因组，核苷酸相似性可以低至50%
- 不需要人工调整比对打分参数
- 可以比对pan-genome
- 准确度更高
progressiveMauve算法的劣势： - 很慢
- 很耗内存
2.1 安装
#下载地址：http://darlinglab.org/mauve/download.html
wget http://darlinglab.org/mauve/snapshots/2015/2015-02-13/linux-x64/mauve_linux_snapshot_2015-02-13.tar.gz
tar -zxvf mauve_linux_snapshot_2015-02-13.tar.gz
rm mauve_linux_snapshot_2015-02-13.tar.gz
mv mauve_snapshot_2015-02-13/ mauve_2015_02_13
#把路径加到环境变量里。因为要用命令行程序，所以路径指定到安装路径下的linux-x64
export PATH="xxxxxx/mauve/linux-x64:$PATH"
source ~/.bashrc
2.2 运行
2.2.1 输入
文件格式
输入文件可以是FastA，Multi-FastA或 GenBank。如果一个文件中有多个序列，即Multi-FastA，Mauve会先将它们合并起来，再将合并后的序列和其他序列比对。
选项
--output：输出文件名称
--collinear：假设输入文件是共线性，即没有重排
重要参数的默认设置
--island-gap-size=<number> Alignment gaps above this size in nucleotides are considered to be islands [20]
似乎不允许换行？因为程序太长，我换行了。然后就报错了···弄成一行后，又正常了
2.2.2 输出
比对完成后，Mauve会产生多个比对相关文件。下面介绍一下每个文件包含的信息及对应的格式。
2.2.2.1 .alignment文件和XMFA格式(.xmfa)
.alignment文件以 eXtended Multi-FastA格式存储Mauve产生的比对数据。
XMFA格式中存储了多个共线性块的比对情况，不同共线性块以=分割。
每个共线性块中，一个基因组有一条对应的fasta格式的序列。其中，定义行给出这条序列位于基因组的位置和方向（正负链）。
这些共线性块共同组成了基因组比对结果。
>seq_num:start1-end1 ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
> seq_num:startN-endN ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
= comments, and optional field-value pairs, i.e. score=12345
> seq_num:start1-end1 ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
> seq_num:startN-endN ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
= comments, and optional field-value pairs, i.e. score=12345
Mauve略微改造了下xmfa格式，要求基因组序列中的核苷酸能且仅能出现一次。
- 对于未对齐且在LCB之外的基因组片段，以XMFA中无gap的单条目形式出现。
- XMFA中的第一个LCB条目会列出所有输入的基因组序列。如果某些基因组未参与该lCB的比对，其对应的坐标范围标记为0-0。
2.2.2.2 Islands文件(.bbcols)
.island文件中存储了比对中发现的genomic islands，以tab键分割。
Island指的是比对中一部分基因组有，而另一部分基因组没有的区域。
一个island由一个序列的基因组坐标定义，其中另一个基因组在比对的那部分中包含长度为n或更长的缺口。缺口的长度n可以人为定义。是不是很拗口，看下例子就明白了。
Genome 0: ACACGTTCGCTTCGAAA
Genome 1: ACAC------TTCGAA-
Genome 2: ATACGATCGCTTCGTAA
设定n=5，则称基因组0和2在5-10位置上有个岛
对应的.islands文件中，一行记录一个岛。格式：基因组A➡️岛最左侧在A中的位置➡️岛最右侧在A中的位置➡️基因组B➡️岛最左侧在B中的位置➡️岛最右侧在B中的位置。
0 4 11 1 4 5
1 4 5 2 4 11
第一行记录在基因组0中，核苷酸4至11与基因组1中的核苷酸4至5对齐。
岛的长度=｜(rightA - leftA) - (rightB - leftB)｜。
2.2.2.4 backbone文件(.backbone)
什么是骨架，支撑性的东西。
2.2.2.4.1 The original Mauve backbone file
.backbone文件存储了在所有基因组中都是保守的比对区域。
保守区域定义为长度大于等于x，期间gap不超过y个核苷酸的alignment片段。
22256 22371 20147 20299 22255 22370
来自第一个基因组的片段[22,256-22,371]分别在第二个[20,147-20,299]和第三个基因组[22,255-22,370]中保守。
2.2.2.4.2 The Progressive Mauve backbone file
跟original mauve不同的是，progressive mauve的backbone文件不再要求在所有基因组中都是保守的，align regions conserved among two or more genomes。backbone文件按照seq_0_leftend列排序。
The Progressive Mauve backbone file
可以从backbone文件中推出island位置
2.2.2.5 相似度矩阵文件
行和列都是基因组文件，按照输入顺序排序。
0代表没有一个相同的核苷酸，1代表每个位置的核苷酸都是相同的。
2.2.2.6 SNP文件
记录了SNP模式（按照输入顺序排序）以及SNP在每个基因组中的位置。
SNP pattern sequence_1 sequence_2 sequence_3
AAT 5276590 5246627 394
TTC 5276784 5246821 588
AAC 5277418 5247455 1222
MAA 5278225 5248262 2030
AAC 5282804 5252841 6609
2.2.2.7 orthologs文件
0:Z03:2818-3750 1:c04:3512-4444 2::2801-3733 3:ECSE_03:2800-3732
0:Z04:3751-5037 1:c05:4445-5731 2::3734-5020 3:ECSE_04:3733-5019
0:Z05:5251-5547 1:c07:5945-6241 1:c08:6021-6269 2::5234-5530 3:ECSE_05:5233-5529
0:Z06:5700-6476 1:c10:6301-7077 3:ECSE_06:5682-6458
第一行列出了4个同源基因，每个基因分别来自一个基因组。0:Z03:2818-3750中：
- 0代表基因组index，index跟输入顺序一致。
- Z03代表注释到的基因的locus_tag。
- 2818-3750代表在基因组上的位置范围。
第2行中，2::2801-3733代表该区域没有注释到的基因，因此没有locus_tag，但是是同源的。
第3行中，基因组1上注释到了两个基因。
第4行中，基因组2上没有注释到基因。
2.2.3 例子
# 比对1, 2 , 3这3个基因组，将输出保存到threeway.xmfa文件中
progressiveMauve --output=threeway.xmfa genome_1.gbk genome_2.gbk genome_3.gbk
相同的测试文件，Mauve用了大约45分钟。这明明比Mugsy快好吧···
