R可视化——曼哈顿图

曼哈顿图(Manhattan plot)是一种GWAS分析中常用的展示基因组数据的散点图。今天，小编就带大家看一下如何使用R语言中的qqman、CMplot及ggplot2三个包绘制曼哈顿图。

基于qqman包绘制

1、安装、加载R包

#安装包#install.package("qqman")#加载包library(qqman)

2、数据

#加载数据——以qqman包自带示例数据gwasresult数据为例df1 <- gwasResults#数据中SNP为SNP名称，CHR为染色体编号，BP为碱基位置，P为p值head(df1)#预览数据

3、绘图

manhattan(df1,#数据 col = c('#30A9DE','#EFDC05','#E53A40','#090707'),#交替使用颜色展示 suggestiveline = -log10(1e-05),#－log10(1e－5)处添加"suggestive"横线 genomewideline = -log10(5e-08),#－log10(5e－10)处添加"genome-wide sigificant"横线 highlight = snpsOfInterest,#内置高亮的snp数据， 也可以对snpOfInterest进行设置 annotatePval = 0.05,#标记p值小于0.05的点 annotateTop = T,#如果为T，则仅批注低于注解阈值的每个染色体上的顶部点，为F则标记所有小于注解阈值的点。 main = "XXXXXXXX"#标题 )##可以使用??manhattan查看具体参数

基于CMplot包绘制

1、安装、加载包

#安装包#install.package("CMplot")#加载包library(CMplot)

2、数据

#加载数据——以CMplot包自带示例数据pig60k数据为例data(pig60K)#预览数据

3、绘图‍

#在工作目录下会生成一系列图CMplot(pig60K,#示例数据 chr.den.col=c("black","green","red"),#SNP密度展示 file="jpg",#绘制图片类型 memo="",#输出文件名中添加一个字符 dpi = 600)#绘制图片的分辨率 #可以使用??CMplot查看具体参数,根据需要进行设置CMplot(Pmap, col=c("#4197d8", "#f8c120", "#413496", "#495226", "#d60b6f", "#e66519", "#d581b7", "#83d3ad", "#7c162c", "#26755d"), bin.size=1e6, bin.range=NULL, bin.legend.num=10, pch=19, type="p", band=1, H=1.5, ylim=NULL, cex.axis=1, lwd.axis=1.5, cex.lab=1.5, plot.type="b", multracks=FALSE, points.alpha=100L, cex=c(0.5,1,1), r=0.3, outward=FALSE, ylab=expression(-log[10](italic(p))), ylab.pos=3, xticks.pos=1, mar = c(3,6,3,3), threshold = NULL, threshold.col="red", threshold.lwd=1, threshold.lty=2, amplify= TRUE, signal.cex = 1.5, signal.pch = 19, signal.col=NULL, signal.line=2, highlight=NULL, highlight.cex=1, highlight.pch=19, highlight.type="p", highlight.col="red", highlight.text=NULL, highlight.text.col="black", highlight.text.cex=1, highlight.text.xadj=NULL, highlight.text.yadj=NULL, highlight.text.font=3, chr.labels=NULL, chr.border=FALSE, chr.labels.angle=0, chr.den.col="black", chr.pos.max=FALSE, cir.band=1, cir.chr=TRUE, cir.chr.h=1.5, cir.legend=TRUE, cir.legend.cex=0.6, cir.legend.col="black", LOG10=TRUE, box=FALSE, conf.int=TRUE, conf.int.col=NULL, file.output=TRUE, file=c("jpg","pdf","tiff"), dpi=300, height=NULL, width=NULL, memo="", main="", main.cex=1.5, main.font=2, trait.legend.ncol=NULL, verbose=TRUE)

基于ggplot2包进行绘制

1、安装、加载包

#安装包# install.packages("ggplot2")# install.packages("tidyverse")# install.packages("ggforce")# install.packages("ggprism")#加载包library(ggplot2)library(tidyverse)library(ggforce)library(ggprism)

2、数据及其处理

#数据，同样以qqman包自带数据gwasresult为例df <- gwasResults###数据处理df %>% group_by(CHR) %>% summarise(df_chr_len=max(BP)) %>% #计算染色体长度 mutate(total = cumsum(df_chr_len) - df_chr_len) %>% select(-df_chr_len) %>% #计算染色体初始位置 left_join(df, ., by="CHR") %>% arrange(CHR, BP) %>% mutate( BPcum = BP + total)->df_SNP_position#计算累计SNP的位置
head(df_SNP_position)#预览数据

3、绘图

#X轴标签位置X_axis <- df_SNP_position %>% group_by(CHR) %>% summarize(center=( max(BPcum) +min(BPcum) ) / 2 )#添加高亮和注释信息：snpsOfInterest中的rs编号和P值大于10的点data <- df_SNP_position %>% mutate( is_highlight=ifelse(SNP %in% snpsOfInterest, "yes", "no")) %>% mutate( is_annotate=ifelse(-log10(P)>10, "yes", "no"))#绘图ggplot(data, aes(x=BPcum, y=-log10(P))) + geom_point(aes(color=as.factor(CHR)),alpha=0.8, size=1.5)+ scale_color_manual(values = rep(c('#30A9DE','#EFDC05','#E53A40','#090707'), 22 ))+#颜色设置 scale_x_continuous(label = X_axis$CHR, breaks= X_axis$center)+#设定X轴 scale_y_continuous(expand = c(0, 0) ) +#去除绘图区和X轴之间的gap geom_hline(yintercept = c(6, -log10(0.05/nrow(df_SNP_position))), #添加阈值线 color = c('green', 'red'),size = 1.2, linetype = c("dotted", "twodash")) + geom_point(data=subset(data, is_highlight=="yes"), color="green", size=2)+facet_zoom(x = BPcum >= 3000 & BPcum <=3500)+#展示某一区域的p值情况 theme_prism(palette = "flames",#使用ggprism包自带主题 base_fontface = "plain", base_family = "serif", base_size = 16, base_line_size = 0.8, axis_text_angle = 45)+ theme(legend.position = "none",#去除图例 panel.grid = element_blank(), panel.border = element_blank(), axis.line.x = element_line(), axis.line.y = element_line())

参考：

1）https://blog.csdn.net/ddxygq/article/details/103555955

2）https://mp.weixin.qq.com/s/IO5rqz8sXaYD5vc5EbkkPw

我就知道你“在看”