在R中对分组SNP数据执行Cochran-Armitage趋势测试的问题排查

问题背景 #

拥有812222行数据，结构如下：

X            snp_id         is_severe encoding_1 encoding_2 encoding_0
1     0  GL000191.1-37698         0          0          1          7
2     1  GL000191.1-37698         1          0          2         11
3     3 GL000192.1-100085         0          5          3          0
4     4 GL000192.1-100085         1          3         10          0

已编写Cochran-Armitage趋势测试函数：

catt <-
  function(y, x0,x1,x2, score = c(0, 1, 2)) {
    h=data.frame(is_severe=c(0,1),
                 encod_0=c(0,0),
                 encod_1=c(0,0),
                 encod_2=c(0,0)
    )
    h$is_severe=y
    h$encod_0=x0
    h$encod_1=x1
    h$encod_2=x2
    my_mat <- as.matrix(h)
    A <- matrix(0, 2, 3)
    colnames(A) <- c(0, 1, 2)
    rownames(A) <- c(0, 1)
    A[1,] =my_mat[1,2:4]
    A[2, ]=my_mat[2,2:4]
    ptt <- prop.trend.test(A[1, ],colSums(A), score = score) 
    
    res <- list(
      chisq = as.numeric(ptt$statistic), 
      p.value = as.numeric(ptt$p.value)
    )
    return(res)
  }

遇到的问题：

1. 使用by()分组执行能得到正确结果，但返回的不是数据框格式：

by(union,union$snp_id, function (union) {catt(union$is_severe,union$encoding_0,union$encoding_1,union$encoding_2)})

输出示例：

union$snp_id: chr17-3343519
$`2x3-table`
  0 1 2
0 1 6 1
1 8 4 1

$chisq
[1] 3.47929

$p.value
[1] 0.0621414

------------------------------------------------------------------------ 
union$snp_id: chr17-33445549
$`2x3-table`
   0 1 2
0 11 2 0
1  8 0 0

$chisq
[1] 1.360324

$p.value
[1] 0.2434814

2. 使用data.table分组时出现错误：

library("data.table")
setDT(union)
union[,.(chisq=catt(union$is_severe,union$encoding_0,union$encoding_1,union$encoding_2)[1],p.value=catt(union$is_severe,union$encoding_0,union$encoding_1,union$encoding_2)[2]),snp_id]

错误信息：

Error in `$<-.data.frame`(`*tmp*`, "is_severe", value = c(0L, 1L, 0L,  : 
  replacement has 812222 rows, data has 2

解决方案 #

方法1：将by()的结果转换为数据框 #

by()返回的是by类对象，可通过do.call(rbind, ...)结合lapply把每个分组结果转为数据框行，再合并：

# 执行by操作并处理每个分组结果
by_result <- by(union, union$snp_id, function(sub_df) {
  res <- catt(sub_df$is_severe, sub_df$encoding_0, sub_df$encoding_1, sub_df$encoding_2)
  # 将列表转为单行数据框
  data.frame(chisq = res$chisq, p.value = res$p.value)
})

# 合并所有分组结果
final_df <- do.call(rbind, by_result)
# 将行名转为snp_id列并重置行名
final_df$snp_id <- rownames(final_df)
rownames(final_df) <- NULL
# 调整列顺序（可选）
final_df <- final_df[, c("snp_id", "chisq", "p.value")]

方法2：修复data.table的用法错误 #

错误原因是在data.table分组表达式中，使用union$引用了整表数据，而非当前分组的子数据。直接引用列名即可，无需加union$：

library(data.table)
setDT(union)

# 分组执行测试并生成数据框格式结果
final_dt <- union[, {
  res <- catt(is_severe, encoding_0, encoding_1, encoding_2)
  .(chisq = res$chisq, p.value = res$p.value)
}, by = snp_id]

# 如需转为普通数据框（可选）
final_df <- as.data.frame(final_dt)

额外优化：简化catt函数（可选） #

原函数中创建数据框再转矩阵的步骤可简化，直接构造测试所需的2x3矩阵：

catt <- function(y, x0, x1, x2, score = c(0, 1, 2)) {
  # 直接构造趋势测试用矩阵
  A <- matrix(c(x0, x1, x2), nrow = 2, byrow = TRUE)
  rownames(A) <- c(0, 1)
  colnames(A) <- c(0, 1, 2)
  
  ptt <- prop.trend.test(A[1, ], colSums(A), score = score)
  
  list(
    chisq = as.numeric(ptt$statistic),
    p.value = as.numeric(ptt$p.value)
  )
}

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