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

anchorwave对存在染色体融合，易位，全基因组复制的基因组进行比对，本文以斑马鱼和金鱼为例进行基因组比对。

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

wget ftp://ftp.ensembl.org/pub/release-102/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-102/gff3/danio_rerio/Danio_rerio.GRCz11.102.chr.gff3.gz
wget https://research.nhgri.nih.gov/goldfish/download/carAur01.sm.fa
gunzip *gz

2.查看基因组间的共线性和全基因组复制情况

2.1使用samtools概览两个基因组

通过使用samtools的faidx命令，我们可以得到fai格式的文件，知道斑马鱼有25条染色体，金鱼有50条常染色体，推测金鱼之间经历了一次全基因组复制，并且第15923983行是金鱼的最后一条染色体，进行重新命名。

samtools faidx carAur01.sm.fa
samtools faidx Danio_rerio.GRCz11.dna.primary_assembly.fa

head -15923983 carAur01.sm.fa > goldfish.fa

2.2提取参考基因组的中的保守序列

这一步需要使用必须使用anchorwave，这样可以过滤一些不必要的anchors，最小化的减少下一步中minimap2的限制。-r代表要输入的参考基因组，-i代表他的注释文件，-o代表他输出的所有的保守序列的文件。

anchorwave gff2seq -r Danio_rerio.GRCz11.dna.primary_assembly.fa -i Danio_rerio.GRCz11.102.chr.gff3 -o cds.fa

2.3将保守序列（本例中为蛋白质编码序列）分别映射到两个基因组上

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

minimap2 -x splice -t 11 -k 12 -a -p 0.4 -N 20 Danio_rerio.GRCz11.dna.primary_assembly.fa cds.fa > ref.sam
minimap2 -x splice -t 11 -k 12 -a -p 0.4 -N 20 goldfish.fa cds.fa > carAur.sam

2.4使用anchorwave进行共线性的识别

proali 代表对有染色体结构变异如易位，染色体融合，全基因组复制的基因组进行比对。-i代表输入参考基因组的注释文件，-r 和-s 参考和查询基因组文件，-a 和-ar映射到两个基因组的结果，-R代表参考基因组覆盖度是2，而查询基因组经历一次全基因组复制事件，所以-Q为1，-ns代表不产生新的anchor。

之后就可以利用输出的align1.anchors进行共线性区域可视化了。

anchorwave proali -i Danio_rerio.GRCz11.102.chr.gff3 -r Danio_rerio.GRCz11.dna.primary_assembly.fa -a carAur.sam -as cds.fa -ar ref.sam -s goldfish.fa -n align1.anchors -R 2 -Q 1 -ns

使用R语言中的ggplot2进行绘图，代码如下

library(ggplot2)
changetoM <- function ( position ){
position=position/1000000;
paste(position, "M", sep="")
}

data =read.table("align1.anchors", head=TRUE)
data = data[which(data$refChr %in% c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20",
"21", "22", "23", "24", "25", "MT" )),]

data = data[which(data$queryChr %in% c("LG1","LG2","LG3", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10", "LG11", "LG12", "LG13", "LG14", "LG15",
"LG16", "LG17", "LG18", "LG19", "LG20", "LG21", "LG22", "LG23", "LG24", "LG25", "LG26", "LG27", "LG28","LG28B", "LG29", "LG30", "LG30F", "LG31", "LG32", "LG33", "LG34", "LG35", "LG36", "LG36F", "LG37", "LG37M","LG38", "LG39", "LG40", "LG41", "LG42", "LG42F", "LG43", "LG44", "LG44F", "LG45", "LG45M", "LG46", "LG47",
"LG48", "LG48F", "LG49", "LG49B", "LG50")),]

data$refChr = factor(data$refChr, levels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19",
"20", "21", "22", "23", "24", "25", "MT" ))

data$queryChr = factor(data$queryChr, levels=c("LG1","LG2","LG3", "LG4", "LG5", "LG6", "LG7", "LG8", "LG9", "LG10", "LG11", "LG12", "LG13", "LG14",
"LG15", "LG16", "LG17", "LG18", "LG19", "LG20", "LG21", "LG22", "LG23", "LG24", "LG25", "LG26","LG27", "LG28", "LG28B", "LG29", "LG30", "LG30F", "LG31", "LG32", "LG33", "LG34", "LG35","LG36","LG36F", "LG37", "LG37M", "LG38", "LG39", "LG40", "LG41", "LG42", "LG42F", "LG43", "LG44", "LG44F","LG45", "LG45M", "LG46", "LG47", "LG48", "LG48F", "LG49", "LG49B", "LG50" ))

ggplot(data=data, aes(x=queryStart, y=referenceStart))+geom_point(size=0.5, aes(color=strand))+
facet_grid(refChr~queryChr, scales="free", space="free" )+ theme_grey(base_size = 120) +
labs(x="goldfish", y="zebrafish")+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=1, colour = "black") )

绘图如下，共线性证明了两个染色体之间存在全基因组复制，染色体重排等，

 

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

proali 代表对有染色体结构变异如易位，染色体融合，全基因组复制的基因组进行比对。-i代表输入参考基因组的注释文件，-r 和-s 参考和查询基因组文件，-a 和-ar映射到两个基因组的结果，-R代表参考基因组覆盖度是2，而查询基因组经历一次全基因组复制事件，所以-Q为1，不写-ns代表产生新的anchor。-t 代表线程为10，-o和-f均为输出的文件。

anchorwave proali -i Danio_rerio.GRCz11.102.chr.gff3 -r Danio_rerio.GRCz11.dna.primary_assembly.fa -a carAur.sam -as cds.fa -ar ref.sam -s goldfish.fa -n align1.anchors -o align1.maf -t 10 -R 2 -Q 1 -f align1.f.maf  > fishlog1 2>&1