anchorwave进行复杂基因组比对（2）

本文采用两个不同品种的拟南芥进行全基因组比对和变异检测。这种比对方法使得每个相对应的染色体名称都一样。并且对于两个相同物种之间存在倒位等染色体变异，它的全基因组比对过程也是类似的。两个基因组进行从头到尾的碱基水平上的全基因组比对。

1.下载基因组的序列文件和参考基因组的注释文件

使用gean工具进行从sdi格式转化为fasta格式。

#如果没有安装GEAN，可以通过以下方式进行安装，这只是文件格式的转化，
如果直接下载的格式是fasta格式这个过程是不必要的。
git clone https://github.com/baoxingsong/GEAN.git
cd GEAN
cmake CMakeLists.txt
make

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命令进行合成一个。
cat TAIR10_chr1.fas TAIR10_chr2.fas TAIR10_chr3.fas TAIR10_chr4.fas TAIR10_chr5.fas > tair10.fa

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

如果得到的两个基因组的名称不一样则需要进行，名称规范化。

2. 提取参考基因组的中的保守序列 和 将保守序列分别匹配到两个基因组上

-x splice 代表将保守序列映射到基因组，-t 代表线程数，-k 代表kmer， -a代表输出sam格式的文件，-p和-N代表如果它们的比对得分大于初次比对得分的0.4将会保留最多保留20个结果，还要输入他们的基因组和保守序列。

anchorwave gff2seq -r tair10.fa -i TAIR10_GFF3_genes.gff -o cds.fa

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

3.对两个品种拟南芥基因组进行变异检测和全基因组比对

在这个例子中我们使用“genoAli”命令而不是上一篇文章中的“proali”-v输出变异检测结果的vcf文件，-i代表输入参考基因组的注释文件，-r 和-s 参考和查询基因组文件，-o和-f 输出变异比对结果的maf格式文件，-n文件可以用于可视化结果。如果过程中消耗的内存过大，可以减小滑动窗口的大小。此外，-m设置最小的外显子长度，不要随意修改，并且和上个步骤一样。

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 

对输出的共线性进行作图，代码如下

library(ggplot2)
#设置一个坐标函数
changetoM <- function ( position ){
  position=position/1000000;
  paste(position, "M", sep="")
}
#读取数据
data =read.table("bur_0.anchors", header=TRUE)
#筛选染色体并设置为因子
data = data[which(data$refChr %in% c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5")),]
data = data[which(data$queryChr %in% c("Chr1", "Chr2", "Chr3","Chr4","Chr5")),]
data$refChr = factor(data$refChr, levels=c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5"))
data$queryCh = factor(data$queryChr, levels=c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5" ))
#使用ggplot2包画图并美化它.
ggplot(data=data, aes(x=queryStart, y=referenceStart))+
  geom_point(size=0.5, aes(color=strand)) +
  facet_grid(refChr~queryChr, scales="free", space="free") +
  labs(x="Query_bur_0", y="Reference_tair10")+scale_x_continuous(labels=changetoM) +
  scale_y_continuous(labels=changetoM) +
  theme(axis.line = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(fill =NA,color="black", size=0.5, linetype="solid"),
        axis.text.y = element_text( colour = "black"),
        legend.position='none',
        axis.text.x = element_text(angle=300,hjust=0, vjust=0.5, colour = "black") )

 两基因组间不存在染色体重排全基因加倍等现象，输出的MAF文件说明两基因组间有snp，indels等变异现象.