论文是 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()函数保存
上半部分具体的数据格式可以参考之前的推文
跟着Nature Genetics学画图：R语言ggbio包画基因结构图
下半部分的数据格式
image.png
这个原图中有7个品种，我这边就不全部准备了，我这边只准备3个
- 第一列是品种的名字
- 第二列是snp的位置
- 第三列是snp在图上的y轴位置，从-1开始，每多一个品种就减一
- 第四列是碱基类型
- 第五列是碱基的分类 A代表 变异的碱基，R是参考序列的碱基
第一步是加载需要用到的R包
代码语言：javascript复制library(ggh4x)
library(ggplot2)
library(ggbio)
library(GenomicRanges)
第二步是准备作图数据
代码语言：javascript复制df<-read.csv("NG/waxy.gtf",
header=F,
sep="t")
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")
这个数据的内容用文字表达可能得写好多内容，后面争取出视频内容进行介绍
最后是画图代码
代码语言：javascript复制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="Nipponbaren(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
好了，今天的内容就先到这里了
注视频号 小明的数据分析笔记本，下载链接就在下面视频的评论区
