RNA-seq经RUVg归一化产生负值导致plotRLE报错的解决方案问询

调用RUVg时输入未预处理的原始整数计数矩阵，RUVg的线性校正逻辑会调整原始计数，原始计数中的0值校正后易出现负值，而plotRLE默认对输入矩阵做对数转换，遇到负值就会触发报错。代码中还有两处语法问题：RUVg参数中的k <- k应为k = k，最后一行绘图的dev.off缺少括号，应为dev.off()。

方案1：预处理原始计数再输入RUVg（更推荐） #

先将原始计数转换为带偏移量的对数CPM矩阵，从源头避免0值校正后出现负值的问题，且更符合差异表达分析的常规预处理逻辑：

library(edgeR)
library(RUVSeq)

runRUVg <- function(dds, symbol=TRUE, k, dirn='.', ...) {
  # 原始计数转log2(CPM+1)，消除0值影响
  dds_mat <- data.matrix(dds)
  dds_norm <- cpm(dds_mat, log=TRUE, prior.count = 1)
  seq <- newSeqExpressionSet(dds_norm)
  
  a <- rownames(seq)
  if (symbol) {
    UsedHKG <- intersect(a,hkg_symbols)  # 同时存在于数据集和管家基因列表的基因
    UnusedHKG <- setdiff(a,hkg_symbols)  # 数据集中不存在的管家基因
  } else { 
    UsedHKG <- intersect(a,hkg_ids) 
    UnusedHKG <- setdiff(a,hkg_ids)
  }
  # 修正k的传参写法
  seqRUVg <- RUVg(seq, UsedHKG, k = k, round = FALSE)
  
  # 导出绘图
  png(filename="unnorm.png",width=2000,height=1600, res=300)
  # 输入已经是对数转换后的值，关闭plotRLE的默认对数转换
  plotRLE(seq, outline=F, main='未归一化', ylim=c(-2.5,2.5), log=FALSE)
  dev.off()
  
  png(filename="RUVgnorm.png",width=2000,height=1600, res=300)
  plotRLE(seqRUVg, outline=F, main='RUVg归一化后', ylim=c(-2.5,2.5), log=FALSE)
  dev.off()
  
  return(seqRUVg)
}

方案2：校正后处理负值 #

如果必须保留原始计数作为RUVg输入，可对校正后的负值做截断处理：

library(RUVSeq)

runRUVg <- function(dds, symbol=TRUE, k, dirn='.', ...) {
  dds <- data.matrix(dds)
  seq <- newSeqExpressionSet(dds)
  a <- rownames(seq)
  if (symbol) {
    UsedHKG <- intersect(a,hkg_symbols)
    UnusedHKG <- setdiff(a,hkg_symbols)
  } else { 
    UsedHKG <- intersect(a,hkg_ids)
    UnusedHKG <- setdiff(a,hkg_ids)
  }
  seqRUVg <- RUVg(seq, UsedHKG, k = k, round = TRUE)
  # 处理负值：所有小于0的计数替换为0
  norm_counts <- counts(seqRUVg)
  norm_counts[norm_counts < 0] <- 0
  seqRUVg_fixed <- newSeqExpressionSet(norm_counts)
  
  # 导出绘图
  png(filename="unnorm.png",width=2000,height=1600, res=300)
  plotRLE(seq, outline=F, main='未归一化', ylim=c(-2.5,2.5))
  dev.off()
  
  png(filename="RUVgnorm.png",width=2000,height=1600, res=300)
  plotRLE(seqRUVg_fixed, outline=F, main='RUVg归一化后', ylim=c(-2.5,2.5))
  dev.off()
  
  return(seqRUVg_fixed)
}

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