R语言里可视化多序列比对（paf格式）的R包：pafr

pafr包的参考链接

https://cran.r-project.org/web/packages/pafr/vignettes/Introduction_to_pafr.html

首先用minimap2比对两个基因组

这里我用NCBI下载的两个拟南芥基因组做演示

下载两个基因组

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

解压重命名

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

conda install minimap2

比对

minimap2 -ax asm5 target.fna query.fna > arabidopsis_aln.paf

这个最终的比对结果有900多兆，自己的电脑R语言读取应该很吃力，下面的操作还是使用这个R包自带的数据吧

这里修正一个错误 -ax参数最终结果是sam格式 -cx参数才是paf格式

minimap2 -cx asm5 target.fna query.fna > arabidopsis_aln.paf

接下来是R语言里的操作

安装pafr包

install.packages("pafr")

加载需要用到的R包

library(pafr)
library(tidyverse)
library(ggplot2)

读取数据

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()

共线性的点图

dotplot(fungi.paf.2)

image.png

指定染色体的共线性

plot_synteny(fungi.paf.2, 
             q_chrom="Q_chr3", 
             t_chrom="T_chr4", 
             centre=TRUE) +
  theme_bw()

image.png

这里有个问题好像是只能1对1，研究研究他的代码，看看能不能改成可以多对一

覆盖度

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

今天推文的示例数据和代码可以在公众号后台回复20220317获取

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！