如何在Matrix_eQTL_engine开展的eQTL分析中应用Bonferroni校正

前置说明 #

你提到的MatrixEQTL默认使用FDR校正属于常见误解：默认pvOutputThreshold参数对应的是原始未校正p值的过滤阈值，工具本身不会自动对输出结果做FDR校正，多重检验校正通常需要自己计算阈值或者后处理结果实现。

操作步骤 #

步骤1：计算Bonferroni校正阈值 #

Bonferroni校正的核心逻辑是将预设的显著性水平α（通常取0.05）除以总检验次数，总检验次数为你纳入分析的SNP数量乘以基因数量。R代码如下：

# 计算总检验量：SNP数 × 基因数
n_total_tests <- nrow(snps) * nrow(gene)
# 设定全局显著性水平
alpha <- 0.05
# 计算Bonferroni校正后的p值阈值
bonferroni_cutoff <- alpha / n_total_tests

步骤2：替换分析代码中的阈值 #

将Matrix_eQTL_engine函数中的pvOutputThreshold参数替换为上一步计算得到的bonferroni_cutoff即可，修改后的代码如下：

me <- Matrix_eQTL_engine(
    snps = snps,
    gene = gene,
    cvrt = cvrt,
    output_file_name = output_file_name,
    pvOutputThreshold = bonferroni_cutoff,
    useModel = modelLINEAR, 
    errorCovariance = errorCovariance, 
    verbose = TRUE,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE)

补充后处理方案 #

如果你希望先保留所有关联结果，再手动做Bonferroni校正筛选，可以对输出的原始p值调用R原生p.adjust函数处理：

# 提取所有关联的原始p值
raw_pvalues <- me$all$pvalue
# 计算Bonferroni校正后p值
me$all$p_bonferroni <- p.adjust(raw_pvalues, method = "bonferroni")
# 筛选校正后p值小于0.05的显著关联
sig_results <- me$all[me$all$p_bonferroni < 0.05, ]

注意事项 #

Bonferroni校正属于非常严格的多重检验校正方法，当检验量极大时阈值会极低，容易出现假阴性。如果需要控制家族错误率（FWER）同时保留更高的检验效能，可以替换为Holm校正，只需将p.adjust的method参数改为"holm"即可。

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