之前分享过一篇推文介绍过这个内容 R语言ggplot2包画曼哈顿图的一个简单小例子,但是当时自己不太懂曼哈顿图,实现是直接借助ggplot2的geom_jitter()
这个函数实现的。这个函数并不会考虑每个变异位点的位置,而实际的曼哈顿图是需要根据变异位点的位置来画的。今天的推文重新介绍一下ggplot2绘制曼哈顿图的代码。数据集就使用之前的推文中用到的数据跟着Nature Genetics学GWAS分析:emmax软件gwas分析/qqman包展示结果,这个数据太大,出图有些慢,只随机选取了其中1%的数据 (这个数据我自己的存储路径population.genomics/gwas/NG.tomato/at/)。
R语言中也有现成的包和函数可以直接画曼哈顿图,我这里选择用ggplot2来画是因为出图后可以非常方便的组合其他的图,比如可以叠加一个基因结构的图,然后再拼一个展示不同基因型表型差异的图。这些如果是用ggplot2来做,都可以用代码实现,省去了手动拼图的过程。
image.png
然后还需要准备一个染色体长度的文件
image.png
library(tidyverse)
library(readxl)
library(ggrastr)
gwas.results<-read_delim("D:/R_4_1_0_working_directory/env001/data/20230912/at_gwas01.txt",
delim = "\t")
gwas.results %>% dim()
chr.len<-read_excel("D:/R_4_1_0_working_directory/env001/data/20230912/at_gwas.xlsx")
chr.len %>% colnames()
gwas.results %>% colnames()
接下来的思路是给每个变异位点的位置加上染色体的长度,如果是1号染色体就加0,2号染色体的变异位点给加上1号染色体的长度,3号染色体的变异位点位置加上(1+2)号染色体的长度,以此类推
chr.len %>% pull(LEN) %>% cumsum() -> x1
chr.len %>% pull(LEN) -> x2
head(x1,-1)
temp.df<-data.frame(chromo=chr.len %>% pull(CHR),
chr_len=c(0,head(x1,-1)))
gwas.results %>%
left_join(temp.df,by=c("CHR"="chromo")) %>%
mutate(new_pos=BP+chr_len) -> new.gwas.results
构造好数据就可以画图了
cols<-c("#bdbadb","#fdb363","#f98177","#80b3d4","gray")
pdf(file = "gwas_man.pdf",width = 8,height = 5)
ggplot(data=new.gwas.results,
aes(x=new_pos,y=-log10(P)))+
geom_point_rast(aes(color=CHR),show.legend = FALSE)+
scale_x_continuous(breaks =c(0,head(x1,-1)) + x2/2 ,
labels =temp.df %>% pull(chromo))+
theme_bw()+
theme(panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line())+
scale_color_manual(values = cols)+
scale_y_continuous(expand = expansion(mult = c(0,0)),
limits = c(0,5),
breaks = c(0,1,2,3,4,5))+
labs(x=NULL)
dev.off()
image.png
如果觉得每条染色体之间比较拥挤就给每条染色体的长度多加一些,我这里的示例数据是每条染色体的长度多加了10M
chr.len<-read_excel("D:/R_4_1_0_working_directory/env001/data/20230912/at_gwas.xlsx")
chr.len %>% colnames()
chr.len %>%
mutate(LEN=LEN+10000000) -> chr.len
chr.len %>% pull(LEN) %>% cumsum() -> x1
chr.len %>% pull(LEN) -> x2
temp.df<-data.frame(chromo=chr.len %>% pull(CHR),
chr_len=c(0,head(x1,-1)))
gwas.results %>%
left_join(temp.df,by=c("CHR"="chromo")) %>%
mutate(new_pos=BP+chr_len) -> new.gwas.results
cols<-c("#bdbadb","#fdb363","#f98177","#80b3d4","gray")
pdf(file = "gwas_man01.pdf",width = 8,height = 5)
ggplot(data=new.gwas.results,
aes(x=new_pos,y=-log10(P)))+
geom_point_rast(aes(color=CHR),show.legend = FALSE)+
scale_x_continuous(breaks =c(0,head(x1,-1)) + x2/3 ,
labels =temp.df %>% pull(chromo))+
theme_bw()+
theme(panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line())+
scale_color_manual(values = cols)+
scale_y_continuous(expand = expansion(mult = c(0,0)),
limits = c(0,5),
breaks = c(0,1,2,3,4,5))+
labs(x=NULL) -> p2
print(p2)
dev.off()
image.png
这样看起来好像美观一点