写在前面：在上亿年的进化历程中，基因组经历了大大小小的改变。从小的核苷酸突变、插入、缺失到大的基因缺失、重复、基因组重排和水平基因转移。基因组比对可以通过比对序列中的同源区域，找到DNA中的改变。理解每种改变的速率和模式将有助于了解多种多样的生物学过程。
经过google，高分推荐的multiple whole genome alignment软件有Mavue, Mugsy和Multiz三个。这三个软件都可以通过conda安装。
但相比于conda，源码附带的readme是我们了解这个软件的重要参考，因此我更偏爱源码。
补充：
过了两个月，发现mauve和mugsy比对都有一些小问题。一些blast可以找出的homologous gene ，在比对结果中却找不到···。我打算换方法，先Lastz比对，然后Multiz合并。
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列排序。
可以从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快好吧···
3. 先pairwise alignment，后合并
有LASTZ+ Multiz和LAST+ Multiz两种方案。鉴于UCSC用lastz，我就先试试第一种方案。
- 受网友Mr_我爱读文献启发，我才知道还有这一种办法，感谢感谢
- 与LASTZ相比，LAST的输入基因组不用提前屏蔽掉重复序列。
What distinguishes LAST from BLAST and similar tools (e.g. BLAT, LASTZ, YASS)? The main difference is that it copes more efficiently with repeat-rich sequences (e.g. genomes). For example: it can align reads to genomes without repeat-masking, without becoming overwhelmed by repetitive hits.- Multiz是内嵌在TBA中的一个reference-biased 多基因组比对工具。UCSC用它来生成多基因组比对。通常，先让其他序列和reference genome做两两比对，再以两两比对的结果为输入，供Multiz处理。
- 参考链接：
很大一部分参考了ucsc genomewiki，写得很清楚。
https://darencard.net/blog/2019-11-01-whole-genome-alignment-tutorial/
http://genomewiki.ucsc.edu/index.php/Whole_genome_alignment_howto
https://www.biostars.org/p/43225/
http://seqanswers.com/forums/showthread.php?t=29383
http://www.bx.psu.edu/miller_lab/dist/tba_howto.pdf
https://github.com/mcfrith/last-genome-alignments/issues/7
3.1 安装
因为用的是UCSC的自动化流程-doBlastzChainNet.pl，所以需要配置一些东西
- lastz
- parasol
- kent command line binaries
- 其他不在库里的一些小程序
除了lastz，UCSC上有后三种的详细安装教程。安装的时候很容易搜到。
3.2 pipeline
3.2.1 调整lastz的score文件
lastz_D target[multiple] query --inferonly=lastzd.control --infscores=q
--inferonly: 表明该程序只用于推断参数文件
--infscores: 存储生成的score文件
control文件选的默认的文件：见https://lastz.github.io/lastz/
cat lastzd.control
# base the inference on alignments in the middle half by identity
min_identity = 25.0% # 25th percentile
max_identity = 75.0% # 75th percentile
# scale scores so max substitution score will be 100, and only use
# alignments scoring at least as well as 20 ideal matches
inference_scale = 100 # max substitution score
hsp_threshold = 20*inference_scale
gapped_threshold = hsp_threshold
# allow substitution score inference to iterate at most 20 times;
# don't perform gap penalty inference -- instead hardwire gap penalties
# relative to max substitution
max_sub_iterations = 20
max_gap_iterations = 0
gap_open_penalty = 4*inference_scale
gap_extend_penalty = 0.3*inference_scale
# use all seedword positions (don't sample)
step = 1
# adjust for entropy when qualifying HSPs
entropy = on
3.2.2 UCSC-doBlastzChainNet.pl
ucsc上的doBlastzChainNet.pl可以自动完成pairwise alignement的所有步骤，包括partition, blastz, cat, chainRun, chainMerge, net, load, download, cleanup, syntenicNet。最后直接得到可用于Multiz的结果文件。
3.2.2.1数据准备
基因组文件需要先压缩成gz格式，然后再用faToTwoBit将基因组文件转为.2bit格式，用twoBitInfo获取各染色体的信息。
# example
gzip ps128.fna
faToTwoBit ps128.fna.gz ps128.2bit
twoBitInfo ps128.2bit stdout | sort -k2,2nr > ps128.chrom.sizes
3.2.2.2 参数文件准备
创建一个名为DEF的文件，里面为lastz所用参数。
TMPDIR目录要提前创建
参数意义详见：http://genomewiki.ucsc.edu/index.php/DoBlastzChainNet.pl
# ps128 vs 299v
PATH=/data/bin:/data/scripts # kent command line binaries存储位置
BLASTZ=/data/bin/lastz-1.04.00 # lastz
BLASTZ_H=2000
BLASTZ_Y=3400
BLASTZ_L=4000
BLASTZ_K=2200
BLASTZ_Q=/data/lastz/HoxD55.q # score文件，这里要换成我们调整好的score文件
# TARGET: ps128
SEQ1_DIR=path/ps128/ps128.2bit # target genome
SEQ1_LEN=path/ps128/ps128.chrom.sizes # target chromosome size
SEQ1_CHUNK=32100000 #将target genome切成多大一份
SEQ1_LAP=10000 # target chunk之间的重叠程度
SEQ1_LIMIT=18
# QUERY: 299v
SEQ2_DIR=path/ps128/trackData/299v/299v.2bit
SEQ2_LEN=path/ps128/trackData/299v/299v.chrom.sizes
SEQ2_CHUNK=1000000
SEQ2_LIMIT=2000
SEQ2_LAP=0
BASE=/path/ps128/trackData/299v
TMPDIR=path/dev/shm
3.2.2.3 比对
为了缩短比对时间，将genome切成了chunk。为了充分利用多核系统，UCSC开发了parasol系统，用于管理任务。
a. 比对前，要先开启parasol
initParasol start
b. 比对中
nohup bash -c 'time (doBlastzChainNet.pl `pwd`/DEF -verbose=2 -noDbNameCheck \
-workhorse=localhost -bigClusterHub=localhost -skipDownload \
-dbHost=localhost -smallClusterHub=localhost \
-trackHub -fileServer=localhost -syntenicNet)' > do.log 2>&1 &
c. 比对后
initParasol stop
最后 BASE/axtChain/target.query.synNet.maf为结果文件，作为Multiz的输入文件。
3.3.3 Multiz
合并pairwise的结果：
multiz ps128_299v.maf.gz ps128_ATG-K2.maf.gz 1 out1 out2 > multipleAlign.maf
写在最后：
- 第三部分安装写的很简略，实际步骤还是挺多的。
- 第三部分的parasol只能配置在组里的服务器上，ParasolLSF可以将它配置在slurm上，但我还没学会···
-
除了上面两个缺点，doBlastzChainNet.pl还是很好用哒。
