R：LEfSe(线性判别分析)

rm(list=ls()) #清空工作环境
setwd("C:\\Users\\Administrator\\Desktop\\LDA") #设置工作目录
library(tidyverse) #包含了一系列与数据分析和可视化相关的包
library(microeco) #生态学分析的包
library(magrittr) #提供了用于简化代码的管道操作符 %>%

feature_table <- read.table('feature_table.txt', header = TRUE, row.names = 1)
sample_table <- read.table('sample_table.txt', header = TRUE, row.names = 1)
tax_table <- read.table('tax_table.txt', header = TRUE, row.names = 1)
#head(feature_table)[1:6,1:6]
#head(sample_table)[1:6, ]
#head(tax_table)[,1:6]

dataset <- microtable$new(sample_table = sample_table,
                          otu_table = feature_table, 
                          tax_table = tax_table)
#用于创建 microtable 对象。这个对象用于存储微生物组数据，通常包括样本表（sample_table）、OTU（操作性分类单元）表（feature_table）、和分类信息表（tax_table）。
dataset
#用于执行差异丰度分析（Differential Abundance Analysis）
lefse <- trans_diff$new(dataset = dataset,  #用于执行差异丰度分析，其中包含了不同组间的比较
                        method = "lefse", #选择差异丰度分析的方法为 LEfSe（LDA Effect Size）
                        group = "Group",  #指定了用于分组的变量
                        alpha = 0.01, # 是显著性水平，用于确定哪些差异是显著的
                        lefse_subgroup = "Subgroup", #指定子组别的变量名,NULL表示在 LEfSe 中不使用子组别。
                        taxa_level = "Genus", #分析的分类水平
                        #nresam = 0.9, #抽取原始样本的比例
                        lefse_min_subsam = 8
                        ) 
# 查看分析结果
head(lefse$res_diff[lefse$res_diff$LDA > 2.5, ])

# 将分析结果保存为TXT文件
write.table(lefse$res_diff[lefse$res_diff$LDA > 2.5, ], file = "lefse_result.txt", sep = "\t", quote = FALSE, row.names = TRUE)

# 找出LDA大于3的结果的索引
index <- which(lefse$res_diff$LDA > 2.5)
# 定义一个颜色映射
color_map <- c("B73" = "#00BFFF", "Mo17" = "#FF4500")
# 在你的图形代码中使用这个颜色映射
plot <- lefse$plot_diff_bar(use_number = index, 
                            width = 0.8, 
                            group_order = c("B73", "Mo17"))

plot <- plot + scale_color_manual(values = color_map) #颜色映射
plot <- plot + scale_fill_manual(values = color_map)
# 添加标题，并调整位置
plot <- plot + ggtitle("Biomarker (Genus)") +
  theme(plot.title = element_text(hjust = 0,vjust = 0, size = 23))
#横坐标的标题内容
plot <- plot + ylab("LDA SCORE")
plot <- plot + theme(axis.title.x = element_text(hjust = 0.46, vjust = -1))
#y轴刻度标签的字体大小
plot <- plot + theme(axis.text.y = element_text(size = 14))
# 调整图例的位置、大小、外观
plot <- plot + theme(text = element_text(size = 18, color = "darkblue")) ## 改变字体大小和颜色
# 使用legend.justification和legend.direction可以更灵活地定义图例的位置
plot <- plot + theme(legend.position = c(0.7, 1.07), 
                     legend.justification = c(0.7, 1.07), 
                     legend.direction = "horizontal", 
                     legend.text = element_text(size = 18), 
                     legend.key.size = unit(0.9, "cm"))
# 添加注释
#plot <- plot + annotate("text", x = 1, y = 1, label = "Your annotation here")
# 添加一个框把整个图框起来，可以自定义框的颜色和线宽
#plot <- plot + theme(plot.background = element_rect(color = "black", size = 1))
# 调整图形边缘与内容之间的距离 上、右、下、左
plot <- plot + theme(plot.margin = margin(2, 1, 1, 1, "cm"))
#使用的字体是新罗马
plot <- plot + theme(text = element_text(family = "Times New Roman"))
# 打印图形
print(plot)

# 保存图形为PNG格式，可以自定义大小、分辨率等
ggsave(filename = "your_plot2.png", plot = plot, width = 10, height = 10, dpi = 600)