論文是 Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice
今天重復(fù)的圖是來自于論文中的figure3 b
之前的推文已經(jīng)介紹過 上半部分的基因結(jié)果的畫法羞延，
今天的推文介紹下半部分SNP位點(diǎn)的堿基類型的實(shí)現(xiàn)辦法，背景顏色這里借助的是ggplot2包中的geom_tile()函數(shù)导饲；表示堿基的文本借助的是geom_text()函數(shù)
這里最開始的思路是借助aplot這個(gè)包的拼圖功能實(shí)現(xiàn)的，但是上下兩個(gè)部分拼接的時(shí)候遇到了報(bào)錯(cuò)广匙，使用patchwork拼接的時(shí)候也遇到了報(bào)錯(cuò)，報(bào)錯(cuò)的內(nèi)容忘記保存了，暫時(shí)不知道如何解決，使用ggbio這個(gè)包做的圖可以繼續(xù)使用ggplot2的函數(shù)疊加幻捏，但是如果使用ggplot2的拼圖方式卻不行。使用ggbio這個(gè)包做的圖也不能使用ggsave()函數(shù)保存
上半部分具體的數(shù)據(jù)格式可以參考之前的推文
下半部分的數(shù)據(jù)格式
這個(gè)原圖中有7個(gè)品種命咐，我這邊就不全部準(zhǔn)備了篡九，我這邊只準(zhǔn)備3個(gè)
- 第一列是品種的名字
- 第二列是snp的位置
- 第三列是snp在圖上的y軸位置，從-1開始醋奠，每多一個(gè)品種就減一
- 第四列是堿基類型
- 第五列是堿基的分類 A代表 變異的堿基榛臼，R是參考序列的堿基
第一步是加載需要用到的R包
library(ggh4x)
library(ggplot2)
library(ggbio)
library(GenomicRanges)
第二步是準(zhǔn)備作圖數(shù)據(jù)
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")
這個(gè)數(shù)據(jù)的內(nèi)容用文字表達(dá)可能得寫好多內(nèi)容，后面爭取出視頻內(nèi)容進(jìn)行介紹
最后是畫圖代碼
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()
這個(gè)是最終的結(jié)果
好了窜司，今天的內(nèi)容就先到這里了
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1沛善、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子；2塞祈、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)金刁、基因組學(xué)、群體遺傳學(xué)文獻(xiàn)閱讀筆記织咧；3胀葱、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記漠秋！
示例數(shù)據(jù)和代碼獲取步驟
公眾號(hào)的推文寫了獲取示例數(shù)據(jù)和代碼的步驟笙蒙，可以到公眾號(hào)查看對(duì)應(yīng)推文
