论文
A saturated map of common genetic variants associated with human height
https://www.nature.com/articles/s41586-022-05275-y
s41586-022-05275-y.pdf
代码没有公开，但是作图数据基本都公开了，争取把每个图都重复一遍
今天的推文重复论文中的Figure1
代码
setwd("data/20221014")
library(readxl)
fig1<-read_excel("Figure1.xlsx")
colnames(fig1)
library(tidyverse)
library(stringr)
str_replace_all("[0,5e-100]","\\(|5e-|\\]|\\[","") %>%
str_split_fixed(",",n=2) %>%
as.data.frame() %>%
pull(V1) %>% as.numeric()
str_replace_all("[0,5e-100]","\\(|5e-|\\]|\\[","") %>%
str_split_fixed(",",n=2) %>%
as.data.frame() %>%
pull(V2) %>% as.numeric()
fig1 %>%
mutate(max_value=str_replace_all(`P-value Caregory`,"\\(|5e-|\\]|\\[","") %>%
str_split_fixed(",",n=2) %>%
as.data.frame() %>%
pull(V1) %>% as.numeric(),
min_value=str_replace_all(`P-value Caregory`,"\\(|5e-|\\]|\\[","") %>%
str_split_fixed(",",n=2) %>%
as.data.frame() %>%
pull(V2) %>% as.numeric()) %>%
mutate(group=case_when(
min_value == 100 & max_value == 0 ~ "group01",
min_value == 50 & max_value == 100 ~ "group02",
min_value == 20 & max_value == 50 ~ "group03",
min_value == 10 & max_value == 20 ~ "group04",
min_value == 8 & max_value == 10 ~ "group05",
)) -> new.fig1
table(new.fig1$group)
library(ggplot2)
library(ggh4x)
library(cowplot)
ggplot(data=new.fig1,
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`,
color=group))+
geom_point( key_glyph = rectangle_key_glyph(color=color,
fill=color,
padding = margin(3, 3, 3, 3)))+
scale_color_manual(values = c("group01"="#ee82ee",
"group02"="#2e8b57",
"group03"="#1e90ff",
"group04"="#daa520",
"group05"="#cdc673"),
name="",
labels=c("group01"="P < 5 × 10–100 (672 SNPs)",
"group02"="5 × 10–50 > P > 5 × 10–100 (1,110 SNPs)",
"group03"="5 × 10–20 > P > 5 × 10–50 (3,513 SNPs)",
"group04"="5 × 10–10 > P > 5 × 10–20 (5,192 SNPs)",
"group05"="5 × 10–8 > P > 5 × 10–10 (1,624 SNPs)"))+
theme_bw()+
theme(panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line(),
legend.position = c(0.7,0.8))+
scale_x_continuous(breaks = c(0.01,0.05,0.1,0.2,0.3,0.4,0.5),
labels = c(1,5,10,20,30,40,50))+
scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1,0.2,0.3),
limits = c(-0.3,0.3))+
guides(x=guide_axis_truncated(trunc_lower = 0.01,
trunc_upper = 0.5),
y=guide_axis_truncated(trunc_lower = -0.3,
trunc_upper = 0.3))+
labs(x="MAF (%) in cross-ancestry meta-analysis",
y="Joint effect sizes (s.d.) of minor alleles\nin cross-ancestry meta-analysis")+
geom_hline(yintercept = 0,color="gray")+
geom_smooth(data = new.fig1 %>%
filter(group=="group01") %>%
filter(`Join Effect of Minor Allele`<0),
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`),
method = 'loess',
formula = 'y~x',
se=FALSE,color="gray",
show.legend = FALSE)+
geom_smooth(data = new.fig1 %>%
filter(group=="group01") %>%
filter(`Join Effect of Minor Allele`>0),
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`),
method = 'loess',
formula = 'y~x',
se=FALSE,color="gray",
show.legend = FALSE)
关于曲线不太清楚是用什么数据做的，这里直接自动添加拟合曲线
图例里的文本上下标 出图后再编辑吧
关于散点图今天还新学到一个知识点是：散点图的点如果非常多，如果输出pdf文件的话，pdf文件会非常大，比如GWAS里常用的曼哈顿图，这个pdf文件如果非常大后续如果想要编辑这个pdf文件会比较麻烦。
关于如何解决这个问题，看到一个讨论群里有人讨论，他们提到一个办法是可以把散点栅格化 （栅格化是什么意思暂时不太明白）可以借助R包ggrastr
对应的github主页是
https://github.com/VPetukhov/ggrastr
正好我们今天的推文内容是数据量比较多的散点图，我们可以按照这个做法试试，这里参考微信公众号推文 https://mp.weixin.qq.com/s/ou0cjD8dLMNaDLk588KSwQ
安装ggrastr这个R包
install.packages('ggrastr')如果要把点栅格化，只需要把对应的散点图函数geom_point()换成geom_point_rast()
library(ggrastr)
p2<-ggplot(data=new.fig1,
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`,
color=group))+
geom_point_rast( key_glyph = rectangle_key_glyph(color=color,
fill=color,
padding = margin(3, 3, 3, 3)),
size=0.1,
raster.dpi = getOption("ggrastr.default.dpi", 300))+
scale_color_manual(values = c("group01"="#ee82ee",
"group02"="#2e8b57",
"group03"="#1e90ff",
"group04"="#daa520",
"group05"="#cdc673"),
name="",
labels=c("group01"="P < 5 × 10–100 (672 SNPs)$)",
"group02"="5 × 10–50 > P > 5 × 10–100 (1,110 SNPs)",
"group03"="5 × 10–20 > P > 5 × 10–50 (3,513 SNPs)",
"group04"="5 × 10–10 > P > 5 × 10–20 (5,192 SNPs)",
"group05"="5 × 10–8 > P > 5 × 10–10 (1,624 SNPs)"))+
theme_bw()+
theme(panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line(),
legend.position = c(0.7,0.8))+
scale_x_continuous(breaks = c(0.01,0.05,0.1,0.2,0.3,0.4,0.5),
labels = c(1,5,10,20,30,40,50))+
scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1,0.2,0.3),
limits = c(-0.3,0.3))+
guides(x=guide_axis_truncated(trunc_lower = 0.01,
trunc_upper = 0.5),
y=guide_axis_truncated(trunc_lower = -0.3,
trunc_upper = 0.3))+
labs(x="MAF (%) in cross-ancestry meta-analysis",
y="Joint effect sizes (s.d.) of minor alleles\nin cross-ancestry meta-analysis")+
geom_hline(yintercept = 0,color="gray")+
geom_smooth(data = new.fig1 %>%
filter(group=="group01") %>%
filter(`Join Effect of Minor Allele`<0),
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`),
method = 'loess',
formula = 'y~x',
se=FALSE,color="gray",
show.legend = FALSE)+
geom_smooth(data = new.fig1 %>%
filter(group=="group01") %>%
filter(`Join Effect of Minor Allele`>0),
aes(x=`Minor Allele Frequency`,
y=`Join Effect of Minor Allele`),
method = 'loess',
formula = 'y~x',
se=FALSE,color="gray",
show.legend = FALSE)
pdf("p1.pdf",width = 6,height = 6)
p1
dev.off()
pdf("p2.pdf",width = 6,height = 6)
p2
dev.off()输出的p2如果放大 点是会变模糊的
两个文件的大小也不一样，栅格化之前是700k，栅格化之后只有200k
示例数据和代码可以给公众号推文点赞，点击在看，最后留言获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
