論文是 Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice

image.png

今天重復(fù)的圖是來自于論文中的figure3 b

image.png

之前的推文已經(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ù)格式

image.png

這個(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é)果

image.png

好了，今天的內(nèi)容就先到這里了

歡迎大家關(guān)注我的公眾號

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：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ù)和代碼獲取步驟

公眾號的推文寫了獲取示例數(shù)據(jù)和代碼的步驟，可以到公眾號查看對應(yīng)推文