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

anchorwave以玉米B73和Mo17为例对存在倒位的基因组

进行全基因组比对。

玉米基因组大约有2300M，在操作过程中请注意自己的内存使用情况。

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

下载玉米玉米B73品种的基因组和注释文件和Mo17品种的基因组，并解压他们。

wget ftp://ftp.ensemblgenomes.org/pub/plants/release-34/fasta/\
zea_mays/dna/Zea_mays.AGPv4.dna.toplevel.fa.gz
gunzip Zea_mays.AGPv4.dna.toplevel.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-34/gff3/\
zea_mays/Zea_mays.AGPv4.34.gff3.gz
gunzip Zea_mays.AGPv4.34.gff3.gz
wget https://download.maizegdb.org/Zm-Mo17-REFERENCE-CAU-1.0\
/Zm-Mo17-REFERENCE-CAU-1.0.fa.gz
gunzip Zm-Mo17-REFERENCE-CAU-1.0.fa.gz

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

使用之前确保这两个基因组的命名方法是一致的。使用grep命令查看使用染色体名称是否相同，使用sed命令进行重新命名。使用minimap2过程中消耗的内存太大，可以降低线程数（-t）。

grep ">" Zm-Mo17-REFERENCE-CAU-1.0.fa
grep ">" Zea_mays.AGPv4.dna.toplevel.fa
sed -i 's/>chr/>/g' Zm-Mo17-REFERENCE-CAU-1.0.fa
anchorwave gff2seq -i Zea_mays.AGPv4.34.gff3 \
-r Zea_mays.AGPv4.dna.toplevel.fa -o cds.fa
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 \
Zm-Mo17-REFERENCE-CAU-1.0.fa cds.fa > cds.sam
minimap2 -x splice -t 10 -k 12 -a -p 0.4 -N 20 \
Zea_mays.AGPv4.dna.toplevel.fa cds.fa > ref.sam

3.使用anchorwave进行共线性的识别

-IV代表不识别新的anchor，识别结果输出一个align.anchors文件，此文件可以用于共线性作图。

以下使用ggplot2进行作图。

anchorwave genoAli -i Zea_mays.AGPv4.34.gff3 -as cds.fa \
-r Zea_mays.AGPv4.dna.toplevel.fa -a cds.sam -ar ref.sam \
-s Zm-Mo17-REFERENCE-CAU-1.0.fa -n align.anchors -IV

#使用R画一个点图
#加载ggplot2包，为坐标名设置一个规范的刻度值。
library(ggplot2)
changetoM <- function ( position ){
  position=position/1000000;
  paste(position, "M", sep="")
}
#读取共显性区以及它们属于哪个染色体
data =read.table("align.anchors", header=TRUE)
#将选取的染色体设置为因子
data = data[which(data$refChr %in% c("1", "2", "3", "4",
                                     "5", "6", "7", "8", "9", "10")),]
data = data[which(data$queryChr %in% c("1", "2", "3", "4",
                                       "5", "6", "7", "8", "9", "10")),]
data$refChr = factor(data$refChr, levels=c("1", "2", "3", "4",
                                           "5", "6", "7", "8", "9", "10"))
data$queryCh = factor(data$queryChr, levels=c("1", "2", "3", "4",
                                              "5", "6", "7", "8", "9", "10"))
#使用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="qurymaize", y="referencemaize")+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( size=5,colour = "black"),
        legend.position='none',
        axis.text.x = element_text(angle=300,size=6,hjust=0, vjust=0.5, colour = "black") )

蓝色代表正常匹配结果，红色发生了倒位等染色体重排现象 

4.使用anchorwave进行全基因组比对

此过程anchorwave使用最快速的WFA算法和a 2-piece affine gap cost 策略进行单碱基分辨率的全基因组比对。此过程消耗的时间和内存可能较大，建议大型计算机上运行。

```

使用genoAli请注意，ref和qurey的染色体的名字要相同

```

anchorwave genoAli -i Zea_mays.AGPv4.34.gff3 -as cds.fa\
-r Zea_mays.AGPv4.dna.toplevel.fa -a cds.sam -ar ref.sam\
-s Zm-Mo17-REFERENCE-CAU-1.0.fa -n anchors\ 
-o anchorwave.maf -f anchorwave.f.maf 

更多详细内容：

相关论文如下：https://www.pnas.org/doi/10.1073/pnas.2113075119

其他代码如下：https://github.com/baoxingsong/anchorwave