pafr包的参考链接
https://cran.r-project.org/web/packages/pafr/vignettes/Introduction_to_pafr.html
首先用minimap2比对两个基因组
这里我用NCBI下载的两个拟南芥基因组做演示
下载两个基因组
代码语言：javascript复制wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_genomic.fna.gz
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/019/202/795/GCA_019202795.1_ASM1920279v1/GCA_019202795.1_ASM1920279v1_genomic.fna.gz
解压重命名
代码语言：javascript复制gunzip GCA_019202795.1_ASM1920279v1_genomic.fna.gz
mv GCA_019202795.1_ASM1920279v1_genomic.fna query.fna
grep ">" query.fna | wc -l # 这个里有218 条序列
gunzip GCF_000001735.4_TAIR10.1_genomic.fna.gz
mv GCF_000001735.4_TAIR10.1_genomic.fna target.fna
grep ">" target.fna | wc -l ## 这个里有7条序列
minimap2的安装
直接使用conda安装 https://github.com/lh3/minimap2
代码语言：javascript复制conda install minimap2
比对
代码语言：javascript复制minimap2 -ax asm5 target.fna query.fna > arabidopsis_aln.paf
这个最终的比对结果有900多兆，自己的电脑R语言读取应该很吃力，下面的操作还是使用这个R包自带的数据吧
接下来是R语言里的操作
安装pafr包
代码语言：javascript复制install.packages("pafr")
加载需要用到的R包
代码语言：javascript复制library(pafr)
library(tidyverse)
library(ggplot2)
读取数据
代码语言：javascript复制fungi.paf.1<-read_paf("fungi.paf")
fungi.paf.1 %>% as.data.frame() %>% head()
fungi.paf.2<-read_paf("fungi.paf",include_tags = FALSE)
fungi.paf.2 %>% as.data.frame() %>% head()
共线性的点图
代码语言：javascript复制dotplot(fungi.paf.2)
image.png
指定染色体的共线性
代码语言：javascript复制plot_synteny(fungi.paf.2,
q_chrom="Q_chr3",
t_chrom="T_chr4",
centre=TRUE)
theme_bw()
image.png
这里有个问题好像是只能1对1，研究研究他的代码，看看能不能改成可以多对一
覆盖度
代码语言：javascript复制plot_coverage(fungi.paf.2) -> p1
plot_coverage(fungi.paf.2,fill='qname') -> p2
plot_coverage(fungi.paf.2,fill='qname')
scale_fill_brewer(palette = "Set1") -> p3
library(patchwork)
p1
p2 theme(legend.position = "none")
p3 theme(legend.position = "none")
image.png
