封装edgeR标准化代码为函数时出现`$<-.data.frame`错误的原因排查

问题背景 #

单独执行edgeR标准化代码可正常运行：

d0 = DGEList(raw)
group = as.factor(metadata$Benefit)
d0$samples$group = group
keep.exprs = edgeR::filterByExpr(d0, group = group)
d0 = d0[keep.exprs,, keep.lib.sizes = FALSE]
d0 = calcNormFactors(d0 , method = 'TMM' )

但封装为limmaNorm函数后，调用d0 = limmaNorm(metadata, raw, 'Benefit')时出现报错：

Error in `$<-.data.frame`(`*tmp*`, group, value = integer(0)) : 
  replacement has 0 rows, data has 49 

封装的函数代码：

limmaNorm = function(meta_data, dataset, column){
  
  d = DGEList(dataset)
  group = as.factor(meta_data$column)
  d$samples$group = group
  keep.exprs = edgeR::filterByExpr(d, group = group)
  d = d[keep.exprs,, keep.lib.sizes = FALSE]
  d = calcNormFactors(d , method = 'TMM' )
  return(d)
  
}

附数据结构示例：
raw数据框：

structure(list(Pt1 = c(436L, 0L, 565L, 67L, 6824L, 239L, 1045L
), Pt10 = c(306L, 0L, 980L, 338L, 21827L, 192L, 706L), Pt101 = c(279L, 
0L, 1630L, 228L, 2485L, 74L, 1189L), Pt103 = c(1124L, 0L, 1069L, 
23L, 1337L, 129L, 484L), Pt106 = c(290L, 3L, 800L, 3108L, 9424L, 
159L, 1184L), Pt11 = c(108L, 0L, 429L, 69L, 5105L, 182L, 907L
), Pt17 = c(178L, 0L, 244L, 36L, 1481L, 151L, 1049L)), row.names = c("A1BG", 
"NAT2", "ADA", "CDH2", "AKT3", "ZBTB11-AS1", "MED6"), class = "data.frame")

metadata数据框：

structure(list(Cohort = c("NIV3-PROG", "NIV3-NAIVE", "NIV3-PROG", 
"NIV3-PROG", "NIV3-PROG", "NIV3-NAIVE", "NIV3-PROG"), Response = c("PD", 
"SD", "PR", "PD", "PD", "PD", "PD"), `Dead/Alive
(Dead = True)` = c(TRUE, 
TRUE, FALSE, FALSE, TRUE, TRUE, TRUE), `Time to Death
(weeks)` = c(22.85714286, 
36.57142857, 119.1428571, 69.14285714, 13, 119.5714286, 8.142857143
), Subtype = c("CUTANEOUS", "CUTANEOUS", "CUTANEOUS", "CUTANEOUS", 
"MUCOSAL", "CUTANEOUS", "CUTANEOUS"), `Mutational
 Subtype` = c("NA", 
"NF1", "TripleWt", "TripleWt", "BRAF", "BRAF", "TripleWt"), `M Stage` = c("M1C", 
"M1A", "M1A", "M1B", "M1C", "NA", "M1C"), `Mutation Load` = c("NA", 
"75", "10", "21", "700", "106", "20"), `Neo-antigen Load` = c("NA", 
"33", "5", "5", "219", "67", "13"), `Neo-peptide Load` = c("NA", 
"56", "6", "11", "273", "187", "32"), `Cytolytic Score` = c("977.86911190000001", 
"65.840716889999996", "1392.1422339999999", "1108.8620289999999", 
"645.54163300000005", "602.6740413", "20.904544959999999"), Benefit = c("NoResponse", 
"NoResponse", "Response", "NoResponse", "NoResponse", "NoResponse", 
"NoResponse")), row.names = c("Pt1", "Pt10", "Pt101", "Pt103", 
"Pt106", "Pt11", "Pt17"), class = "data.frame")

报错原因 #

问题出在函数中取分组列的代码：

group = as.factor(meta_data$column)

$符号是按字面列名提取数据框列，这里它会尝试寻找名为column的列，但你的metadata中并没有这个列，因此得到的是NULL，转成因子后就是长度为0的空向量。

当把这个空向量赋值给d$samples$group时，空向量的行数（0）和样本数据的行数（49）不匹配，就触发了报错。

解决方法 #

改用[[ ]]来提取列，它可以解析字符串变量，正确获取参数column指定的列：

limmaNorm = function(meta_data, dataset, column){
  
  d = DGEList(dataset)
  # 用[[ ]]替代$，解析字符串变量column
  group = as.factor(meta_data[[column]])
  d$samples$group = group
  keep.exprs = edgeR::filterByExpr(d, group = group)
  d = d[keep.exprs,, keep.lib.sizes = FALSE]
  d = calcNormFactors(d , method = 'TMM' )
  return(d)
  
}

调用修正后的函数即可正常执行。

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