R语言~组间差异分析并标注abcd

加载工具包

library(dplyr)
library(gghalves)
library(agricolae)

导入数据并整理

dat = read.table('..GlobalAtlasv2_scores_v5_Manu.txt', header = T, sep = '\t',check.names = FALSE)
dat = subset(dat, Eco_corrected !='Moss_heath')
names(dat) = gsub('_wc2','',names(dat))
names(dat) = gsub('_c$','',names(dat))
names(dat) = gsub('_18s_richness','', names(dat))
names(dat) = gsub('Richness_bacteria','Bacteria', names(dat))
names(dat) = gsub('WaterHoldingCapacity','WHC',names(dat))
names(dat) = gsub('Proportion_decomposers','Saprobe',names(dat))
names(dat) = gsub('Proportion_symbiotic_fungi','Mycorrhizal', names(dat))
names(dat) = gsub('Pathogen_control','Pathogenctrl',names(dat))
names(dat) = gsub('Soil_respiration','Respiration',names(dat))
names(dat) = gsub('Soil_salinity','Salinity',names(dat))
names(dat) = gsub('Soil_pH','pH', names(dat))
names(dat) = gsub('Clay_silt','Clay',names(dat))
names(dat) = gsub('Soil_ORC_g_kg','SOC', names(dat))
names(dat) = gsub('Plant_cover','Pcov', names(dat))
names(dat) = gsub('Aridity_Index_v3','AI', names(dat))
names(dat) = gsub('Eco_corrected','vegetation', names(dat))
dat$Forest_corrected = NULL
dat = dat %>% filter(AI<2) %>%
  mutate(Order_ID = paste0('S', Order_ID),
         Latitude = abs(Latitude),
         Longitude = sin(Longitude),
         NPP = NPP/10000,
         Pathogenctrl = -(log(abs(Pathogenctrl)+1)),
         Mycorrhizal = log(Mycorrhizal+1),
         Saprobe = log(Saprobe +1),
         NO3 = log(NO3+1),
         NH4 = log(NH4+1),
         PO4 = log(PO4+1),
         SOC = log(SOC+1),
         NAG = log(NAG+1),
         # pH = abs(pH-7),
         Respiration = log(Respiration+1),
         Salinity = log(Salinity+1))

with(dat, table(vegetation))

简历变量组合对象，方便后续调用

env_var = c("Latitude", "MAT", 'TSEA', "PSEA", 'AI',  "Pcov", "pH", 'Salinity', 'Clay')

查看数据分布的正态性

qq = dat[c(div_var, env_var)] %>% scale() %>% data.frame() %>% melt() %>% 
  ggplot(aes(sample = value, color = variable)) +
  facet_wrap(.~variable,scales = 'free') +
  geom_abline()+
  stat_qq(show.legend = F) 

数据整理和变量重命名

env_dat = dat[c(env_var,'vegetation')] %>% reshape2::melt()
env_dat$vegetation = factor(env_dat$vegetation, levels = c("Forest",'Grassland','Shrubland'),labels = c("F",'G','S'))

对多个变量分别进行组间差异分析

env_hsd = lapply(split(env_dat, f = env_dat$variable),
                 function(x) {aov(value ~ vegetation, data = x)%>% HSD.test(trt = 'vegetation',group = T)}
) %>%
  lapply('[[', 'groups') %>%
  do.call('rbind',.) %>%
  rownames_to_column() %>%
  rename(letter = groups)

整理结果变量名称，并规定显著性标记的坐标位置

env_hsd$variable = str_split_fixed(env_hsd$rowname,'[.]',2)[,1]
env_hsd$vegetation = str_split_fixed(env_hsd$rowname,'[.]',2)[,2]
env_hsd = env_dat %>% 
  group_by(variable) %>% 
  mutate(yaxis = max(value)) %>% 
  select(variable, yaxis) %>% 
  unique.data.frame() %>% merge(env_hsd, by = 'variable')

设定变量的出图顺序，绘图

env_dat$variable = factor(env_dat$variable, 
                          levels = c('Latitude','MAT','AI','PSEA','TSEA','Pcov','Clay','pH','Salinity'),
                          labels = c(expression(paste("Latitude (abs)")),
                                     expression(paste('MAT (℃)')),
                                     expression(paste('Aridity index')),
                                     expression(paste('TSEA (SD x 100)')),
                                     expression(paste('PSEA (CV)')),
                                     expression(paste('Plant coverage (%)')),
                                     expression(paste('Clay + silt (%)')),
                                     expression(paste('Soil pH')), 
                                     expression(paste('Soil salinity (log, μS cm'^-1,')'))))

env_hsd$variable = factor(env_hsd$variable,
                          levels = c('Latitude','MAT','AI','PSEA','TSEA','Pcov','Clay','pH','Salinity'),
                          labels = c(expression(paste("Latitude (abs)")),
                                     expression(paste('MAT (℃)')),
                                     expression(paste('Aridity index')),
                                     expression(paste('TSEA (SD x 100)')),
                                     expression(paste('PSEA (CV)')),
                                     expression(paste('Plant coverage (%)')),
                                     expression(paste('Clay + silt (%)')),
                                     expression(paste('Soil pH')),
                                     expression(paste('Soil salinity (log, μS cm'^-1,')'))))

env_plot = ggplot(env_dat, aes(vegetation, value))+
  geom_half_violin(aes(fill = vegetation, color = vegetation),
                   position=position_nudge(x=0.1,y=0),
                   side='R',adjust=1.2,trim=F,color=NA,alpha=0.8,show.legend = F)+
  # geom_boxplot(aes(fill = vegetation), width = 0.1, alpha = 0.7,outlier.alpha = 0,show.legend = F)+
  geom_point(aes(x = as.numeric(vegetation)-0.1,y = value, color = vegetation),
             position = position_jitter(width =0.1),size =0.2, shape = 20, show.legend = F)+
  geom_text(data = env_hsd, aes(x=vegetation, y = yaxis, label =letter), 
            size = 4,position = position_nudge(x = 0.3))+
  facet_wrap(.~variable, scales = 'free', strip.position = 'left',labeller = label_parsed)+
  theme_classic()+
  stat_compare_means(label = 'p.signif', label.x.npc = 0.5, label.y.npc = 1,  
                     method = 'anova', color = 'black',
                     symnum.args = list(cutpoints = c(0,0.001,0.01,0.05,Inf), 
                                        symbols = c("***","**","*",'ns')))+
  labs(x = NULL, y = NULL)+
  theme(strip.placement = 'outside',
        strip.background = element_blank(),
        strip.text = element_text(color = 'black'),
        axis.text = element_text(color = 'black'),
        panel.spacing.y = unit(1,'lines'))
ggsave(env_plot,filename = 'env.png', width = 8, height = 6, dpi = 300, bg = 'white')