论文
Pangenome analysis reveals genomic variations associated with domestication traits in broomcorn millet
https://doi.org/10.1038/s41588-023-01571-z
论文中提供大部分图的原始作图数据，我们可以试着用论文中提供的原始数据来复现一下论文中的图
今天的推文来复现一下论文中的figure1b 和figure1c
image.png
部分示例数据截图
image.png
image.png
读取数据
代码语言：javascript复制library(readxl)
library(tidyverse)
dat<-read_excel("data/20231201/41588_2023_1571_MOESM3_ESM.xlsx",
sheet = "Supplementary Table 1",
skip = 1)
dat
dat<-dat %>%
mutate(Accession=str_replace(`Alternative sample name`,"Clean_",""))
fig1b.dat<-read_excel("data/20231201/41588_2023_1571_MOESM4_ESM.xlsx",
sheet = "Fig. 1b",
skip = 1)
fig1b.dat
fig1c.dat<-read_excel("data/20231201/41588_2023_1571_MOESM4_ESM.xlsx",
sheet = "Fig. 1c",
skip = 1)
fig1c.dat
论文中给landrace还分为 Chinese landraces (landrace C) from the European and Central Asian landraces (landrace O)
生成这个分组
代码语言：javascript复制dat %>%
filter(str_count(Location,"[A-Z]")>=2) %>%
pull(Location) %>% table() %>%
as.data.frame() %>%
rename("province"=".") %>%
filter(!(province=="USA"|province=="Czech Republic"|province=="South Korea")) %>%
pull(province) %>%
as.character() -> province.X
这个代码的作用是把中国的省份摘出来，思路是中国省份的拼音是最少2个大写字母，然后再去除其他
figure1b作图代码
代码语言：javascript复制fig1b.dat %>%
left_join(dat,by=c("Accession"="Accession")) %>%
mutate(new.group=case_when(
Location %in% province.X ~ "C",
TRUE ~ "O"
)) %>%
mutate(new.group02=case_when(
Group == "Landrace" ~ paste(Group,new.group,sep=" "),
TRUE ~ Group
)) %>%
ggplot(aes(x=PC1,y=PC2))
geom_point(aes(color=new.group02),
size=2)
scale_color_manual(values = c("Cultivar"="#e86a10",
"Landrace C"="#3a78c4",
"Landrace O"="#56a4aa",
"Wild"="#f1ab00"),
name=NULL)
theme_bw(base_size = 15)
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line(),
legend.position = c(0.2,0.2))
labs(x="PC1 (27.9%)",y="PC2 (16.9%)")
guides(color=guide_legend(override.aes = list(size=5)))
image.png
figure1c的作图代码，思路和1b基本一致
代码语言：javascript复制fig1c.dat %>%
left_join(dat,by=c("Accession"="Accession")) %>%
mutate(new.group=case_when(
Location %in% province.X ~ "C",
TRUE ~ "O"
)) %>%
mutate(new.group02=case_when(
Group == "Landrace" ~ paste(Group,new.group,sep=" "),
TRUE ~ Group
)) %>%
ggplot(aes(x=PC1,y=PC3))
geom_point(aes(color=new.group02),
size=2)
scale_color_manual(values = c("Cultivar"="#e86a10",
"Landrace C"="#3a78c4",
"Landrace O"="#56a4aa",
"Wild"="#f1ab00"),
name=NULL)
theme_bw(base_size = 15)
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line(),
legend.position = c(0.2,0.2))
labs(x="PC1 (27.9%)",y="PC2 (12.8%)")
guides(color=guide_legend(override.aes = list(size=5)))
image.png
组合图
代码语言：javascript复制p1b<-fig1b.dat %>%
left_join(dat,by=c("Accession"="Accession")) %>%
mutate(new.group=case_when(
Location %in% province.X ~ "C",
TRUE ~ "O"
)) %>%
mutate(new.group02=case_when(
Group == "Landrace" ~ paste(Group,new.group,sep=" "),
TRUE ~ Group
)) %>%
ggplot(aes(x=PC1,y=PC2))
geom_point(aes(color=new.group02),
size=2)
scale_color_manual(values = c("Cultivar"="#e86a10",
"Landrace C"="#3a78c4",
"Landrace O"="#56a4aa",
"Wild"="#f1ab00"),
name=NULL)
theme_bw(base_size = 15)
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line(),
legend.position = c(0.2,0.2))
labs(x="PC1 (27.9%)",y="PC2 (16.9%)")
guides(color=guide_legend(override.aes = list(size=5)))
p1b
p1c<-fig1c.dat %>%
left_join(dat,by=c("Accession"="Accession")) %>%
mutate(new.group=case_when(
Location %in% province.X ~ "C",
TRUE ~ "O"
)) %>%
mutate(new.group02=case_when(
Group == "Landrace" ~ paste(Group,new.group,sep=" "),
TRUE ~ Group
)) %>%
ggplot(aes(x=PC1,y=PC3))
geom_point(aes(color=new.group02),
size=2)
scale_color_manual(values = c("Cultivar"="#e86a10",
"Landrace C"="#3a78c4",
"Landrace O"="#56a4aa",
"Wild"="#f1ab00"),
name=NULL)
theme_bw(base_size = 15)
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line(),
legend.position = c(0.2,0.2))
labs(x="PC1 (27.9%)",y="PC2 (12.8%)")
guides(color=guide_legend(override.aes = list(size=5)))
p1c
library(patchwork)
p1c p1b
image.png
示例数据可以到论文中下载，或者给推文打赏1元获取我整理的示例数据和代码
