跟着Nature Genetics学作图：R语言ggplot2散点图展示主成分分析（PCA）的结果

用户7010445

发布于 2023-12-01 14:50:54

2580

发布于 2023-12-01 14:50:54

文章被收录于专栏：小明的数据分析笔记本小明的数据分析笔记本

论文

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

读取数据

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)

生成这个分组

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作图代码

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基本一致

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

组合图

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