使用tidyverse提取线性模型斜率p值并整合至绘图

问题背景 #

拥有5个连续发育阶段的多基因表达数据，希望通过线性模型检验基因表达随发育的变化趋势（上升、下降或稳定）。当前代码中以Group（分类因子）拟合线性模型时，无法直接得到斜率的p值（用于判断趋势是否显著，即斜率是否为0），同时需要将每个基因的斜率p值整合到可视化图中。

解决步骤 #

1. 转换Group为连续变量，获取斜率的p值 #

原代码中Group是分类因子，拟合的线性模型生成的是各阶段与参考阶段的对比项，而非整体趋势的斜率。需将Group转换为连续数值（提取Stage后的数字），再拟合线性模型，此时模型系数即为斜率，对应的p值就是检验斜率是否为0的结果。

修改后的代码：

# 加载所需包
library(tidyverse)
library(magrittr)
library(broom)

# 生成示例数据
set.seed(123)
num_genes <- 5
num_groups <- 5
exp <- data.frame()
for (gene_id in 1:num_genes) {
  gene_name <- paste("Gene", gene_id, sep = "_")
  for (group_id in 1:num_groups) {
    group_name <- paste("Stage", group_id, sep = "_")
    expression_values <- rnorm(10, mean = 10, sd = 2)
    group_data <- data.frame(Gene = gene_name, Group = group_name, Expression = expression_values)
    exp <- rbind(exp, group_data)
  }
}
exp$Group <- factor(exp$Group, c('Stage_1',  'Stage_2', 'Stage_3', 'Stage_4', 'Stage_5'))

# 将Group转换为连续数值变量
exp <- exp %>%
  mutate(Group_num = as.integer(str_extract(Group, "\\d+")))

# 拟合线性模型，提取斜率的p值
dat.lm <- exp %>%
  group_by(Gene) %>%
  group_modify(~ broom::tidy(lm(Expression ~ Group_num, data = .x))) %>%
  filter(term == "Group_num") # 仅保留斜率相关结果

head(dat.lm)

此时dat.lm中的p.value列即为每个基因斜率的显著性p值，可用于判断表达趋势是否显著。

2. 将p值整合到绘图中 #

可通过geom_text()将p值标注在对应基因的拟合线附近，或根据p值设置线条样式区分显著性。以下示例将p值标注在拟合线末端：

# 合并表达数据与p值数据
exp_with_p <- exp %>%
  left_join(dat.lm %>% select(Gene, p.value) %>% rename(slope_p = p.value), by = "Gene")

# 绘制带p值标注的趋势图
ggplot(exp_with_p, aes(x = Group_num, y = Expression, color = Gene, group = Gene)) +
  geom_jitter(alpha = 0.3) + # 展示原始散点数据
  geom_smooth(method = lm, se = T, alpha = 0.1, aes(fill = Gene)) +
  # 在拟合线末端标注p值
  geom_text(data = exp_with_p %>% group_by(Gene) %>% slice(1),
            aes(x = 5.2, y = predict(lm(Expression ~ Group_num, data = .)), 
                label = paste("p =", round(slope_p, 4))),
            hjust = 0) +
  scale_x_continuous(breaks = 1:5, labels = paste("Stage", 1:5)) +
  labs(x = "发育阶段", y = "基因表达量") +
  theme_bw()

若需添加显著性标记（如*代表p<0.05），可修改label参数：

label = case_when(
  slope_p < 0.001 ~ "***",
  slope_p < 0.01 ~ "**",
  slope_p < 0.05 ~ "*",
  TRUE ~ "ns"
)

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