论文是 Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice
image.png
今天重复的图是来自于论文中的figure3 b
image.png
之前的推文已经介绍过 上半部分的基因结果的画法，
今天的推文介绍下半部分SNP位点的碱基类型的实现办法，背景颜色这里借助的是ggplot2包中的geom_tile()函数；表示碱基的文本借助的是geom_text()函数
这里最开始的思路是借助aplot这个包的拼图功能实现的，但是上下两个部分拼接的时候遇到了报错，使用patchwork拼接的时候也遇到了报错，报错的内容忘记保存了，暂时不知道如何解决，使用ggbio这个包做的图可以继续使用ggplot2的函数叠加，但是如果使用ggplot2的拼图方式却不行。使用ggbio这个包做的图也不能使用ggsave()函数保存
上半部分具体的数据格式可以参考之前的推文
下半部分的数据格式
image.png
这个原图中有7个品种，我这边就不全部准备了，我这边只准备3个
- 第一列是品种的名字
- 第二列是snp的位置
- 第三列是snp在图上的y轴位置，从-1开始，每多一个品种就减一
- 第四列是碱基类型
- 第五列是碱基的分类 A代表 变异的碱基，R是参考序列的碱基
第一步是加载需要用到的R包
library(ggh4x)
library(ggplot2)
library(ggbio)
library(GenomicRanges)
第二步是准备作图数据
waxy<-GRanges("chr06",IRanges(df$V4,end=df$V5,group=df$V3))
waxy_snp<-read.csv("NG/waxy_snp.csv")
head(waxy_snp)
cultivar<-data.frame(x=1765000,
y=unique(waxy_snp$y_location),
label=unique(waxy_snp$cultivars))
snp<-data.frame(xmin=unique(waxy_snp$x_location),
xmax=unique(waxy_snp$x_location),
ymin=0.6,
ymax=1.4)
snp_segment<-data.frame(xmin=unique(waxy_snp$x_location),
xmax=unique(waxy_snp$x_location),
ymin=-0.5,
ymax=0.6)
atg<-data.frame(x=1766921,y=1.5,label="ATG")
这个数据的内容用文字表达可能得写好多内容，后面争取出视频内容进行介绍
最后是画图代码
pdf(file = "NG/waxy-2.pdf",width = 12,height = 4)
autoplot(waxy,aes(fill=group),geom="alignment")+
theme_bw()+
scale_x_continuous(limits = c(1764000,1771000),
breaks = c(seq(1764000,1771000,by=1000)),
position = "top")+
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.line.x = element_line(),
axis.ticks.length = unit(0.2,'cm'))+
guides(x=guide_axis_truncated(trunc_lower = 1764000,
trunc_upper = 1771000))+
scale_fill_manual(values = c("#92d050","#a6a6a6","#a6a6a6"))+
#theme(aspect.ratio = 0.1)+
#scale_y_continuous(limits = c(0,3))+
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank())+
annotate(geom = "text",x=1764500,y=1,
label="Nipponbare\n(waxy:Os06g0133000)")+
coord_cartesian(clip="off")+
ggnewscale::new_scale_fill()+
geom_tile(data=waxy_snp,aes(x=x_location,
y=y_location,
fill=type),
width=100)+
scale_fill_manual(values = c("#ffccff","#99ccff"),
labels=c("Alternative type",
"Reference type"))+
geom_text(data = waxy_snp,aes(x=x_location,
y=y_location,
label=label))+
geom_text(data=cultivar,aes(x=x,y=y,label=label))+
geom_segment(data=snp,aes(x=xmin,xend=xmax,y=ymin,yend=ymax),
color="red")+
geom_segment(data=snp_segment,aes(x=xmin,
xend=xmax,
y=ymin,
yend=ymax),
lty="dashed",color="grey")+
geom_label(data=atg,aes(x=x,y=y,label=label),
fill="blue",
color="white",
label.r = unit(0,'mm'),
nudge_y = 0.3)
dev.off()
这个是最终的结果
image.png
好了，今天的内容就先到这里了
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
示例数据和代码获取步骤
公众号的推文写了获取示例数据和代码的步骤，可以到公众号查看对应推文
