R语言数据处理那些事（二十七）计算蛋白组的差异倍数和显著性

缘由

生信之路，道阻且长，记录经验，但渡有缘人。

碎碎念

及时当勉励，岁月不饶人 ——陶渊明

代码

在数据分析中，我们经常需要 “按组计算” ——
比如蛋白组学中：

对每个蛋白 (Protein)，再按每个肽段 (Peptide) 分组，
计算每个实验组 (group) 的平均表达量。

在组学数据分析（蛋白质组、代谢组）中，我们经常需要：

• 按某个描述列（如蛋白 ID、代谢物 ID）分组
• 对每一组内的两个实验组（control vs treat）做 t 检验
• 自动判断是否方差齐性（Levene test）
• 自动选择 t 检验类型（方差齐 / 非方差齐）
• 输出 log2FC、FC、p 值等结果

1、模拟示例数据

library(dplyr)
library(rstatix)
 
set.seed(123)
 
# 示例数据：每个 ID 有 6 个样本（3 control + 3 treat）
df3 <- data.frame(
  ID = rep(paste0("Protein_", 1:5), each = 6),
  group = rep(rep(c("Control", "Treat"), each = 3), 5),
  value = rnorm(30, mean = rep(c(10, 12), each = 3, times = 5), sd = 1)
)
 
head(df3, 10)

---

2 、函数 1：自动检验方差齐性 + t 检验

cal_test <- function(data, value_col = "value"){
    
    # 创建 formula，例如 value ~ group
    fml <- as.formula(paste(value_col, "~ group"))
    
    # Levene 检验：判断方差是否齐性
    var_res <- data %>% levene_test(fml)
    var_equal <- var_res$p > 0.05  # p > 0.05 → Variance equal
    
    # 根据方差齐性自动选择 t 检验
    if (var_equal) {
        stat_res <- t_test(data, fml, var.equal = TRUE, detailed = TRUE) %>%
            mutate(method = "t-test equal var")
    } else {
        stat_res <- t_test(data, fml, var.equal = FALSE, detailed = TRUE) %>%
            mutate(method = "t-test no equal var")
    }
    
    return(stat_res)
}

✔ 解释

• levene_test()：用于判断两组是否方差齐性
• 如果 p > 0.05 → 方差齐 → 使用经典 Student t-test
• 如果 p ≤ 0.05 → 非方差齐 → 使用 Welch t-test
• 输出包含均值、标准差、统计值、p 值等

---

3. 函数 2：计算 log2FC、FC，整理结果

get_rstatix_mean <- function(stat_res, desc_col, control, treat){
    
    stat_res1 <- stat_res %>%
        mutate(
            # log2 fold change
            log2FC = case_when(
                group1 == treat & group2 == control ~ estimate1 - estimate2,
                group1 == control & group2 == treat ~ estimate2 - estimate1,
                TRUE ~ NA_real_
            ),
            FC = 2^log2FC,
            compare_group = paste0(treat, "_vs_", control)
        ) %>%
        dplyr::select(all_of(desc_col), FC, p, compare_group, method)
    
    return(stat_res1)
}

✔ 解释

• estimate1、estimate2 是 rstatix t_test 输出的组均值
• log2FC = Treat − Control
• FC = 2
• 添加比较组名和检验方法

---

4. 对所有 ID 批量进行 t 检验

假设：

desc_col <- "ID"
control <- "Control"
treat   <- "Treat"

批量分析：

stat_list <- df3 %>%
    filter(group %in% c(control, treat)) %>%
    mutate(group = factor(group)) %>%
    group_by(!!sym(desc_col)) %>%
    group_map(
        ~ {
            desc_val <- as.character(.y[[1]])
            
            # 每个 ID 的 t 检验
            stat_res <- cal_test(.x, value_col = "value") %>%
                mutate(!!sym(desc_col) := desc_val)
            
            # 计算 log2FC、FC 等
            get_rstatix_mean(stat_res, desc_col, control, treat)
        }
    )

合并结果：

final_df <- bind_rows(stat_list)
final_df

---

5、最终结果示例

输出格式：

ID	FC	p	compare_group	method
Protein_1	5.21	0.0488	Treat_vs_Control	t-test equal var
Protein_2	7.34	0.0147	Treat_vs_Control	t-test no equal var