加载R包

library(stringr)
library(RColorBrewer)
library(scales)
library(reshape2)
library(ggpubr)
library(dplyr)
library(ggsci)

加载数据

load('../../Result/Bacteria/SILVA138/Filter_Outliers/00_load_data/.RData')

rm(physeq, result, envs, idx)

确认topN优势物种

otutab = data.frame(t(t(otutab)/colSums(otutab)))

top_bac = c('Alphaproteobacteria','Gammaproteobacteria','Actinobacteriota','Cyanobacteria','Firmicutes','Bacteroidota','Acidobacteriota','Chloroflexi','Planctomycetota','Nitrospirota')

phylum = merge(otutab, taxa['Phylum'], by = 'row.names')
phylum$Row.names = NULL
phylum = subset(phylum, Phylum %in% top_bac)
phylum = aggregate(.~Phylum, sum, data = phylum)
rownames(phylum) = phylum$Phylum
phylum$Phylum = NULL

sub_phylum = merge(otutab, taxa['Class'], by = 'row.names')
sub_phylum$Row.names = NULL
sub_phylum = subset(sub_phylum, Class %in% top_bac)
sub_phylum = aggregate(.~Class, sum, data = sub_phylum)
rownames(sub_phylum) = sub_phylum$Class
sub_phylum$Class = NULL
sub_phylum = sub_phylum[c(1,3),]

第一张图

mixed = rbind(phylum, sub_phylum)
mixed = merge(t(mixed), map[c('MAT','MAP','HII')], by = 'row.names')
mixed$Row.names = NULL

mixed_env = melt(mixed, measure.vars = c('MAT','MAP','HII'), variable.name = 'ENV', value.name = 'Number') %>%
  melt(id.vars = c('ENV','Number'), variable.name = 'taxa', value.name = 'value')

mixed_env$ENV = factor(mixed_env$ENV, levels = c('MAP','MAT','HII'), labels = c('MAP','MAT','Human influence'))
mixed_env$taxa = gsub('Gamma','γ-',mixed_env$taxa)
mixed_env$taxa = gsub('Alpha','α-',mixed_env$taxa)
mixed_env_reg = ggplot(subset(mixed_env, value !=0), aes(x = Number, y = value, color = taxa)) +
  # geom_point() +
  facet_grid(.~ENV, scales = "free",switch = 'both') +
  scale_color_brewer(palette = 'Paired') +
  geom_smooth(method = 'lm', formula = y~x) +
  scale_y_continuous(labels = percent, expand = c(0,0), breaks = c(0,.1,.2,.3,.4)) +
  guides(color = guide_legend(override.aes = list(size = 2, fill = NA)),
         shape = guide_legend(override.aes = list(size = 5))) +
  labs(x = NULL, y = 'Proportion', color = NULL) +
  theme_classic()+
  theme(strip.text = element_text(color = 'black'),
        axis.title = element_text(color = 'black'),
        axis.text = element_text(color = 'black'),
        strip.background = element_blank(),
        strip.placement = 'outside',
        legend.key.size = unit(0.3,'cm'),
        legend.text = element_text(color = 'black'),
        legend.title = element_text(color = 'black'))
# ggsave(mixed_env_reg, filename = 'mixed_env_reg.jpg', width = 8, height = 3, dpi = 600)

同理，第二张图

load('../../Result/Fungi/00_closed/00_load_data/.RData')

rm(physeq, result, envs, idx)

# 确认topN优势物种
otutab = otutab[rownames(map)]
otutab = otutab[rowSums(otutab)>0,]
otutab = data.frame(t(t(otutab)/colSums(otutab)))
taxa = taxa[rownames(otutab),]

class = merge(otutab, taxa['Class'], by = 'row.names')
class$Row.names = NULL

Ascomycota = c("Sordariomycetes","Dothideomycetes","Eurotiomycetes", "Lecanoromycetes", "Saccharomycetes", "Leotiomycetes")
Basidiomycota = c("Agaricomycetes","Malasseziomycetes","Microbotryomycetes", "Tremellomycetes")
Mortierellomycota = "Mortierellomycetes"

class$Class = ifelse(class$Class %in% c(Ascomycota, Basidiomycota, Mortierellomycota), as.character(class$Class), 'Others')
class = aggregate(.~Class, sum, data = class)
rownames(class) = class$Class
class$Class = NULL
class_stat = merge(t(class), map['Material'], by = 'row.names')
class_stat$Row.names = NULL
class_stat = aggregate(.~Material, mean, data = class_stat) %>% 
  melt(variable.name = 'Class', value.name = 'value')

class_stat$Class = factor(class_stat$Class, levels = c(Ascomycota, "Others", Basidiomycota, Mortierellomycota))
class_stat$Material = factor(class_stat$Material, 
                             levels = c('brick','mural','dolomite','soil','limestone','canvas','wood',
                                        'granite','basalt','sandstone','marble','calcite'))

p_taxa = ggplot(class_stat, aes(x = Material, y = value, fill = Class)) +
  # facet_grid(.~Material,scales = 'free',switch = 'x')+
  geom_bar(stat = "identity", position = 'stack',width = 0.7) +
  scale_fill_manual(values = c(brewer.pal(9, 'Reds')[7:2], 'grey', brewer.pal(9,'Greens')[7:4], 'Steelblue')) +
  # scale_fill_brewer('Reds') +
  scale_y_continuous(labels = percent, expand = c(0,0)) +
  guides(fill = guide_legend(ncol = 1)) +
  theme_classic() +
  labs(x = NULL, y = 'Proportion',fill = NULL) +
  theme(axis.title = element_text(color = 'black'),
        axis.text = element_text(color = 'black'),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        legend.key.size = unit(0.3,'cm'),
        legend.text = element_text(color = 'black'))
# ggsave(class_plot_facet, filename = 'class_facet1.jpg', width = 10, height = 4, dpi = 600)

第三张图

taxa$Genus = gsub(' ',"",taxa$Genus)
trait = read.table('../../Result/Fungi/FungalTraits.txt', header = T, sep = '\t')
trait$primary_lifestyle = ifelse(trait$primary_lifestyle == "", 'unknown', as.character(trait$primary_lifestyle))

pathogens = unique(c(grep("parasite", trait$primary_lifestyle, value = T),grep('patho',trait$primary_lifestyle, value = T), 'sooty_mold'))
decomposers = unique(grep("saprotroph",trait$primary_lifestyle, value =T))
symbionts = c(unique(grep("symbio", trait$primary_lifestyle, value =T)),'ectomycorrhizal','lichenized','arbuscular_mycorrhizal','ericoid_mycorrhizal')
setdiff(unique(trait$primary_lifestyle), c(pathogens, decomposers, symbionts))

trait$Guild = ifelse(trait$primary_lifestyle %in% decomposers | trait$Secondary_lifestyle %in% decomposers, 'Decomposers',
                     ifelse(trait$primary_lifestyle %in% pathogens  | trait$Secondary_lifestyle %in% pathogens, 'Pathogens',
                            ifelse(trait$primary_lifestyle %in% symbionts  | trait$Secondary_lifestyle %in% symbionts, 'Symbionts',"Others")))


# trait$Guild = ifelse(trait$primary_lifestyle %in% c(""))
taxa = merge(taxa, trait[c('GENUS','Guild')], by.x = 'Genus', by.y = 'GENUS', all.x = T)
taxa$Guild = ifelse(is.na(taxa$Guild), 'Others',as.character(taxa$Guild))

func_table = merge(taxa[c('otuid','Guild')], otutab, by.x = 'otuid', by.y = 'row.names')

rownames(func_table) = func_table$otuid
func_table$otuid = NULL
func = aggregate(.~Guild, func_table, sum)
rownames(func) = func$Guild
func$Guild = NULL
sort(rowMeans(func), decreasing = T)

# for raw otu tables
func_t = merge(data.frame(t(func)), map['Material'], by = 'row.names')
rownames(func_t) = func_t$Row.names
func_t$Row.names = NULL

func_t_df = melt(func_t, variable.name = 'Guilds', value.name = 'proportion') %>% 
  aggregate(proportion~Material*Guilds, mean, data = .)
func_t_df$Material = factor(func_t_df$Material, 
                            levels = c('brick','mural','dolomite','soil','limestone','canvas','wood',
                                       'granite','basalt','sandstone','marble','calcite'))

func_t_df$Guilds = gsub('Decomposers','Saprotrophs', func_t_df$Guilds)
p_function = ggplot(func_t_df, aes(Material, proportion, fill = Guilds)) +
  # facet_grid(.~Material,scales = 'free',switch = 'x')+
  geom_bar(stat = "identity", position = 'stack',width = 0.7) +
  scale_fill_jama() +
  # scale_fill_brewer('Reds') +
  scale_y_continuous(labels = percent, expand = c(0,0)) +
  guides(fill = guide_legend(ncol = 1)) +
  theme_classic() +
  labs(x = NULL, y = 'Proportion',fill = NULL) +
  theme(axis.title = element_text(color = 'black'),
        axis.text = element_text(color = 'black'),
        axis.text.x = element_text(angle = 30, hjust = 1),
        # axis.ticks.x = element_blank(),
        panel.grid = element_blank(),
        legend.text = element_text(color = 'black'),
        legend.key.size = unit(0.3,'cm'),
        plot.margin = margin(l = 6, b = 6, t = 6, r = 36.6))

合并

p = ggarrange(mixed_env_reg, p_taxa,  p_function, nrow = 3, heights = c(1,0.8,0.8), labels = c('A','B','C'))

Figure3.jpg