如何循环实现Tukey事后检验并绘制带显著性标记的基因表达图

一、前置准备（确认数据结构与加载依赖包） #

首先，假设你已经完成了数据整理和双因素ANOVA，你的统计数据框（比如叫gene_stats）应该包含以下列：

• gene：基因名称
• Time：时间分组（比如0h、24h、48h）
• Genotype：基因型（WT、mut1）
• mean_expr：分组均值
• se：标准误

先加载需要的R包，没装的话先执行安装：

# install.packages(c("dplyr", "lsmeans", "ggplot2", "multcompView"))

library(dplyr)
library(lsmeans)  # 用于Tukey多重比较
library(ggplot2)  # 绘图
library(multcompView)  # 生成紧凑字母标记

二、批量执行Tukey多重比较并获取紧凑字母标记 #

我们用dplyr的分组操作，给每个基因单独拟合模型并做Tukey检验，最后把显著性字母合并到统计数据中：

# 按基因分组，批量执行Tukey多重比较
tukey_results <- gene_stats %>%
  group_by(gene) %>%
  do(
    # 拟合双因素交互模型
    model = lm(mean_expr ~ Genotype * Time, data = .),
    # 用lsmeans做Tukey多重比较并生成紧凑字母
    tukey = lsmeans(., ~ Genotype * Time, adjust = "tukey") %>%
      cld(adjust = "tukey", Letters = letters, sort = FALSE)
  ) %>%
  unnest(tukey)  # 展开嵌套的结果

# 把字母标记合并回原始统计数据框（匹配基因、基因型、时间）
gene_stats_with_letters <- gene_stats %>%
  left_join(
    tukey_results %>% select(gene, Genotype, Time, .group),
    by = c("gene", "Genotype", "Time")
  ) %>%
  mutate(
    # 清理字母标记（去掉前面的分组编号和空格）
    sig_letters = gsub("^\\d+ ", "", .group)
  )

说明：cld()函数生成的.group列会带分组编号，我们用gsub()把它清理成纯字母，方便绘图时直接标注。

三、循环绘制每个基因的柱状图 #

这里用purrr的map()函数批量生成图，也可以换成for循环，每个基因输出一张带误差棒和显著性字母的柱状图：

library(purrr)

# 获取所有唯一的基因名称
gene_list <- unique(gene_stats_with_letters$gene)

# 批量绘图并保存（或者直接在RStudio中查看）
map(gene_list, function(g) {
  # 筛选当前基因的数据
  plot_data <- filter(gene_stats_with_letters, gene == g)
  
  # 绘制柱状图
  p <- ggplot(plot_data, aes(x = Time, y = mean_expr, fill = Genotype)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.7) +
    # 添加标准误误差棒
    geom_errorbar(
      aes(ymin = mean_expr - se, ymax = mean_expr + se),
      position = position_dodge(width = 0.8),
      width = 0.2
    ) +
    # 标注显著性字母（位置在误差棒上方）
    geom_text(
      aes(y = mean_expr + se + 0.05*max(mean_expr), label = sig_letters),
      position = position_dodge(width = 0.8),
      size = 4
    ) +
    # 自定义基因型填充色
    scale_fill_manual(values = c("WT" = "#1f77b4", "mut1" = "#ff7f0e")) +
    # 标题和坐标轴标签
    labs(title = paste("Gene Expression:", g),
         x = "Time", y = "Mean Expression", fill = "Genotype") +
    # 主题优化
    theme_bw() +
    theme(
      plot.title = element_text(hjust = 0.5, size = 14),
      axis.title = element_text(size = 12),
      legend.position = "top"
    )
  
  # 按基因名称保存图片（可选）
  ggsave(paste0(g, "_expr_plot.png"), p, width = 6, height = 4, dpi = 300)
  
  # 返回绘图对象，方便在RStudio中查看
  return(p)
})

提示：如果不想保存图片，删掉ggsave()那一行即可。你可以根据需求调整误差棒高度、字母位置、颜色方案等细节。

四、注意事项 #

• 如果你的ANOVA结果显示某基因的交互项不显著，可以把模型改成mean_expr ~ Genotype + Time，这样Tukey检验会更合理。
• 确保Time列是因子类型，否则ggplot会把它当成连续变量处理：gene_stats$Time <- as.factor(gene_stats$Time)
• 如果部分基因的Tukey检验没有显著性差异，sig_letters会显示相同的字母，这是正常结果。

内容的提问来源于stack exchange，提问作者Dendrobium