使用missMethyl排序Mc行名关联度时遇xtfrm报错及结果异常

需求 #

基于FDR校正p值，对Mc的行名按其与目标因子的关联度排序，识别关联度最低的行名。

执行代码 #

library(missMethyl)

df <- rbind(clin$subtype, meth)
rownames(df)[rownames(df) == "1"] <- "subtype"
rownames(clin.info) <- clin.info$Sample.ID
type <- as.character(head(df, n=1))
type.df <- head(df, n=1)
group <- factor(type, levels=c("1a", "1b", "1c", "2a", "2b", "2c"))
mat <- apply(as.matrix(df), 2, as.numeric)
rownames(mat) <- rownames(df)
rownames(mat) <- make.names(rownames(mat), unique=TRUE) # 保留唯一索引
mat <- mat[-1,]

# 提取Illumina阴性对照数据
INCs <- getINCs(rgSet)
head(INCs)

# 将阴性对照数据添加到M值矩阵
Mc <- rbind(mat,INCs)
# 创建标记阴性对照的向量
ctl1 <- rownames(Mc) %in% rownames(INCs)
table(ctl1)

rfit1 <- RUVfit(Y = Mc, X = group, ctl = ctl1) # 第一阶段分析
rfit2 <- RUVadj(Y = Mc, fit = rfit1)

top1 <- topRUV(rfit2, num=Inf, p.BH = 1)

missMethyl包中topRUV函数定义 #

topRUV <- function (fitsum, number = 10, sort.by = c("p","F.p"), p.BH = 1){
  
  tab <- fitsum$C
  
  if (p.BH < 1) {
    
    sig <- (tab[,grepl("p.BH_", colnames(tab))] < p.BH)
    
    if (any(is.na(sig)))
      sig[is.na(sig)] <- FALSE
    
    if (all(!sig))
      return(data.frame())
    
    tab <- tab[sig,]
  }
  
  sort.by <- match.arg(sort.by)
  
  ord <- switch(sort.by, p = order(tab[,grepl("p_", colnames(tab))], 
                                   decreasing = FALSE), 
                F.p = order(tab$F.p, decreasing=FALSE))
  
  if (nrow(tab) < number)
    number <- nrow(tab)
  
  if (number < 1)
    return(data.frame())
  
  top <- ord[1:number]
  
  tab[top,]
}

运行异常 #

• 警告信息：In xtfrm.data.frame(x) : cannot xtfrm data frames
• 输出结果：top1前6行所有值均为NA，行名为"NA"、"NA.1"等

数据结构示例 #

> dput(meth[1:20,1:20])
structure(c(5973.04545454545, 7167, ..., dimnames = list(c("cg00000957", ...), c("TCGA.Y8.A8S1.01A", ...)))

> dput(clin[1:20,])
structure(list(subtype = c("2a", "1a", ...), age = c("61", "58", ...)), row.names = c("TCGA.Y8.A8S1.01A", ...), class = "data.frame")

问题分析与修正方案 #

核心问题原因 #

1. 数据类型混乱：原代码将字符型的临床分组变量clin$subtype与数值型的甲基化数据meth行绑定，导致整个df变为字符型。后续转数值时，分组行的字符会被转为NA，污染了整个矩阵，最终导致RUV分析结果异常。
2. 排序逻辑错误触发警告：topRUV中order()函数传入了数据框而非数值向量，因为前面的错误导致tab[,grepl("p_", colnames(tab))]返回的是包含NA的异常数据框，而非p值向量，从而触发cannot xtfrm data frames警告。

修正步骤 #

1. 正确构造分组变量

直接从clin数据框提取分组，避免数据类型混合：

# 确保样本顺序与甲基化数据完全匹配
stopifnot(all(colnames(meth) == rownames(clin)))
# 构造分组因子，指定水平顺序
group <- factor(clin$subtype, levels=c("1a", "1b", "1c", "2a", "2b", "2c"))

2. 正确处理甲基化数据矩阵

无需绑定分组变量，直接转换甲基化数据为数值矩阵：

mat <- as.matrix(meth)
# 确保行名唯一（若存在重复探针名）
rownames(mat) <- make.names(rownames(mat), unique=TRUE)

3. 重新执行RUV分析与结果提取

library(missMethyl)

# 提取Illumina阴性对照数据
INCs <- getINCs(rgSet)
# 合并甲基化数据与阴性对照
Mc <- rbind(mat, INCs)
# 标记阴性对照行
ctl1 <- rownames(Mc) %in% rownames(INCs)

# 重新拟合RUV模型
rfit1 <- RUVfit(Y = Mc, X = group, ctl = ctl1)
rfit2 <- RUVadj(Y = Mc, fit = rfit1)

# 按FDR校正p值升序排序（最显著在前，最不显著在后）
top1 <- topRUV(rfit2, num=Inf, p.BH = 1, sort.by = "p")

4. 识别关联度最低的行名

关联度最低对应FDR校正p值最大（最不显著），直接取排序结果的最后一行：

# 获取关联度最低的行名
lowest_assoc_row <- rownames(top1)[nrow(top1)]
# 查看末尾5行确认
tail(top1)

额外注意事项 #

• 原代码中rownames(clin.info) <- clin.info$Sample.ID需确保clin.info对象已定义，且其Sample.ID与甲基化数据的样本ID匹配，否则会导致后续分析样本不对应。
• 若仍存在NA值，可检查甲基化数据是否有缺失值，或阴性对照数据是否正常加载。

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