论文
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、生物信息学入门学习资料及自己的学习笔记！
网友评论