跟著Nature Communications學(xué)作圖：R語言ggplot2散點組合誤差線展示響應(yīng)比（Response ratio）

論文

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í)素材

這篇論文是公眾號的一位讀者留言，說這篇論文提供了數(shù)據(jù)和代碼，但是代碼比較長，看起來比較吃力。我看了論文中提供的代碼，大體上能夠看懂，爭取抽時間把論文中提供的代碼都復(fù)現(xiàn)一下。因為論文中的圖都對應(yīng)著提供了作圖數(shù)據(jù)，我們想復(fù)現(xiàn)論文中的圖。關(guān)于用原始數(shù)據(jù)分析的部分后續(xù)有時間在單獨介紹。

論文中提供的代碼鏈接

https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-020-16881-7/MediaObjects/41467_2020_16881_MOESM8_ESM.txt

今天的推文我們復(fù)現(xiàn)論文中的figure1

image.png

論文中提供的作圖數(shù)據(jù)如下，excel存儲

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)

首先是第一個小圖a

論文中的代碼是用RR作為Y軸，GCFs作為X軸，然后再通過coord_flip()函數(shù)整體旋轉(zhuǎn)；論文中關(guān)于字體上小標是用expression函數(shù)實現(xiàn)的，這里我們使用latex2exp這個R包

代碼我們參考論文中的代碼，但是不完全按照他的寫

數(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這個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這個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ù)可以給公眾號推文點贊，點擊在看，最后留言獲取

查rma()函數(shù)找到了這個鏈接

http://www.simonqueenborough.info/R/specialist/meta-analysis#:~:text=The%20function%20rma()%20(random,compute%20effect%20sizes%20before%20modelling.&text=Random%20effect%20model%20can%20be,%2D%2D%2DFixed%20effect%20model%20cannot.

http://www.simonqueenborough.info/R/intro/index.html

歡迎大家關(guān)注我的公眾號

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)、群體遺傳學(xué)文獻閱讀筆記；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記！

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容