基于置换的多比较t检验：蛋白质组学数据FDR校正求助

数据背景 #

处理的大规模蛋白质组学数据结构如下（tibble格式）：

> shark_long_AN_E1_E2
# A tibble: 11,154 × 9
   MatchGroup_AUTO   genename    `Unique peptides` name          value Extraction Tissue replicate imputation 
   <fct>             <chr>       <fct>             <chr>         <dbl> <fct>      <fct>  <fct>     <chr>      
 1 XM_038815272.1.p1 UNANNOTATED 8                 E1_AN_1_Imp   1.10  E1         AN     1         E1_AN_1_Imp
 2 XM_038815272.1.p1 UNANNOTATED 8                 E1_AN_2_Imp   1.96  E1         AN     2         E1_AN_2_Imp
 3 XM_038815272.1.p1 UNANNOTATED 8                 E1_AN_3_Imp   0.994 E1         AN     3         E1_AN_3_Imp
 4 XM_038815272.1.p1 UNANNOTATED 8                 E2_AN_1_Imp1  7.05  E2         AN     1         Imp1       
 5 XM_038815272.1.p1 UNANNOTATED 8                 E2_AN_2_Imp1  7.45  E2         AN     2         Imp1       
 6 XM_038815272.1.p1 UNANNOTATED 8                 E2_AN_3_Imp1  8.07  E2         AN     3         Imp1       

需求：针对每个MatchGroup_AUTO，在E1和E2（Extraction分组）间做t检验，用置换法校正多假设检验的FDR。

现有步骤 #

已完成原始均值差计算和全局value置换：

1. 计算原始均值差 #

mean_diff <- shark_filtered %>%
  group_by(MatchGroup_AUTO) %>%
  summarise(mean_diff = mean(value[Extraction == 'E1']) - mean(value[Extraction == 'E2']))

2. 全局置换value #

set.seed(1979)  

n <- nrow(shark_filtered)  

P <- 100 

variable <- shark_filtered$value  

PermSamples <- matrix(0, nrow = n, ncol = P)

for (i in 1:P){
  PermSamples[, i] <- sample(variable, size = n, replace = FALSE)
}

后续操作指导 #

步骤1：关联置换后value与分组信息 #

将每一列置换后的value与原始数据的分组（MatchGroup_AUTO、Extraction等）绑定，生成置换数据集，再计算每个MatchGroup_AUTO的均值差。用purrr包批量处理的代码如下：

library(dplyr)
library(purrr)

# 提取原始数据的分组信息（剔除value列）
group_info <- shark_filtered %>% select(-value)

# 遍历每个置换样本列，生成置换数据集并计算均值差
perm_mean_diffs <- map_dfr(1:P, function(i) {
  # 绑定当前置换列的value到分组信息
  perm_data <- bind_cols(group_info, value = PermSamples[, i])
  # 计算当前置换下每个MatchGroup_AUTO的均值差
  perm_data %>%
    group_by(MatchGroup_AUTO) %>%
    summarise(perm_mean_diff = mean(value[Extraction == 'E1']) - mean(value[Extraction == 'E2']),
              perm_id = i)
})

步骤2：计算置换法p值 #

基于原始均值差和置换后的均值差，统计每个MatchGroup_AUTO的置换p值（原假设：E1与E2无差异）：

# 合并原始均值差与置换结果
combined_results <- mean_diff %>%
  left_join(perm_mean_diffs, by = "MatchGroup_AUTO")

# 计算置换p值：统计置换中绝对值≥原始绝对值的次数占比
perm_p_values <- combined_results %>%
  group_by(MatchGroup_AUTO, mean_diff) %>%
  summarise(p_value = mean(abs(perm_mean_diff) >= abs(mean_diff)), .groups = "drop")

步骤3：FDR校正 #

用置换得到的p值进行FDR校正（以BH方法为例）：

perm_p_values <- perm_p_values %>%
  mutate(fdr = p.adjust(p_value, method = "BH"))

优化提示：更合理的置换方式 #

当前使用的全局value置换，不如每个MatchGroup_AUTO内部置换Extraction标签贴合t检验原假设——即同一蛋白的E1/E2样本混合后重新分配分组标签，更符合“同一蛋白的E1和E2无差异”的假设逻辑，代码示例：

set.seed(1979)
P <- 100

# 每个置换循环中，对每个MatchGroup_AUTO内部置换Extraction标签
perm_mean_diffs_optimized <- map_dfr(1:P, function(i) {
  shark_filtered %>%
    group_by(MatchGroup_AUTO) %>%
    mutate(perm_extraction = sample(Extraction, size = n(), replace = FALSE)) %>%
    ungroup() %>%
    group_by(MatchGroup_AUTO) %>%
    summarise(perm_mean_diff = mean(value[perm_extraction == 'E1']) - mean(value[perm_extraction == 'E2']),
              perm_id = i)
})

后续计算p值和FDR的步骤与之前一致。

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