在R中获取40种蛋白对应2个分类临床变量的Anova检验各p值的方法

你当前仅得到1个p值的核心原因是aov()函数单次调用默认仅处理1个响应变量（即单种蛋白），没有批量遍历所有40种蛋白分别拟合模型，你可以按照下方方法批量提取所有蛋白对应的p值。

前提说明 #

假设你的数据框命名为df，结构符合以下规则：前2列为需要纳入分析的分类型临床变量，后续列均为40种蛋白的定量数据。
请将后续代码中的「临床变量1名」「临床变量2名」替换为你数据中实际的两个临床变量的列名，蛋白列的范围也要根据你的实际数据调整

批量运行代码 #

# 1. 提取所有蛋白对应的列名，自行调整列的范围匹配你的数据
protein_cols <- colnames(df)[3:42]

# 2. 遍历所有蛋白拟合ANOVA模型并提取p值
p_value_list <- lapply(protein_cols, function(pro_name) {
  # 构建模型公式，若需加入两个临床变量的交互效应，可将+改为*
  model_form <- as.formula(paste0(pro_name, " ~ 临床变量1名 + 临床变量2名"))
  aov_fit <- aov(model_form, data = df)
  # 提取两个临床变量对应的p值
  anova_table <- anova(aov_fit)
  p_vals <- anova_table$`Pr(>F)`[1:2]
  # 返回结果
  return(c(蛋白名称 = pro_name, 变量1_p值 = p_vals[1], 变量2_p值 = p_vals[2]))
})

# 3. 将结果转换为数据框方便查看和导出
p_value_df <- do.call(rbind, p_value_list)
p_value_df <- as.data.frame(p_value_df, stringsAsFactors = FALSE)
p_value_df[, 2:3] <- apply(p_value_df[, 2:3], 2, as.numeric)

# 4. 查看所有蛋白的p值结果
View(p_value_df)

补充操作 #

• 若需要做多重检验校正，可使用p.adjust()函数对p值列直接校正，示例如下：

# 对第一个临床变量的p值做FDR校正
p_value_df$变量1_校正后p值 <- p.adjust(p_value_df$变量1_p值, method = "fdr")

• 如需导出结果到本地，可使用write.csv()函数：

write.csv(p_value_df, "所有蛋白ANOVA_p值结果.csv", row.names = FALSE, fileEncoding = "GBK")

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