【来自读者投稿】
随着地球生物基因组计划的开展，将有越来越多的上百万个真核生物复杂基因组被测序、拼装。基因组比对是对这些基因组的最核心分析之一。
2021年12月PNAS发表了一款名为anchorwave的全基因比对工具。由于目前的全基因组比对软件主要针对人或者动物，而anchorwawve依据植物基因组的特点开发，对植物基因组非常敏感和准确，同时也适用于动物和人，能够同时识别短的和长的插入或删除。
anchorwave对存在染色体融合，易位，全基因组复制的基因组比对十分敏感准确，对发生全基因组复制的基因间进行比对时需要手动设置WGD次数，所以重复序列基本不会错误。此外它还可以进行变异检测等功能，目前正在开发多序列比对，系统发育关系等功能。
本文以拟南芥的两个不同品种，对不存在染色体变异的基因组间进行全基因组的比对。
具体步骤
将基因组的fasta及注释gff文件和查询基因组的fasta文件作为输入。首先将提取序列中的保守序列(一般为外显子或者蛋白编码基因)作为靶点生成一个文件，然后将这个文件中的保守序列分别比对到参考和查询基因组分别生成两个文件。之后选择genoAli或proali功能进行全基因组比对，proali用于存在染色体融合，易位，全基因组复制的基因组比对。本例中利用软件中的genoAli命令提取基因组的共线性区域，将以上的文件作为输入作为用genoAli命令执行比对输出maf文件。
anchorwave流程图：
安装与使用
anchorwave安装与下载：
#源代码安装
git clone https://github.com/baoxingsong/anchorwave.git
cd anchorwave
cmake ./
make
#conda安装anchorwave
conda install -c bioconda -c conda-forge anchorwave
#docker安装anchorwave
docker build -f docker/Dockerfile -t anchorwave ./
docker run -it anchorwave anchorwave
下载两个基因组的序列文件和参考基因组的注释文件
#拟南芥参考基因组下载的每条染色体，使用cat命令进行合成一个完整的染色体。
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr1.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr2.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr3.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr4.fas
wget ftp://ftp.arabidopsis.org/home/tair/Sequences/whole_chromosomes/TAIR10_chr5.fas
cat TAIR10_chr1.fas TAIR10_chr2.fas TAIR10_chr3.fas TAIR10_chr4.fas TAIR10_chr5.fas > tair10.fa
下载的查询基因为SDI格式需要通过软件GEAN转化为FASTA文件格式。
wget http://mtweb.cs.ucl.ac.uk/mus/www/19genomes/variants.SDI/bur_0.v7c.sdi
gean pseudogeno -r tair10.fa -v bur_0.v7c.sdi -o bur_0.fa
提取参考基因组的中的保守序列并将保守序列分别匹配到两个基因组上
这一步需要使用必须使用anchorwave，这样可以过滤一些不必要的anchors，最小化的减少下一步中minimap2的限制。-r代表要输入的参考基因组，-i代表他的注释文件，-o代表他输出的所有的保守序列的文件。
anchorwave gff2seq -r tair10.fa -i TAIR10_GFF3_genes.gff -o cds.fa
-x splice 代表将保守序列映射到基因组，-t 代表线程数，-k 代表kmer， -a代表输出sam格式的文件，-p和-N代表如果它们的比对得分大于初次比对得分的0.4将会保留最多保留20个结果，还要输入他们的基因组和保守序列。
minimap2 -x splice -t 6 -k 12 -a -p 0.4 -N 20 tair10.fa cds.fa > ref.sam
minimap2 -x splice -t 6 -k 12 -a -p 0.4 -N 20 bur_0.fa cds.fa > bur_0.sam
使用anchorwave进行全基因组比对
在这个例子中我们使用genoAli命令，它一般用于不存在染色体重排或存在少量倒位的基因组间的全基因组比对，-i代表输入参考基因组的注释文件，-r 和-s 参考和查询基因组文件，-o和-f 输出变异比对结果的maf格式文件，-n输出共线性文件可以用于可视化结果。如果过程中消耗的内存过大，可以减小滑动窗口的大小。
anchorwave genoAli -i TAIR10_GFF3_genes.gff -as cds.fa -r tair10.fa -a bur_0.sam -ar ref.sam -s bur_0.fa -v bur_0.vcf -n bur_0.anchors -o bur_0.maf -f bur_0.f.maf -t 1
相关论文如下：https://www.pnas.org/doi/10.1073/pnas.2113075119
其他代码如下：https://github.com/baoxingsong/anchorwave
