跟着Science学作图：R语言ggplot2画箭头展示变量对主成分的贡献

用户7010445

发布于 2022-05-23 07:38:31

80400

代码可运行

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

运行总次数：0

代码可运行

论文

https://www.science.org/doi/10.1126/science.abk0989

image.png

最近朋友圈好多人都在转这个论文，我也找来看了看，论文研究的内容看的还是一知半解。

论文用到的数据代码都是公开的，我们可以学习一下其中的代码

代码链接

https://github.com/James-S-Santangelo/glue_pc

今天的图文重复论文中的Figure 2B

image.png

这个图的图注写的是The eigenvectors for environmental variables, colored according to their contribution to PC2

这里为什么只展示对PC2的贡献暂时还不明白。主要是论文的研究内容看不明白

本篇推文只记录画图代码了

还是先做主成分分析

library(readr)
dat01<-read_csv("phenotypic-analyses/sciencefig2A.csv")
dim(dat01)
colnames(dat01)
dat02<-read_csv("phenotypic-analyses/sciencefig2A_group_info.csv")
dim(dat02)
colnames(dat02)

library(vegan)
enviroPCA <- rda(dat01, 
                 scale = TRUE, na.action = "na.omit")

eig <- enviroPCA$CA$eig
percent_var <- eig * 100 / sum(eig)
PC1_varEx <- round(percent_var[1], 1)  # Percent variance explained by PC1
PC2_varEx <- round(percent_var[2], 1)  # Percent variance explained by PC2

计算物种贡献百分比

论文里提供的代码里放了一个参考链接 https://stackoverflow.com/questions/50177409/how-to-calculate-species-contribution-percentages-for-vegan-rda-cca-objects

contrib <- round(100*scores(enviroPCA, display = "sp", scaling = 0)[,2]^2, 3)

生成作图数据

library(tidyverse)
enviroPCA_vars  <- scores(enviroPCA, display = 'species', choices = c(1, 2), scaling = 2) %>% 
  as.data.frame() %>% 
  rownames_to_column(., var = 'var') %>% 
  mutate(contrib = contrib)

准备配色

library(wesanderson)
pal <- wes_palette("Darjeeling1", 3, type = "continuous")

作图主题的一些设置

library(ggplot2)
ng1 <- theme(aspect.ratio=0.7,panel.background = element_blank(),
             panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             panel.border=element_blank(),
             axis.line.x = element_line(color="black",size=1),
             axis.line.y = element_line(color="black",size=1),
             axis.ticks=element_line(size = 1, color="black"),
             axis.ticks.length=unit(0.25, 'cm'),
             axis.text=element_text(color="black",size=15),
             axis.title=element_text(color="black",size=1),
             axis.title.y=element_text(vjust=2,size=17),
             axis.title.x=element_text(vjust=0.1,size=17),
             axis.text.x=element_text(size=15),
             axis.text.y=element_text(size=15),
             strip.text.x = element_text(size = 10, colour = "black",face = "bold"),
             strip.background = element_rect(colour="black"),
             legend.position = "top", legend.direction="vertical",
             legend.text=element_text(size=17), legend.key = element_rect(fill = "white"),
             legend.title = element_text(size=17),legend.key.size = unit(1.0, "cm"))

最后的作图代码

enviroPCA_variableContrib <- ggplot() +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(data = enviroPCA_vars, aes(x = 0, xend = PC1, y=0, yend = PC2, color = contrib), 
               size = 2, arrow = arrow(length = unit(0.02, "npc")), alpha = 1) +
  geom_text(data = enviroPCA_vars,
            aes(x = PC1, y = PC2, label = var,
                hjust = "inward", vjust =  0.5 * (1 - sign(PC1))),
            color = "black", size = 3.5) + 
  xlab(sprintf("PC1 (%.1f%%)", PC1_varEx)) + ylab(sprintf("PC2 (%.1f%%)", PC2_varEx)) +
  scale_colour_gradientn(colours = rev(pal), breaks = seq(from = 5, to = 25, by = 5)) +
  # scale_x_continuous(breaks = seq(from = -1, to = 1, by = 0.25)) +
  # scale_y_continuous(breaks = seq(from = -1, to = 1, by = 0.25)) +
  ng1 + theme(legend.position = "top",
              legend.direction="horizontal",
              # legend.title = element_blank(),
              legend.key.size = unit(0.5, "cm"),
              legend.spacing.x = unit(0.1, "cm"),
              legend.text = element_text(size=10)) +
  guides(color = guide_colourbar(barwidth = 10, barheight = 0.5))
enviroPCA_variableContrib