批量处理mRNA通路数据时遇类型转换错误的技术求助

需求说明 #

需要对数据框列表mrna.pathways中的每个数据框，迭代执行针对单个通路（如ADIPOGENESIS）的处理代码。

单通路数据框处理代码示例 #

library(biomaRt)
library(dplyr)
library(tidyr)
library(tibble)
library(ELMER)

exp <- mrna.pathways$ADIPOGENESIS
genes <- rownames(exp)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", GRCh = 37, verbose = T)
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
g_list <- getBM(attributes = c('ensembl_gene_id','hgnc_symbol'), 
                filters='hgnc_symbol',values=genes,
                mart=ensembl, useCache = FALSE)
g_list <- g_list[!duplicated(g_list$ensembl_gene_id) | duplicated(g_list$ensembl_gene_id, fromLast = TRUE), ]

idx <- exp %>%
  rownames_to_column(var = 'symbol') %>%
  pull(symbol) %>%
  match(table = g_list$hgnc_symbol)
exp <- exp %>%
  mutate(row_name = g_list$ensembl_gene_id[idx]) %>%
  group_by(row_name) %>%
  mutate(across(everything(), mean)) %>%
  drop_na() %>%
  ungroup() %>%
  rownames_to_column() %>%
  column_to_rownames('row_name')

exp$rowname <- NULL
exp <- exp[, colnames(exp) %in% colnames(meth)]

exp <- na.omit(exp)
meth <- na.omit(meth)

sample.info <- S4Vectors::DataFrame(primary=colnames(meth))
TSS <- getTSS(genome="hg38", TSS=list(upstream=2000, downstream=2000))

promoter.features <- get.feature.probe(feature = NULL,TSS,genome = "hg38",met.platform = "450K",TSS.range = list(upstream = 2000, downstream = 2000),promoter = T,rm.chr = NULL)

mae <- createMAE(exp = exp, met = meth, met.platform ="450K",colData = sample.info,filter.probes = promoter.features,genome = "hg38") 

meth.adipogenesis <- meth[rownames(meth) %in% mae@ExperimentList@listData$`DNA methylation`@rowRanges@ranges@NAMES,]
meth.adipogenesis <- meth.adipogenesis[ , order(names(meth.adipogenesis))]

尝试的迭代代码 #

ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", GRCh = 37, verbose = T)
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)
g_list <- getBM(attributes = c('ensembl_gene_id','hgnc_symbol'), 
                filters='hgnc_symbol',values=genes,
                mart=ensembl, useCache = FALSE)
g_list <- g_list[!duplicated(g_list$ensembl_gene_id) | duplicated(g_list$ensembl_gene_id, fromLast = TRUE), ]

lapply(mrna.pathways, function(x) {
  genes <- rownames(x)
  idx <- x %>%
    rownames_to_column(var = 'symbol') %>%
    pull(symbol) %>%
    match(table = g_list$hgnc_symbol)
  x <- x %>%
    mutate(row_name = g_list$ensembl_gene_id[idx]) %>%
    group_by(row_name) %>%
    mutate(across(everything(), mean)) %>%
    drop_na() %>%
    ungroup() %>%
    rownames_to_column() %>%
    column_to_rownames('row_name')
  
  x$rowname <- NULL
  x <- x[, colnames(x) %in% colnames(meth)]
  
  x <- na.omit(x)
  meth <- na.omit(meth)
  
  sample.info <- S4Vectors::DataFrame(primary=colnames(meth))
  TSS <- getTSS(genome="hg38", TSS=list(upstream=2000, downstream=2000))
  
  promoter.features <- get.feature.probe(feature = NULL,TSS,genome = "hg38",met.platform = "450K",TSS.range = list(upstream = 2000, downstream = 2000),promoter = T,rm.chr = NULL)
  
  mae <- createMAE(exp = x, met = meth, met.platform ="450K",colData = sample.info,filter.probes = promoter.features,genome = "hg38") 
  
  # methylation.pathways <- meth[rownames(meth) %in% mae@ExperimentList@listData$`DNA methylation`@rowRanges@ranges@NAMES,]
  # methylation.pathways <- methylation.pathways[ , order(names(methylation.pathways))]
  x} )

报错信息 #

Error in makeSummarizedExperimentFromDataFrame(aux[, !grepl("external_gene_name|ensembl_gene_id|entrezgene",  : 
  failed to coerce non-range columns to 'numeric'
In addition: Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :

# List of dataframes
methylation.pathways <- list(HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_WNT_BETA_CATENIN_SIGNALING)

输入数据 #

甲基化数据示例（meth） #

> dput(meth[1:5,1:5])
structure(list(TCGA.Y8.A8S1.01 = c(2.18703449627152, 2.19469161550175, 
-4.08698324454124, -4.1561371941108, -3.45287644116425), TCGA.Y8.A8RZ.01 = c(2.32778952882647, 
2.34278400744779, -3.75343850250278, -4.30774465586903, -4.32830755550417
), TCGA.Y8.A8RY.01 = c(2.52979099210225, 3.14802506855705, -2.14602292062711, 
-4.36197762548339, -4.40610672766836), TCGA.Y8.A897.01 = c(2.14168822194255, 
1.80634363823259, -3.47041057815363, -4.29818228422328, -4.41784030526389
), TCGA.Y8.A896.01 = c(2.65073982842331, 2.61095380298199, -1.89501392291243, 
-4.14543318132608, -4.11547721741994)), row.names = c("cg00000957", 
"cg00001349", "cg00001583", "cg00002028", "cg00002719"), class = "data.frame")

通路基因数据示例（hallmark） #

> hallmark
structure(list(V3 = c("JUNB", "PGK1", "FDPS", "ARHGEF2", "CD74"
), V4 = c("CXCL2", "PDK1", "CYP51A1", "CLASP1", "CTNNB1"), V5 = c("ATF3", 
"GBE1", "IDI1", "KIF11", "JAG2"), V6 = c("NFKBIA", "PFKL", "FDFT1", 
"KIF23", "NOTCH1"), V7 = c("ALDOA", "JUNB", "DHCR7", "ALS2", 
"DLL1")), row.names = c("HALLMARK_TNFA_SIGNALING_VIA_NFKB", "HALLMARK_HYPOXIA", 
"HALLMARK_CHOLESTEROL_HOMEOSTASIS", "HALLMARK_MITOTIC_SPINDLE", 
"HALLMARK_WNT_BETA_CATENIN_SIGNALING"), class = "data.frame")

参考文献内容翻译 #

矩阵B的行为样本，列为某一通路的基因。通过主成分分析（PCA），将矩阵B分解为不相关成分，得到Gᵢₚ∈Rⁿ×ᵠ，其中q为主成分（PCs）的数量。该通路pᵢ的处理任务也在CNV矩阵（C）和DNA甲基化矩阵（M）上执行，分别得到Cᵢₚ、Mᵢₚ∈Rⁿ×ᵠ。对全部146条通路重复此过程，最终得到合并矩阵Gₚ、Cₚ和Mₚ∈Rⁿ×¹⁴⁶ᵠ。

问题分析与解决建议 #

1. 报错根源：createMAE函数要求表达矩阵（exp）的所有列必须是数值型，但当前处理后的x可能包含非数值列，或列类型转换失败。
2. 关键问题点：
   • 迭代代码中，g_list提前用未定义的全局genes变量生成，导致基因匹配错误，应在每个迭代内部根据当前数据框行名获取基因列表后再调用getBM。
   • meth <- na.omit(meth)放在迭代函数内部，会修改全局甲基化数据，导致后续迭代样本不匹配，应改为使用局部副本。
   • 处理后的x可能残留非数值列或未完全清除NA，导致createMAE无法转换为数值矩阵。
3. 修正步骤：
   • 将getBM相关代码移到lapply函数内部，确保每个通路用自身基因列表获取ensembl ID映射。
   • 全局预处理甲基化数据，迭代内部使用局部副本避免修改全局对象。
   • 强制转换表达矩阵所有列为数值型，清除非数值列。
   • 添加空矩阵判断，避免无效的createMAE调用。

修正后的迭代代码示例：

library(biomaRt)
library(dplyr)
library(tidyr)
library(tibble)
library(ELMER)

# 初始化biomaRt连接（只需一次）
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", GRCh = 37, verbose = T)
ensembl <- useDataset(dataset = "hsapiens_gene_ensembl", mart = ensembl)

# 全局预处理甲基化数据
meth_clean <- na.omit(meth)
sample.info <- S4Vectors::DataFrame(primary=colnames(meth_clean))
TSS <- getTSS(genome="hg38", TSS=list(upstream=2000, downstream=2000))
promoter.features <- get.feature.probe(feature = NULL,TSS,genome = "hg38",met.platform = "450K",TSS.range = list(upstream = 2000, downstream = 2000),promoter = T,rm.chr = NULL)

# 迭代处理每个通路
methylation.pathways <- lapply(mrna.pathways, function(x) {
  genes <- rownames(x)
  # 针对当前通路基因获取ensembl ID映射
  g_list <- getBM(attributes = c('ensembl_gene_id','hgnc_symbol'), 
                  filters='hgnc_symbol',values=genes,
                  mart=ensembl, useCache = FALSE)
  g_list <- g_list[!duplicated(g_list$ensembl_gene_id) | duplicated(g_list$ensembl_gene_id, fromLast = TRUE), ]
  
  idx <- x %>%
    rownames_to_column(var = 'symbol') %>%
    pull(symbol) %>%
    match(table = g_list$hgnc_symbol)
  
  x_processed <- x %>%
    mutate(row_name = g_list$ensembl_gene_id[idx]) %>%
    drop_na(row_name) %>%
    group_by(row_name) %>%
    mutate(across(everything(), mean)) %>%
    ungroup() %>%
    column_to_rownames('row_name') %>%
    mutate(across(everything(), as.numeric)) %>%
    select(intersect(colnames(.), colnames(meth_clean))) %>%
    na.omit()
  
  if(nrow(x_processed) > 0 && ncol(x_processed) > 0){
    mae <- createMAE(exp = x_processed, met = meth_clean, met.platform ="450K",colData = sample.info,filter.probes = promoter.features,genome = "hg38") 
    meth_filtered <- meth_clean[rownames(meth_clean) %in% rownames(mae[["DNA methylation"]]), ]
    meth_filtered <- meth_filtered[, order(colnames(meth_filtered))]
    return(meth_filtered)
  } else {
    warning("Empty expression matrix after processing, returning NULL")
    return(NULL)
  }
})

# 移除空结果
methylation.pathways <- Filter(Negate(is.null), methylation.pathways)

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