使用 tidyverse 从线性模型中提取斜率的 pValue 并将其添加到绘图中

Question

这涉及到之前提出的问题，关于拟合基因表达数据的线性或倾斜模型。 我有各种基因在 5 个连续发育阶段的表达值，并且想要创建一个线性模型来测试基因在发育过程中的增加、减少或保持稳定。 由于我的数据点相当少，因此线性模型是最合适的。我资助了一种从线性模型中提取一些数据的方法，但无法获得斜率的 pValue，它可以告诉我斜率实际上是水平的（因此在开发过程中没有变化）。 下面，找到示例数据以及我设法得到的内容。

模型中的数据仅提供截距 p 值，但不提供斜率 p 值。

作为奖励，如果能够整合图中每个基因的斜率 p 值，那就太好了。

# packages
library('tidyverse') # ggplot & dplyr
library('magrittr') # pipe operations
library(broom) # tidy models

## generate data
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)  # Change 10 to your desired sample size
    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'))
head(exp)


# get data from linear model: missing: pvalue of the slope
dat.lm <- exp %>%
  group_by(Gene) %>%
  group_modify(~ broom::tidy(lm(Expression ~ Group, data = .x)))
head(dat.lm)

# plot
exp %>%
  ggplot(aes(x=Group, y=Expression, color = Gene, group = Gene)) +
  geom_smooth(method = lm, se = T, alpha = 0.1, aes(fill = Gene))

Answer 1

这里有几点需要注意。 首先，

geom_smooth(method = "lm")

不会产生与

broom::tidy(lm())

相同的输出。 lm 模型的拟合估计每个组/阶段的系数。因此，这里实际上没有斜率，而只是增加或减少截距表达式的系数。

geom_smooth()

将因子变量组视为连续变量，因此能够估计斜率。

因此，为了使模型拟合的结果与

geom_smooth()

一致并能够使用相应的 p 值计算斜率，我们需要首先将因子转换为数字。

然后我们可以通过

tableGrob()

包中的

gridExtra

添加估计值和 p 值。

library(tidyverse) 
library(magrittr) 
library(broom) 
library(moderndive)
library(ggplot2)
library(gridExtra)

## generate data
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'))
head(exp)

exp$Group <- as.numeric(exp$Group)

dat.lm <- exp %>%
  group_by(Gene) %>%
  group_modify(~ broom::tidy(lm(Expression ~ Group, data = .x)))
head(dat.lm)

exp %>%
  ggplot(aes(x = Group, y = Expression, color = Gene, group = Gene)) +
  geom_smooth(method = "lm", se = T, alpha = 0.1, aes(fill = Gene)) + 
  annotation_custom(tableGrob(dat.lm[,c(1,2,3,6)],
                              rows = NULL,
                              theme = ttheme_minimal(base_size = 6)),
                    xmin = 1.0, xmax = 2, ymin = 10.5, ymax = 11.3)

结果：