論文
Meta-analysis of the impacts of global change factors on soil microbial diversity and functionality
https://www.nature.com/articles/s41467-020-16881-7#Sec15
論文里提供了數(shù)據(jù)和代碼栓始，很好的學(xué)習(xí)素材
這篇論文是公眾號(hào)的一位讀者留言岔擂，說(shuō)這篇論文提供了數(shù)據(jù)和代碼，但是代碼比較長(zhǎng)裁替，看起來(lái)比較吃力阀溶。我看了論文中提供的代碼抽莱，大體上能夠看懂葵陵，爭(zhēng)取抽時(shí)間把論文中提供的代碼都復(fù)現(xiàn)一下。因?yàn)檎撐闹械膱D都對(duì)應(yīng)著提供了作圖數(shù)據(jù)蘸泻，我們想復(fù)現(xiàn)論文中的圖琉苇。關(guān)于用原始數(shù)據(jù)分析的部分后續(xù)有時(shí)間在單獨(dú)介紹。
論文中提供的代碼鏈接
今天的推文我們復(fù)現(xiàn)論文中的figure1
image.png
論文中提供的作圖數(shù)據(jù)如下悦施，excel存儲(chǔ)
image.png
加載需要用到的R包
library(readxl)
library(tidyverse)
library(latex2exp)
library(ggplot2)
讀取數(shù)據(jù)
metaresult<-read_excel("data/20221129/41467_2020_16881_MOESM9_ESM.xlsx",
sheet = 'Fig1')
colnames(metaresult)
首先是第一個(gè)小圖a
論文中的代碼是用RR作為Y軸并扇，GCFs作為X軸，然后再通過(guò)coord_flip()函數(shù)整體旋轉(zhuǎn)抡诞；論文中關(guān)于字體上小標(biāo)是用expression函數(shù)實(shí)現(xiàn)的穷蛹，這里我們使用latex2exp這個(gè)R包
代碼我們參考論文中的代碼，但是不完全按照他的寫(xiě)
數(shù)據(jù)整理和作圖代碼
data1<-metaresult %>%
filter(Variables=="Richness"|Variables=="Shannon")
data1$GCFs
data1<-data1 %>%
mutate(GCFs=factor(GCFs,
levels = c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"))
)
data1 %>% colnames()
ggplot(data = data1,
aes(x=`Weighted means of RR`,
y=`GCFs`,
shape=Variables))+
geom_vline(xintercept=0,linetype = "dashed",size=0.2)+
geom_errorbarh(aes(xmin=`Lower confidence intervals`,
xmax=`Upper confidence intervals`),
position=position_dodge(0.8),
height=0.2)+
geom_point(position=position_dodge(0.8),
size=3, stroke = 0.3,
aes(fill=GCFs),
show.legend = FALSE)+
geom_text(aes(y =`GCFs` , x = `Upper confidence intervals`+0.015,
label = `Sample sizes`),
position = position_dodge(width = 0.8),
vjust = 0.4, hjust=0.4, size = 4,
check_overlap = FALSE)+
geom_segment(y = 11.6, x = -Inf,
yend = 11.6, xend = Inf,
colour = "black",size=0.4)+
scale_shape_manual(values=c("Richness"=21,"Shannon"=22))+
scale_x_continuous(limits=c(-0.17,0.17),
breaks = c(-0.16,-0.08,0,0.08,0.16))+
scale_y_discrete(breaks=c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"),
labels=c(TeX(r"($N \times P \times K$)"),
TeX(r"($N \times P$)"),
TeX(r"($N \times PPT$+)"),
TeX(r"($W \times eCO_2$)"),
"LUC","N","P","PPT+","PPT-",
TeX(r"($eCO_2$)"),
"W"))+
labs(x = "Global change factors ", y = "RR of alpha diversity",colour = 'black')+
theme(legend.title = element_blank(),
legend.position=c(0.2,0.94),
legend.key = element_rect(fill = "white",size = 2),
legend.key.width = unit(0.5,"lines"),
legend.key.height= unit(0.8,"lines"),
legend.background = element_blank(),
legend.text=element_text(size=6),
panel.background = element_rect(fill = 'white', colour = 'white'),
axis.title=element_text(size=9),
axis.text.y = element_text(colour = 'black', size = 8),
axis.text.x = element_text(colour = 'black', size = 8),
axis.line = element_line(colour = 'black',size=0.4),
axis.line.y = element_blank(),
axis.ticks = element_line(colour = 'black',size=0.4),
axis.ticks.y = element_blank())
輸出結(jié)果
image.png
小圖b
data2<-metaresult %>%
filter(Variables=="Beta Diversity")
data2$GCFs
data2<-data2 %>%
mutate(GCFs=factor(GCFs,
levels = c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"))
)
data2 %>% colnames()
ggplot(data = data2,
aes(x=`Weighted means of RR`,
y=`GCFs`))+
geom_vline(xintercept=0,linetype = "dashed",size=0.2)+
geom_errorbarh(aes(xmin=`Lower confidence intervals`,
xmax=`Upper confidence intervals`),
height=0.2)+
geom_point(size=3, stroke = 0.3,
shape=21,
aes(fill=GCFs),
show.legend = FALSE)+
geom_text(aes(y =`GCFs` , x = `Upper confidence intervals`+0.1,
label = `Sample sizes`),
#position = position_dodge(width = 0.8),
vjust = 0.4, hjust=0.4, size = 4,
check_overlap = FALSE)+
geom_segment(y = 11.6, x = -Inf,
yend = 11.6, xend = Inf,
colour = "black",size=0.4)+
scale_x_continuous(limits=c(-0.6,1.1),breaks = c(-0.5,0,0.5,1))+
scale_y_discrete(breaks=c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"),
labels=c(TeX(r"($N \times P \times K$)"),
TeX(r"($N \times P$)"),
TeX(r"($N \times PPT$+)"),
TeX(r"($W \times eCO_2$)"),
"LUC","N","P","PPT+","PPT-",
TeX(r"($eCO_2$)"),
"W"))+
labs(y = "Global change factors ",
x = "RR of alpha diversity",
colour = 'black')+
theme(legend.title = element_blank(),
legend.position=c(0.2,0.9),
legend.key = element_rect(fill = "white",size = 2),
legend.key.width = unit(0.5,"lines"),
legend.key.height= unit(0.8,"lines"),
legend.background = element_blank(),
legend.text=element_text(size=6),
panel.background = element_rect(fill = 'white', colour = 'white'),
axis.title=element_text(size=9),
axis.text.y = element_text(colour = 'black', size = 8),
axis.text.x = element_text(colour = 'black', size = 8),
axis.line = element_line(colour = 'black',size=0.4),
axis.line.y = element_blank(),
axis.ticks = element_line(colour = 'black',size=0.4),
axis.ticks.y = element_blank())
image.png
小圖c
data3<-metaresult %>%
filter(Variables=="Community structure")
data3$GCFs
data3<-data3 %>%
mutate(GCFs=factor(GCFs,
levels = c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"))
)
data3 %>% colnames()
ggplot(data = data3,
aes(x=`Weighted means of RR`,
y=`GCFs`))+
geom_vline(xintercept=0,linetype = "dashed",size=0.2)+
geom_errorbarh(aes(xmin=`Lower confidence intervals`,
xmax=`Upper confidence intervals`),
height=0.2)+
geom_point(size=3, stroke = 0.3,
shape=21,
aes(fill=GCFs),
show.legend = FALSE)+
geom_text(aes(y =`GCFs` , x = `Upper confidence intervals`+0.1,
label = `Sample sizes`),
#position = position_dodge(width = 0.8),
vjust = 0.4, hjust=0.4, size = 4,
check_overlap = FALSE)+
geom_segment(y = 11.6, x = -Inf,
yend = 11.6, xend = Inf,
colour = "black",size=0.4)+
scale_x_continuous(limits=c(-0.6,2.0),breaks = c(-0.5,0,0.5,1,1.5,2.0))+
scale_y_discrete(breaks=c("N_P_K","N_P","N_PPT+",
"W_eCO2","LUC","N","P",
"PPT+","PPT-","eCO2","W"),
labels=c(TeX(r"($N \times P \times K$)"),
TeX(r"($N \times P$)"),
TeX(r"($N \times PPT$+)"),
TeX(r"($W \times eCO_2$)"),
"LUC","N","P","PPT+","PPT-",
TeX(r"($eCO_2$)"),
"W"))+
labs(y = "Global change factors ",
x = "RR of community structure",
colour = 'black')+
theme(legend.title = element_blank(),
legend.position=c(0.2,0.9),
legend.key = element_rect(fill = "white",size = 2),
legend.key.width = unit(0.5,"lines"),
legend.key.height= unit(0.8,"lines"),
legend.background = element_blank(),
legend.text=element_text(size=6),
panel.background = element_rect(fill = 'white', colour = 'white'),
axis.title=element_text(size=9),
axis.text.y = element_text(colour = 'black', size = 8),
axis.text.x = element_text(colour = 'black', size = 8),
axis.line = element_line(colour = 'black',size=0.4),
axis.line.y = element_blank(),
axis.ticks = element_line(colour = 'black',size=0.4),
axis.ticks.y = element_blank())
圖b和圖c是一樣的
最后是拼圖
論文中提供的拼圖代碼是用ggpubr這個(gè)R包做的
ggpubr::ggarrange(p1, p2, p3,
widths = c(7, 5, 5),
ncol = 3, nrow = 1,
labels = c("a", "b", "c"),
font.label=list(size=12),
hjust = 0, vjust = 1)
我自己更習(xí)慣使用patchwork這個(gè)R包
library(patchwork)
p1+
p2+theme(axis.text.y = element_blank(),
axis.title.y = element_blank())+
p3+theme(axis.text.y = element_blank(),
axis.title.y = element_blank())+
plot_annotation(tag_levels = "a")+
plot_layout(widths = c(7, 5, 5))
最終結(jié)果
image.png
示例數(shù)據(jù)和代碼可以自己到論文中下載昼汗，如果需要我推文中的代碼和數(shù)據(jù)可以給公眾號(hào)推文點(diǎn)贊肴熏，點(diǎn)擊在看，最后留言獲取
查rma()函數(shù)找到了這個(gè)鏈接
http://www.simonqueenborough.info/R/intro/index.html
歡迎大家關(guān)注我的公眾號(hào)
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號(hào) 主要分享：1顷窒、R語(yǔ)言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡(jiǎn)單小例子蛙吏；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)鞋吉、基因組學(xué)出刷、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3坯辩、生物信息學(xué)入門(mén)學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！
