7

我有一个大约 35,000 行，7 列的数据框。它看起来像这样：

头(nuc)

  chr feature    start      end   gene_id    pctAT    pctGC length
1   1     CDS 67000042 67000051 NM_032291 0.600000 0.400000     10
2   1     CDS 67091530 67091593 NM_032291 0.609375 0.390625     64
3   1     CDS 67098753 67098777 NM_032291 0.600000 0.400000     25
4   1     CDS 67101627 67101698 NM_032291 0.472222 0.527778     72
5   1     CDS 67105460 67105516 NM_032291 0.631579 0.368421     57
6   1     CDS 67108493 67108547 NM_032291 0.436364 0.563636     55

gene_id 是一个因子，它有大约 3,500 个独特的级别。我想，对于每个级别的gene_id 得到min(start), max(end), mean(pctAT), mean(pctGC), 和sum(length).

我尝试使用 lapply 和 do.call ，但它需要永远 +30 分钟才能运行。我正在使用的代码是：

nuc_prof = lapply(levels(nuc$gene_id), function(gene){
  t = nuc[nuc$gene_id==gene, ]
  return(list(gene_id=gene, start=min(t$start), end=max(t$end), pctGC =
              mean(t$pctGC), pct = mean(t$pctAT), cdslength = sum(t$length))) 
})
nuc_prof = do.call(rbind, nuc_prof)

我确定我做错了什么来减慢速度。我还没有等待它完成，因为我确信它可以更快。有任何想法吗？

4

3 回答 3

15

由于我处于传福音的心情……这就是快速data.table解决方案的样子：

library(data.table)
dt <- data.table(nuc, key="gene_id")

dt[,list(A=min(start),
         B=max(end),
         C=mean(pctAT),
         D=mean(pctGC),
         E=sum(length)), by=key(dt)]
#      gene_id        A        B         C         D   E
# 1: NM_032291 67000042 67108547 0.5582567 0.4417433 283
# 2:       ZZZ 67000042 67108547 0.5582567 0.4417433 283
于 2012-06-15T16:14:22.693 回答
8

do.call在大型物体上可能会非常慢。我认为这是由于它如何构建调用，但我不确定。一个更快的选择是data.table包。或者，正如@Andrie 在评论中建议的那样，tapply用于每个计算和cbind结果。

关于您当前实现的注释：您可以使用该split函数将您的 data.frame 分解为您可以循环的 data.frames 列表，而不是在您的函数中进行子集设置。

g <- function(tnuc) {
  list(gene_id=tnuc$gene_id[1], start=min(tnuc$start), end=max(tnuc$end),
       pctGC=mean(tnuc$pctGC), pct=mean(tnuc$pctAT), cdslength=sum(tnuc$length))
}
nuc_prof <- lapply(split(nuc, nuc$gene_id), g)
于 2012-06-15T16:05:36.813 回答
5

正如其他人所提到的 -do.call大型对象存在问题，我最近发现它在大型数据集上到底有多慢。为了说明这个问题，这里有一个使用带有大型回归对象的简单摘要调用的基准测试（使用 rms 包的 cox 回归）：

> model <- cph(Surv(Time, Status == "Cardiovascular") ~  
+              Group + rcs(Age, 3) + cluster(match_group),
+              data=full_df, 
+              x=TRUE, y=TRUE)

> system.time(s_reg <- summary(object = model))
   user  system elapsed 
   0.00    0.02    0.03 
> system.time(s_dc <- do.call(summary, list(object = model)))
   user  system elapsed 
 282.27    0.08  282.43 
> nrow(full_df)
[1] 436305

虽然该data.table解决方案是一个很好的方法，但它不包含完整的功能，do.call因此我想我会分享我的fastDoCall功能 - Hadley Wickhams 的修改建议在 R-mailing 列表中进行hack 。它在 Gmisc-package 1.0-version 中可用（尚未在 CRAN 上发布，但您可以在此处找到它）。基准是：

> system.time(s_fc <- fastDoCall(summary, list(object = model)))
   user  system elapsed 
   0.03    0.00    0.06

该函数的完整代码如下：

fastDoCall <- function(what, args, quote = FALSE, envir = parent.frame()){
  if (quote)
    args <- lapply(args, enquote)

  if (is.null(names(args))){
    argn <- args
    args <- list()
  }else{
    # Add all the named arguments
    argn <- lapply(names(args)[names(args) != ""], as.name)
    names(argn) <- names(args)[names(args) != ""]
    # Add the unnamed arguments
    argn <- c(argn, args[names(args) == ""])
    args <- args[names(args) != ""]
  }

  if (class(what) == "character"){
    if(is.character(what)){
      fn <- strsplit(what, "[:]{2,3}")[[1]]
      what <- if(length(fn)==1) {
        get(fn[[1]], envir=envir, mode="function")
      } else {
        get(fn[[2]], envir=asNamespace(fn[[1]]), mode="function")
      }
    }
    call <- as.call(c(list(what), argn))
  }else if (class(what) == "function"){
    f_name <- deparse(substitute(what))
    call <- as.call(c(list(as.name(f_name)), argn))
    args[[f_name]] <- what
  }else if (class(what) == "name"){
    call <- as.call(c(list(what, argn)))
  }

  eval(call,
       envir = args,
       enclos = envir)
}
于 2014-08-18T20:27:49.347 回答