R：microtable包线性判别分析LEfSe

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

feature_table <- read.table('Bac_species.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_s.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() #计算每个分类等级的分类群丰度

t1 <- trans_diff$new(dataset = dataset, method = "lefse", group = "Group", taxa_level = "Species", alpha = 0.01, lefse_subgroup = NULL)

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

library(ggplot2)
library(gridExtra)

# 定义分组的颜色
group_colors <- c("B73" = "#8FC9E2", "Mo17" = "#ECC97F")  # 可以根据需要选择颜色

# 生成条形图并设置填充色和边框色
g1 <- t1$plot_diff_bar(threshold = 2.5, width = 0.8, group_order = c("B73", "Mo17")) +
  scale_fill_manual(values = group_colors) +  # 设置条形填充颜色
  scale_color_manual(values = group_colors) +  # 设置条形边框颜色
  theme(legend.position = "none",axis.text.x = element_text(size = 12, face = "bold")
        ,text = element_text(family = "Times New Roman"))

# 生成丰度图并设置填充色和边框色
g2 <- t1$plot_diff_abund(threshold = 2.5, group_order = c("B73", "Mo17"), select_taxa = t1$plot_diff_bar_taxa) +
  scale_fill_manual(values = group_colors) +  # 设置丰度图填充颜色
  scale_color_manual(values = group_colors) +  # 设置丰度图边框颜色
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),axis.text.x = element_text(size = 12, face = "bold")
        ,text = element_text(family = "Times New Roman"))

# 生成并排显示的图形对象
combined_plot <- gridExtra::grid.arrange(g1, g2, ncol = 2, nrow = 1, widths = c(2, 1.7))
# 保存图形为 PNG 格式
ggsave("Species_LEfSe.png", plot = combined_plot, width = 10, height = 6, dpi = 800)