加载工具包

library(phyloseq)
library(dplyr)
library(ggpubr)
library(ggridges)
library(scales)
library(RColorBrewer)
library(vegan)
library(reshape2)
library(rfPermute)
library(cowplot)

定义简单函数以合并距离矩阵

mtrx2cols = function(m1,m2,m3,m4,val1,val2,val3,val4){
  lt = lower.tri(m1)
  res = data.frame(row = row(m1,as.factor = T)[lt], 
                   col = col(m1,as.factor = T)[lt],
                   val1 = m1[lt], val2= m2[lt], val3 = m3[lt], val4 = m4[lt])
  names(res)[3:6] = c(val1,val2, val3,val4)
  return(res)
}

加载数据

load('../../Result/Bacteria/SILVA138/Filter_Outliers/00_load_data/.RData')
rm(result, otutab)
bac_map = map
bac_physeq = physeq 

计算距离矩阵

bac_bray = distance(bac_physeq, method = 'bray', type = 'sample')
bac_bray_df = as.data.frame(as.matrix(bac_bray))

bac_MAT_df = as.data.frame(as.matrix(dist(bac_map['MAT'],method = 'euclidean')))
bac_MAP_df = as.data.frame(as.matrix(dist(bac_map['MAP'],method = 'euclidean')))
bac_HII_df = as.data.frame(as.matrix(dist(bac_map['HII'],method = 'euclidean')))
names(bac_MAT_df) == names(bac_bray_df)
bac_map$Exposed = ifelse(bac_map$State == 'Exposed', 1, 0)
bac_map = subset(bac_map, 
                 select = -c(BioProject, Center.Name, geo_loc_country, geo_loc_continent, 
                             geo_loc_name, Instrument, Longitude, Sample_Name, Organism, 
                             Primers, FP_sequence, RP_sequence,Site, Source))

bac_ddr = mtrx2cols(bac_bray_df,bac_MAP_df,bac_MAT_df,bac_HII_df,'Bray','MAP','MAT','HII')
head(bac_ddr)
bac_ddr = melt(bac_ddr, id.vars = c('row','col','Bray'), variable.name = 'climate')
bac_ddr = merge(bac_ddr, bac_map['State'], by.x = 'row', by.y = 'row.names') %>%
  merge(bac_map['State'], by.x = 'col', by.y = 'row.names')
bac_ddr$State = ifelse(bac_ddr$State.x == 'Exposed' & bac_ddr$State.y == 'Exposed', 'Exposed',
                   ifelse(bac_ddr$State.x == 'Closed' & bac_ddr$State.y == 'Closed','Closed','Others'))

bac_ddr = subset(bac_ddr, State != 'Others', select = -c(State.x, State.y))

bac_ddr$State = factor(bac_ddr$State, levels = c('Exposed','Closed'))
bac_ddr$climate = factor(bac_ddr$climate, 
                         levels = c('MAP', 'MAT','HII'),
                         labels = c(bquote(Delta~MAP),bquote(Delta~MAT),bquote(Delta~HII)))

绘图距离衰减曲线

bac_ddr_plt = ggplot(bac_ddr, aes(value, 1-Bray, color = State)) +
  facet_wrap(.~climate, scales = 'free', strip.position = 'bottom', labeller = label_parsed, nrow = 3)+
  ggpubr::color_palette("jco")+
  geom_smooth(method = 'lm', formula = y ~ x, show.legend = F)+
  labs(x = NULL, y = 'Community similarity') +
  theme_classic()+
  theme(legend.position = c(0.8, 0.8),
        legend.text = element_text(color = 'black'),
        panel.spacing.y = unit(-.5,'lines'),
        strip.placement = 'outside',
        strip.background = element_blank(),
        strip.text = element_text(color = 'black',vjust = 3),
        axis.text = element_text(color = 'black'),
        axis.title = element_text(color = 'black'))

非度量多维尺度分析 NMDS

bac_NMDS = ordinate(bac_physeq, method="NMDS", distance='bray')
bac_NMDS_vec = bac_NMDS$points %>% merge(bac_map, by = 'row.names')

bac_otutab = data.frame(otu_table(bac_physeq))
bac_otutab = bac_otutab[,rownames(bac_map)]
bac_otutab = bac_otutab[rowSums(bac_otutab)>0,]

方差分解 ANOVA并提取结果用于绘图

bac_permanova = adonis2(t(bac_otutab) ~ Exposed+Category+Material+HII+Climate, data = bac_map, permutations = 999, by = 'margin')
bac_permanova_stat = data.frame(bac_permanova)[1:5,] %>% arrange(desc(F))
bac_labl = data.frame(MDS1 = -0.55, MDS2 = 0.35, 
                       label = paste0(paste0(rownames(bac_permanova_stat),': ', 
                                      round(bac_permanova_stat$F,2), collapse = '\n')))
bac_labl$label =gsub('Exposed','Storage',bac_labl$label)

pal = colorRampPalette(c("#3E5CC5","#E6EB00","#65B48E", "#E64E00"))

绘图

library(cowplot) 
bac_NMDS_vec$State = factor(bac_NMDS_vec$State, levels = c('Exposed','Closed'))

bac_margin = ggplot(bac_NMDS_vec, aes(MDS1, fill = State)) +
  geom_density(alpha = .7, show.legend = F) + 
  theme_nothing()+
  ggpubr::fill_palette("jco")+
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  theme(plot.margin = margin(t = 1, l = 13, b = 1, r = 1, unit = 'mm'))

bac_pmain <- ggplot(bac_NMDS_vec, aes(x = MDS1, y = MDS2))+
  theme_classic()+
  geom_point(aes(color = MAP, shape = Category), size = 2, show.legend = F)+
  geom_text(data = bac_labl, aes(label = label),hjust = 0) +
  annotate('text', x = -0.3, y = 0.56, label = 'bold (PERMANOVA~(F~value))', parse = T) +
  labs(color = 'MAP') +
  scale_color_gradient2(low = "white",mid = 'red',midpoint = 1200,high = "purple", na.value = NA) +
  theme(legend.position = 'left',
        legend.key.height = unit(0.3,'cm'),
        legend.background = element_blank(),
        legend.text = element_text(color = 'black', size = 10),
        axis.title = element_text(color = 'black', size = 10),
        axis.text = element_text(color = 'black', size = 10))

合并细菌图

p1 <- ggarrange(bac_margin, bac_pmain, nrow = 2, heights = c(0.2, 1))

同理对于真菌

load('../../Result/Fungi/00_closed/00_load_data/.RData')
rm(result, otutab)
fungi_map = map
fungi_physeq = physeq

# (5) beta diversity, DDR
fungi_bray = distance(fungi_physeq, method = 'bray', type = 'sample')
fungi_bray_df = as.data.frame(as.matrix(fungi_bray))

fungi_MAT_df = as.data.frame(as.matrix(dist(fungi_map['MAT'],method = 'euclidean')))
fungi_MAP_df = as.data.frame(as.matrix(dist(fungi_map['MAP'],method = 'euclidean')))
fungi_HII_df = as.data.frame(as.matrix(dist(fungi_map['HII'],method = 'euclidean')))
names(fungi_MAT_df) == names(fungi_bray_df)
fungi_map$Exposed = ifelse(fungi_map$State == 'Exposed', 1, 0)
fungi_map = subset(fungi_map, 
                 select = -c(BioProject, Center.Name, geo_loc_country, geo_loc_continent, 
                             geo_loc_name, Instrument, Longitude, Sample_Name, Organism, 
                             Primers, FP_sequence, RP_sequence,Site, Source))

fungi_ddr = mtrx2cols(fungi_bray_df,fungi_MAP_df,fungi_MAT_df,fungi_HII_df,'Bray','MAP','MAT','HII')
head(fungi_ddr)
fungi_ddr = melt(fungi_ddr, id.vars = c('row','col','Bray'), variable.name = 'climate')
fungi_ddr = merge(fungi_ddr, fungi_map['State'], by.x = 'row', by.y = 'row.names') %>%
  merge(fungi_map['State'], by.x = 'col', by.y = 'row.names')
fungi_ddr$State = ifelse(fungi_ddr$State.x == 'Exposed' & fungi_ddr$State.y == 'Exposed', 'Exposed',
                       ifelse(fungi_ddr$State.x == 'Closed' & fungi_ddr$State.y == 'Closed','Closed','Others'))

fungi_ddr = subset(fungi_ddr, State != 'Others', select = -c(State.x, State.y))

fungi_ddr$State = factor(fungi_ddr$State, levels = c('Exposed','Closed'))
fungi_ddr$climate = factor(fungi_ddr$climate, 
                           levels = c('MAP', 'MAT','HII'),
                           labels = c(expression(Delta~MAP),expression(Delta~MAT),expression(Delta~HII)))

fungi_ddr_plt = ggplot(fungi_ddr, aes(value, 1-Bray, color = State)) +
  facet_wrap(.~climate, scales = 'free', strip.position = 'bottom', labeller = label_parsed, nrow = 3)+
  ggpubr::color_palette("jco")+
  # geom_point()+
  geom_smooth(method = 'lm', formula = y ~ x, show.legend = F)+
  labs(x = NULL, y = 'Community similarity') +
  theme_classic()+
  theme(legend.position = c(0.8, 0.8),
        legend.text = element_text(color = 'black'),
        panel.spacing.y = unit(-.5,'lines'),
        strip.placement = 'outside',
        strip.background = element_blank(),
        strip.text = element_text(color = 'black',vjust = 3),
        axis.text = element_text(color = 'black'),
        axis.title = element_text(color = 'black'))

fungi_NMDS = ordinate(fungi_physeq, method="NMDS", distance='bray')
fungi_NMDS_vec = fungi_NMDS$points %>% merge(fungi_map, by = 'row.names')

fungi_otutab = data.frame(otu_table(fungi_physeq))
fungi_otutab = fungi_otutab[,rownames(fungi_map)]
fungi_otutab = fungi_otutab[rowSums(fungi_otutab)>0,]

fungi_permanova = adonis2(t(fungi_otutab) ~ Exposed+Category+Material+HII+Climate, data = fungi_map, permutations = 999, by = 'margin')
fungi_permanova_stat = data.frame(fungi_permanova)[1:5,] %>% arrange(desc(F))
fungi_labl = data.frame(MDS1 = 0.8, MDS2 = 0.4, 
                      label = paste0(paste0(rownames(fungi_permanova_stat),': ', 
                                            round(fungi_permanova_stat$F,2), collapse = '\n')))
fungi_labl$label =gsub('Exposed','Storage',fungi_labl$label)

fungi_NMDS_vec$State = factor(fungi_NMDS_vec$State, levels = c('Exposed','Closed'))
fungi_margin = ggplot(fungi_NMDS_vec, aes(MDS1, fill = State)) +
  geom_density(alpha = .7) + 
  theme_nothing()+
  labs(fill = 'Storage') +
  ggpubr::fill_palette("jco")+
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  theme(legend.position = 'left',
        legend.key.size = unit(0.3,'cm'),
        plot.margin = margin(t = 1, l = 16, b = 1, r = 1, unit = 'mm'))


fungi_pmain <- ggplot(fungi_NMDS_vec, aes(x = MDS1, y = MDS2))+
  theme_classic()+
  geom_point(aes(color = MAP, shape = Category), size = 2)+
  geom_text(data = fungi_labl, aes(label = label), hjust = 1) +
  annotate('text', x = 0.5, y = .65, label = 'bold (PERMANOVA~(F~value))', parse = T) +
  labs(color = 'MAP') +
  scale_color_gradient2(low = "white",mid = 'red',midpoint = 1200,high = "purple", na.value = NA) +
  guides(fill = guide_legend(order = 1), shape = guide_legend(override.aes = list(size = 3), order = 2)) +
  
  theme(legend.position = 'left',
        legend.key.height = unit(0.3,'cm'),
        legend.background = element_blank(),
        legend.text = element_text(color = 'black', size = 10),
        axis.title = element_text(color = 'black', size = 10),
        axis.text = element_text(color = 'black', size = 10))
p2 <- ggarrange(fungi_margin, fungi_pmain, nrow = 2, heights = c(0.2, 1))

合并真菌和细菌图

p = ggarrange(p1, p2, ncol = 2, widths = c(0.75,1),labels = c('A','                 B'))
ggsave(p, filename = 'Figure2.jpg', width = 10, height = 5, dpi = 600)

Figure2.jpg

ddr.jpg