ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
103 stars 30 forks source link

跟着Nature Genetics学作图:R语言ggplot2散点图展示主成分分析(PCA)的结果 #4179

Closed ixxmu closed 9 months ago

ixxmu commented 9 months ago

https://mp.weixin.qq.com/s/jHVlmMZHAn_kWLHcO1oknQ

ixxmu commented 9 months ago

跟着Nature Genetics学作图:R语言ggplot2散点图展示主成分分析(PCA)的结果 by 小明的数据分析笔记本

论文

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

image.png

示例数据可以到论文中下载,或者给推文打赏1元获取我整理的示例数据和代码

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享:1、R语言和python做数据分析和数据可视化的简单小例子;2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记;3、生物信息学入门学习资料及自己的学习笔记!