论文

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

image.png

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

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

代码链接

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

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

image.png

这幅图展示的是城乡环境之间的差异，做的是主成分分析，首先看下他的数据集的样子

inpath<-"paper_dat/glue_pc/phenotypic-analyses/data/clean/popMeans_allCities_withEnviro/"
csv.files<-list.files(path = inpath,pattern = "*.csv")
csv.files

总共有160个csv文件，其中一个部分截图如下

image.png

原论文提供的数据并没有这个数据集，需要通过运行论文的一系列代码才能得到这个数据，这个论文的代码量是真的大，非常值得好好学习

首先将github链接上所有文件都下载下来，然后将Rstudio的工作目录设置到

getwd()
setwd("phenotypic-analyses/")

目录下

然后打开main.R这个代码文件
一直运行到如下行代码

source("scripts/r/data-processing/popMeans_addEnviroData.R")

运行这行代码的时候会遇到报错 （这里报错界面忘记截图了）

我的处理方式是将 utilityFunctions.R这个文件里的涉及到noquote()函数都删掉，运行完就得到了作图需要用到的数据集。

这里的noquote函数具体起到什么作用我暂时还没太理解

接下来是运行这行代码

source("scripts/r/analyses/enviroAnalyses.R")

这行代码如果运行成功会得到pca分析的结果，但是我运行的时候遇到报错

Error in rlm.default(x, y, weights, method = method, wt.method = wt.method,  : 
  'x' is singular: singular fits are not implemented in 'rlm'
Called from: rlm.default(x, y, weights, method = method, wt.method = wt.method, 
    ...)

暂时没看懂这个报错是什么意思

最后通过反复看 enviroAnalyses.R这个文件里的代码，通过去掉三个城市的数据集的方式得到了pca分析的输出数据

以上是运行论文代码排查错误的过程，以下是做pca分析并且作图的代码

现在是两个输入文件

• sciencefig2A.csv
• sciencefig2A_group_info.csv

首先是读取数据集

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)

主成分分析

这个论文里选择的是vegan这个包里的rda函数

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

image.png

计算PC1 PC2的解释率

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

提取PC1 PC2的作图数据

library(tidyverse)

enviroPCA_sites  <- scores(enviroPCA, display = 'sites', choices = c(1, 2), scaling = 1) %>%
  as.data.frame() %>% 
  dplyr::select(PC1, PC2) %>% 
  mutate(habitat = dat02$habitat,
         city = dat02$cities)

准备作图配色

library(wesanderson)
pal <- wes_palette('Darjeeling1', 5, type = 'discrete')
urban_col <- pal[4]
rural_col <- pal[2]
cols <- c(urban_col, rural_col)
cols

ggplot2作图主题设置

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_plot <- ggplot(enviroPCA_sites, aes(x = PC1, y = PC2)) + 
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_line(aes(group = city), alpha = 0.7) +
  geom_point(size = 2.75, shape = 21, colour = "black", aes(fill =  habitat)) +
  stat_ellipse(aes(colour = habitat), size = 1.5, 
               level = 0.95,
               show.legend = F) +
  xlab(sprintf("PC1 (%.1f%%)", PC1_varEx)) + 
  ylab(sprintf("PC2 (%.1f%%)", PC2_varEx)) +
  scale_colour_manual(values = rev(cols)) +
  scale_fill_manual(values = rev(cols)) +
  scale_x_continuous(breaks = seq(from = -0.6, to = 0.6, by = 0.2), 
                     labels = scales::comma) +
  scale_y_continuous(breaks = seq(from = -0.45, to = 0.45, by = 0.15), 
                     labels = scales::comma) +
  ng1 + theme(legend.position = "top", 
              legend.direction="horizontal",
              legend.text = element_text(size=15), 
              legend.key = element_rect(fill = "white"),
              legend.title = element_blank(),
              legend.key.size = unit(0.5, "cm"),
              legend.spacing.x = unit(0.1, "cm"))
enviroPCA_plot

image.png

今天推文的示例数据和代码可以在推文下留言20220422获取

欢迎大家关注我的公众号

小明的数据分析笔记本

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