如何为含GLM交互项的分面ggplot自动添加Tukey检验注释

针对模型包含Treatment*Age交互项的情况，可借助emmeans和multcompView包批量处理每个基因的Tukey检验，并将分组注释自动添加到分面箱线图中，步骤如下：

1. 加载所需R包 #

library(ggplot2)
library(dplyr)
library(forcats)
library(emmeans)
library(multcompView)

2. 批量生成Tukey检验分组字母 #

对每个基因单独拟合GLM模型，计算Treatment与Age交互项的边际均值，通过Tukey检验生成组间差异的分组字母：

tukey_letters <- Raw_data %>%
  group_by(Gene) %>%
  do({
    # 拟合与绘图y轴一致的模型（log2(FC)）
    model <- glm(log2(FC) ~ Treatment * Age, data = .)
    # 计算交互项的边际均值
    emm <- emmeans(model, ~ Treatment * Age)
    # 执行Tukey多重比较
    tukey <- pairs(emm, adjust = "tukey")
    # 生成分组字母（相同字母表示组间无显著差异，p>0.05）
    letters <- multcompLetters(emmeans::contrast(tukey, method = "pairwise")$p.value,
                               threshold = 0.05)$Letters
    # 整理为绘图可用的数据框
    as.data.frame(emm) %>%
      mutate(Letter = as.character(letters),
             Gene = unique(.$Gene))
  }) %>%
  ungroup()

3. 调整注释的y轴位置 #

由于分面采用自由缩放的y轴，需为每个基因的每个分组计算注释的位置（置于分组最大值上方）：

annotation_pos <- Raw_data %>%
  group_by(Gene, Treatment, Age) %>%
  summarise(max_val = max(log2(FC)), .groups = "drop") %>%
  left_join(tukey_letters, by = c("Gene", "Treatment", "Age")) %>%
  mutate(y_pos = max_val + 0.1 * (max(max_val) - min(max_val))) # 可根据图的比例调整偏移系数

4. 将注释添加到分面图中 #

修改原绘图代码，加入geom_text绘制分组字母，通过position_dodge与箱线图对齐：

Test_plot <- ggplot(data = Raw_data,
                    aes(x = Age, y = log2(FC), fill = Treatment)) +
  geom_boxplot() +
  stat_boxplot(geom ='errorbar', width = 0.4,
               position = position_dodge(width = 0.75)) +
  geom_point(aes(group = Treatment),
             position = position_dodge(width = 0.75),
             pch = 1, size = 4, alpha = 0.5) +
  stat_summary(aes(group = Treatment),
               position = position_dodge(width = 0.75),
               geom = "point",
               fun = "mean",
               colour = "black",
               size = 2,
               pch =24,
               fill = "red") +
  # 添加Tukey分组注释
  geom_text(data = annotation_pos,
            aes(x = Age, y = y_pos, label = Letter, group = Treatment),
            position = position_dodge(width = 0.75),
            size = 4,
            fontface = "bold") +
  theme_bw() +
  theme(axis.title.x=element_blank(), 
        text = element_text(size=15),  
        plot.title = element_text(size=15, hjust = 0.5), 
        axis.text.x = element_text(size=10), 
        axis.text.y = element_text(size=8))+ 
  facet_wrap(Gene ~ ., ncol = 4, scales = "free") +
  ylab(expression(Log[2] ~ Fold ~ Change))

Test_plot

关键说明 #

• 模型响应变量需与绘图y轴保持一致（均使用log2(FC)），避免检验结果与可视化不匹配
• emmeans可灵活处理GLM的交互项边际均值，解决了传统aov+TukeyHSD在复杂模型下的局限性
• 分组字母由multcompView自动生成，无需手动整理检验结果

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