ggplot2画图——将筛选到的突变位点（基因）画到染色体上

今天需要画一个图，把筛选到的SNP画到染色体上，可惜在网上逛了一圈一直找不到心仪的r包，只好自己动手了。成图如下所示：

首先需要一个关于染色体总长的文件，类似下图，具体怎么获得大家可以各显神通，这边推荐用tbtools。

 

另外需要存有SNP位点信息的文件，至少需要包含位于哪条染色体，位于染色体上的哪个位置两个信息，类似下图：

 处理一下染色体长度文件：

chr <- read.table("chrlength.txt") #读取文件
colnames(chr) <- c("Chr","End") #命名列名
chr <- chr[order(chr$Chr),] #根据染色体名对数据进行排序
chr$x <- seq(1,56) + rep(seq(0,26,2),each=4) #计算染色体在图上的横坐标。这里我有14*4条的四倍型染色体，每条染色体宽为1，共56条染色体所以用seq(1,56)生成每条染色体的横坐标。且我希望同源染色体紧挨着，非同源染色体之间相隔两倍的染色体宽度，所以只要生成一个类似于：0000111122223333的数组然后和seq(1,56)相加就可以实现，这里用的方法是：rep(seq(0,26,2),each=4)

head(chr)
      Chr      End x
41 Chr01A 47186193 1
42 Chr01B 45237200 2
43 Chr01C 45157756 3
44 Chr01D 39476532 4
17 Chr02A 37619394 7
18 Chr02B 35749136 8

这样就可以用ggplot2画染色体了，试试看：

 ggplot()+
     geom_rect(chr,mapping=aes(ymin=0,ymax=End,xmin=x-0.5,xmax=x+0.5),
color="black",fill="white")
#用了geom_rect画矩形，柱状图没准也可以，有兴趣的可以去试试。
#geom_rect四个参数控制矩形的四条边ymin,ymax,xmin,xmax

 效果还可以，再处理一下SNP的文件：

#我有4组突变文件所以要读取四次，分别是与温度\降水相关的snp\indel
snptem <- read.table("tem_snp_uniq.txt")
snpwet <- read.table("wet_snp_uniq.txt")
indeltem <- read.table("tem_indel_uniq.txt")
indelwet <- read.table("wet_indel_uniq.txt")

#给SNP上个标记叫SNP，indel叫Indel,再给温度相关的标记为tem，降水相关叫wet,最后用rbind()合并为一个数据框
snptem$type <- snpwet$type <- "SNP" ;indeltem$type <- indelwet$type <- "Indel"
snptem$type2 <- indeltem$type2 <- "tem" ;indelwet$type2 <- snpwet$type2 <- "wet"
variation <- rbind(snptem,snpwet,indeltem,indelwet)
#由于我的SNP位置信息是Chr01A:100这种形式的，由:分割，前者是染色体，后者是在染色体上的位置
#用tidyr包的separate函数分割一下
library(tidyr)
variation <- separate(variation,col = V1,into = c("chr","pos"),sep = ":")
variation$pos <- as.numeric(variation$pos) #转化为数字格式

head(variation)
     chr      pos type type2
1 Chr01A  8203552  SNP   tem
2 Chr01B 30258368  SNP   tem
3 Chr01B 30258404  SNP   tem
4 Chr01B 30258416  SNP   tem
5 Chr01B 30258460  SNP   tem
6 Chr01B 30258484  SNP   tem

#由于有些SNP所在位置不在核染色体上，所以要去除，这里用%in%判断
variation <- variation[variation$chr %in% chr$Chr,]
#再根据上一个数据框chr里算好的染色体位置得到SNP的横坐标
variation$x <- chr$x[match(variation$chr , chr$Chr)] 

head(variation)
     chr      pos type type2 x
1 Chr01A  8203552  SNP   tem 1
2 Chr01B 30258368  SNP   tem 2
3 Chr01B 30258404  SNP   tem 2
4 Chr01B 30258416  SNP   tem 2
5 Chr01B 30258460  SNP   tem 2
6 Chr01B 30258484  SNP   tem 2

数字已经处理完成了，开始画图：

ggplot()+
   geom_rect(chr,mapping=aes(ymin=0,ymax=End,xmin=x-0.5,xmax=x+0.5),
color="black",fill="white")+
   geom_segment(data = variation, mapping=aes(x = x-0.5, y = pos, xend = x+0.5, yend = pos, colour = type2),
size=1.2,alpha=0.7)
#在原先geom_rect的基础上加了geom_segment用来绘制线段。
#geom_segment也是四个参数控制线段起点和终点，x,y,xend,yend。x和xend保持和geom_rect一致，y用pos行控制

 效果不错，美化一下：

install.packages("ggchicklet", repos = "https://cinc.rud.is")
library(ggchicklet)
ggplot()+
   ggchicklet:::geom_rrect(chr,mapping=aes(ymin=0,ymax=End,xmin=x-0.5,xmax=x+0.5),color="black",fill="white",r = unit(0.5, 'npc'))+
   geom_segment(data = variation, mapping=aes(x = x-0.5, y = pos, xend = x+0.5, yend = pos, colour = type2),size=1.2,alpha=0.7)+
   geom_point(data = variation[variation$type=="Indel",], mapping=aes(x = x, y = pos, fill = type2 ),colour=rgb(0,0,0,0),shape=24,size=2)+
   scale_color_manual(values = c("#A52A2A", "#0073C2"))+
   scale_fill_manual(values = c("#FF0000", "#1E90FF"))+
   theme(legend.position = 'none',
         panel.grid = element_blank(),panel.background = element_blank(),
         axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())

这里用到了ggchicklet这个包，主要是想要让染色体边缘变成圆角的，具体看参考。

成图如下所示：

 我非常满意，大功告成！

参考：

1、ggplot2画图——点图想要描边但是又需要去除边框

2、ggplot2绘制精美圆角矩形图加点图