如何构建数据集以运行计数比时序二项式GLM并按基因输出检验结果

核心逻辑 #

你现有的数据格式无需转长格式，R的glm原生支持两列矩阵作为二项分布的响应变量，第一列为等位基因1的计数，第二列为等位基因2的计数，完全匹配你需要的(A1:A2)作为响应的模型结构。我们通过分组拟合模型+似然比检验的方式批量获取每个基因的Day整体效应p值。

完整代码 #

1. 依赖包加载 #

# 用到dplyr做分组处理，broom做模型结果提取，都属于tidyverse生态
library(dplyr)
library(broom)

2. 数据预处理 #

# 首先将Day转为因子类型，符合你不把Day作为数值变量的要求
df <- df %>%
  mutate(Day = as.factor(Day))

3. 批量拟合模型+提取p值 #

gene_p_result <- df %>%
  # 按基因分组，自动遍历所有基因
  group_by(Gene) %>%
  summarise(
    # 拟合含Day因子的二项式GLM
    full_model = list(glm(cbind(A.count_1, A.count_2) ~ Day, family = binomial)),
    # 拟合仅含截距的零模型，用于似然比检验
    null_model = list(glm(cbind(A.count_1, A.count_2) ~ 1, family = binomial)),
    # 似然比检验获取Day的整体效应p值
    raw_pval = anova(full_model[[1]], null_model[[1]], test = "Chisq")$`Pr(>Chi)`[2]
  ) %>%
  # 保留所需列，也可以自行保留模型对象用于后续查看细节
  select(Gene, raw_pval) %>%
  # 可选：添加多重检验校正后的p值，上百个基因检验推荐做校正
  mutate(adj_pval = p.adjust(raw_pval, method = "fdr"))

结果说明 #

运行上述代码后得到的gene_p_result数据框包含三列：

• Gene：基因编号
• raw_pval：每个基因对应的Day整体效应原始p值
• adj_pval：FDR校正后的p值，可根据需求选择校正方法

注意：多水平分类变量的整体显著性不能直接取模型summary中的单系数p值，上述方案通过似然比检验对比含Day的完整模型和截距零模型，得到的是Day对计数比的全局效应p值，完全匹配你的分析需求。如果需要查看单个Day水平与参照水平的差异显著性，可以提取对应基因的full_model对象用summary()查看即可。

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