R语言～多图合并
2022-07-22 本文已影响0人
Oodelay
加载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
