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
這個最終的比對結(jié)果有900多兆，自己的電腦R語言讀取應(yīng)該很吃力爱榕，下面的操作還是使用這個R包自帶的數(shù)據(jù)吧
這里修正一個錯誤 -ax參數(shù)最終結(jié)果是sam格式 -cx參數(shù)才是paf格式
minimap2 -cx asm5 target.fna query.fna > arabidopsis_aln.paf
接下來是R語言里的操作
安裝pafr包
install.packages("pafr")
加載需要用到的R包
library(pafr)
library(tidyverse)
library(ggplot2)
讀取數(shù)據(jù)
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
今天推文的示例數(shù)據(jù)和代碼可以在公眾號后臺回復(fù)20220317獲取
歡迎大家關(guān)注我的公眾號
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子昨凡；2土匀、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)就轧、基因組學(xué)田度、群體遺傳學(xué)文獻閱讀筆記镇饺；3乎莉、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！
