R：microtable包计算相对丰度堆叠柱状图

rm(list = ls())
setwd("C:\\Users\\Administrator\\Desktop\\New_microtable") #设置工作目录
library(microeco)
library(magrittr)
library(dplyr)
library(tibble)

feature_table <- read.table('Bac_genus.txt', header = TRUE, row.names = 1, sep = "\t", fill = TRUE) #特征表
# 检查并处理缺失值
if (any(is.na(feature_table))) {
  feature_table[is.na(feature_table)] <- 0  # 将缺失值填充为零
}
sample_table <- read.table('sample_table.txt', header = TRUE, row.names = 1, sep = "\t") #样本表
tax_table <- read.table('tax_table_g.txt', header = TRUE, row.names = 1, sep = "\t", fill = TRUE) #分类表

dataset <- microtable$new(sample_table = sample_table,
                          otu_table = feature_table, 
                          tax_table = tax_table)

dataset$tidy_dataset() #整理和预处理数据集
#数据清洗：移除或填补缺失值、异常值等。
#数据标准化：确保数据符合一定的格式，比如统一的数据类型。
#数据整合：如果有多个表格，确保它们之间的链接正确无误。

dataset$sample_sums() %>% range #计算并查看样本总数的范围

dataset$rarefy_samples(sample.size = 1000000) #执行重采样，标准化样本中的测序深度

dataset$sample_sums() %>% range #计算并查看标准化后样本总数的范围


dataset$cal_abund() #计算每个分类等级的分类群丰度
#class(dataset$taxa_abund)
dataset$taxa_abund$Genus[1:5, 1:5]

t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10, groupmean = "Subgroup")
#t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10, groupmean = "Group")
#t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 10)

write.table(t1$data_abund, file = "trans_abund_data_DAS.txt", sep = "\t", row.names = TRUE, col.names = NA, quote = FALSE)

########################################################################################

# 清空当前环境中的所有对象
rm (list = ls ())
setwd("C:\\Users\\Administrator\\Desktop\\New_microtable") #设置工作目录
library(ggplot2)   # 用于绘图
library(ggalluvial) # 用于绘制柱状图背后的条带
library(grid) # 用于访问unit函数
ra <- as.matrix(read.table("Genus_top15.txt", row.names =1, 
                           header = F, sep = "\t")) # 读入相对丰度数据并转换为矩阵方便后续数据整理

group <- c("B73", "Mo17")  # 品种变量
code <- c("B73_DAS28","B73_DAS42","B73_DAS56","B73_DAS70","Mo17_DAS28","Mo17_DAS42","Mo17_DAS56","Mo17_DAS70") # 8组处理变量
code_rep <- rep(code, each = 15)
taxa_rep <- rep(rownames(ra), 8)
cultivar_rep <- rep(group, each = 60)
abundance_vec <- as.vector(ra)
dat <- data.frame(code = code_rep,
                  Phylum = taxa_rep,
                  cultivar = cultivar_rep,
                  abundance = abundance_vec)

dat$Phylum <- factor(dat$Phylum, levels = c("Streptomyces","Sphingobium","Sphingomonas","Bradyrhizobium",
                                            "Variovorax","Lysobacter","Pseudomonas","Acidovorax","Nocardioides",
                                            "Burkholderia","Stenotrophomonas","Mesorhizobium","Microbacterium","Rhizobium","Roseateles"))
ggplot(dat, aes(x = code, y = abundance, fill = Phylum))+
  geom_bar(stat = "identity", width = 0.6, alpha = 1)+ 
  geom_flow(aes(alluvium = Phylum), alpha = 0.4) + 
  scale_fill_manual(values = c(
    "#F2A2C7", "#E9C46A", "#88C695", "#F4E3A1", "#76B4BD", "#C9A8C5", "#6D6DB5", "#5BA8A0", 
    "#9CAD60", "#C0C0C0", "#BC8F8F", "#76B4BD", "#6AB47B", "#89CFF0", "#8A2BE2"
  ))+
  theme_bw()+ 
  facet_grid(.~cultivar, scales = "free_x", space = "free_x")+ 
  xlab("")+ 
  ylab("Relative abundance (%)")+ 
  theme(panel.grid.major.x = element_blank(),
        axis.text.x = element_text(size = rel(1.3), angle = 0, face = "plain", color = "black", family = "Times New Roman"),
        axis.text.y = element_text(size=rel(1.8), face = "plain", color = "black", family = "Times New Roman"),
        legend.text = element_text(size = rel(1), family = "Times New Roman"),
        strip.text.x = element_text(size = rel(2), face = "bold", family = "Times New Roman"),
        strip.background = element_rect(fill = "lightblue", size = 1),
        panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
        axis.title.y = element_text(size = 18, family = "Times New Roman"),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.key.size = unit(0.6, "cm"),
        legend.spacing.x = unit(0.1, "cm"),
        legend.box.margin = margin(t = -20, unit = "pt")) + 
  guides(fill = guide_legend(ncol = 7)) + 
  scale_y_continuous(expand = c(0, 0.01))

ggsave("Abundance2.png", width = 10, height = 7, dpi = 1200) # 图片导出，导出为pdf文件，设置图片长和宽