R：计算均值和标准差并进行多重比较

rm (list = ls ()) #清除所有变量
setwd("C:\\Users\\Administrator\\Desktop\\新建文件夹\\FAPROTAX微生物功能预测") #设置工作目录

# 加载所需的包
library(tidyverse)
library(agricolae)
library(dplyr)

# 读取数据
abundance_file <- "vegan.txt" # 功能丰度文件路径
group_file <- "metadata.txt" # 样本分组文件路径

abundance_data <- read.csv(abundance_file, sep = "\t", row.names = 1, check.names = FALSE)
group_data <- read.csv(group_file, sep = "\t", header = FALSE, col.names = c("Sample", "Group"))

# 整理分组信息
# 假设分组文件有两列，第一列为样本名，第二列为对应的分组
colnames(group_data) <- c("Sample", "Group")
# 转换丰度数据为长格式
abundance_long <- as.data.frame(t(abundance_data))
abundance_long$Sample <- rownames(abundance_long)
melted_data <- reshape2::melt(abundance_long, id.vars = "Sample")
# 合并分组信息
merged_data <- merge(melted_data, group_data, by = "Sample")

# 计算每个功能每一组的均值和标准差
group_summary <- merged_data %>%
  group_by(variable, Group) %>%
  summarise(Mean = mean(value), SD = sd(value))

# 指定分组的顺序
groups_order <- c("B73_DAS28", "B73_DAS42", "B73_DAS56", "B73_DAS70", 
                  "Mo17_DAS28", "Mo17_DAS42", "Mo17_DAS56", "Mo17_DAS70")

# 确保merged_data中的Group列是因子类型，并设置正确的水平
merged_data$Group <- factor(merged_data$Group, levels = groups_order)

# 对每个功能进行ANOVA分析和Tukey HSD测试
analysis_results <- list()
for(feature in unique(merged_data$variable)) {
  feature_data <- subset(merged_data, variable == feature)
  model <- aov(value ~ Group, data = feature_data)
  tukey <- HSD.test(model, "Group", console = TRUE)
  analysis_results[[feature]] <- tukey
}

# 准备输出数据框
output_data <- data.frame(Feature = character(),
                          Group = factor(levels = groups_order),
                          Mean = numeric(),
                          SD = numeric(),
                          Significance = character())

# 填充输出数据
for(feature in names(analysis_results)) {
  groups_info <- analysis_results[[feature]]$groups
  for(group_name in groups_order) { # 按指定顺序处理分组
    if(group_name %in% rownames(groups_info)){
      group_data <- subset(merged_data, variable == feature & Group == group_name)
      mean_val <- mean(group_data$value)
      sd_val <- sd(group_data$value)
      sig_mark <- as.character(groups_info[group_name, "groups"])
      output_data <- rbind(output_data, data.frame(Feature = feature,
                                                   Group = group_name,
                                                   Mean = mean_val,
                                                   SD = sd_val,
                                                   Significance = sig_mark))
    }
  }
}

# 输出结果到TXT文件
write.table(output_data, "output_results.txt", row.names = FALSE, sep = "\t", quote = FALSE)