R语言ggplot2绘制曼哈顿图展示GWAS分析的结果

用户7010445

发布于 2023-09-21 08:25:02

7880

发布于 2023-09-21 08:25:02

之前分享过一篇推文介绍过这个内容 R语言ggplot2包画曼哈顿图的一个简单小例子，但是当时自己不太懂曼哈顿图，实现是直接借助ggplot2的geom_jitter()这个函数实现的。这个函数并不会考虑每个变异位点的位置，而实际的曼哈顿图是需要根据变异位点的位置来画的。今天的推文重新介绍一下ggplot2绘制曼哈顿图的代码。数据集就使用之前的推文中用到的数据跟着Nature Genetics学GWAS分析：emmax软件gwas分析/qqman包展示结果，这个数据太大，出图有些慢，只随机选取了其中1%的数据 (这个数据我自己的存储路径population.genomics/gwas/NG.tomato/at/)。

R语言中也有现成的包和函数可以直接画曼哈顿图，我这里选择用ggplot2来画是因为出图后可以非常方便的组合其他的图，比如可以叠加一个基因结构的图，然后再拼一个展示不同基因型表型差异的图。这些如果是用ggplot2来做，都可以用代码实现，省去了手动拼图的过程。

首先是gwas结果的部分截图

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()