使用ComBat_seq进行RNA-seq批次校正时遇pnbinom函数错误求助

问题场景 #

尝试对原始RNA-Seq计数数据进行批次校正，运行ComBat_seq时触发错误。

执行代码 #

library("sva")
batch <- c(rep(1, 9),rep(2, 19))

adjusted <- ComBat_seq(matr_red_gene %>% select(-GENE), batch=batch, group=FALSE)

错误输出 #

Found 2 batches
Using null model in ComBat-seq.
Adjusting
for 0 covariate(s) or covariate level(s)
Estimating
dispersions
Fitting the GLM model
Shrinkage off - using GLM
estimates for parameters
Adjusting the data
Error in
pnbinom(counts_sub[a, b] - 1, mu = old_mu[a, b], size = 1/old_phi[a])
: Non-numeric argument to mathematical function

解决建议 #

• 检查输入数据类型：ComBat_seq要求输入的计数矩阵必须是纯数值型（整数或数值）。执行str(matr_red_gene %>% select(-GENE))查看列类型，若存在字符型列，需确认是否误保留了非计数数据，或把字符型的缺失值（如"NA"）转为真实的NA后再处理。
• 转换为矩阵格式：dplyr的select返回的是data.frame，部分版本的ComBat_seq对矩阵输入兼容性更好。修改代码将输入转为矩阵：

  count_matrix <- as.matrix(matr_red_gene %>% select(-GENE))
  adjusted <- ComBat_seq(count_matrix, batch=batch, group=FALSE)
  

• 排查非数值/缺失值：检查计数数据中是否存在非数值元素或缺失值。执行any(is.na(count_matrix))和any(!sapply(count_matrix, is.numeric))确认，若有缺失值可考虑用na.omit()过滤（需谨慎评估影响），或替换为合理值；若有非数值元素则需清理数据来源。
• 验证离散度估计环节：报错发生在数据校正阶段，可能是离散度（old_phi）计算出现非数值结果。可先手动估计离散度，检查是否有异常值：

  library(edgeR)
  dge <- DGEList(counts=count_matrix)
  dge <- estimateDisp(dge, design=model.matrix(~batch))
  summary(dge$common.dispersion)
  

  若离散度结果异常，需检查样本是否存在极端低表达或异常样本，考虑过滤低表达基因或异常样本后再尝试。

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