如何对分两类的计数数据各列执行Fisher精确检验？

问题背景 #

现有样本数据，前8列为分类元数据（示例中仅前2列：sampleID、category），后续列是各基因的存在/缺失标记（1/0）。样本按category分为两类，需要对每个基因列执行Fisher精确检验，将结果整理为包含对应基因名的表格。此前尝试的for循环、apply方法因列引用错误无法正常运行。

示例数据：

mydata <- data.frame(
  sampleID = c("A", "B", "C", "D", "E", "F", "G"),
  category = c("high", "low", "high", "high", "low", "high", "low"),
  Gene1 = c(1, 1, 0, 0, 0, 1, 1),
  Gene2 = c(0, 1, 1, 1, 1, 1, 0),
  Gene3 = c(0, 0, 0, 1, 1, 1, 1)
)

单独对某一列执行检验的可行代码：

# 构建列联表
contingency_table <- with(mydata, table(category, Gene1))  
# 执行Fisher检验
fisher.test(contingency_table)

此前尝试的错误代码问题：

• for循环中直接用x（列索引）构建表，table(category, x)实际是将索引值与category关联，而非基因列数据
• apply版本中错误引用Site（应为category），且列引用逻辑错误

解决方案 #

方法1：for循环实现（易理解，适合新手） #

# 创建空数据框存储结果，预设列名
result_df <- data.frame(
  Gene = character(),
  P_value = numeric(),
  Odds_Ratio = numeric(),
  Confidence_Lower = numeric(),
  Confidence_Upper = numeric(),
  stringsAsFactors = FALSE
)

# 遍历所有基因列：示例中跳过前2列，实际使用时替换为9:ncol(mydata)
for (x in 3:ncol(mydata)) {
  # 获取当前列的基因名
  gene_name <- colnames(mydata)[x]
  # 提取当前基因列的数据
  gene_data <- mydata[[x]]
  # 构建列联表：按category分组，统计基因的1/0数量
  contingency_table <- table(mydata$category, gene_data)
  
  # 执行Fisher精确检验，设置conf.int=TRUE计算置信区间
  test_result <- fisher.test(contingency_table, conf.int = TRUE)
  
  # 将结果添加到数据框中
  result_df <- rbind(result_df, data.frame(
    Gene = gene_name,
    P_value = test_result$p.value,
    Odds_Ratio = test_result$estimate,
    Confidence_Lower = test_result$conf.int[1],
    Confidence_Upper = test_result$conf.int[2],
    stringsAsFactors = FALSE
  ))
}

# 查看结果
print(result_df)
# 导出结果到CSV文件
write.csv(result_df, "fisher_test_results.csv", row.names = FALSE)

代码解释

• 空结果数据框：提前定义存储结构，避免循环中重复创建带来的性能问题
• 列遍历逻辑：3:ncol(mydata)指定从第3列开始遍历，实际场景替换为9:ncol(mydata)跳过前8列元数据
• 列引用方式：mydata[[x]]通过列索引提取对应列的向量，colnames(mydata)[x]获取基因名，这是解决之前错误的核心
• Fisher检验参数：conf.int=TRUE开启置信区间计算，方便后续分析
• 结果存储：用rbind将每一行结果追加到数据框，最后导出为CSV文件

方法2：apply函数实现（更简洁） #

# 定义处理单个基因列的函数
fisher_fun <- function(gene_col) {
  # 构建列联表：mydata$category是分组变量，gene_col是当前基因列的向量
  contingency_table <- table(mydata$category, gene_col)
  # 执行检验并返回结果列表
  test_result <- fisher.test(contingency_table, conf.int = TRUE)
  # 返回需要的结果项
  return(c(
    P_value = test_result$p.value,
    Odds_Ratio = test_result$estimate,
    Conf_Lower = test_result$conf.int[1],
    Conf_Upper = test_result$conf.int[2]
  ))
}

# 对基因列批量应用函数：示例中取第3到最后一列，实际替换为9:ncol(mydata)
gene_results <- apply(mydata[, 3:ncol(mydata)], MARGIN = 2, FUN = fisher_fun)

# 将结果转置为数据框，并添加基因名列
result_df <- as.data.frame(t(gene_results))
result_df$Gene <- rownames(result_df)
# 调整列顺序，把Gene列放前面
result_df <- result_df[, c("Gene", "P_value", "Odds_Ratio", "Conf_Lower", "Conf_Upper")]

# 查看结果
print(result_df)
# 导出结果
write.csv(result_df, "fisher_test_results_apply.csv", row.names = FALSE)

代码解释

• 自定义函数fisher_fun：接收单个基因列的向量，完成列联表构建、检验执行，并返回需要的统计量
• apply参数：MARGIN=2表示按列处理数据，mydata[,3:ncol(mydata)]筛选出所有基因列
• 结果转置：apply默认返回的结果是按行存储统计量，用t()转置后更符合表格格式
• 列名调整：将行名（原基因列名）转为单独的Gene列，优化结果可读性

关键注意事项 #

• 确保category列是分类变量（因子或字符型均可），基因列是0/1的数值型或整数型
• 若列联表中出现0值，Fisher检验仍可正常执行，但需注意结果的解释合理性
• 批量检验后建议添加多重检验校正（如p.adjust(result_df$P_value, method = "fdr")），避免假阳性

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